Divide by zero error #13
Description
At certain bottom conditions, source depths, and receiver depths, I get a divide by zero error
fn = 'FCTD_SoundSpeed_Colosi_NESCAFE.mat'
fctd_mat = io.loadmat(fn)
fctd = xr.DataArray(
np.flip(fctd_mat['cc'], axis=0),
dims=['range','depth'],
coords={
'range':fctd_mat['r'].flatten()*1000,
'depth':fctd_mat['z'].flatten(),
'latitude':('range', fctd_mat['lat'].flatten()),
'longitude':('range', fctd_mat['lon'].flatten())
}
).sortby('range')
fctd_pilot = fctd.assign_coords({'range':fctd.range.values - (201222-153500)})
# Add 10000 meter
fctd_pilot = xr.concat(
(
fctd_pilot,
xr.ones_like(fctd_pilot[:,-2:-1]).assign_coords({'depth':[10000]})*1640
),
dim='depth'
)
# %%
fctd_grid = fctd_pilot.interp({'range':np.arange(0,153500,500), 'depth':np.arange(0,6000)})
bz, az = signal.butter(4, 1/100, fs=1)
br, ar = signal.butter(4, 1/10e3, fs=1/500)
fctd_smooth = xrs.filtfilt(xrs.filtfilt(fctd_grid, dim='range', b=br, a=ar), dim='depth', b=bz, a=az)
# %%
bathy_depths = np.linspace(5050, 4950, 100)
#bathy_depths = np.linspace(5000, 5000, 100)
#bathy_depths = np.linspace(4908, 4908, 100)
bathy_ranges = np.linspace(0,155e3,100)
bathy = xr.DataArray(bathy_depths, dims=['range'], coords={'range':bathy_ranges})
# %%
environment = pr.OceanEnvironment2D(sound_speed=fctd_smooth, bathymetry=bathy)
# %%
source_depth = 2350
source_range = 0
launch_angles = np.linspace(-30,30,6000)
reciever_range = 153.5e3
rays = pr.shoot_rays(
source_depth,
source_range,
launch_angles,
reciever_range,
1000,
environment,
)
# %%
rays.plot_ray_fan(alpha=0.005)
# %%
rays.plot_time_front()
# %%
eigs = pr.find_eigenrays(
rays,
[900, 925, 950, 975],
source_depth,
source_range,
reciever_range,
10000,
environment,
ztol=1,
max_iter=20,
)
# %%
eigs.plot(alpha=0.3)
# %%
eigs.save_mat('eigenrays2350_900mrec.mat')---------------------------------------------------------------------------
RemoteTraceback Traceback (most recent call last)
RemoteTraceback:
"""
Traceback (most recent call last):
File "/home/khoekstra/miniconda3/envs/NESCAFE/lib/python3.12/multiprocessing/pool.py", line 125, in worker
result = (True, func(*args, **kwds))
^^^^^^^^^^^^^^^^^^^
File "/mnt/c/users/kmhoe/documents/research/python/pygenray/src/pygenray/eigenrays.py", line 113, in _find_single_eigenray
ray = pr.shoot_ray(source_depth, source_range, -bisection_theta, receiver_range, num_range_save, environment)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/mnt/c/users/kmhoe/documents/research/python/pygenray/src/pygenray/launch_rays.py", line 229, in shoot_ray
sols, full_ray, n_bottom, n_surface = _shoot_ray_array(
^^^^^^^^^^^^^^^^^
File "/mnt/c/users/kmhoe/documents/research/python/pygenray/src/pygenray/launch_rays.py", line 326, in _shoot_ray_array
sol = _shoot_ray_segment(
^^^^^^^^^^^^^^^^^^^
File "/mnt/c/users/kmhoe/documents/research/python/pygenray/src/pygenray/launch_rays.py", line 574, in _shoot_ray_segment
sol = scipy.integrate.solve_ivp(
^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/khoekstra/miniconda3/envs/NESCAFE/lib/python3.12/site-packages/scipy/integrate/_ivp/ivp.py", line 655, in solve_ivp
message = solver.step()
^^^^^^^^^^^^^
File "/home/khoekstra/miniconda3/envs/NESCAFE/lib/python3.12/site-packages/scipy/integrate/_ivp/base.py", line 197, in step
success, message = self._step_impl()
^^^^^^^^^^^^^^^^^
File "/home/khoekstra/miniconda3/envs/NESCAFE/lib/python3.12/site-packages/scipy/integrate/_ivp/rk.py", line 144, in _step_impl
y_new, f_new = rk_step(self.fun, t, y, self.f, h, self.A,
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/khoekstra/miniconda3/envs/NESCAFE/lib/python3.12/site-packages/scipy/integrate/_ivp/rk.py", line 64, in rk_step
K[s] = fun(t + c * h, y + dy)
^^^^^^^^^^^^^^^^^^^^^^
File "/home/khoekstra/miniconda3/envs/NESCAFE/lib/python3.12/site-packages/scipy/integrate/_ivp/base.py", line 154, in fun
return self.fun_single(t, y)
^^^^^^^^^^^^^^^^^^^^^
File "/home/khoekstra/miniconda3/envs/NESCAFE/lib/python3.12/site-packages/scipy/integrate/_ivp/base.py", line 23, in fun_wrapped
return np.asarray(fun(t, y), dtype=dtype)
^^^^^^^^^
File "/home/khoekstra/miniconda3/envs/NESCAFE/lib/python3.12/site-packages/scipy/integrate/_ivp/ivp.py", line 593, in fun
return fun(t, x, *args)
^^^^^^^^^^^^^^^^
ZeroDivisionError: division by zero
"""
The above exception was the direct cause of the following exception:
ZeroDivisionError Traceback (most recent call last)
Cell In[31], line 1
----> 1 eigs = pr.find_eigenrays(
2 rays,
3 [900, 925, 950, 975],
4 source_depth,
5 source_range,
6 reciever_range,
7 10000,
8 environment,
9 ztol=1,
10 max_iter=20,
11 )
File /mnt/c/users/kmhoe/documents/research/python/pygenray/src/pygenray/eigenrays.py:80, in find_eigenrays(rays, receiver_depths, source_depth, source_range, receiver_range, num_range_save, environment, ztol, max_iter, num_workers)
78 num_workers = mp.cpu_count()
79 with mp.Pool(num_workers) as pool:
---> 80 results = list(tqdm(pool.imap(_find_single_eigenray, args_list), total=len(args_list), desc="Finding eigenrays"))
82 # Filter out None results and add successful rays
83 for result in results:
File ~/miniconda3/envs/NESCAFE/lib/python3.12/site-packages/tqdm/std.py:1181, in tqdm.__iter__(self)
1178 time = self._time
1180 try:
-> 1181 for obj in iterable:
1182 yield obj
1183 # Update and possibly print the progressbar.
1184 # Note: does not call self.update(1) for speed optimisation.
File ~/miniconda3/envs/NESCAFE/lib/python3.12/multiprocessing/pool.py:873, in IMapIterator.next(self, timeout)
871 if success:
872 return value
--> 873 raise value
ZeroDivisionError: division by zero