Repository: OceanParcels/parcels
Pull request #340 (2018-04-11 13:22:00)

Spatially varying Brownian diffusion
Added Kernel for Spatially Varying Brownian motion, following the ideas by Ross & Sharples 2004 and Spagnol et al. 2002

Also added a unit test, although this one is not very strict at the moment

This commit addresses #68


Patch:
diff --git a/parcels/kernels/diffusion.py b/parcels/kernels/diffusion.py
index 9957915a..405791d1 100644
--- a/parcels/kernels/diffusion.py
+++ b/parcels/kernels/diffusion.py
@@ -2,7 +2,7 @@ from parcels import rng as random
 import math
 
 
-__all__ = ['BrownianMotion2D']
+__all__ = ['BrownianMotion2D', 'SpatiallyVaryingBrownianMotion2D']
 
 
 def BrownianMotion2D(particle, fieldset, time, dt):
@@ -14,3 +14,26 @@ def BrownianMotion2D(particle, fieldset, time, dt):
     particle.lat += random.uniform(-1., 1.) * math.sqrt(2*math.fabs(dt)*kh_meridional/r)
     kh_zonal = fieldset.Kh_zonal[time, particle.lon, particle.lat, particle.depth]
     particle.lon += random.uniform(-1., 1.) * math.sqrt(2*math.fabs(dt)*kh_zonal/r)
+
+
+def SpatiallyVaryingBrownianMotion2D(particle, fieldset, time, dt):
+    # Diffusion equations for particles in non-uniform diffusivity fields
+    # from Ross & Sharples (2004, doi:10.4319/lom.2004.2.289)
+    # and Spagnol et al. (2002, doi:10.3354/meps235299)
+
+    # regular Brownian motion step
+    r = 1/3.
+    kh_meridional = fieldset.Kh_meridional[time, particle.lon, particle.lat, particle.depth]
+    Ry = random.uniform(-1., 1.) * math.sqrt(2*math.fabs(dt)*kh_meridional/r)
+    kh_zonal = fieldset.Kh_zonal[time, particle.lon, particle.lat, particle.depth]
+    Rx = random.uniform(-1., 1.) * math.sqrt(2*math.fabs(dt)*kh_zonal/r)
+
+    # Deterministic 'boost' out of areas of low diffusivity
+    dKdx = fieldset.dKh_zonal_dx[time, particle.lon, particle.lat, particle.depth]
+    dKdy = fieldset.dKh_meridional_dy[time, particle.lon, particle.lat, particle.depth]
+    CorrectionX = dKdx * math.fabs(dt)
+    CorrectionY = dKdy * math.fabs(dt)
+
+    # diffuse particle as sum of Brownian motion and deterministic 'boost'
+    particle.lon += Rx + CorrectionX
+    particle.lat += Ry + CorrectionY


Test patch:
diff --git a/tests/test_diffusion.py b/tests/test_diffusion.py
index a0071f69..a0fb7c53 100644
--- a/tests/test_diffusion.py
+++ b/tests/test_diffusion.py
@@ -1,9 +1,11 @@
 from parcels import (FieldSet, Field, RectilinearZGrid, ParticleSet, BrownianMotion2D,
-                     JITParticle, ScipyParticle)
+                     SpatiallyVaryingBrownianMotion2D, JITParticle, ScipyParticle,
+                     Geographic, GeographicPolar)
 from parcels import rng as random
 from datetime import timedelta as delta
 import numpy as np
 import pytest
+from scipy import stats
 
 ptype = {'scipy': ScipyParticle, 'jit': JITParticle}
 
@@ -53,6 +55,47 @@ def test_fieldKh_Brownian(mesh, mode, xdim=200, ydim=100, kh_zonal=100, kh_merid
     assert np.allclose(np.mean(lats), 0, atol=tol)
 
 
+@pytest.mark.parametrize('mesh', ['spherical', 'flat'])
+@pytest.mark.parametrize('mode', ['scipy', 'jit'])
+def test_fieldKh_SpatiallyVaryingBrownianMotion(mesh, mode, xdim=200, ydim=100):
+    """Test SpatiallyVaryingDiffusion on a non-uniform diffusivity field
+    with a linear gradient in one direction"""
+    mesh_conversion = 1/1852./60 if mesh is 'spherical' else 1
+    fieldset = zeros_fieldset(mesh=mesh, xdim=xdim, ydim=ydim, mesh_conversion=mesh_conversion)
+
+    Kh = np.zeros((ydim, xdim), dtype=np.float32)
+    for x in range(xdim):
+        Kh[:, x] = np.tanh(fieldset.U.lon[x]/fieldset.U.lon[-1]*10.)*xdim/2.+xdim/2. + 100.
+
+    grid = RectilinearZGrid(lon=fieldset.U.lon, lat=fieldset.U.lat, mesh=mesh)
+    fieldset.add_field(Field('Kh_zonal', Kh, grid=grid))
+    fieldset.add_field(Field('Kh_meridional', Kh, grid=grid))
+
+    dKh_zonal_dx, _ = fieldset.Kh_zonal.gradient()
+    _, dKh_meridional_dy = fieldset.Kh_meridional.gradient()
+    fieldset.add_field(Field('dKh_zonal_dx', dKh_zonal_dx, grid=grid))
+    fieldset.add_field(Field('dKh_meridional_dy', dKh_meridional_dy, grid=grid))
+    if mesh is 'spherical':
+        fieldset.dKh_zonal_dx.units = GeographicPolar()
+        fieldset.dKh_meridional_dy.units = Geographic()
+
+    npart = 100
+    runtime = delta(days=1)
+
+    random.seed(1234)
+    pset = ParticleSet(fieldset=fieldset, pclass=ptype[mode],
+                       lon=np.zeros(npart), lat=np.zeros(npart))
+    pset.execute(pset.Kernel(SpatiallyVaryingBrownianMotion2D),
+                 runtime=runtime, dt=delta(hours=1))
+
+    lats = np.array([p.lat for p in pset])
+    lons = np.array([p.lon for p in pset])
+    tol = 2000*mesh_conversion  # effectively 2000 m errors (because of low numbers of particles)
+    assert np.allclose(np.mean(lons), 0, atol=tol)
+    assert np.allclose(np.mean(lats), 0, atol=tol)
+    assert(stats.skew(lons) > stats.skew(lats))
+
+
 @pytest.mark.parametrize('mode', ['scipy', 'jit'])
 @pytest.mark.parametrize('lambd', [1, 5])
 def test_randomexponential(mode, lambd, npart=1000):