-
Notifications
You must be signed in to change notification settings - Fork 120
/
advectiondiffusion.py
112 lines (82 loc) · 5.23 KB
/
advectiondiffusion.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
"""Collection of pre-built advection-diffusion kernels
See `this tutorial <https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_diffusion.ipynb>`_ for a detailed explanation"""
import math
import parcels.rng as ParcelsRandom
__all__ = ['DiffusionUniformKh', 'AdvectionDiffusionM1', 'AdvectionDiffusionEM', ]
def AdvectionDiffusionM1(particle, fieldset, time):
"""Kernel for 2D advection-diffusion, solved using the Milstein scheme
at first order (M1).
Assumes that fieldset has fields `Kh_zonal` and `Kh_meridional`
and variable `fieldset.dres`, setting the resolution for the central
difference gradient approximation. This should be (of the order of) the
local gridsize.
This Milstein scheme is of strong and weak order 1, which is higher than the
Euler-Maruyama scheme. It experiences less spurious diffusivity by
including extra correction terms that are computationally cheap.
The Wiener increment `dW` is normally distributed with zero
mean and a standard deviation of sqrt(dt).
"""
# Wiener increment with zero mean and std of sqrt(dt)
dWx = ParcelsRandom.normalvariate(0, math.sqrt(math.fabs(particle.dt)))
dWy = ParcelsRandom.normalvariate(0, math.sqrt(math.fabs(particle.dt)))
Kxp1 = fieldset.Kh_zonal[time, particle.depth, particle.lat, particle.lon + fieldset.dres]
Kxm1 = fieldset.Kh_zonal[time, particle.depth, particle.lat, particle.lon - fieldset.dres]
dKdx = (Kxp1 - Kxm1) / (2 * fieldset.dres)
u = fieldset.U[time, particle.depth, particle.lat, particle.lon]
bx = math.sqrt(2 * fieldset.Kh_zonal[time, particle.depth, particle.lat, particle.lon])
Kyp1 = fieldset.Kh_meridional[time, particle.depth, particle.lat + fieldset.dres, particle.lon]
Kym1 = fieldset.Kh_meridional[time, particle.depth, particle.lat - fieldset.dres, particle.lon]
dKdy = (Kyp1 - Kym1) / (2 * fieldset.dres)
v = fieldset.V[time, particle.depth, particle.lat, particle.lon]
by = math.sqrt(2 * fieldset.Kh_meridional[time, particle.depth, particle.lat, particle.lon])
# Particle positions are updated only after evaluating all terms.
particle.lon += u * particle.dt + 0.5 * dKdx * (dWx**2 + particle.dt) + bx * dWx
particle.lat += v * particle.dt + 0.5 * dKdy * (dWy**2 + particle.dt) + by * dWy
def AdvectionDiffusionEM(particle, fieldset, time):
"""Kernel for 2D advection-diffusion, solved using the Euler-Maruyama
scheme (EM).
Assumes that fieldset has fields `Kh_zonal` and `Kh_meridional`
and variable `fieldset.dres`, setting the resolution for the central
difference gradient approximation. This should be (of the order of) the
local gridsize.
The Euler-Maruyama scheme is of strong order 0.5 and weak order 1.
The Wiener increment `dW` is normally distributed with zero
mean and a standard deviation of sqrt(dt).
"""
# Wiener increment with zero mean and std of sqrt(dt)
dWx = ParcelsRandom.normalvariate(0, math.sqrt(math.fabs(particle.dt)))
dWy = ParcelsRandom.normalvariate(0, math.sqrt(math.fabs(particle.dt)))
Kxp1 = fieldset.Kh_zonal[time, particle.depth, particle.lat, particle.lon + fieldset.dres]
Kxm1 = fieldset.Kh_zonal[time, particle.depth, particle.lat, particle.lon - fieldset.dres]
dKdx = (Kxp1 - Kxm1) / (2 * fieldset.dres)
ax = fieldset.U[time, particle.depth, particle.lat, particle.lon] + dKdx
bx = math.sqrt(2 * fieldset.Kh_zonal[time, particle.depth, particle.lat, particle.lon])
Kyp1 = fieldset.Kh_meridional[time, particle.depth, particle.lat + fieldset.dres, particle.lon]
Kym1 = fieldset.Kh_meridional[time, particle.depth, particle.lat - fieldset.dres, particle.lon]
dKdy = (Kyp1 - Kym1) / (2 * fieldset.dres)
ay = fieldset.V[time, particle.depth, particle.lat, particle.lon] + dKdy
by = math.sqrt(2 * fieldset.Kh_meridional[time, particle.depth, particle.lat, particle.lon])
# Particle positions are updated only after evaluating all terms.
particle.lon += ax * particle.dt + bx * dWx
particle.lat += ay * particle.dt + by * dWy
def DiffusionUniformKh(particle, fieldset, time):
"""Kernel for simple 2D diffusion where diffusivity (Kh) is assumed uniform.
Assumes that fieldset has constant fields `Kh_zonal` and `Kh_meridional`.
These can be added via e.g.
fieldset.add_constant_field("Kh_zonal", kh_zonal, mesh=mesh)
fieldset.add_constant_field("Kh_meridional", kh_meridional, mesh=mesh)
where mesh is either 'flat' or 'spherical'
This kernel assumes diffusivity gradients are zero and is therefore more efficient.
Since the perturbation due to diffusion is in this case isotropic independent, this
kernel contains no advection and can be used in combination with a seperate
advection kernel.
The Wiener increment `dW` is normally distributed with zero
mean and a standard deviation of sqrt(dt).
"""
# Wiener increment with zero mean and std of sqrt(dt)
dWx = ParcelsRandom.normalvariate(0, math.sqrt(math.fabs(particle.dt)))
dWy = ParcelsRandom.normalvariate(0, math.sqrt(math.fabs(particle.dt)))
bx = math.sqrt(2 * fieldset.Kh_zonal[particle])
by = math.sqrt(2 * fieldset.Kh_meridional[particle])
particle.lon += bx * dWx
particle.lat += by * dWy