OceanParcels / Parcels

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

Kernel Died Automatically (Ocean Parcels) #1367

Closed BenOwusu084 closed 1 year ago

BenOwusu084 commented 1 year ago

This code from ocean parcels gives me issues with kernels when I ran the particles for 20 days and more but works fine if the simulation days are less than 20

#This script shows the trajectory of a single particle through the ocean
#************************************************************************

import numpy as np
import matplotlib.pyplot as plt
import netCDF4 as nc4
from mpl_toolkits.basemap import Basemap
import matplotlib.colors as colors
from parcels import FieldSet, ParticleSet, JITParticle, AdvectionRK4_3D, ErrorCode, ScipyParticle, plotTrajectoriesFile
from glob import glob
from datetime import timedelta as delta
#from parcels import logger, XarrayDecodedFilter
#logger.addFilter(XarrayDecodedFilter())  # Add a filter for the xarray decoding warning

location  = '/home/bejamino/Desktop/DataFiles/'

evel_Data = nc4.Dataset(location + 'EVELMASS/Monthly_Avg_Files/EVEL_01_avg.nc')
wvel_Data = nc4.Dataset(location + 'WVELMASS/Monthly_Avg_Files/WVEL_01_avg.nc')
nvel_Data = nc4.Dataset(location + 'NVELMASS/Monthly_Avg_Files/NVEL_01_avg.nc')

#loading the coordinates of the data
lonGrid       = evel_Data.variables['longitude']
latGrid       = evel_Data.variables['latitude'] 

#Making the streamplot at the surface in January
EVEL      = evel_Data.variables['EVELMASS'][0,0,:,:]
NVEL      = nvel_Data.variables['NVELMASS'][0,0,:,:]

ufiles = sorted(glob(location + 'EVELMASS/Monthly_Avg_Files/EVEL_*.nc'))
vfiles = sorted(glob(location + 'NVELMASS/Monthly_Avg_Files/NVEL_*.nc'))
wfiles = sorted(glob(location + 'WVELMASS/Monthly_Avg_Files/WVEL_*.nc'))

filenames = {'U': ufiles, 'V': vfiles,'W': wfiles}

variables = {'U': 'EVELMASS', 'V': 'NVELMASS','W': 'WVELMASS'}

dimensions = {'U': {'lon': 'longitude', 'lat': 'latitude', 'depth': 'Z',  'time': 'time'},
              'V': {'lon': 'longitude', 'lat': 'latitude', 'depth': 'Z',  'time': 'time'},
              'W': {'lon': 'longitude', 'lat': 'latitude', 'depth': 'Zl', 'time': 'time'}}

fset = FieldSet.from_netcdf(filenames, variables, dimensions, time_periodic=True)

pset = ParticleSet.from_line(fieldset=fset, pclass=JITParticle,
                             size   = 1, #the number of particles that are released
                             start  = (-30, -40), #defining the starting location (lon,lat)
                             finish = (-30, 30), #defining the finishing location (lon,lat)
                             depth  =100) #the depth 

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

#defining the particle locations for 8 particle sets (somewhere in the Atlantic)
#p_lons=np.array([-30,-30,-30,-30,-30,-30,-30,-30])
#p_lats=np.array([-40,-30,-20,-10,0,10,20,30])
#p_dep=np.array([100,100,100,100,100,100,100,100])

#pset = ParticleSet(fieldset=fset, pclass=JITParticle, lon=p_lons, lat=p_lats, depth=p_dep)

#kernels = pset.Kernel(AdvectionRK4_3D)
#pset.execute(kernels, runtime=delta(days=30), dt=delta(hours=6))

depth_level = 8
print("Level[%d] depth is: [%g %g]" % (depth_level, fset.W.grid.depth[depth_level], fset.W.grid.depth[depth_level+1]))

def colorGradient(val, minval, maxval):
   # min_depth = 100
    #max_depth = 1000

    if val < minval:
        return (0,0,0) #returns black particle

    if val > maxval:
        return (1,0,0) #returns red particle

#interpolating continuously the depth
    v = (val - minval) / (maxval - minval)
    return (v,0,0)

#Generating the trajectories for the particles
nParticles = len(pset)
nDays = 20
nSteps = 20 #this show how often the particle is propagated for nDays (if it goes in total beyong 360 is crushes)
kernels = pset.Kernel(AdvectionRK4_3D)
plocations = [[] for _ in range(nParticles)]
for i in range(nSteps):
    pset.execute(kernels, runtime=delta(days=nDays), dt=delta(hours=6)
                 ,recovery={ErrorCode.ErrorOutOfBounds: DeleteParticle})
    for k, p in enumerate(pset):
        if k == 0:
            print(f"Step {i}, particle {k}:")
            print(p)
            print(f"lon: {p.lon}, lat: {p.lat}, depth: {p.depth}")
        plocations[k].append((p.lon, p.lat, p.depth))

#making the plot
plt.ion()
fig = plt.figure(figsize=(12,8))

# miller projection with Basemaps
map = Basemap(projection='cyl',llcrnrlat=-90,urcrnrlat=90,llcrnrlon=-180,urcrnrlon=180)

# plot coastlines, draw label meridians and parallels.
map.drawmeridians(np.arange(0,360,60),labels=[False,False,False,True],linewidth=(0.0))
map.drawparallels(np.arange(-80,80,20),labels=[True,False,False,False],linewidth=(0.0))

# Fill the continents 
map.fillcontinents(color='0.7',lake_color='0.7')
map.drawmapboundary(fill_color='0.7')
map.fillcontinents(color='#ffe2ab',lake_color='0.7')
map.drawcoastlines()

magnitude  = (EVEL ** 2 + NVEL ** 2) ** 0.5       

#Plotting the stream pattern
map.streamplot(x=lonGrid, y=latGrid, u=EVEL, v=NVEL, linewidth=2, density=2, color=magnitude)

#plt.title('Put the title here', fontsize=8)

plt.xlabel('Longitude [$^\circ$]', labelpad=20, fontsize=12)
plt.ylabel('Latitude [$^\circ$]', labelpad=40, fontsize=12)

# Plot trajectories    
min_depth, max_depth = 100, 1000

#colorBy = "depth"
colorBy = "step"
for j, pos in enumerate(plocations):
    for i in range(nSteps):
        lon, lat, depth = pos[i]
        if colorBy == "depth":
            map.scatter(lon, lat, s=4, color=colorGradient(depth, min_depth, max_depth))
        elif colorBy == "step":
            map.scatter(lon, lat, s=4, color=colorGradient(i, 0, nSteps))
erikvansebille commented 1 year ago

Can I ask, why did you close this? Did you fix/solve the problem?

BenOwusu084 commented 1 year ago

Hi,

Yes it was solved but I realised it was only partially solved. I only had to include the function (allow_time_extrapolation= True) and the particle simulate beyond 360 days however is I run it for 100 days it gives an OutOfBounds Error. I will try and fix it but if I’m not then I’ll open it in GitHub. It also appears that particles that are very shallow hit the boundary and produces an error message.

Thanks for the enquiry anyway.

Regards, Benjamin


Benjamin Owusu
Doctoral Researcher - Biogeochemistry Ocean Modeling
Institute of Chemistry and biology of the Marine Environment
Carl-von-Ossietzky University of Oldenburg
Carl-von-Ossietzky Str. 9-11
Room: W16-0-010
26129 Oldenburg
+49 163 906 9586

From: Erik van Sebille @.***> Sent: 25 May 2023 17:14:02 To: OceanParcels/parcels Cc: Benjamin Owusu; State change Subject: Re: [OceanParcels/parcels] Kernel Died Automatically (Ocean Parcels) (Issue #1367)

ACHTUNG! Diese E-Mail kommt von Extern! WARNING! This email originated off-campus.

Can I ask, why did you close this? Did you fix/solve the problem?

— Reply to this email directly, view it on GitHubhttps://github.com/OceanParcels/parcels/issues/1367#issuecomment-1563084081, or unsubscribehttps://github.com/notifications/unsubscribe-auth/A6B3RXSGKKZQZLXNZXY5I5TXH5ZLVANCNFSM6AAAAAAYOUKTOA. You are receiving this because you modified the open/close state.Message ID: @.***>