New Advection Diffusion Algorithms (Bug in Milstein + template for 3D Advection)

[ElizaJayne11]

Hi

This issue is continuing from a discussion my student and I started in issue #849, which is related to the new Advection + Diffusion kernels implemented in Parcels version 2.2.

As noted earlier, we are wanting to run a simulation with 3D advection + 3D spatially-varying horizontal diffusion (i.e. Kh varies with lat/lon and depth). I am very familiar with both RK4 and Milstein methods, and since the current Parcels code only works for 2D, I am developing a custom kernel that extends what you did to also advect in the vertical, and to use the local Kh values at the particle's current depth. Below I list 3 points related to this effort. I am very happy to continue interacting with you to address these points, as well as aiding with future improvements to Parcels numerical methods.

1. I discovered an errror in the formula you used for Milstein. Specifically, in your terms involving dKdx and dKdy, your term in parentheses is wrong. It should be the random number ^2 + 1 and then all of that multiplied by delta t/2, which is not the same as what you have. The correct formula is show in the code below using your notation.

2. Below is my draft of the custom kernel for the 3D Advection with the Milstein method for 2D diffusion at each depth. I welcome you feedback about this draft code. I should note that I am a Python newbie (my expertise is MATLAB) and am struggling with figuring out Python's syntax and data-types in order to make things readable and efficient. This is NOT how I would write this in MATLAB.

3. I am confused about whether my units are correct. For example, it appears that in Parcels, velocities are already converted to degrees per second in fieldset, so we can use u,v,w directly when calculating displacements in lat lon in the custom kernel. In contrast, it appears that in the fieldset, kH remains in m^2/ second, so the units need to be adjusted in the custom kernel. This is not obvious, given that when I run your 2D Adv-Dif with the builtin Milstein (and your earlier built-in kernels BrownianMotion2D and SpatiallyVaryingBrownianMotion), it appears that Kh has been converted to degrees^2 / second. Can you please confirm this?

Best wishes
Wendy

---

from parcels import FieldSet, Field, ParticleSet, JITParticle, ParticleFile, ErrorCode
import numpy as np
from datetime import timedelta as delta
import math
from netCDF4 import Dataset
from parcels import rng as random

def DeleteParticle(particle, fieldset, time):
    particle.delete()

def AdvectionRK43D_DiffusionM12D(particle, fieldset, time):
    # Kernel for 3D advection-solved using fourth order Runge-Kutta (RK4)
    # and 3D Horizontal diffusion -- solved using the Milstein scheme at first order (M1). 
    # Note that this approach does not consider the interaction of advection and diffusion within the time step
    # Milstein requires computation of the local gradient in Kh in both directions
    # These gradients are estimated using a centered finite different approach of O(h^2)
    # Where h = dres (set in the main code). this should be at least an order of magnitude
    # less than the typical grid resolution, but not so small as to introduce catestrophic cancellation
    # Parcels 2.2 has used a uniform random number to determine the random kick
    # Mathematically, a normal distribution is more accurate for a given time step as noted in Parcels 2.2
    # Parcels 2.2 argue that uniform is more efficient, which is generally true when comparing individual calls
    # However, since the uniform only converges on the correct distribution after about 5 steps, 
    # it is not overall more efficient than a normal distribution. 
    # Parcel 2.2. also argue that the random increments remain bounded. 
    # It is not immediately apparent to me why "boundedness" is an advantage. 
    # Despite that,I have kept the code below using a uniform distribution
    # in part b/c horizontal diffusion is already very approximate (our diffusivities are not precise!)
    # There is also a long history of Ocean Modelers using a uniformly distributed random number (e.g. Visser, Ross & Sharples etc)

    # Current space-time location of particle
    lonp = particle.lon
    latp = particle.lat
    depthp = particle.depth
    dt = particle.dt

    # RK4 for Advection
    # Changed variable naming from Parcels 2.2 to 
    # (1) improve efficiency as we were instructed that it is faster to rename variables that are "." than to recall and 
    # (2) to make the code easier to read
    # Note: I have followed Parcels 2.2. and left the RK4 to be computing slopes within the time step instead of displacements
    # Both approaches are used in the literature 

    
    #start at current space-time point using those values to get slope
    (u1, v1, w1) = fieldset.UVW[time, depthp, latp, lonp]
    
    #Get estimate for slope at midpoint using previous estimate for slope
    lon1, lat1, dep1 = (lonp + u1*.5*dt, latp + v1*.5*dt, depthp + w1*.5*dt) 
    (u2, v2, w2) = fieldset.UVW[time + .5*dt, dep1, lat1, lon1]
    
    #Get improved estimate for slope at midpoint using previous estimate for slope
    lon2, lat2, dep2 = (lonp + u2*.5*dt, latp + v2 *.5*dt, depthp + w2*.5*dt)
    (u3, v3,w3) = fieldset.UVW[time + .5 * dt, dep2, lat2, lon2]
    
    #Get estimate for slope at endppoint using previous estimate for slope at midpoint
    lon3, lat3, dep3 = (lonp + u3*dt, latp + v3*dt, depthp+w3*dt)
    (u4, v4, w4) = fieldset.UVW[time + dt, dep3, lat3, lon3]
    
    #Calculate particle displacement due to local advection
    #This assumes that fieldset has already converted u,v,w to degrees
    Advect_lon = ((u1 + 2 * u2 + 2 * u3 + u4) / 6.) * dt 
    Advect_lat = ((v1 + 2 * v2 + 2 * v3 + v4) / 6.) * dt 
    Advect_dep = ((w1 + 2 * w2 + 2 * w3 + w4) / 6.) * dt 

    #Milstein for Horizontal Diffusion
    Khdifferencedist = fieldset.dres #in degrees
    #Note Kh is in m^2/s here, unlike built-in that has already converted to degrees
    #To save on repeated conversions, only done after compute displace ent

    #Sample random number (Wiener increment) with zero mean and std of sqrt(dt) from uniform distribution
    dWx = random.uniform(-1., 1.) * math.sqrt(math.fabs(dt) * 3)
    dWy = random.uniform(-1., 1.) * math.sqrt(math.fabs(dt) * 3)
    # dWx = random.normalvariate(0., 1.) * math.sqrt(math.fabs(dt))
    # dWy = random.normalvariate(0., 1.) * math.sqrt(math.fabs(dt))
  
    #Get estimate of random kick in x-direction based on local diffuvisity (neglects spatial variation)
    bx = math.sqrt(2* fieldset.Kh_zonal[time, depthp, latp, lonp])
    #Get estimate of randome kick in y-direction based on local diffusivity (neglects spatial variation)
    by = math.sqrt(2 * fieldset.Kh_meridional[time, depthp, latp, lonp])
    
    #Get estimate of zonal diffusivity gradient at current location using finite centered finite differenece
    #This derivative is used to correct basic random kick due to variable diffusivity
    Kxp1 = fieldset.Kh_zonal[time, depthp, latp, lonp + Khdifferencedist]
    Kxm1 = fieldset.Kh_zonal[time, depthp, latp, lonp - Khdifferencedist]
    dKdx = (Kxp1 - Kxm1) / (2 * Khdifferencedist)
 
    #Get estimate of meridional gradient at current location using fininte centered difference
    #This derivative is used to correct basic random kick due to variable diffusivity 
    Kyp1 = fieldset.Kh_meridional[time, depthp, latp + Khdifferencedist, lonp]
    Kym1 = fieldset.Kh_meridional[time, depthp, latp - Khdifferencedist, lonp]
    dKdy = (Kyp1 - Kym1) / (2 * Khdifferencedist) 
    
    #Calculate particle horizontal displacement due to local diffusiona and diffusivity gradient
    DiffH_lon = bx*dWx + 0.5 * dKdx * (dWx**2 + 1)*dt #Wonky units
    DiffH_lat = by*dWy + 0.5 * dKdy * (dWy**2 + 1)*dt #Wonky units
    to_lat = 1/1000./1.852/60.
    to_lon = to_lat/math.cos(particle.lat*math.pi/180)
    DiffH_lon=DiffH_lon*to_lon
    DiffH_lat=DiffH_lat*to_lat

    #Particle positions are updated only after evaluating all terms (i.e. Advection + Diffusion simulataneous)
    #Note this approach does not consider the interaction of adv and diff on estimated particle displacements within the time step)
    particle.lon += Advect_lon + DiffH_lon
    particle.lat += Advect_lat + DiffH_lat
    particle.depth += Advect_dep


#Read in physical fields

data_path = '../..//BNAM/PhysicalFields/'

ufiles = data_path + '3D_U_'+ 'September' +'_BNAM.nc'
vfiles = data_path + '3D_V_'+ 'September' +'_BNAM.nc'
wfiles = data_path + '3D_W_'+ 'September' +'_BNAM.nc'
mesh_mask = data_path + '3D_Mask_'+ 'September' +'_BNAM.nc'
        
filenames = {'U': {'lon': mesh_mask, 'lat': mesh_mask, 'depth': wfiles, 'data': ufiles},
                     'V': {'lon': mesh_mask, 'lat': mesh_mask, 'depth': wfiles, 'data': vfiles},
                     'W': {'lon': mesh_mask, 'lat': mesh_mask, 'depth': wfiles, 'data': wfiles}}
        
variables = {'U': 'U',
                     'V': 'V',
                     'W': 'W'}
        
dimensions = {'U': {'lon': 'glamf', 'lat': 'gphif', 'depth': 'deptht', 'time': 'time_counter'},
                      'V': {'lon': 'glamf', 'lat': 'gphif', 'depth': 'deptht', 'time': 'time_counter'},
                      'W': {'lon': 'glamf', 'lat': 'gphif', 'depth': 'deptht', 'time': 'time_counter'}}



fieldset = FieldSet.from_nemo(filenames, variables, dimensions, allow_time_extrapolation=True)


dataset = Dataset(data_path +"Mixing3D_BNAM.nc")
kh_zonal = dataset.variables['kh_zonal'][:]
kh_meridional = dataset.variables['kh_meridional'][:]
fieldset.add_field(Field('Kh_zonal', kh_zonal, grid=fieldset.U.grid))
fieldset.add_field(Field('Kh_meridional', kh_meridional, grid=fieldset.U.grid))
fieldset.add_constant('dres', 0.001) #Changed to avoid potential for catestrophic cancellation. Our grid resolution is ??
# Set coordinates of partiles
a = np.loadtxt('../../SCALLOPS/SCALLOP DATA/ICs/'+ 'ICs1.txt')
#Lon and lat coordinates respectively
lonp = a[:,0] 
latp = a[:,1]
depthp = np.ones(len(lonp))


pset = ParticleSet.from_list(fieldset, JITParticle, lon=lonp, lat=latp, depth=depthp)
kernels = pset.Kernel(AdvectionRK43D_DiffusionM12D) 
pset.execute(kernels, runtime=delta(days=45), dt=delta(hours=1),
             output_file = ParticleFile("Particles_Testing", pset, outputdt=delta(days=1)),
             recovery={ErrorCode.ErrorOutOfBounds: DeleteParticle})

[daanreijnders]

Dear Wendy,

Thanks for working with and contributing to Parcels, and for going through this kernel so extensively. I'll try to answer your questions/remarks on a point-by-point basis.

1. Note that dW is the Wiener increment of which the standard deviation is already scaled with the square root of dt. The kernel implemented in Parcels 2.2 is consistent with Eq (12) from Gräwe (2011) and with a rearranged form of (3.1) from Kloeden & Platen (1999, see p 345). The term of 0.5 * dKdx * (particle.dt) represents the deterministic drift caused by the gradient in diffusivity. The term 0.5 * dldx * (dWx**2) + bx * dWx corresponds to the stochastic part. The latter term + bx * dWx seems to be missing in your implementation.
   
   Could you share the literature source for your implementation of the Milstein kernel? I'd like to find out what causes the discrepancy between our implementations.

2. Replying to your comments:

   • The rationale for using a uniform over a normal distribution is that a normal distribution is computationally more costly, but perhaps this is tinkering in the margins of efficiency, since other Parcels routines are much more costly. Still, simulations are often carried out only taking into account ensemble distributions over more than 5 time steps (and usually with a smaller dt). The 'boundedness' argument is related to the probability of overshooting boundaries (land or pycnocline-like structures) when the random step lies in the tails of the normal distribution. When using a uniform distribution, this is easier to account for (since the random step is bounded sqrt(particle.dt * 3)). Feel free to change the Wiener increment to use the random numbers that you find fitting; that's why it's kept on a separate line.
   • I believe accessing variables with . is only more expensive when it implies making field interpolations, e.g. fieldset.UV[..], but not when it concerns a variable like fieldset.dres. @erikvansebille or @CKehl, could you elaborate on this? Otherwise it's a good change to implement.
   • You're making a good point about the interaction between the advection and diffusion. I've recently seen implementations where the diffusivity is sampled at the midpoint of the advective tilmestep. What do you think is a best practice?
     
     I think that apart from our different implementation of the Milstein kernel, your steps are correct. There is a problem with the units in your implementation though (see next point).

3. Parcels automatically takes cares of the unit conversion for u, v, Kh_zonal and Kh_meridional; the variables in the unit converter_map here. This is implemented for the diffusivities here. I did a dimensional analysis of the particle displacement in the Milstein kernel, and if Kh is in deg^2/s, the displacement is in deg (note also that [dKdx]=deg/s, [dW]=sqrt(s), [dW^2*dKdx]=deg, [bx] = sqrt([Kh])=deg/sqrt(s)). With your implementation, the units are inconsistent (deg/s * (s+1)*s).

Thanks again. Parcels certainly benefits from another set of eyes sifting through the code.

Best,
Daan

[erikvansebille]

I believe accessing variables with . is only more expensive when it implies making field interpolations, e.g. fieldset.UV[..], but not when it concerns a variable like fieldset.dres. @erikvansebille or @CKehl, could you elaborate on this? Otherwise it's a good change to implement.

The conversion is not done in fieldset.dres, but instead when fieldset.Kh_zonal[...] or fieldset.Kh_meridional[...] are called. So there should indeed be no measurable overhead with using fieldset.dres multiple times

[ElizaJayne11]

Hello again

Thank you for your prompt and thorough reply to my message. More importantly, thank you for your continued development and maintenance of Ocean Parcels. Your contributions enable the rest of us to get on with the science because we have such a useful tool.

Below I follow up on your points above.

1. I did have the bx * dWx term (it just came before the stochastic part, see above). However where I went wrong was in recognizing that your dW term did not correspond to the stochastic term I was using in my MATLAB Milstein code, which just the random number and therefore different by the factor sqrt(dt). When I converted my code to Python, I used your dW in lieu of my term, which of course made the equations look and function differently. Now that I make that distinction, I can see that your and my implementations of Milstein are in fact identical. I sincerely apologize for any concern or inconvenience I caused. Thanks again for your checking and corrections.

2a. I do appreciate that simulations are run over a series of time steps, so the uniform will have a better chance of "convergence", depending of course on how variable your Kh is over the displacement occurring in a time step. Still, from a technical viewpoint, each individual step is more accurate with Gauss. So you can get away with a bigger time step (similar to Euler vs. RK4, in that RK4 is more costly per time step, but you need fewer time steps for the same accuracy). Thus my own personal preference is for Gauss. I appreciate that this may lead to more "out of bounds" scenarios. Are those boundary checks done on every particle at every time step? (i.e. would having that affect the efficiency?)

2b. As you likely know, the Visser, 97 algorithm had the diffusivity evaluated at the midpoint of the deterministic part of the time step, which is the only difference between that algorithm and the Euler method. In my (unpublished) investigations that compared different random walk schemes, Visser was virtually indistinguishable from Euler, and both were no where near as good as Milstein. Thus, I don't think that such a mid-point interpolation would increase the accuracy. Moreover, I think that such an interpolation might reduce the efficiency, since you now have to get Kh at another spatial point. My own view, expressed in the earlier issue ticket, is that accurate soIution of the diffusivity part is not critical (at least in my work). However, when the gradients in your currents are large enough to warrant an RK4 approach, you may well want to account for the additional displacements due to diffusion in the RK4 intra-time-step interpolations. I have never tried this, but think it would be worth checking into.

3. I had changed the units because my simulations appeared to be having diffusive kicks that were wrong, and I incorrectly thought it was due to the units. Those errors had to do with my incorrect dW term. Now that I have fixed my Milstein implementation, I was able to use the Kh's as they come out of fieldset, which is clearly in degrees^2/s as you note.

Thanks again
Wendy

[daanreijnders]

Hi Wendy,

I'm glad we've got to the bottom of this, and to know we've both been using the right numerical scheme. If you're interested in making your kernel available to the rest of the Parcels userbase, you're invited to open a pull request in the oceanparcels/parcels_contributions repository. We are trying to only include a small amount of kernels by default, but we've created this space for users to share variations of kernels (or entirely new ones) with each other under the MIT License.

To answer your remaining question about the out-of-bounds check: this is done in each JIT-loop, the length of which is determined by the frequency of writing output data. If a particle goes out-of-bounds within that loop, it will not be evaluated anymore during the rest of the loop. Then, at the end of the loop, the recovery kernel is applied to the out-of-bounds particles.

And with regards to your answer about where to evaluate the diffusivity: I agree that right now implementing a sampling at the mid-point is not worth it. In case an RK4 advection-diffusion algorithm is properly implemented (instead of the current advectionRK4diffusionM1), diffusivities should indeed be sampled at intermediate points, but RK4 is future work. Good to have received your take on it.

Best,
Daan

[SaAlrabeei]

I believe accessing variables with . is only more expensive when it implies making field interpolations, e.g. fieldset.UV[..], but not when it concerns a variable like fieldset.dres. @erikvansebille or @CKehl, could you elaborate on this? Otherwise, it's a good change to implement.

The conversion is not done in fieldset.dres, but instead when fieldset.Kh_zonal[...] or fieldset.Kh_meridional[...] are called. So there should indeed be no measurable overhead with using fieldset.dres multiple times

Hello all , thank you all for this nice discussion, I want to ask about two more points.
When having 3D model ( U,V and W) velocities,
What is the conversion constant of the W from m/s to m/degree! Furthermore,
I am trying to add a vertical diffusivity ( eddy diffusivity which is varying spatially), I don't know also what would be the conversion.

On the other hand, the vertical grid ( given by the depth levels) is not regular, which makes the "dres" constant is not unique although the particle vertical diffusion is small within each time step.

My question is that, do you recommend not to use vertical diffusion in PARCELS so far!

[daanreijnders]

In 3D data, vertical velocities are kept in m/s. Only horizontal velocities are converted, only in the case the grid you use is spherical. So you won't have to worry about a conversion of w.

Good point about dres. In your kernel, you can implement a different dres for the horizontal and vertical directions. There is no best value of dres. The larger its value, the more it smoothens derivatives in the data, and it's up to you to determine how desirable this is. But as a rule of thumb, so far I've used 1/10th of the typical grid resolution.

Whether you should or shouldn't use vertical diffusion is a question inherent to your research question. It may be very insightful for you. So far, we haven't implemented it in a separate kernel, since our own applications have chiefly focused on horizontal diffusion. But vertical diffusion is on the agenda.

[SaAlrabeei]

In 3D data, vertical velocities are kept in m/s. Only horizontal velocities are converted, only in the case the grid you use is spherical. So you won't have to worry about a conversion of w.

Good point about dres. In your kernel, you can implement a different dres for the horizontal and vertical directions. There is no best value of dres. The larger its value, the more it smoothens derivatives in the data, and it's up to you to determine how desirable this is. But as a rule of thumb, so far I've used 1/10th of the typical grid resolution.

Whether you should or shouldn't use vertical diffusion is a question inherent to your research question. It may be very insightful for you. So far, we haven't implemented it in a separate kernel, since our own applications have chiefly focused on horizontal diffusion. But vertical diffusion is on the agenda.

Fair enough,
Thank you for the quick feedback.

[ElizaJayne11]

Hi Daan I am very happy to share my 3D RK4 Advection + Horizontal Milstein Diffusion kernel with the other Parcels users. However, I have never used GitHub and as such I find the pull request explanations totally confusing and full of jargon I don't understand. For example, when I click on you link above, and then on "Pull Requests/New Pull Request" the "Create pull request" link is greyed out. It seems to want me to compare branches and forks in order to upload my code, but I don't know what those things are or how they relate to my code. The process seems very complicated when all I want to do is upload a variant on an existing kernel. Can you please walk me through the steps to make my contribution? Wendy

…

On Fri, 12 Jun 2020 at 05:26, Daan Reijnders ***@***.***> wrote: Hi Wendy, I'm glad we've got to the bottom of this, and to know we've both been using the right numerical scheme. If you're interested in making your kernel available to the rest of the Parcels userbase, you're invited to open a pull request in the oceanparcels/parcels_contributions <https://github.com/OceanParcels/parcels_contributions> repository. We are trying to only include a small amount of kernels by default, but we've created this space for users to share variations of kernels (or entirely new ones) with each other under the MIT License. To answer your remaining question about the out-of-bounds check: this is done in each *JIT-loop*, the length of which is determined by the frequency of writing output data. If a particle goes out-of-bounds within that loop, it will not be evaluated anymore during the rest of the loop. Then, at the end of the loop, the recovery kernel is applied to the out-of-bounds particles. And with regards to your answer about where to evaluate the diffusivity: I agree that right now implementing a sampling at the mid-point is not worth it. In case an RK4 advection-diffusion algorithm is properly implemented (instead of the current advectionRK4diffusionM1), diffusivities should indeed be sampled at intermediate points, but this also requires a cubic interpolation method, which is not yet implemented. So this is future work. Good to have received your take on it. Best, Daan — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#858 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AN7BD2NASPXZMSUVE2AZPC3RWHRE7ANCNFSM4N2XG5YA> .

[daanreijnders]

No problem. Since Github has a bit of a learning curve for new users, you can simply send me your code via email and I can publish it for you (with attribution). Otherwise I think this guide could serve as a walkthrough (without the first step, since you can use the existing repository we have). I can also recommend checking out these resources.

[ElizaJayne11]

Hi Daan I greatly appreciate the links to the GitHub help pages. I will investigate them when I have more time. For now I am going to take you up on your offer to post it for me. Please see attached. It would also be great, if you have the opportunity, to give it a quick proofread before posting. Best wishes Wendy

…

On Mon, 29 Jun 2020 at 09:28, Daan Reijnders ***@***.***> wrote: No problem. Since Github has a bit of a learning curve for new users, you can simply send me your code via email and I can publish it for you (with attribution). Otherwise I think this guide <https://guides.github.com/activities/hello-world/> could serve as a walkthrough (without the first step, since you can use the existing repository we have). I can also recommend checking out these resources <https://try.github.io>. — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#858 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AN7BD2IYFSTS7W7C5WWH3NTRZCCFXANCNFSM4N2XG5YA> .

[mollyjames2]

Hi, did the updated code for this 3D RK4 Advection + Horizontal Milstein Diffusion kernel ever get shared? I've had a look in parcels contributions and cannot see it but i would be very keen to use it.

Many thanks

[ElizaJayne11]

Hi Molly I thought the code had been shared in 2022. I am inexperienced with GitHub and so had Daan Rejiinders submit it for me. If you can't access it, I'm happy to send you my own code, which you can modify to suit your needs. Wendy

…

On Tue, 16 Jan 2024 at 07:29, Mollyjames ***@***.***> wrote: Hi, did the updated code for this 3D RK4 Advection + Horizontal Milstein Diffusion kernel ever get shared? I've had a look in parcels contributions and cannot see it but i would be very keen to use it. Many thanks — Reply to this email directly, view it on GitHub <#858 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AN7BD2PGIOLPD6SUIK3JZ6TYOZQBXAVCNFSM4N2XG5YKU5DIOJSWCZC7NNSXTN2JONZXKZKDN5WW2ZLOOQ5TCOBZGM2TKOBVGI2Q> . You are receiving this because you authored the thread.Message ID: ***@***.***>

[mollyjames2]

Hi Wendy, that would be fantastic if you're happy to do that? my email is moja@pml.ac.uk - Thank you very much