OceanParcels / Parcels

Main code for Parcels (Probably A Really Computationally Efficient Lagrangian Simulator)
https://www.oceanparcels.org
MIT License
291 stars 133 forks source link

ThroughSurfaceError using AdvectionRK4_3D #1210

Closed gacuervol closed 2 years ago

gacuervol commented 2 years ago

Hello,

First I wanted to say thank you for this wonderful Python module. I am using the following data:

!python -m pip install motuclient==1.8.4 --no-cache-dir  
!python -m motuclient --motu https://my.cmems-du.eu/motu-web/Motu --service-id MULTIOBS_GLO_PHY_W_3D_REP_015_007-TDS --product-id dataset-omega-3d-rep-weekly --longitude-min 30 --longitude-max 46 --latitude-min -32 --latitude-max -24 --date-min "2013-03-01" --date-max "2013-06-30" --depth-min 2.5 --depth-max 1482.5 --user 'YOUR_USR' --pwd 'YOUR_USR'
# Open the product with motuclient
ds_3d = xr.open_dataset('/YOUR/PATH/data.nc')

I am currently trying to simulate the 3D advection of 6 particles located in the axis of a cyclonic eddy at different depths. Using the following code:

u_name = ds_3d['uo'].attrs['long_name']
v_name = ds_3d['vo'].attrs['long_name']
w_name = ds_3d['wo'].attrs['long_name']
filenames = {'U': 'madagascar_3d.nc',
             'V': 'madagascar_3d.nc',
             'W': 'madagascar_3d.nc'}
variables = {'U': 'uo', 'V': 'vo', 'W': 'wo'}
dimensions = {'U': {'lat': 'latitude', 'lon': 'longitude', 'depth':'depth', 'time': 'time'},
              'V': {'lat': 'latitude', 'lon': 'longitude', 'depth':'depth', 'time': 'time'},
              'W': {'lat': 'latitude', 'lon': 'longitude', 'depth':'depth', 'time': 'time'}}
fieldset = FieldSet.from_netcdf(filenames, variables, dimensions)
# Define a new particle class
class DistSampleParticle(JITParticle):
    # Variable que mide la distancia de viaje
    distance = Variable('distance', initial=0., dtype=np.float32)  # the distance travelled
    prev_lon = Variable('prev_lon', dtype=np.float32, to_write=False,
                        initial=attrgetter('lon'))  # the previous longitude
    prev_lat = Variable('prev_lat', dtype=np.float32, to_write=False,
                        initial=attrgetter('lat'))  # the previous latitude.
# 2. Defining the particles type and initial conditions in a ParticleSet object
first_day= 2 # Indice del dia que se suelta la particula
plon, plat = 41.375, -26.375
depths_def = [-2.50, -100.25, -306.50, -506.50, -600.50, -1024.25, -1482.50]
pset = ParticleSet(fieldset=fieldset,   # the fields on which the particles are advected
                   pclass=DistSampleParticle,  # the type of particles (JITParticle or ScipyParticle)
                   lon=[plon]*len(depths_def),              # release longitudes 
                   lat=[plat]*len(depths_def),
                   depth=depths_def,
                   time= ds_3d.time.values[first_day])             # release latitudes
print(f'first particles: \n {pset}')   

# Extra kernel
def TotalDistance(particle, fieldset, time):
    # Calculate the distance in latitudinal direction (using 1.11e2 kilometer per degree latitude)
    lat_dist = (particle.lat - particle.prev_lat) * 1.11e2
    # Calculate the distance in longitudinal direction, using cosine(latitude) - spherical earth
    lon_dist = (particle.lon - particle.prev_lon) * 1.11e2 * math.cos(particle.lat * math.pi / 180)
    # Calculate the total Euclidean distance travelled by the particle
    particle.distance += math.sqrt(math.pow(lon_dist, 2) + math.pow(lat_dist, 2))

    particle.prev_lon = particle.lon  # Set the stored values for next iteration.
    particle.prev_lat = particle.lat
# Casting the TotalDistance function to a kernel.
k_dist = pset.Kernel(TotalDistance) 

# 3. Executing an advection kernel on the given fieldset
dias = 67
time_step_hour= 24

output_file = pset.ParticleFile(name="TrackParticles_3d.nc", 
                                outputdt=timedelta(days=1)) # the file name and the time step of the outputs
pset.execute(AdvectionRK4_3D + k_dist,                 # the kernel (which defines how particles move)
             runtime=timedelta(days=dias),              # the total length of the run
             dt=timedelta(days=1),                       # the timestep of the kernel
             output_file=output_file, 
             )
print(f'last particle: {pset}')
# 4. Exporting the simulation output to a netcdf file
output_file.export()

However, when I run the simulation with AdvectionRK4_3D, I get the following error:

first particles: 
 P[70](lon=41.375000, lat=-26.375000, depth=-2.500000, distance=0.000000, time=1209600.000000)
P[71](lon=41.375000, lat=-26.375000, depth=-100.250000, distance=0.000000, time=1209600.000000)
P[72](lon=41.375000, lat=-26.375000, depth=-306.500000, distance=0.000000, time=1209600.000000)
P[73](lon=41.375000, lat=-26.375000, depth=-506.500000, distance=0.000000, time=1209600.000000)
P[74](lon=41.375000, lat=-26.375000, depth=-600.500000, distance=0.000000, time=1209600.000000)
P[75](lon=41.375000, lat=-26.375000, depth=-1024.250000, distance=0.000000, time=1209600.000000)
P[76](lon=41.375000, lat=-26.375000, depth=-1482.500000, distance=0.000000, time=1209600.000000)
Exception ignored in: <function ParticleFileSOA.__del__ at 0x7fc7ddf65290>
Traceback (most recent call last):
  File "/usr/local/lib/python3.7/dist-packages/parcels/particlefile/particlefilesoa.py", line 40, in __del__
  File "/usr/local/lib/python3.7/dist-packages/parcels/particlefile/baseparticlefile.py", line 235, in __del__
  File "/usr/local/lib/python3.7/dist-packages/parcels/particlefile/baseparticlefile.py", line 240, in close
  File "/usr/local/lib/python3.7/dist-packages/parcels/particlefile/particlefilesoa.py", line 125, in export
  File "/usr/local/lib/python3.7/dist-packages/numpy/lib/npyio.py", line 417, in load
FileNotFoundError: [Errno 2] No such file or directory: 'out-BGMHACAX/0/pset_info.npy'
INFO: Compiled ArrayDistSampleParticleAdvectionRK4_3DTotalDistance ==> /tmp/parcels-0/lib01973388a1f5be34894ec40fc35244f7_0.so
---------------------------------------------------------------------------
ThroughSurfaceError                       Traceback (most recent call last)
[<ipython-input-49-d9a868f2affc>](https://localhost:8080/#) in <module>()
     35              dt=timedelta(days=1),                       # the timestep of the kernel
     36              output_file=output_file,
---> 37              recovery={#ErrorCode.ErrorOutOfBounds: DeleteParticle,
     38                        #ErrorCode.ErrorThroughSurface: SubmergeParticle
     39                        }

2 frames
[/usr/local/lib/python3.7/dist-packages/parcels/tools/statuscodes.py](https://localhost:8080/#) in recovery_kernel_through_surface(particle, fieldset, time)
    205         # TODO: JIT does not yet provide the context that created
    206         # the exception. We need to pass that info back from C.
--> 207         raise ThroughSurfaceError(particle, fieldset)
    208     else:
    209         error = particle.exception

ThroughSurfaceError: 0
Particle P[70](lon=41.375000, lat=-26.375000, depth=-2.500000, distance=0.000000, time=1209600.000000)
Time: 2013-03-20T00:00:00.000000000,    timestep dt: 86400.000000
Through-surface sampling by particle at (41.375000, -26.375000, -2.500000)

I must say that no error occurs with the AdvectionRK4 advection kernel.

I really appreciate any help.

axnsantana commented 2 years ago

Hi,

As the message error indicates, this is happening because the particle (at location 41.375000, -26.375000) tried to get through the surface, thus out of water.

You can fix that for example in two ways:

  1. Adding a recovery action for the specific error ThroughSurfaceError, like that:
def DeleteParticle(particle, fieldset, time):
    particle.delete()

pset.execute(AdvectionRK4_3D + k_dist,                 # the kernel (which defines how particles move)
             runtime=timedelta(days=dias),              # the total length of the run
             dt=timedelta(days=1),                       # the timestep of the kernel
             recovery={ErrorCode.ErrorThroughSurface: DeleteParticle},
             output_file=output_file, 
             )

In that case, you will need to import the ErrorCode class and a recovery kernel. The example shows a kernel that deletes the particle above the surface. from parcels import ErrorCode

  1. Another way is to create your own AdvectionRK4_3D kernel limiting particles to the surface (example below). Since a negative depth means a particle above the surface, we limited the depth up to 0 meters using the max function.

    def AdvectionRK4_3D_surface_limited(particle, fieldset, time):
    """Run AdvectionRK4_3D kernel keeping particles on 0-depth, avoiding Through-surface error."""
    (u1, v1, w1) = fieldset.UVW[particle]
    
    lon1 = particle.lon + u1 * 0.5 * particle.dt
    lat1 = particle.lat + v1 * 0.5 * particle.dt
    dep1 = particle.depth + w1 * 0.5 * particle.dt
    dep1 = max(0, dep1)
    (u2, v2, w2) = fieldset.UVW[time + 0.5 * particle.dt, dep1, lat1, lon1, particle]
    
    lon2 = particle.lon + u2 * 0.5 * particle.dt
    lat2 = particle.lat + v2 * 0.5 * particle.dt
    dep2 = particle.depth + w2 * 0.5 * particle.dt
    dep2 = max(0, dep2)
    (u3, v3, w3) = fieldset.UVW[time + 0.5 * particle.dt, dep2, lat2, lon2, particle]
    
    lon3 = particle.lon + u3 * particle.dt
    lat3 = particle.lat + v3 * particle.dt
    dep3 = particle.depth + w3 * particle.dt
    dep3 = max(0, dep3)
    (u4, v4, w4) = fieldset.UVW[time + particle.dt, dep3, lat3, lon3, particle]
    
    particle.lon += (u1 + 2 * u2 + 2 * u3 + u4) / 6.0 * particle.dt
    particle.lat += (v1 + 2 * v2 + 2 * v3 + v4) / 6.0 * particle.dt
    particle.depth += (w1 + 2 * w2 + 2 * w3 + w4) / 6.0 * particle.dt
    
    particle.depth = max(0, particle.depth)
erikvansebille commented 2 years ago

Thanks for picking up this answer, @axnsantana! You are completely right! I'll close the issue now