Extracting sub-domains with indices, but how to handle the longitude periodicity ?

slawchune commented 6 years ago

Hello Erik and Phillipe,

I run Parcels with 1/12° global currents and my velocity data are so big that I feel convenient to extract a subdomain of velocity fields to work with. But some of the drifters that I model cross the [-180°,180°] boundary, that are in fact at the edges of my grid (centred at 0° of longitude).

I tried to set negative indexes to my extractions (for instance : indices_extract= {'depth': [0], 'lat': range(863, 975), 'lon': range(-67, 67)}), which seems to work when plotting the initial state of the particles : parcels_longitude_periodicity_0

But finally after some hours of drift, the particle doesn't want to cross the western boundary to the reach the eastern one.

parcels_longitude_periodicity_50

Instead, the particle falls out of the domain and is removed thanks the DeleteParticle kernel given as example in the tutorial.

Am I doing things wrong ? Sorry if the question was already post, I didn't find anything about this.

Best regards, Stéphane

slawchune commented 6 years ago

Momentarily closing the issue because the answer might be there : https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_periodic_boundaries.ipynb

slawchune commented 6 years ago

It actually works with a periodic halo. I only add the halo if the -180/180° longitude axis is inside my extraction region :

    # add a zonal periodic halo at the (180/-180) longitude axis. 
    # We only perform this if the longitude axis (180/-180) is contained inside the subgrid
    if ( fieldset.U.grid.lon[0] > 0  and fieldset.U.grid.lon[-1] < 0 ):
        fieldset.add_constant('halo_west', fieldset.U.grid.lon[fieldset.U.grid.lon.argmin()])
        fieldset.add_constant('halo_east', fieldset.U.grid.lon[fieldset.U.grid.lon.argmax()])
        fieldset.add_periodic_halo(zonal=True)

parcels_longitude_periodicity_72

Maybe a more elegant solution exists...

delandmeterp commented 6 years ago

Hi Stéphane,

At least you can replace fieldset.U.grid.lon[fieldset.U.grid.lon.argmin()] by np.min(fieldset.U.grid.lon)

But I'm not sure of what you want to do. You want for example to have a subset from 30E to 60E and it works, and in some run you want a subset from 170E to 160W and it does not work? Periodic boundaries cannot currently reconnect the second example. You can reconnect when the left and right longitude points do match (e.g 180W and 180E, but this is not the case). So is it what you need to do?

slawchune commented 6 years ago

Hi Philippe,

Thanks for the suggestion !

Well in the first example, ie a subset from 30E to 60E, I don't need the right/left periodicity (I assume that the particle must stay in my sub domain during all the advection. If not, it is declared out of boundaries and suppressed)

In the second example, due to the extraction indexes that I have chosen (-67, 67 for longitude), I have the current longitude grid :

In [239]: fieldset.U.grid.lon Out[239]: array([ 173.99998, 174.08331, 174.16664, 174.24998, 174.33331, 174.41666, 174.5 , 174.58334, 174.66666, 174.75 , 174.83334, 174.91666, 175. , 175.08334, 175.16666, 175.25 , 175.33334, 175.41666, 175.5 , 175.58334, 175.66666, 175.75 , 175.83334, 175.91666, 176. , 176.08334, 176.16666, 176.25 , 176.33334, 176.41666, 176.5 , 176.58334, 176.66666, 176.75 , 176.83334, 176.91666, 177. , 177.08334, 177.16666, 177.25 , 177.33334, 177.41666, 177.5 , 177.58334, 177.66666, 177.75 , 177.83334, 177.91666, 178. , 178.08334, 178.16666, 178.25 , 178.33334, 178.41666, 178.5 , 178.58334, 178.66666, 178.75 , 178.83334, 178.91666,

, 179.08334, 179.16666, 179.25 , 179.33334, 179.41666, 179.5 , 179.58334, 179.66666, 179.75 , 179.83334, 179.91666, -180. , -179.91667, -179.83333, -179.75 , -179.66667, -179.58333, -179.5 , -179.41667, -179.33333, -179.25 , -179.16667, -179.08333, -179. , -178.91667, -178.83333, -178.75 , -178.66667, -178.58333, -178.5 , -178.41667, -178.33333, -178.25 , -178.16667, -178.08333, -178. , -177.91667, -177.83333, -177.75 , -177.66667, -177.58333, -177.5 , -177.41667, -177.33333, -177.25 , -177.16667, -177.08333, -177. , -176.91667, -176.83333, -176.75 , -176.66667, -176.58333, -176.5 , -176.41667, -176.33333, -176.25 , -176.16667, -176.08333, -176. , -175.91667, -175.83333, -175.75 , -175.66667, -175.58333, -175.5 , -175.41667, -175.33333, -175.25 , -175.16667, -175.08333, -175. , -174.91667, -174.83333, -174.75 , -174.66667, -174.58333, -174.5 , -174.41666, -174.33331, -174.24997, -174.16666, -174.08331], dtype=float32)

1/ If I don't put any periodic boundary, my particle reaches the -180 of longitude point and is declared out of bounds by Parcels.

2/ If I put my periodic boundary as I explained in my last post, with [halo_east=180, halo_west=-179.91667], there is no error and the particles correctly travels from west to east as shown in the images (from post 1 to post 3). The weird thing is that the periodic boundary is now at the centre of the fieldset grid. And as you can see from my longitude grid, I don't have two points of connection like [...,-180,180,...].

Is this working for bad reasons ?

I think I just need a recommendation of how to extract a subdomain that contains the grid periodicity. Sorry if this is unclear.

Cheers, Stéphane

delandmeterp commented 6 years ago

The mesh you specify is not a periodic boundary actually. But yet it has the discontinuity.

We treat such mesh for the case of curvilinear grids (such as in NEMO) but you have a rectilinear grid, and we assume the vector is increasing. Why don't you simply re-set your lon vector? So you set fieldset.U.grid.lon = np.where(fieldset.U.grid.lon < 0, fieldset.U.grid.lon+360, fieldset.U.grid.lon) Such that it will be reset between 0 and 360.

slawchune commented 6 years ago

Do only the grids need to be shifted and fieldset.U.lon stays untouched ? I kept the initialization of my particle at -179.09873.

Applying this produces the same error than previously (out of bounds at 51th hour of drift). The plots only covers the right side of the domain ( -> 174,180 ).

delandmeterp commented 6 years ago

fieldset.U.lon should be modified as well for safety, but fieldset.U.grid.lon is the main vector.

Could you send an example code with a link to the necessary data ro run it (as small as possible) ?

slawchune commented 6 years ago

I Phillipe, you will find the main launch script that I use and one input file example here :

https://atlas.mercator-ocean.fr/s/5MxNKPnfJZny88A

As there is only one daily input files, you will need to modify manually the time window, sorry : /. My problem is the initialisation of one particle very close to the -180/180 meridian.

Thanks a lot for your expertise and advice. Stéphane

delandmeterp commented 6 years ago

Hi Stéphane,

I have reduced your example to a simpler one. See attached (it is a py file, but I added the txt extension to put it here) test_run.py.txt

So, this works easily with no indices. We have never tried it before, but I have checked that you can add non-contiguous indices (see in the file, with option useIndices=True So using the beginning and the end of the lon grid. (@erikvansebille, have you seen that?) BUT:

It is way slower with the indices. There is something super slow when we load the non-contiguous array. So it is to be avoided, at least for now
Why do you want to use indices in your case? It is not a super big mesh, and it fits easily in memory. By the way, do you know that we don't load all time steps at once, but only 3 by 3, at the moment which is needed, specifically to save memory. So you don't need to go for sub-indices for the application you do
Watch out, that with the non-contiguous sub-indices, you have a mesh like [-180, -179, -178, 178, 179], (plus the periodic halo if you add it. This means that a particle located at 0°E will be simply interpolated naturally between -178 and 178. Since you don't want it, you have to add a kernel that checks if -178 < p.lon < 178 and in this case to delete the particle.

delandmeterp commented 6 years ago

A few more comments:

I don't get how you could have a lon vector as the one you showed a few days ago. I guess you did something weird with your indices. As you see in my examples, lon is always a increasing vector, since it was the case with your original data.
You can see that I have ignored the depth indices, but I also ignored depth in the dimensions. This is the equivalent to only taking the [0] index.
Since I am not really sure what was happening in your case and what was your problem, could you tell me if it is now fixed?
If you need further help, please provide a simple example as the one I've sent you, when you reduce the python file to the strictly necessary lines to reproduce the bug, thanks!

slawchune commented 6 years ago

Hi Phillipe !

I have just come back at the office. Thank you for having worked this out ! I guess I tried to optimize unnecessarily and clumsily. I'll avoid the sub-extraction indices if it slows things then. I didn't know that only a fraction of the time window was loaded in memory, this reassures me for the large number of simulation that I plan to launch.

Many thanks, Stéphane

Ps : for the weird non-contiguous longitude vector that I had previously, all I did was to set up negative sub-indexes for the longitude extraction as indices_extract['lon']=[-67, 67]

delandmeterp commented 6 years ago

Good that it works! Yes, we can't use negative indices. To see it, try this piece of code:

import numpy as np
a = np.arange(0,10)
print a[-2:4]

and you'll see it does not work

I added this comment to indices documentation (see #483)

OceanParcels / Parcels

Extracting sub-domains with indices, but how to handle the longitude periodicity ? #468