Double precision

parcels/codegenerator.py

@@ -396,7 +396,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("type_coord", ", ".join(funcvars)))
         if len(transformer.tmp_vars) > 0:
             self.ccode.body.insert(0, c.Value("float", ", ".join(transformer.tmp_vars)))
 

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(float 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(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(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(float x, float y, float 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(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(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(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(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)
 {

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(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(float x, float y, float 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(float x, float y, float 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(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(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(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(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(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(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)

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-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.],

parcels/particle.py

@@ -194,6 +194,12 @@ def __repr__(self):
     def delete(self):
         self.state = ErrorCode.Delete
 
+    @classmethod
+    def set_lonlatdepth_dtype(cls, dtype):
+        cls.lon.dtype = dtype
+        cls.lat.dtype = dtype
+        cls.depth.dtype = dtype
+
 
 class JITParticle(ScipyParticle):
     """Particle class for JIT-based (Just-In-Time) Particle objects

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
@@ -28,10 +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 lonlatdepth_dtype: Floating precision for lon, lat, depth particle coordinates.
+           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')
     """
 
-    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, lonlatdepth_dtype=None, **kwargs):
         self.fieldset = fieldset
         self.fieldset.check_complete()
 
@@ -82,6 +86,14 @@ def convert_to_list(var):
             self.repeatpclass = pclass
             self.repeatkwargs = kwargs
 
+        if lonlatdepth_dtype is None:
+            self.lonlatdepth_dtype = self.lonlatdepth_dtype_from_field_interp_method(fieldset.U)
+        else:
+            self.lonlatdepth_dtype = 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)
         self.ptype = pclass.getPType()
@@ -143,8 +155,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.float32)
-        lat = np.linspace(start[1], finish[1], size, dtype=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]:
             depth = [depth] * size
         return cls(fieldset=fieldset, pclass=pclass, lon=lon, lat=lat, depth=depth, time=time, repeatdt=repeatdt)
@@ -201,6 +215,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 np.float64
+        else:
+            if field.interp_method == 'cgrid_velocity':
+                return np.float64
+        return np.float32
+
     @property
     def size(self):
         return self.particles.size
@@ -290,7 +315,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 = ['-DDOUBLE_COORD_VARIABLES'] if self.lonlatdepth_dtype == np.float64 else None
+                self.kernel.compile(compiler=GNUCompiler(cppargs=cppargs))
                 self.kernel.load_lib()
 
         # Convert all time variables to seconds

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, float *dt)
 {
   float (*data)[f->xdim] = (float (*)[f->xdim]) f->data;
   float u = data[2][1];

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))

tests/test_fieldset.py

@@ -436,7 +436,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):