From 2a9c691e7b9243ad168bf2612069aa8b9f50ed8c Mon Sep 17 00:00:00 2001 From: delandmeter Date: Tue, 26 Feb 2019 18:51:22 +0100 Subject: [PATCH 01/10] change float particle variable to double --- parcels/include/index_search.h | 10 +++++----- parcels/include/parcels.h | 12 ++++++------ parcels/particle.py | 8 ++++---- parcels/particleset.py | 4 ++-- 4 files changed, 17 insertions(+), 17 deletions(-) diff --git a/parcels/include/index_search.h b/parcels/include/index_search.h index d5a48181b2..8e49ee6258 100644 --- a/parcels/include/index_search.h +++ b/parcels/include/index_search.h @@ -38,7 +38,7 @@ typedef enum RECTILINEAR_Z_GRID=0, RECTILINEAR_S_GRID=1, CURVILINEAR_Z_GRID=2, CURVILINEAR_S_GRID=3 } GridCode; -static inline ErrorCode search_indices_vertical_z(float z, int zdim, float *zvals, int *zi, double *zeta) +static inline ErrorCode search_indices_vertical_z(double z, int zdim, float *zvals, int *zi, double *zeta) { if (z < zvals[0] || z > zvals[zdim-1]) {return ERROR_OUT_OF_BOUNDS;} while (*zi < zdim-1 && z > zvals[*zi+1]) ++(*zi); @@ -49,7 +49,7 @@ static inline ErrorCode search_indices_vertical_z(float z, int zdim, float *zval return SUCCESS; } -static inline ErrorCode search_indices_vertical_s(float z, int xdim, int ydim, int zdim, float *zvals, +static inline ErrorCode search_indices_vertical_s(double z, int xdim, int ydim, int zdim, float *zvals, int xi, int yi, int *zi, double xsi, double eta, double *zeta, int z4d, int ti, int tdim, double time, double t0, double t1) { @@ -120,7 +120,7 @@ static inline void reconnect_bnd_indices(int *xi, int *yi, int xdim, int ydim, i } -static inline ErrorCode search_indices_rectilinear(float x, float y, float z, CStructuredGrid *grid, GridCode gcode, +static inline ErrorCode search_indices_rectilinear(double x, double y, double z, CStructuredGrid *grid, GridCode gcode, int *xi, int *yi, int *zi, double *xsi, double *eta, double *zeta, int ti, double time, double t0, double t1) { @@ -213,7 +213,7 @@ static inline ErrorCode search_indices_rectilinear(float x, float y, float z, CS } -static inline ErrorCode search_indices_curvilinear(float x, float y, float z, CStructuredGrid *grid, GridCode gcode, +static inline ErrorCode search_indices_curvilinear(double x, double y, double z, CStructuredGrid *grid, GridCode gcode, int *xi, int *yi, int *zi, double *xsi, double *eta, double *zeta, int ti, double time, double t0, double t1) { @@ -351,7 +351,7 @@ static inline ErrorCode search_indices_curvilinear(float x, float y, float z, CS /* Local linear search to update grid index * params ti, sizeT, time. t0, t1 are only used for 4D S grids * */ -static inline ErrorCode search_indices(float x, float y, float z, CStructuredGrid *grid, +static inline ErrorCode search_indices(double x, double y, double z, CStructuredGrid *grid, int *xi, int *yi, int *zi, double *xsi, double *eta, double *zeta, GridCode gcode, int ti, double time, double t0, double t1) { diff --git a/parcels/include/parcels.h b/parcels/include/parcels.h index 43fd5def8b..c6ef9ea65d 100644 --- a/parcels/include/parcels.h +++ b/parcels/include/parcels.h @@ -101,7 +101,7 @@ static inline ErrorCode spatial_interpolation_tracer_c_grid_3D(int xi, int yi, i } /* Linear interpolation along the time axis */ -static inline ErrorCode temporal_interpolation_structured_grid(float x, float y, float z, double time, CField *f, +static inline ErrorCode temporal_interpolation_structured_grid(double x, double y, double z, double time, CField *f, GridCode gcode, int *xi, int *yi, int *zi, int *ti, float *value, int interp_method) { @@ -303,7 +303,7 @@ static inline ErrorCode spatial_interpolation_UV_c_grid(double xsi, double eta, -static inline ErrorCode temporal_interpolationUV_c_grid(float x, float y, float z, double time, CField *U, CField *V, +static inline ErrorCode temporal_interpolationUV_c_grid(double x, double y, double z, double time, CField *U, CField *V, GridCode gcode, int *xi, int *yi, int *zi, int *ti, float *u, float *v) { @@ -490,7 +490,7 @@ static inline ErrorCode spatial_interpolation_UVW_c_grid(double xsi, double eta, return SUCCESS; } -static inline ErrorCode temporal_interpolationUVW_c_grid(float x, float y, float z, double time, CField *U, CField *V, CField *W, +static inline ErrorCode temporal_interpolationUVW_c_grid(double x, double y, double z, double time, CField *U, CField *V, CField *W, GridCode gcode, int *xi, int *yi, int *zi, int *ti, float *u, float *v, float *w) { @@ -540,7 +540,7 @@ static inline ErrorCode temporal_interpolationUVW_c_grid(float x, float y, float } -static inline ErrorCode temporal_interpolation(float x, float y, float z, double time, CField *f, +static inline ErrorCode temporal_interpolation(double x, double y, double z, double time, CField *f, void *vxi, void *vyi, void *vzi, void *vti, float *value, int interp_method) { @@ -559,7 +559,7 @@ static inline ErrorCode temporal_interpolation(float x, float y, float z, double } } -static inline ErrorCode temporal_interpolationUV(float x, float y, float z, double time, +static inline ErrorCode temporal_interpolationUV(double x, double y, double z, double time, CField *U, CField *V, void *vxi, void *vyi, void *vzi, void *vti, float *valueU, float *valueV, int interp_method) @@ -582,7 +582,7 @@ static inline ErrorCode temporal_interpolationUV(float x, float y, float z, doub } } -static inline ErrorCode temporal_interpolationUVW(float x, float y, float z, double time, +static inline ErrorCode temporal_interpolationUVW(double x, double y, double z, double time, CField *U, CField *V, CField *W, void *vxi, void *vyi, void *vzi, void *vti, float *valueU, float *valueV, float *valueW, int interp_method) diff --git a/parcels/particle.py b/parcels/particle.py index 4fd61f0e94..d5cb3b13dd 100644 --- a/parcels/particle.py +++ b/parcels/particle.py @@ -160,13 +160,13 @@ class ScipyParticle(_Particle): Additional Variables can be added via the :Class Variable: objects """ - lon = Variable('lon', dtype=np.float32) - lat = Variable('lat', dtype=np.float32) - depth = Variable('depth', dtype=np.float32) + lon = Variable('lon', dtype=np.float64) + lat = Variable('lat', dtype=np.float64) + depth = Variable('depth', dtype=np.float64) time = Variable('time', dtype=np.float64) id = Variable('id', dtype=np.int32) fileid = Variable('fileid', dtype=np.int32, to_write=False) - dt = Variable('dt', dtype=np.float32, to_write=False) + dt = Variable('dt', dtype=np.float64, to_write=False) state = Variable('state', dtype=np.int32, initial=ErrorCode.Success, to_write=False) def __init__(self, lon, lat, fieldset, depth=0., time=0., cptr=None): diff --git a/parcels/particleset.py b/parcels/particleset.py index 19d72fbc11..a84fc97fb1 100644 --- a/parcels/particleset.py +++ b/parcels/particleset.py @@ -143,8 +143,8 @@ def from_line(cls, fieldset, pclass, start, finish, size, depth=None, time=None, :param time: Optional start time value for particles. Default is fieldset.U.time[0] :param repeatdt: Optional interval (in seconds) on which to repeat the release of the ParticleSet """ - lon = np.linspace(start[0], finish[0], size, dtype=np.float32) - lat = np.linspace(start[1], finish[1], size, dtype=np.float32) + lon = np.linspace(start[0], finish[0], size, dtype=np.float64) + lat = np.linspace(start[1], finish[1], size, dtype=np.float64) if type(depth) in [int, float]: depth = [depth] * size return cls(fieldset=fieldset, pclass=pclass, lon=lon, lat=lat, depth=depth, time=time, repeatdt=repeatdt) From 1982556e9ff9095077c0a4fe3e9e9c08ea79bf36 Mon Sep 17 00:00:00 2001 From: delandmeter Date: Wed, 27 Feb 2019 15:42:19 +0100 Subject: [PATCH 02/10] fixing broken tests --- tests/customed_header.h | 2 +- tests/test_fieldset.py | 2 +- tests/test_kernel_language.py | 2 +- tests/test_particle_sets.py | 4 ++-- 4 files changed, 5 insertions(+), 5 deletions(-) diff --git a/tests/customed_header.h b/tests/customed_header.h index e861c36c49..d53afbfb4b 100644 --- a/tests/customed_header.h +++ b/tests/customed_header.h @@ -1,5 +1,5 @@ -static inline void func(CField *f, float *lon, float *dt) +static inline void func(CField *f, double *lon, double *dt) { float (*data)[f->xdim] = (float (*)[f->xdim]) f->data; float u = data[2][1]; diff --git a/tests/test_fieldset.py b/tests/test_fieldset.py index ae414f320e..9b30210c13 100644 --- a/tests/test_fieldset.py +++ b/tests/test_fieldset.py @@ -416,7 +416,7 @@ def generate_dataset(xdim, ydim, zdim=1, tdim=1): pset = ParticleSet(fieldset, JITParticle, 0, 0) pset.execute(AdvectionRK4, dt=1) - assert pset[0].lon == 4.5 and pset[0].lat == 10 + assert np.allclose(pset[0].lon, 4.5) and np.allclose(pset[0].lat, 10) def test_fieldset_from_data_gridtypes(xdim=20, ydim=10, zdim=4): diff --git a/tests/test_kernel_language.py b/tests/test_kernel_language.py index a87f4d567b..ce0db24894 100644 --- a/tests/test_kernel_language.py +++ b/tests/test_kernel_language.py @@ -208,7 +208,7 @@ def func(U, lon, dt): if c_inc == 'str': c_include = """ - static inline void func(CField *f, float *lon, float *dt) + static inline void func(CField *f, double *lon, double *dt) { float (*data)[f->xdim] = (float (*)[f->xdim]) f->data; float u = data[2][1]; diff --git a/tests/test_particle_sets.py b/tests/test_particle_sets.py index 3ad1fa2ed5..04865a88e2 100644 --- a/tests/test_particle_sets.py +++ b/tests/test_particle_sets.py @@ -31,8 +31,8 @@ def test_pset_create_lon_lat(fieldset, mode, npart=100): @pytest.mark.parametrize('mode', ['scipy', 'jit']) def test_pset_create_line(fieldset, mode, npart=100): - lon = np.linspace(0, 1, npart, dtype=np.float32) - lat = np.linspace(1, 0, npart, dtype=np.float32) + lon = np.linspace(0, 1, npart) + lat = np.linspace(1, 0, npart) pset = ParticleSet.from_line(fieldset, size=npart, start=(0, 1), finish=(1, 0), pclass=ptype[mode]) 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) From b356a98eb5db23532c4d2a300d1708c46d083404 Mon Sep 17 00:00:00 2001 From: delandmeter Date: Wed, 27 Feb 2019 16:27:24 +0100 Subject: [PATCH 03/10] remove all explicit settings of plon plat to float32 --- tests/test_advection.py | 42 ++++++++--------- tests/test_fieldset_sampling.py | 24 +++++----- tests/test_kernel_execution.py | 32 ++++++------- tests/test_kernel_language.py | 16 +++---- tests/test_particle_sets.py | 80 ++++++++++++++++----------------- tests/test_particles.py | 18 ++++---- tests/test_scripts.py | 2 +- 7 files changed, 107 insertions(+), 107 deletions(-) diff --git a/tests/test_advection.py b/tests/test_advection.py index 4897c902d4..216350f36f 100644 --- a/tests/test_advection.py +++ b/tests/test_advection.py @@ -47,17 +47,17 @@ def test_advection_zonal(lon, lat, depth, mode, npart=10): assert fieldset2D.U.creation_log == 'from_data' pset2D = ParticleSet(fieldset2D, pclass=ptype[mode], - lon=np.zeros(npart, dtype=np.float32) + 20., - lat=np.linspace(0, 80, npart, dtype=np.float32)) + lon=np.zeros(npart) + 20., + lat=np.linspace(0, 80, npart)) pset2D.execute(AdvectionRK4, runtime=delta(hours=2), dt=delta(seconds=30)) assert (np.diff(np.array([p.lon for p in pset2D])) > 1.e-4).all() dimensions['depth'] = depth fieldset3D = FieldSet.from_data(data3D, dimensions, mesh='spherical', transpose=True) pset3D = ParticleSet(fieldset3D, pclass=ptype[mode], - lon=np.zeros(npart, dtype=np.float32) + 20., - lat=np.linspace(0, 80, npart, dtype=np.float32), - depth=np.zeros(npart, dtype=np.float32) + 10.) + lon=np.zeros(npart) + 20., + lat=np.linspace(0, 80, npart), + depth=np.zeros(npart) + 10.) pset3D.execute(AdvectionRK4, runtime=delta(hours=2), dt=delta(seconds=30)) assert (np.diff(np.array([p.lon for p in pset3D])) > 1.e-4).all() @@ -73,8 +73,8 @@ def test_advection_meridional(lon, lat, mode, npart=10): fieldset = FieldSet.from_data(data, dimensions, mesh='spherical', transpose=True) pset = ParticleSet(fieldset, pclass=ptype[mode], - lon=np.linspace(-60, 60, npart, dtype=np.float32), - lat=np.linspace(0, 30, npart, dtype=np.float32)) + lon=np.linspace(-60, 60, npart), + lat=np.linspace(0, 30, npart)) delta_lat = np.diff(np.array([p.lat for p in pset])) pset.execute(AdvectionRK4, runtime=delta(hours=2), dt=delta(seconds=30)) assert np.allclose(np.diff(np.array([p.lat for p in pset])), delta_lat, rtol=1.e-4) @@ -94,9 +94,9 @@ def test_advection_3D(mode, npart=11): fieldset = FieldSet.from_data(data, dimensions, mesh='flat', transpose=True) pset = ParticleSet(fieldset, pclass=ptype[mode], - lon=np.zeros(npart, dtype=np.float32), - lat=np.zeros(npart, dtype=np.float32) + 1e2, - depth=np.linspace(0, 1, npart, dtype=np.float32)) + lon=np.zeros(npart), + lat=np.zeros(npart) + 1e2, + depth=np.linspace(0, 1, npart)) time = delta(hours=2).total_seconds() pset.execute(AdvectionRK4, runtime=time, dt=delta(seconds=30)) assert np.allclose([p.depth*time for p in pset], [p.lon for p in pset], atol=1.e-1) @@ -182,8 +182,8 @@ def fieldset_stationary(xdim=100, ydim=100, maxtime=delta(hours=6)): ('RK45', 1e-5)]) def test_stationary_eddy(fieldset_stationary, mode, method, rtol, npart=1): fieldset = fieldset_stationary - lon = np.linspace(12000, 21000, npart, dtype=np.float32) - lat = np.linspace(12500, 12500, npart, dtype=np.float32) + lon = np.linspace(12000, 21000, npart) + lat = np.linspace(12500, 12500, npart) pset = ParticleSet(fieldset, pclass=ptype[mode], lon=lon, lat=lat) endtime = delta(hours=6).total_seconds() pset.execute(kernel[method], dt=delta(minutes=3), endtime=endtime) @@ -195,9 +195,9 @@ def test_stationary_eddy(fieldset_stationary, mode, method, rtol, npart=1): @pytest.mark.parametrize('mode', ['scipy', 'jit']) def test_stationary_eddy_vertical(mode, npart=1): - lon = np.linspace(12000, 21000, npart, dtype=np.float32) - lat = np.linspace(10000, 20000, npart, dtype=np.float32) - depth = np.linspace(12500, 12500, npart, dtype=np.float32) + lon = np.linspace(12000, 21000, npart) + lat = np.linspace(10000, 20000, npart) + depth = np.linspace(12500, 12500, npart) endtime = delta(hours=6).total_seconds() xdim = ydim = 100 @@ -263,8 +263,8 @@ def fieldset_moving(xdim=100, ydim=100, maxtime=delta(hours=6)): ('RK45', 1e-5)]) def test_moving_eddy(fieldset_moving, mode, method, rtol, npart=1): fieldset = fieldset_moving - lon = np.linspace(12000, 21000, npart, dtype=np.float32) - lat = np.linspace(12500, 12500, npart, dtype=np.float32) + lon = np.linspace(12000, 21000, npart) + lat = np.linspace(12500, 12500, npart) pset = ParticleSet(fieldset, pclass=ptype[mode], lon=lon, lat=lat) endtime = delta(hours=6).total_seconds() pset.execute(kernel[method], dt=delta(minutes=3), endtime=endtime) @@ -307,8 +307,8 @@ def fieldset_decaying(xdim=100, ydim=100, maxtime=delta(hours=6)): ('RK45', 1e-5)]) def test_decaying_eddy(fieldset_decaying, mode, method, rtol, npart=1): fieldset = fieldset_decaying - lon = np.linspace(12000, 21000, npart, dtype=np.float32) - lat = np.linspace(12500, 12500, npart, dtype=np.float32) + lon = np.linspace(12000, 21000, npart) + lat = np.linspace(12500, 12500, npart) pset = ParticleSet(fieldset, pclass=ptype[mode], lon=lon, lat=lat) endtime = delta(hours=6).total_seconds() pset.execute(kernel[method], dt=delta(minutes=3), endtime=endtime) @@ -344,8 +344,8 @@ def test_decaying_eddy(fieldset_decaying, mode, method, rtol, npart=1): npart = args.particles pset = ParticleSet(fieldset, pclass=ptype[args.mode], - lon=np.linspace(4000, 21000, npart, dtype=np.float32), - lat=np.linspace(12500, 12500, npart, dtype=np.float32)) + lon=np.linspace(4000, 21000, npart), + lat=np.linspace(12500, 12500, npart)) if args.verbose: print("Initial particle positions:\n%s" % pset) pset.execute(kernel[args.method], dt=delta(minutes=3), runtime=delta(hours=6)) diff --git a/tests/test_fieldset_sampling.py b/tests/test_fieldset_sampling.py index 415e2d026e..190c6ad54a 100644 --- a/tests/test_fieldset_sampling.py +++ b/tests/test_fieldset_sampling.py @@ -213,13 +213,13 @@ def test_fieldset_sample_particle(mode, k_sample_uv, lat_flip, npart=120): fieldset = FieldSet.from_data(data, dimensions, mesh='flat', transpose=True) - lon = np.linspace(-170, 170, npart, dtype=np.float32) - lat = np.linspace(-80, 80, npart, dtype=np.float32) - pset = ParticleSet(fieldset, pclass=pclass(mode), lon=lon, lat=np.zeros(npart, dtype=np.float32) + 70.) + lon = np.linspace(-170, 170, npart) + lat = np.linspace(-80, 80, npart) + pset = ParticleSet(fieldset, pclass=pclass(mode), lon=lon, lat=np.zeros(npart) + 70.) pset.execute(pset.Kernel(k_sample_uv), endtime=1., dt=1.) assert np.allclose(np.array([p.v for p in pset]), lon, rtol=1e-6) - pset = ParticleSet(fieldset, pclass=pclass(mode), lat=lat, lon=np.zeros(npart, dtype=np.float32) - 45.) + pset = ParticleSet(fieldset, pclass=pclass(mode), lat=lat, lon=np.zeros(npart) - 45.) pset.execute(pset.Kernel(k_sample_uv), endtime=1., dt=1.) assert np.allclose(np.array([p.u for p in pset]), lat, rtol=1e-6) @@ -228,14 +228,14 @@ def test_fieldset_sample_particle(mode, k_sample_uv, lat_flip, npart=120): def test_fieldset_sample_geographic(fieldset_geometric, mode, k_sample_uv, npart=120): """ Sample a fieldset with conversion to geographic units (degrees). """ fieldset = fieldset_geometric - lon = np.linspace(-170, 170, npart, dtype=np.float32) - lat = np.linspace(-80, 80, npart, dtype=np.float32) + lon = np.linspace(-170, 170, npart) + lat = np.linspace(-80, 80, npart) - pset = ParticleSet(fieldset, pclass=pclass(mode), lon=lon, lat=np.zeros(npart, dtype=np.float32) + 70.) + pset = ParticleSet(fieldset, pclass=pclass(mode), lon=lon, lat=np.zeros(npart) + 70.) pset.execute(pset.Kernel(k_sample_uv), endtime=1., dt=1.) assert np.allclose(np.array([p.v for p in pset]), lon, rtol=1e-6) - pset = ParticleSet(fieldset, pclass=pclass(mode), lat=lat, lon=np.zeros(npart, dtype=np.float32) - 45.) + pset = ParticleSet(fieldset, pclass=pclass(mode), lat=lat, lon=np.zeros(npart) - 45.) pset.execute(pset.Kernel(k_sample_uv), endtime=1., dt=1.) assert np.allclose(np.array([p.u for p in pset]), lat, rtol=1e-6) @@ -244,14 +244,14 @@ def test_fieldset_sample_geographic(fieldset_geometric, mode, k_sample_uv, npart def test_fieldset_sample_geographic_polar(fieldset_geometric_polar, mode, k_sample_uv, npart=120): """ Sample a fieldset with conversion to geographic units and a pole correction. """ fieldset = fieldset_geometric_polar - lon = np.linspace(-170, 170, npart, dtype=np.float32) - lat = np.linspace(-80, 80, npart, dtype=np.float32) + lon = np.linspace(-170, 170, npart) + lat = np.linspace(-80, 80, npart) - pset = ParticleSet(fieldset, pclass=pclass(mode), lon=lon, lat=np.zeros(npart, dtype=np.float32) + 70.) + pset = ParticleSet(fieldset, pclass=pclass(mode), lon=lon, lat=np.zeros(npart) + 70.) pset.execute(pset.Kernel(k_sample_uv), endtime=1., dt=1.) assert np.allclose(np.array([p.v for p in pset]), lon, rtol=1e-6) - pset = ParticleSet(fieldset, pclass=pclass(mode), lat=lat, lon=np.zeros(npart, dtype=np.float32) - 45.) + pset = ParticleSet(fieldset, pclass=pclass(mode), lat=lat, lon=np.zeros(npart) - 45.) pset.execute(pset.Kernel(k_sample_uv), endtime=1., dt=1.) # Note: 1.e-2 is a very low rtol, so there seems to be a rather # large sampling error for the JIT correction. diff --git a/tests/test_kernel_execution.py b/tests/test_kernel_execution.py index f586fd217e..5081e3ee3b 100644 --- a/tests/test_kernel_execution.py +++ b/tests/test_kernel_execution.py @@ -35,8 +35,8 @@ def fieldset(xdim=20, ydim=20): ]) def test_execution_endtime(fieldset, mode, start, end, substeps, dt, npart=10): pset = ParticleSet(fieldset, pclass=ptype[mode], time=start, - lon=np.linspace(0, 1, npart, dtype=np.float32), - lat=np.linspace(1, 0, npart, dtype=np.float32)) + lon=np.linspace(0, 1, npart), + lat=np.linspace(1, 0, npart)) pset.execute(DoNothing, endtime=end, dt=dt) assert np.allclose(np.array([p.time for p in pset]), end) @@ -52,8 +52,8 @@ def test_execution_endtime(fieldset, mode, start, end, substeps, dt, npart=10): ]) def test_execution_runtime(fieldset, mode, start, end, substeps, dt, npart=10): pset = ParticleSet(fieldset, pclass=ptype[mode], time=start, - lon=np.linspace(0, 1, npart, dtype=np.float32), - lat=np.linspace(1, 0, npart, dtype=np.float32)) + lon=np.linspace(0, 1, npart), + lat=np.linspace(1, 0, npart)) t_step = abs(end - start) / substeps for _ in range(substeps): pset.execute(DoNothing, runtime=t_step, dt=dt) @@ -66,8 +66,8 @@ def test_execution_runtime(fieldset, mode, start, end, substeps, dt, npart=10): def test_pset_execute_dt_0(fieldset, mode, time, dt, npart=2): def SetLat(particle, fieldset, time): particle.lat = .6 - lon = np.linspace(0, 1, npart, dtype=np.float32) - lat = np.linspace(1, 0, npart, dtype=np.float32) + lon = np.linspace(0, 1, npart) + lat = np.linspace(1, 0, npart) pset = ParticleSet(fieldset, pclass=ptype[mode], lon=lon, lat=lat) pset.execute(SetLat, endtime=time, dt=dt) @@ -91,8 +91,8 @@ def TimedFail(particle, fieldset, time): return ErrorCode.Success pset = ParticleSet(fieldset, pclass=ptype[mode], - lon=np.linspace(0, 1, npart, dtype=np.float32), - lat=np.linspace(1, 0, npart, dtype=np.float32)) + lon=np.linspace(0, 1, npart), + lat=np.linspace(1, 0, npart)) error_thrown = False try: pset.execute(TimedFail, endtime=20., dt=2.) @@ -112,8 +112,8 @@ def PythonFail(particle, fieldset, time): return ErrorCode.Success pset = ParticleSet(fieldset, pclass=ptype[mode], - lon=np.linspace(0, 1, npart, dtype=np.float32), - lat=np.linspace(1, 0, npart, dtype=np.float32)) + lon=np.linspace(0, 1, npart), + lat=np.linspace(1, 0, npart)) error_thrown = False try: pset.execute(PythonFail, endtime=20., dt=2.) @@ -131,8 +131,8 @@ def MoveRight(particle, fieldset, time): particle.lon += 0.1 pset = ParticleSet(fieldset, pclass=ptype[mode], - lon=np.linspace(0, 1, npart, dtype=np.float32), - lat=np.linspace(1, 0, npart, dtype=np.float32)) + lon=np.linspace(0, 1, npart), + lat=np.linspace(1, 0, npart)) error_thrown = False try: pset.execute(MoveRight, endtime=10., dt=1.) @@ -152,8 +152,8 @@ def MoveRight(particle, fieldset, time): def MoveLeft(particle, fieldset, time): particle.lon -= 1. - lon = np.linspace(0.05, 0.95, npart, dtype=np.float32) - lat = np.linspace(1, 0, npart, dtype=np.float32) + lon = np.linspace(0.05, 0.95, npart) + lat = np.linspace(1, 0, npart) pset = ParticleSet(fieldset, pclass=ptype[mode], lon=lon, lat=lat) pset.execute(MoveRight, endtime=10., dt=1., recovery={ErrorCode.ErrorOutOfBounds: MoveLeft}) @@ -171,8 +171,8 @@ def MoveRight(particle, fieldset, time): def DeleteMe(particle, fieldset, time): particle.delete() - lon = np.linspace(0.05, 0.95, npart, dtype=np.float32) - lat = np.linspace(1, 0, npart, dtype=np.float32) + lon = np.linspace(0.05, 0.95, npart) + lat = np.linspace(1, 0, npart) pset = ParticleSet(fieldset, pclass=ptype[mode], lon=lon, lat=lat) pset.execute(MoveRight, endtime=10., dt=1., recovery={ErrorCode.ErrorOutOfBounds: DeleteMe}) diff --git a/tests/test_kernel_language.py b/tests/test_kernel_language.py index ce0db24894..fcdba5e874 100644 --- a/tests/test_kernel_language.py +++ b/tests/test_kernel_language.py @@ -41,8 +41,8 @@ def test_expression_int(fieldset, mode, name, expr, result, npart=10): class TestParticle(ptype[mode]): p = Variable('p', dtype=np.float32) pset = ParticleSet(fieldset, pclass=TestParticle, - lon=np.linspace(0., 1., npart, dtype=np.float32), - lat=np.zeros(npart, dtype=np.float32) + 0.5) + lon=np.linspace(0., 1., npart), + lat=np.zeros(npart) + 0.5) pset.execute(expr_kernel('Test%s' % name, pset, expr), endtime=1., dt=1.) assert(np.array([result == particle.p for particle in pset]).all()) @@ -60,8 +60,8 @@ def test_expression_float(fieldset, mode, name, expr, result, npart=10): class TestParticle(ptype[mode]): p = Variable('p', dtype=np.float32) pset = ParticleSet(fieldset, pclass=TestParticle, - lon=np.linspace(0., 1., npart, dtype=np.float32), - lat=np.zeros(npart, dtype=np.float32) + 0.5) + lon=np.linspace(0., 1., npart), + lat=np.zeros(npart) + 0.5) pset.execute(expr_kernel('Test%s' % name, pset, expr), endtime=1., dt=1.) assert(np.array([result == particle.p for particle in pset]).all()) @@ -83,8 +83,8 @@ def test_expression_bool(fieldset, mode, name, expr, result, npart=10): class TestParticle(ptype[mode]): p = Variable('p', dtype=np.float32) pset = ParticleSet(fieldset, pclass=TestParticle, - lon=np.linspace(0., 1., npart, dtype=np.float32), - lat=np.zeros(npart, dtype=np.float32) + 0.5) + lon=np.linspace(0., 1., npart), + lat=np.zeros(npart) + 0.5) pset.execute(expr_kernel('Test%s' % name, pset, expr), endtime=1., dt=1.) if mode == 'jit': assert(np.array([result == (particle.p == 1) for particle in pset]).all()) @@ -188,8 +188,8 @@ def test_random_float(fieldset, mode, rngfunc, rngargs, npart=10): class TestParticle(ptype[mode]): p = Variable('p', dtype=np.float32 if rngfunc == 'randint' else np.float32) pset = ParticleSet(fieldset, pclass=TestParticle, - lon=np.linspace(0., 1., npart, dtype=np.float32), - lat=np.zeros(npart, dtype=np.float32) + 0.5) + lon=np.linspace(0., 1., npart), + lat=np.zeros(npart) + 0.5) series = random_series(npart, rngfunc, rngargs, mode) kernel = expr_kernel('TestRandom_%s' % rngfunc, pset, 'random.%s(%s)' % (rngfunc, ', '.join([str(a) for a in rngargs]))) diff --git a/tests/test_particle_sets.py b/tests/test_particle_sets.py index 04865a88e2..edf18f3bab 100644 --- a/tests/test_particle_sets.py +++ b/tests/test_particle_sets.py @@ -22,8 +22,8 @@ def fieldset(xdim=40, ydim=100): @pytest.mark.parametrize('mode', ['scipy', 'jit']) def test_pset_create_lon_lat(fieldset, mode, npart=100): - lon = np.linspace(0, 1, npart, dtype=np.float32) - lat = np.linspace(1, 0, npart, dtype=np.float32) + lon = np.linspace(0, 1, npart) + lat = np.linspace(1, 0, npart) pset = ParticleSet(fieldset, lon=lon, lat=lat, pclass=ptype[mode]) 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) @@ -40,8 +40,8 @@ def test_pset_create_line(fieldset, mode, npart=100): @pytest.mark.parametrize('mode', ['scipy', 'jit']) def test_pset_create_list_with_customvariable(fieldset, mode, npart=100): - lon = np.linspace(0, 1, npart, dtype=np.float32) - lat = np.linspace(1, 0, npart, dtype=np.float32) + lon = np.linspace(0, 1, npart) + lat = np.linspace(1, 0, npart) class MyParticle(ptype[mode]): v = Variable('v') @@ -100,8 +100,8 @@ def test_pset_create_field_curvi(npart=100): @pytest.mark.parametrize('mode', ['scipy', 'jit']) def test_pset_create_with_time(fieldset, mode, npart=100): - lon = np.linspace(0, 1, npart, dtype=np.float32) - lat = np.linspace(1, 0, npart, dtype=np.float32) + lon = np.linspace(0, 1, npart) + lat = np.linspace(1, 0, npart) time = 5. pset = ParticleSet(fieldset, lon=lon, lat=lat, pclass=ptype[mode], time=time) assert np.allclose([p.time for p in pset], time, rtol=1e-12) @@ -189,8 +189,8 @@ def DoNothing(particle, fieldset, time): @pytest.mark.parametrize('mode', ['scipy', 'jit']) def test_pset_access(fieldset, mode, npart=100): - lon = np.linspace(0, 1, npart, dtype=np.float32) - lat = np.linspace(1, 0, npart, dtype=np.float32) + lon = np.linspace(0, 1, npart) + lat = np.linspace(1, 0, npart) pset = ParticleSet(fieldset, lon=lon, lat=lat, pclass=ptype[mode]) assert(pset.size == 100) assert np.allclose([pset[i].lon for i in range(pset.size)], lon, rtol=1e-12) @@ -208,8 +208,8 @@ def __init__(self, *args, **kwargs): self.n = 2 pset = ParticleSet(fieldset, pclass=TestParticle, - lon=np.linspace(0, 1, npart, dtype=np.float32), - lat=np.linspace(1, 0, npart, dtype=np.float32)) + lon=np.linspace(0, 1, npart), + lat=np.linspace(1, 0, npart)) assert(pset.size == 100) # FIXME: The float test fails with a conversion error of 1.e-8 # assert np.allclose([p.p - 0.33 for p in pset], np.zeros(npart), rtol=1e-12) @@ -218,8 +218,8 @@ def __init__(self, *args, **kwargs): @pytest.mark.parametrize('mode', ['scipy', 'jit']) def test_pset_add_explicit(fieldset, mode, npart=100): - lon = np.linspace(0, 1, npart, dtype=np.float32) - lat = np.linspace(1, 0, npart, dtype=np.float32) + lon = np.linspace(0, 1, npart) + lat = np.linspace(1, 0, npart) pset = ParticleSet(fieldset, lon=[], lat=[], pclass=ptype[mode]) for i in range(npart): particle = ptype[mode](lon=lon[i], lat=lat[i], fieldset=fieldset) @@ -231,8 +231,8 @@ def test_pset_add_explicit(fieldset, mode, npart=100): @pytest.mark.parametrize('mode', ['scipy', 'jit']) def test_pset_add_shorthand(fieldset, mode, npart=100): - lon = np.linspace(0, 1, npart, dtype=np.float32) - lat = np.linspace(1, 0, npart, dtype=np.float32) + lon = np.linspace(0, 1, npart) + lat = np.linspace(1, 0, npart) pset = ParticleSet(fieldset, lon=[], lat=[], pclass=ptype[mode]) for i in range(npart): pset += ptype[mode](lon=lon[i], lat=lat[i], fieldset=fieldset) @@ -257,11 +257,11 @@ def AddLat(particle, fieldset, time): @pytest.mark.parametrize('mode', ['scipy', 'jit']) def test_pset_merge_inplace(fieldset, mode, npart=100): pset1 = ParticleSet(fieldset, pclass=ptype[mode], - lon=np.linspace(0, 1, npart, dtype=np.float32), - lat=np.linspace(1, 0, npart, dtype=np.float32)) + lon=np.linspace(0, 1, npart), + lat=np.linspace(1, 0, npart)) pset2 = ParticleSet(fieldset, pclass=ptype[mode], - lon=np.linspace(0, 1, npart, dtype=np.float32), - lat=np.linspace(0, 1, npart, dtype=np.float32)) + lon=np.linspace(0, 1, npart), + lat=np.linspace(0, 1, npart)) assert(pset1.size == 100) assert(pset2.size == 100) pset1.add(pset2) @@ -272,11 +272,11 @@ def test_pset_merge_inplace(fieldset, mode, npart=100): @pytest.mark.parametrize('mode', ['scipy', 'jit']) def test_pset_merge_duplicate(fieldset, mode, npart=100): pset1 = ParticleSet(fieldset, pclass=ptype[mode], - lon=np.linspace(0, 1, npart, dtype=np.float32), - lat=np.linspace(1, 0, npart, dtype=np.float32)) + lon=np.linspace(0, 1, npart), + lat=np.linspace(1, 0, npart)) pset2 = ParticleSet(fieldset, pclass=ptype[mode], - lon=np.linspace(0, 1, npart, dtype=np.float32), - lat=np.linspace(0, 1, npart, dtype=np.float32)) + lon=np.linspace(0, 1, npart), + lat=np.linspace(0, 1, npart)) pset3 = pset1 + pset2 assert(pset1.size == 100) assert(pset2.size == 100) @@ -285,8 +285,8 @@ def test_pset_merge_duplicate(fieldset, mode, npart=100): @pytest.mark.parametrize('mode', ['scipy', 'jit']) def test_pset_remove_index(fieldset, mode, npart=100): - lon = np.linspace(0, 1, npart, dtype=np.float32) - lat = np.linspace(1, 0, npart, dtype=np.float32) + lon = np.linspace(0, 1, npart) + lat = np.linspace(1, 0, npart) pset = ParticleSet(fieldset, lon=lon, lat=lat, pclass=ptype[mode]) for ilon, ilat in zip(lon[::-1], lat[::-1]): p = pset.remove(-1) @@ -298,8 +298,8 @@ def test_pset_remove_index(fieldset, mode, npart=100): @pytest.mark.xfail(reason="Particle removal has not been implemented yet") @pytest.mark.parametrize('mode', ['scipy', 'jit']) def test_pset_remove_particle(fieldset, mode, npart=100): - lon = np.linspace(0, 1, npart, dtype=np.float32) - lat = np.linspace(1, 0, npart, dtype=np.float32) + lon = np.linspace(0, 1, npart) + lat = np.linspace(1, 0, npart) pset = ParticleSet(fieldset, lon=lon, lat=lat, pclass=ptype[mode]) for ilon, ilat in zip(lon[::-1], lat[::-1]): p = pset.remove(pset[-1]) @@ -315,8 +315,8 @@ def DeleteKernel(particle, fieldset, time): particle.delete() pset = ParticleSet(fieldset, pclass=ptype[mode], - lon=np.linspace(0, 1, npart, dtype=np.float32), - lat=np.linspace(1, 0, npart, dtype=np.float32)) + lon=np.linspace(0, 1, npart), + lat=np.linspace(1, 0, npart)) pset.execute(pset.Kernel(DeleteKernel), endtime=1., dt=1.0) assert(pset.size == 40) @@ -327,8 +327,8 @@ def AddLat(particle, fieldset, time): particle.lat += 0.1 pset = ParticleSet(fieldset, pclass=ptype[mode], - lon=np.linspace(0, 1, npart, dtype=np.float32), - lat=np.zeros(npart, dtype=np.float32)) + lon=np.linspace(0, 1, npart), + lat=np.zeros(npart)) k_add = pset.Kernel(AddLat) for _ in range(n): pset.execute(k_add, runtime=1., dt=1.0) @@ -341,8 +341,8 @@ def AddLat(particle, fieldset, time): particle.lat += 0.1 pset = ParticleSet(fieldset, pclass=ptype[mode], - lon=np.linspace(0, 1, npart, dtype=np.float32), - lat=np.zeros(npart, dtype=np.float32)) + lon=np.linspace(0, 1, npart), + lat=np.zeros(npart)) k_add = pset.Kernel(AddLat) for _ in range(n): pset.execute(k_add, runtime=1., dt=1.0) @@ -370,8 +370,8 @@ def test_density(fieldset, mode, area_scale): def test_pfile_array_remove_particles(fieldset, mode, tmpdir, npart=10): filepath = tmpdir.join("pfile_array_remove_particles") pset = ParticleSet(fieldset, pclass=ptype[mode], - lon=np.linspace(0, 1, npart, dtype=np.float32), - lat=0.5*np.ones(npart, dtype=np.float32)) + lon=np.linspace(0, 1, npart), + lat=0.5*np.ones(npart)) pfile = pset.ParticleFile(filepath) pfile.write(pset, 0) pset.remove(3) @@ -383,8 +383,8 @@ def test_pfile_array_remove_all_particles(fieldset, mode, tmpdir, npart=10): filepath = tmpdir.join("pfile_array_remove_particles") pset = ParticleSet(fieldset, pclass=ptype[mode], - lon=np.linspace(0, 1, npart, dtype=np.float32), - lat=0.5*np.ones(npart, dtype=np.float32)) + lon=np.linspace(0, 1, npart), + lat=0.5*np.ones(npart)) pfile = pset.ParticleFile(filepath) pfile.write(pset, 0) for _ in range(npart): @@ -405,8 +405,8 @@ def Update_v(particle, fieldset, time): class MyParticle(ptype[mode]): v_once = Variable('v_once', dtype=np.float32, initial=0., to_write='once') age = Variable('age', dtype=np.float32, initial=0.) - lon = np.linspace(0, 1, npart, dtype=np.float32) - lat = np.linspace(1, 0, npart, dtype=np.float32) + lon = np.linspace(0, 1, npart) + lat = np.linspace(1, 0, npart) pset = ParticleSet(fieldset, pclass=MyParticle, lon=lon, lat=lat, repeatdt=0.1) pset.execute(pset.Kernel(Update_v), endtime=1, dt=0.1, output_file=pset.ParticleFile(name=filepath, outputdt=0.1)) @@ -428,8 +428,8 @@ def move_west(particle, fieldset, time): def DeleteP(particle, fieldset, time): particle.delete() - lon = np.linspace(0.05, 0.95, npart, dtype=np.float32) - lat = np.linspace(0.95, 0.05, npart, dtype=np.float32) + lon = np.linspace(0.05, 0.95, npart) + lat = np.linspace(0.95, 0.05, npart) (dt, runtime) = (0.1, 0.8) lon_end = lon - runtime/dt*0.1 diff --git a/tests/test_particles.py b/tests/test_particles.py index 2bdb64fc5e..f4fa639786 100644 --- a/tests/test_particles.py +++ b/tests/test_particles.py @@ -24,8 +24,8 @@ class TestParticle(ptype[mode]): p_double = Variable('p_double', dtype=np.float64, initial=11.) p_int = Variable('p_int', dtype=np.int32, initial=12.) pset = ParticleSet(fieldset, pclass=TestParticle, - lon=np.linspace(0, 1, npart, dtype=np.float32), - lat=np.linspace(1, 0, npart, dtype=np.float32)) + lon=np.linspace(0, 1, npart), + lat=np.linspace(1, 0, npart)) def addOne(particle, fieldset, time): particle.p_float += 1. @@ -68,21 +68,21 @@ class TestParticle(ptype[mode]): def test_variable_init_relative(fieldset, mode, npart=10): """Test that checks relative initialisation of custom variables""" class TestParticle(ptype[mode]): - p_base = Variable('p_base', dtype=np.float32, initial=10.) - p_relative = Variable('p_relative', dtype=np.float32, + p_base = Variable('p_base', dtype=np.float64, initial=10.) + p_relative = Variable('p_relative', dtype=np.float64, initial=attrgetter('p_base')) - p_offset = Variable('p_offset', dtype=np.float32, + p_offset = Variable('p_offset', dtype=np.float64, initial=attrgetter('p_base')) - p_lon = Variable('p_lon', dtype=np.float32, + p_lon = Variable('p_lon', dtype=np.float64, initial=attrgetter('lon')) - p_lat = Variable('p_lat', dtype=np.float32, + p_lat = Variable('p_lat', dtype=np.float64, initial=attrgetter('lat')) def __init__(self, *args, **kwargs): super(TestParticle, self).__init__(*args, **kwargs) self.p_offset += 2. - lon = np.linspace(0, 1, npart, dtype=np.float32) - lat = np.linspace(1, 0, npart, dtype=np.float32) + lon = np.linspace(0, 1, npart) + lat = np.linspace(1, 0, npart) pset = ParticleSet(fieldset, pclass=TestParticle, lon=lon, lat=lat) # Adjust base variable to test for aliasing effects for p in pset: diff --git a/tests/test_scripts.py b/tests/test_scripts.py index 31fc87d99c..7a473b4954 100644 --- a/tests/test_scripts.py +++ b/tests/test_scripts.py @@ -16,7 +16,7 @@ def create_outputfiles(dir): endtime = delta(hours=24) x = 3. * (1. / 1.852 / 60) y = (fieldset.U.lat[0] + x, fieldset.U.lat[-1] - x) - lat = np.linspace(y[0], y[1], npart, dtype=np.float32) + lat = np.linspace(y[0], y[1], npart) fp = dir.join("DelayParticle.nc") output_file = pset.ParticleFile(name=fp, outputdt=delaytime) From 8c10d1ad2dddcc0e5c60a6132b5937aa4196c26f Mon Sep 17 00:00:00 2001 From: delandmeter Date: Thu, 28 Feb 2019 16:52:36 +0100 Subject: [PATCH 04/10] tmp particle position are double --- parcels/codegenerator.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/parcels/codegenerator.py b/parcels/codegenerator.py index 47695ea7af..dd21d15967 100644 --- a/parcels/codegenerator.py +++ b/parcels/codegenerator.py @@ -394,7 +394,7 @@ def generate(self, py_ast, funcvars): funcvars.remove(kvar) self.ccode.body.insert(0, c.Value('ErrorCode', 'err')) if len(funcvars) > 0: - self.ccode.body.insert(0, c.Value("float", ", ".join(funcvars))) + self.ccode.body.insert(0, c.Value("double", ", ".join(funcvars))) if len(transformer.tmp_vars) > 0: self.ccode.body.insert(0, c.Value("float", ", ".join(transformer.tmp_vars))) From 694c8df291227e90aba0303c2ad0ee7236f2044c Mon Sep 17 00:00:00 2001 From: delandmeter Date: Thu, 28 Feb 2019 17:35:31 +0100 Subject: [PATCH 05/10] larger tol in rk45 --- parcels/kernels/advection.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/parcels/kernels/advection.py b/parcels/kernels/advection.py index 0a3278f0cb..4c7edc1cd2 100644 --- a/parcels/kernels/advection.py +++ b/parcels/kernels/advection.py @@ -57,8 +57,8 @@ def AdvectionRK45(particle, fieldset, time): Times-step dt is halved if error is larger than tolerance, and doubled if error is smaller than 1/10th of tolerance, with tolerance set to - 1e-9 * dt by default.""" - tol = [1e-9] + 1e-6 * dt by default.""" + tol = [1e-6] c = [1./4., 3./8., 12./13., 1., 1./2.] A = [[1./4., 0., 0., 0., 0.], [3./32., 9./32., 0., 0., 0.], From 9fc6bd1b5b8dbb25bd162f7a79f14e5affa43c15 Mon Sep 17 00:00:00 2001 From: delandmeter Date: Fri, 1 Mar 2019 10:17:41 +0100 Subject: [PATCH 06/10] tol was too small --- parcels/kernels/advection.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/parcels/kernels/advection.py b/parcels/kernels/advection.py index 4c7edc1cd2..a1b1f2cbeb 100644 --- a/parcels/kernels/advection.py +++ b/parcels/kernels/advection.py @@ -57,8 +57,8 @@ def AdvectionRK45(particle, fieldset, time): Times-step dt is halved if error is larger than tolerance, and doubled if error is smaller than 1/10th of tolerance, with tolerance set to - 1e-6 * dt by default.""" - tol = [1e-6] + 1e-5 * dt by default.""" + tol = [1e-5] c = [1./4., 3./8., 12./13., 1., 1./2.] A = [[1./4., 0., 0., 0., 0.], [3./32., 9./32., 0., 0., 0.], From e4302e9e2a82630895a4e61d44e2a35cb8e60845 Mon Sep 17 00:00:00 2001 From: delandmeter Date: Fri, 1 Mar 2019 18:18:04 +0100 Subject: [PATCH 07/10] coordinates_var_precision is now a param of pset By default it is single precision for A-grid U field and double precision for C-grid --- parcels/codegenerator.py | 2 +- parcels/include/index_search.h | 16 +++++++++++----- parcels/include/parcels.h | 12 ++++++------ parcels/particle.py | 19 +++++++++++++++---- parcels/particleset.py | 27 +++++++++++++++++++++++---- tests/customed_header.h | 2 +- tests/test_kernel_language.py | 6 ++++-- tests/test_particle_sets.py | 24 ++++++++++++------------ tests/test_particles.py | 21 ++++++++++++--------- 9 files changed, 85 insertions(+), 44 deletions(-) diff --git a/parcels/codegenerator.py b/parcels/codegenerator.py index dd21d15967..461cd7c920 100644 --- a/parcels/codegenerator.py +++ b/parcels/codegenerator.py @@ -394,7 +394,7 @@ def generate(self, py_ast, funcvars): funcvars.remove(kvar) self.ccode.body.insert(0, c.Value('ErrorCode', 'err')) if len(funcvars) > 0: - self.ccode.body.insert(0, c.Value("double", ", ".join(funcvars))) + self.ccode.body.insert(0, c.Value("type_coord", ", ".join(funcvars))) if len(transformer.tmp_vars) > 0: self.ccode.body.insert(0, c.Value("float", ", ".join(transformer.tmp_vars))) diff --git a/parcels/include/index_search.h b/parcels/include/index_search.h index 8e49ee6258..d08c0ec3f8 100644 --- a/parcels/include/index_search.h +++ b/parcels/include/index_search.h @@ -10,6 +10,12 @@ extern "C" { #define CHECKERROR(res) do {if (res != SUCCESS) return res;} while (0) +#ifdef DOUBLE_COORD_VARIABLES +typedef double type_coord; +#else +typedef float type_coord; +#endif + typedef struct { int gtype; @@ -38,7 +44,7 @@ typedef enum RECTILINEAR_Z_GRID=0, RECTILINEAR_S_GRID=1, CURVILINEAR_Z_GRID=2, CURVILINEAR_S_GRID=3 } GridCode; -static inline ErrorCode search_indices_vertical_z(double z, int zdim, float *zvals, int *zi, double *zeta) +static inline ErrorCode search_indices_vertical_z(type_coord z, int zdim, float *zvals, int *zi, double *zeta) { if (z < zvals[0] || z > zvals[zdim-1]) {return ERROR_OUT_OF_BOUNDS;} while (*zi < zdim-1 && z > zvals[*zi+1]) ++(*zi); @@ -49,7 +55,7 @@ static inline ErrorCode search_indices_vertical_z(double z, int zdim, float *zva return SUCCESS; } -static inline ErrorCode search_indices_vertical_s(double z, int xdim, int ydim, int zdim, float *zvals, +static inline ErrorCode search_indices_vertical_s(type_coord z, int xdim, int ydim, int zdim, float *zvals, int xi, int yi, int *zi, double xsi, double eta, double *zeta, int z4d, int ti, int tdim, double time, double t0, double t1) { @@ -120,7 +126,7 @@ static inline void reconnect_bnd_indices(int *xi, int *yi, int xdim, int ydim, i } -static inline ErrorCode search_indices_rectilinear(double x, double y, double z, CStructuredGrid *grid, GridCode gcode, +static inline ErrorCode search_indices_rectilinear(type_coord x, type_coord y, type_coord z, CStructuredGrid *grid, GridCode gcode, int *xi, int *yi, int *zi, double *xsi, double *eta, double *zeta, int ti, double time, double t0, double t1) { @@ -213,7 +219,7 @@ static inline ErrorCode search_indices_rectilinear(double x, double y, double z, } -static inline ErrorCode search_indices_curvilinear(double x, double y, double z, CStructuredGrid *grid, GridCode gcode, +static inline ErrorCode search_indices_curvilinear(type_coord x, type_coord y, type_coord z, CStructuredGrid *grid, GridCode gcode, int *xi, int *yi, int *zi, double *xsi, double *eta, double *zeta, int ti, double time, double t0, double t1) { @@ -351,7 +357,7 @@ static inline ErrorCode search_indices_curvilinear(double x, double y, double z, /* Local linear search to update grid index * params ti, sizeT, time. t0, t1 are only used for 4D S grids * */ -static inline ErrorCode search_indices(double x, double y, double z, CStructuredGrid *grid, +static inline ErrorCode search_indices(type_coord x, type_coord y, type_coord z, CStructuredGrid *grid, int *xi, int *yi, int *zi, double *xsi, double *eta, double *zeta, GridCode gcode, int ti, double time, double t0, double t1) { diff --git a/parcels/include/parcels.h b/parcels/include/parcels.h index c6ef9ea65d..a80925fd3b 100644 --- a/parcels/include/parcels.h +++ b/parcels/include/parcels.h @@ -101,7 +101,7 @@ static inline ErrorCode spatial_interpolation_tracer_c_grid_3D(int xi, int yi, i } /* Linear interpolation along the time axis */ -static inline ErrorCode temporal_interpolation_structured_grid(double x, double y, double z, double time, CField *f, +static inline ErrorCode temporal_interpolation_structured_grid(type_coord x, type_coord y, type_coord z, double time, CField *f, GridCode gcode, int *xi, int *yi, int *zi, int *ti, float *value, int interp_method) { @@ -303,7 +303,7 @@ static inline ErrorCode spatial_interpolation_UV_c_grid(double xsi, double eta, -static inline ErrorCode temporal_interpolationUV_c_grid(double x, double y, double z, double time, CField *U, CField *V, +static inline ErrorCode temporal_interpolationUV_c_grid(type_coord x, type_coord y, type_coord z, double time, CField *U, CField *V, GridCode gcode, int *xi, int *yi, int *zi, int *ti, float *u, float *v) { @@ -490,7 +490,7 @@ static inline ErrorCode spatial_interpolation_UVW_c_grid(double xsi, double eta, return SUCCESS; } -static inline ErrorCode temporal_interpolationUVW_c_grid(double x, double y, double z, double time, CField *U, CField *V, CField *W, +static inline ErrorCode temporal_interpolationUVW_c_grid(type_coord x, type_coord y, type_coord z, double time, CField *U, CField *V, CField *W, GridCode gcode, int *xi, int *yi, int *zi, int *ti, float *u, float *v, float *w) { @@ -540,7 +540,7 @@ static inline ErrorCode temporal_interpolationUVW_c_grid(double x, double y, dou } -static inline ErrorCode temporal_interpolation(double x, double y, double z, double time, CField *f, +static inline ErrorCode temporal_interpolation(type_coord x, type_coord y, type_coord z, double time, CField *f, void *vxi, void *vyi, void *vzi, void *vti, float *value, int interp_method) { @@ -559,7 +559,7 @@ static inline ErrorCode temporal_interpolation(double x, double y, double z, dou } } -static inline ErrorCode temporal_interpolationUV(double x, double y, double z, double time, +static inline ErrorCode temporal_interpolationUV(type_coord x, type_coord y, type_coord z, double time, CField *U, CField *V, void *vxi, void *vyi, void *vzi, void *vti, float *valueU, float *valueV, int interp_method) @@ -582,7 +582,7 @@ static inline ErrorCode temporal_interpolationUV(double x, double y, double z, d } } -static inline ErrorCode temporal_interpolationUVW(double x, double y, double z, double time, +static inline ErrorCode temporal_interpolationUVW(type_coord x, type_coord y, type_coord z, double time, CField *U, CField *V, CField *W, void *vxi, void *vyi, void *vzi, void *vti, float *valueU, float *valueV, float *valueW, int interp_method) diff --git a/parcels/particle.py b/parcels/particle.py index d5cb3b13dd..9844c571d5 100644 --- a/parcels/particle.py +++ b/parcels/particle.py @@ -160,13 +160,13 @@ class ScipyParticle(_Particle): Additional Variables can be added via the :Class Variable: objects """ - lon = Variable('lon', dtype=np.float64) - lat = Variable('lat', dtype=np.float64) - depth = Variable('depth', dtype=np.float64) + lon = Variable('lon', dtype=np.float32) + lat = Variable('lat', dtype=np.float32) + depth = Variable('depth', dtype=np.float32) time = Variable('time', dtype=np.float64) id = Variable('id', dtype=np.int32) fileid = Variable('fileid', dtype=np.int32, to_write=False) - dt = Variable('dt', dtype=np.float64, to_write=False) + dt = Variable('dt', dtype=np.float32, to_write=False) state = Variable('state', dtype=np.int32, initial=ErrorCode.Success, to_write=False) def __init__(self, lon, lat, fieldset, depth=0., time=0., cptr=None): @@ -194,6 +194,17 @@ def __repr__(self): def delete(self): self.state = ErrorCode.Delete + @classmethod + def set_coordinate_precision(cls, precision): + if precision == 'double': + cls.lon.dtype = np.float64 + cls.lat.dtype = np.float64 + cls.depth.dtype = np.float64 + else: + cls.lon.dtype = np.float32 + cls.lat.dtype = np.float32 + cls.depth.dtype = np.float32 + class JITParticle(ScipyParticle): """Particle class for JIT-based (Just-In-Time) Particle objects diff --git a/parcels/particleset.py b/parcels/particleset.py index a84fc97fb1..4ab3dff9bf 100644 --- a/parcels/particleset.py +++ b/parcels/particleset.py @@ -5,6 +5,7 @@ from parcels.particlefile import ParticleFile from parcels.tools.loggers import logger from parcels.grid import GridCode +from parcels.field import NestedField, SummedField import numpy as np import progressbar import time as time_module @@ -31,7 +32,7 @@ class ParticleSet(object): Other Variables can be initialised using further arguments (e.g. v=... for a Variable named 'v') """ - def __init__(self, fieldset, pclass=JITParticle, lon=None, lat=None, depth=None, time=None, repeatdt=None, **kwargs): + def __init__(self, fieldset, pclass=JITParticle, lon=None, lat=None, depth=None, time=None, repeatdt=None, coordinates_var_precision=None, **kwargs): self.fieldset = fieldset self.fieldset.check_complete() @@ -82,6 +83,21 @@ def convert_to_list(var): self.repeatpclass = pclass self.repeatkwargs = kwargs + if coordinates_var_precision is None: + if type(fieldset.U) in [SummedField, NestedField]: + self.coordinates_var_precision = 'single' + for f in fieldset.U: + if f.interp_method == 'cgrid_velocity': + self.coordinates_var_precision = 'double' + break + else: + self.coordinates_var_precision = 'double' if fieldset.U.interp_method == 'cgrid_velocity' else 'single' + else: + self.coordinates_var_precision = coordinates_var_precision + assert self.coordinates_var_precision in ['single', 'double'], \ + 'particle coordinate variable precision is either set at single or double' + JITParticle.set_coordinate_precision(self.coordinates_var_precision) + size = len(lon) self.particles = np.empty(size, dtype=pclass) self.ptype = pclass.getPType() @@ -143,8 +159,10 @@ def from_line(cls, fieldset, pclass, start, finish, size, depth=None, time=None, :param time: Optional start time value for particles. Default is fieldset.U.time[0] :param repeatdt: Optional interval (in seconds) on which to repeat the release of the ParticleSet """ - lon = np.linspace(start[0], finish[0], size, dtype=np.float64) - lat = np.linspace(start[1], finish[1], size, dtype=np.float64) + + lonlat_type = np.float64 if fieldset.U.interp_method == 'cgrid_velocity' else np.float32 + lon = np.linspace(start[0], finish[0], size, dtype=lonlat_type) + lat = np.linspace(start[1], finish[1], size, dtype=lonlat_type) if type(depth) in [int, float]: depth = [depth] * size return cls(fieldset=fieldset, pclass=pclass, lon=lon, lat=lat, depth=depth, time=time, repeatdt=repeatdt) @@ -290,7 +308,8 @@ def execute(self, pyfunc=AdvectionRK4, endtime=None, runtime=None, dt=1., # Prepare JIT kernel execution if self.ptype.uses_jit: self.kernel.remove_lib() - self.kernel.compile(compiler=GNUCompiler()) + cppargs = ['-DFLOAT_COORD_VARIABLES'] if self.coordinates_var_precision == 'double' else None + self.kernel.compile(compiler=GNUCompiler(cppargs=cppargs)) self.kernel.load_lib() # Convert all time variables to seconds diff --git a/tests/customed_header.h b/tests/customed_header.h index d53afbfb4b..a9b2072bf3 100644 --- a/tests/customed_header.h +++ b/tests/customed_header.h @@ -1,5 +1,5 @@ -static inline void func(CField *f, double *lon, double *dt) +static inline void func(CField *f, double *lon, float *dt) { float (*data)[f->xdim] = (float (*)[f->xdim]) f->data; float u = data[2][1]; diff --git a/tests/test_kernel_language.py b/tests/test_kernel_language.py index fcdba5e874..9aca269620 100644 --- a/tests/test_kernel_language.py +++ b/tests/test_kernel_language.py @@ -200,7 +200,9 @@ class TestParticle(ptype[mode]): @pytest.mark.parametrize('mode', ['scipy', 'jit']) @pytest.mark.parametrize('c_inc', ['str', 'file']) def test_c_kernel(fieldset, mode, c_inc): - pset = ParticleSet(fieldset, pclass=ptype[mode], lon=[0.5], lat=[0]) + coord_type = 'single' if c_inc == 'str' else 'double' + pset = ParticleSet(fieldset, pclass=ptype[mode], lon=[0.5], lat=[0], + coordinates_var_precision=coord_type) def func(U, lon, dt): u = U.data[0, 2, 1] @@ -208,7 +210,7 @@ def func(U, lon, dt): if c_inc == 'str': c_include = """ - static inline void func(CField *f, double *lon, double *dt) + static inline void func(CField *f, float *lon, float *dt) { float (*data)[f->xdim] = (float (*)[f->xdim]) f->data; float u = data[2][1]; diff --git a/tests/test_particle_sets.py b/tests/test_particle_sets.py index edf18f3bab..9407034d34 100644 --- a/tests/test_particle_sets.py +++ b/tests/test_particle_sets.py @@ -22,8 +22,8 @@ def fieldset(xdim=40, ydim=100): @pytest.mark.parametrize('mode', ['scipy', 'jit']) def test_pset_create_lon_lat(fieldset, mode, npart=100): - lon = np.linspace(0, 1, npart) - lat = np.linspace(1, 0, npart) + 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=ptype[mode]) 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) @@ -31,8 +31,8 @@ def test_pset_create_lon_lat(fieldset, mode, npart=100): @pytest.mark.parametrize('mode', ['scipy', 'jit']) def test_pset_create_line(fieldset, mode, npart=100): - lon = np.linspace(0, 1, npart) - lat = np.linspace(1, 0, npart) + lon = np.linspace(0, 1, npart, dtype=np.float32) + lat = np.linspace(1, 0, npart, dtype=np.float32) pset = ParticleSet.from_line(fieldset, size=npart, start=(0, 1), finish=(1, 0), pclass=ptype[mode]) 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) @@ -40,8 +40,8 @@ def test_pset_create_line(fieldset, mode, npart=100): @pytest.mark.parametrize('mode', ['scipy', 'jit']) def test_pset_create_list_with_customvariable(fieldset, mode, npart=100): - lon = np.linspace(0, 1, npart) - lat = np.linspace(1, 0, npart) + lon = np.linspace(0, 1, npart, dtype=np.float32) + lat = np.linspace(1, 0, npart, dtype=np.float32) class MyParticle(ptype[mode]): v = Variable('v') @@ -189,8 +189,8 @@ def DoNothing(particle, fieldset, time): @pytest.mark.parametrize('mode', ['scipy', 'jit']) def test_pset_access(fieldset, mode, npart=100): - lon = np.linspace(0, 1, npart) - lat = np.linspace(1, 0, npart) + 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=ptype[mode]) assert(pset.size == 100) assert np.allclose([pset[i].lon for i in range(pset.size)], lon, rtol=1e-12) @@ -220,7 +220,7 @@ def __init__(self, *args, **kwargs): def test_pset_add_explicit(fieldset, mode, npart=100): lon = np.linspace(0, 1, npart) lat = np.linspace(1, 0, npart) - pset = ParticleSet(fieldset, lon=[], lat=[], pclass=ptype[mode]) + pset = ParticleSet(fieldset, lon=[], lat=[], pclass=ptype[mode], coordinates_var_precision='double') for i in range(npart): particle = ptype[mode](lon=lon[i], lat=lat[i], fieldset=fieldset) pset.add(particle) @@ -231,8 +231,8 @@ def test_pset_add_explicit(fieldset, mode, npart=100): @pytest.mark.parametrize('mode', ['scipy', 'jit']) def test_pset_add_shorthand(fieldset, mode, npart=100): - lon = np.linspace(0, 1, npart) - lat = np.linspace(1, 0, npart) + lon = np.linspace(0, 1, npart, dtype=np.float32) + lat = np.linspace(1, 0, npart, dtype=np.float32) pset = ParticleSet(fieldset, lon=[], lat=[], pclass=ptype[mode]) for i in range(npart): pset += ptype[mode](lon=lon[i], lat=lat[i], fieldset=fieldset) @@ -287,7 +287,7 @@ def test_pset_merge_duplicate(fieldset, mode, npart=100): def test_pset_remove_index(fieldset, mode, npart=100): lon = np.linspace(0, 1, npart) lat = np.linspace(1, 0, npart) - pset = ParticleSet(fieldset, lon=lon, lat=lat, pclass=ptype[mode]) + pset = ParticleSet(fieldset, lon=lon, lat=lat, pclass=ptype[mode], coordinates_var_precision='double') for ilon, ilat in zip(lon[::-1], lat[::-1]): p = pset.remove(-1) assert(p.lon == ilon) diff --git a/tests/test_particles.py b/tests/test_particles.py index f4fa639786..253cafbe0c 100644 --- a/tests/test_particles.py +++ b/tests/test_particles.py @@ -65,25 +65,28 @@ class TestParticle(ptype[mode]): @pytest.mark.parametrize('mode', ['scipy', 'jit']) -def test_variable_init_relative(fieldset, mode, npart=10): +@pytest.mark.parametrize('coord_type', ['single', 'double']) +def test_variable_init_relative(fieldset, mode, coord_type, npart=10): """Test that checks relative initialisation of custom variables""" + lonlat_type = np.float64 if coord_type == 'double' else np.float32 + class TestParticle(ptype[mode]): - p_base = Variable('p_base', dtype=np.float64, initial=10.) - p_relative = Variable('p_relative', dtype=np.float64, + p_base = Variable('p_base', dtype=lonlat_type, initial=10.) + p_relative = Variable('p_relative', dtype=lonlat_type, initial=attrgetter('p_base')) - p_offset = Variable('p_offset', dtype=np.float64, + p_offset = Variable('p_offset', dtype=lonlat_type, initial=attrgetter('p_base')) - p_lon = Variable('p_lon', dtype=np.float64, + p_lon = Variable('p_lon', dtype=lonlat_type, initial=attrgetter('lon')) - p_lat = Variable('p_lat', dtype=np.float64, + p_lat = Variable('p_lat', dtype=lonlat_type, initial=attrgetter('lat')) def __init__(self, *args, **kwargs): super(TestParticle, self).__init__(*args, **kwargs) self.p_offset += 2. - lon = np.linspace(0, 1, npart) - lat = np.linspace(1, 0, npart) - pset = ParticleSet(fieldset, pclass=TestParticle, lon=lon, lat=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, coordinates_var_precision=coord_type) # Adjust base variable to test for aliasing effects for p in pset: p.p_base += 3. From aa047b1e27afec3856f63fd0c541303b9a1e08f1 Mon Sep 17 00:00:00 2001 From: delandmeter Date: Fri, 1 Mar 2019 18:22:46 +0100 Subject: [PATCH 08/10] add doc string for new parameter --- parcels/particleset.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/parcels/particleset.py b/parcels/particleset.py index 4ab3dff9bf..e70bf9982f 100644 --- a/parcels/particleset.py +++ b/parcels/particleset.py @@ -29,6 +29,9 @@ class ParticleSet(object): :param depth: Optional list of initial depth values for particles. Default is 0m :param time: Optional list of initial time values for particles. Default is fieldset.U.grid.time[0] :param repeatdt: Optional interval (in seconds) on which to repeat the release of the ParticleSet + :param coordinates_var_precision: Floating precision for lon, lat, depth particle coordinates. + It is either 'single' or 'double'. Default is 'single' if fieldset.U.interp_method is 'linear' + and 'double' if the interpolation method is 'cgrid_velocity' Other Variables can be initialised using further arguments (e.g. v=... for a Variable named 'v') """ From dc575cb16a490a892482dc3c1b66338be6ec41c0 Mon Sep 17 00:00:00 2001 From: delandmeter Date: Mon, 4 Mar 2019 11:04:45 +0100 Subject: [PATCH 09/10] renaming coordinates_var_precision + Erik's comments --- parcels/particleset.py | 39 ++++++++++++++++++++--------------- tests/test_kernel_language.py | 2 +- tests/test_particle_sets.py | 4 ++-- tests/test_particles.py | 2 +- 4 files changed, 26 insertions(+), 21 deletions(-) diff --git a/parcels/particleset.py b/parcels/particleset.py index 344225f0c6..971c7687f7 100644 --- a/parcels/particleset.py +++ b/parcels/particleset.py @@ -29,13 +29,13 @@ class ParticleSet(object): :param depth: Optional list of initial depth values for particles. Default is 0m :param time: Optional list of initial time values for particles. Default is fieldset.U.grid.time[0] :param repeatdt: Optional interval (in seconds) on which to repeat the release of the ParticleSet - :param coordinates_var_precision: Floating precision for lon, lat, depth particle coordinates. + :param lonlatdepth_dtype: Floating precision for lon, lat, depth particle coordinates. It is either 'single' or 'double'. Default is 'single' if fieldset.U.interp_method is 'linear' and 'double' if the interpolation method is 'cgrid_velocity' Other Variables can be initialised using further arguments (e.g. v=... for a Variable named 'v') """ - def __init__(self, fieldset, pclass=JITParticle, lon=None, lat=None, depth=None, time=None, repeatdt=None, coordinates_var_precision=None, **kwargs): + def __init__(self, fieldset, pclass=JITParticle, lon=None, lat=None, depth=None, time=None, repeatdt=None, lonlatdepth_dtype=None, **kwargs): self.fieldset = fieldset self.fieldset.check_complete() @@ -86,20 +86,13 @@ def convert_to_list(var): self.repeatpclass = pclass self.repeatkwargs = kwargs - if coordinates_var_precision is None: - if type(fieldset.U) in [SummedField, NestedField]: - self.coordinates_var_precision = 'single' - for f in fieldset.U: - if f.interp_method == 'cgrid_velocity': - self.coordinates_var_precision = 'double' - break - else: - self.coordinates_var_precision = 'double' if fieldset.U.interp_method == 'cgrid_velocity' else 'single' + if lonlatdepth_dtype is None: + self.lonlatdepth_dtype = self.lonlatdepth_dtype_from_field_interp_method(fieldset.U) else: - self.coordinates_var_precision = coordinates_var_precision - assert self.coordinates_var_precision in ['single', 'double'], \ - 'particle coordinate variable precision is either set at single or double' - JITParticle.set_coordinate_precision(self.coordinates_var_precision) + self.lonlatdepth_dtype = lonlatdepth_dtype + assert self.lonlatdepth_dtype in ['single', 'double'], \ + 'lon lat depth precision should be set to either single or double' + JITParticle.set_coordinate_precision(self.lonlatdepth_dtype) size = len(lon) self.particles = np.empty(size, dtype=pclass) @@ -163,7 +156,8 @@ def from_line(cls, fieldset, pclass, start, finish, size, depth=None, time=None, :param repeatdt: Optional interval (in seconds) on which to repeat the release of the ParticleSet """ - lonlat_type = np.float64 if fieldset.U.interp_method == 'cgrid_velocity' else np.float32 + lonlatdepth_dtype = cls.lonlatdepth_dtype_from_field_interp_method(fieldset.U) + lonlat_type = np.float64 if lonlatdepth_dtype == 'double' else np.float32 lon = np.linspace(start[0], finish[0], size, dtype=lonlat_type) lat = np.linspace(start[1], finish[1], size, dtype=lonlat_type) if type(depth) in [int, float]: @@ -222,6 +216,17 @@ def from_field(cls, fieldset, pclass, start_field, size, mode='monte_carlo', dep return cls(fieldset=fieldset, pclass=pclass, lon=lon, lat=lat, depth=depth, time=time, repeatdt=repeatdt) + @staticmethod + def lonlatdepth_dtype_from_field_interp_method(field): + if type(field) in [SummedField, NestedField]: + for f in field: + if f.interp_method == 'cgrid_velocity': + return 'double' + else: + if field.interp_method == 'cgrid_velocity': + return 'double' + return 'single' + @property def size(self): return self.particles.size @@ -311,7 +316,7 @@ def execute(self, pyfunc=AdvectionRK4, endtime=None, runtime=None, dt=1., # Prepare JIT kernel execution if self.ptype.uses_jit: self.kernel.remove_lib() - cppargs = ['-DFLOAT_COORD_VARIABLES'] if self.coordinates_var_precision == 'double' else None + cppargs = ['-DFLOAT_COORD_VARIABLES'] if self.lonlatdepth_dtype == 'double' else None self.kernel.compile(compiler=GNUCompiler(cppargs=cppargs)) self.kernel.load_lib() diff --git a/tests/test_kernel_language.py b/tests/test_kernel_language.py index c83887bffd..4ffebd2557 100644 --- a/tests/test_kernel_language.py +++ b/tests/test_kernel_language.py @@ -235,7 +235,7 @@ class TestParticle(ptype[mode]): def test_c_kernel(fieldset, mode, c_inc): coord_type = 'single' if c_inc == 'str' else 'double' pset = ParticleSet(fieldset, pclass=ptype[mode], lon=[0.5], lat=[0], - coordinates_var_precision=coord_type) + lonlatdepth_dtype=coord_type) def func(U, lon, dt): u = U.data[0, 2, 1] diff --git a/tests/test_particle_sets.py b/tests/test_particle_sets.py index 9407034d34..655befdfa6 100644 --- a/tests/test_particle_sets.py +++ b/tests/test_particle_sets.py @@ -220,7 +220,7 @@ def __init__(self, *args, **kwargs): def test_pset_add_explicit(fieldset, mode, npart=100): lon = np.linspace(0, 1, npart) lat = np.linspace(1, 0, npart) - pset = ParticleSet(fieldset, lon=[], lat=[], pclass=ptype[mode], coordinates_var_precision='double') + pset = ParticleSet(fieldset, lon=[], lat=[], pclass=ptype[mode], lonlatdepth_dtype='double') for i in range(npart): particle = ptype[mode](lon=lon[i], lat=lat[i], fieldset=fieldset) pset.add(particle) @@ -287,7 +287,7 @@ def test_pset_merge_duplicate(fieldset, mode, npart=100): def test_pset_remove_index(fieldset, mode, npart=100): lon = np.linspace(0, 1, npart) lat = np.linspace(1, 0, npart) - pset = ParticleSet(fieldset, lon=lon, lat=lat, pclass=ptype[mode], coordinates_var_precision='double') + pset = ParticleSet(fieldset, lon=lon, lat=lat, pclass=ptype[mode], lonlatdepth_dtype='double') for ilon, ilat in zip(lon[::-1], lat[::-1]): p = pset.remove(-1) assert(p.lon == ilon) diff --git a/tests/test_particles.py b/tests/test_particles.py index 253cafbe0c..0fdcb3987b 100644 --- a/tests/test_particles.py +++ b/tests/test_particles.py @@ -86,7 +86,7 @@ def __init__(self, *args, **kwargs): self.p_offset += 2. 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, coordinates_var_precision=coord_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. From 6c7c3c5ae2a054d4de7e7229f4ba4190521dfc35 Mon Sep 17 00:00:00 2001 From: delandmeter Date: Mon, 4 Mar 2019 12:36:51 +0100 Subject: [PATCH 10/10] pset.lonlatdepth_dtype is now either np.float32 or np.float64 --- parcels/particle.py | 13 ++++--------- parcels/particleset.py | 21 ++++++++++----------- tests/test_kernel_language.py | 2 +- tests/test_particle_sets.py | 4 ++-- tests/test_particles.py | 2 +- 5 files changed, 18 insertions(+), 24 deletions(-) diff --git a/parcels/particle.py b/parcels/particle.py index 9844c571d5..4bef1445c1 100644 --- a/parcels/particle.py +++ b/parcels/particle.py @@ -195,15 +195,10 @@ def delete(self): self.state = ErrorCode.Delete @classmethod - def set_coordinate_precision(cls, precision): - if precision == 'double': - cls.lon.dtype = np.float64 - cls.lat.dtype = np.float64 - cls.depth.dtype = np.float64 - else: - cls.lon.dtype = np.float32 - cls.lat.dtype = np.float32 - cls.depth.dtype = np.float32 + def set_lonlatdepth_dtype(cls, dtype): + cls.lon.dtype = dtype + cls.lat.dtype = dtype + cls.depth.dtype = dtype class JITParticle(ScipyParticle): diff --git a/parcels/particleset.py b/parcels/particleset.py index 971c7687f7..4ae3eb2f0f 100644 --- a/parcels/particleset.py +++ b/parcels/particleset.py @@ -30,8 +30,8 @@ class ParticleSet(object): :param time: Optional list of initial time values for particles. Default is fieldset.U.grid.time[0] :param repeatdt: Optional interval (in seconds) on which to repeat the release of the ParticleSet :param lonlatdepth_dtype: Floating precision for lon, lat, depth particle coordinates. - It is either 'single' or 'double'. Default is 'single' if fieldset.U.interp_method is 'linear' - and 'double' if the interpolation method is 'cgrid_velocity' + It is either np.float32 or np.float64. Default is np.float32 if fieldset.U.interp_method is 'linear' + and np.float64 if the interpolation method is 'cgrid_velocity' Other Variables can be initialised using further arguments (e.g. v=... for a Variable named 'v') """ @@ -90,9 +90,9 @@ def convert_to_list(var): self.lonlatdepth_dtype = self.lonlatdepth_dtype_from_field_interp_method(fieldset.U) else: self.lonlatdepth_dtype = lonlatdepth_dtype - assert self.lonlatdepth_dtype in ['single', 'double'], \ - 'lon lat depth precision should be set to either single or double' - JITParticle.set_coordinate_precision(self.lonlatdepth_dtype) + assert self.lonlatdepth_dtype in [np.float32, np.float64], \ + 'lon lat depth precision should be set to either np.float32 or np.float64' + JITParticle.set_lonlatdepth_dtype(self.lonlatdepth_dtype) size = len(lon) self.particles = np.empty(size, dtype=pclass) @@ -156,8 +156,7 @@ def from_line(cls, fieldset, pclass, start, finish, size, depth=None, time=None, :param repeatdt: Optional interval (in seconds) on which to repeat the release of the ParticleSet """ - lonlatdepth_dtype = cls.lonlatdepth_dtype_from_field_interp_method(fieldset.U) - lonlat_type = np.float64 if lonlatdepth_dtype == 'double' else np.float32 + lonlat_type = cls.lonlatdepth_dtype_from_field_interp_method(fieldset.U) lon = np.linspace(start[0], finish[0], size, dtype=lonlat_type) lat = np.linspace(start[1], finish[1], size, dtype=lonlat_type) if type(depth) in [int, float]: @@ -221,11 +220,11 @@ def lonlatdepth_dtype_from_field_interp_method(field): if type(field) in [SummedField, NestedField]: for f in field: if f.interp_method == 'cgrid_velocity': - return 'double' + return np.float64 else: if field.interp_method == 'cgrid_velocity': - return 'double' - return 'single' + return np.float64 + return np.float32 @property def size(self): @@ -316,7 +315,7 @@ def execute(self, pyfunc=AdvectionRK4, endtime=None, runtime=None, dt=1., # Prepare JIT kernel execution if self.ptype.uses_jit: self.kernel.remove_lib() - cppargs = ['-DFLOAT_COORD_VARIABLES'] if self.lonlatdepth_dtype == 'double' else None + cppargs = ['-DDOUBLE_COORD_VARIABLES'] if self.lonlatdepth_dtype == np.float64 else None self.kernel.compile(compiler=GNUCompiler(cppargs=cppargs)) self.kernel.load_lib() diff --git a/tests/test_kernel_language.py b/tests/test_kernel_language.py index 4ffebd2557..b09e851c4d 100644 --- a/tests/test_kernel_language.py +++ b/tests/test_kernel_language.py @@ -233,7 +233,7 @@ class TestParticle(ptype[mode]): @pytest.mark.parametrize('mode', ['scipy', 'jit']) @pytest.mark.parametrize('c_inc', ['str', 'file']) def test_c_kernel(fieldset, mode, c_inc): - coord_type = 'single' if c_inc == 'str' else 'double' + coord_type = np.float32 if c_inc == 'str' else np.float64 pset = ParticleSet(fieldset, pclass=ptype[mode], lon=[0.5], lat=[0], lonlatdepth_dtype=coord_type) diff --git a/tests/test_particle_sets.py b/tests/test_particle_sets.py index 655befdfa6..ef2d7b88dc 100644 --- a/tests/test_particle_sets.py +++ b/tests/test_particle_sets.py @@ -220,7 +220,7 @@ def __init__(self, *args, **kwargs): def test_pset_add_explicit(fieldset, mode, npart=100): lon = np.linspace(0, 1, npart) lat = np.linspace(1, 0, npart) - pset = ParticleSet(fieldset, lon=[], lat=[], pclass=ptype[mode], lonlatdepth_dtype='double') + pset = ParticleSet(fieldset, lon=[], lat=[], pclass=ptype[mode], lonlatdepth_dtype=np.float64) for i in range(npart): particle = ptype[mode](lon=lon[i], lat=lat[i], fieldset=fieldset) pset.add(particle) @@ -287,7 +287,7 @@ def test_pset_merge_duplicate(fieldset, mode, npart=100): def test_pset_remove_index(fieldset, mode, npart=100): lon = np.linspace(0, 1, npart) lat = np.linspace(1, 0, npart) - pset = ParticleSet(fieldset, lon=lon, lat=lat, pclass=ptype[mode], lonlatdepth_dtype='double') + pset = ParticleSet(fieldset, lon=lon, lat=lat, pclass=ptype[mode], lonlatdepth_dtype=np.float64) for ilon, ilat in zip(lon[::-1], lat[::-1]): p = pset.remove(-1) assert(p.lon == ilon) diff --git a/tests/test_particles.py b/tests/test_particles.py index 0fdcb3987b..7155309fe1 100644 --- a/tests/test_particles.py +++ b/tests/test_particles.py @@ -65,7 +65,7 @@ class TestParticle(ptype[mode]): @pytest.mark.parametrize('mode', ['scipy', 'jit']) -@pytest.mark.parametrize('coord_type', ['single', 'double']) +@pytest.mark.parametrize('coord_type', [np.float32, np.float64]) def test_variable_init_relative(fieldset, mode, coord_type, npart=10): """Test that checks relative initialisation of custom variables""" lonlat_type = np.float64 if coord_type == 'double' else np.float32