Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Break-up analytical function and test #56

Closed
wants to merge 4 commits into from
Closed

Break-up analytical function and test #56

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
45f8c92
wrap up issue#49, see if dask[array] works for pip
MaceKuailv Jun 18, 2023
c00e16f
make it possible to use string as input for oceinterp
MaceKuailv Jun 18, 2023
2671915
break-up analytical
MaceKuailv Jun 18, 2023
93723bc
remove part of analytical step to cross_cell_wall
MaceKuailv Jun 18, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion pyproject.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@ dependencies = [
"numpy",
"pandas",
"scipy",
"dask",
"dask[array]",
"xarray"
]
description = "A python package that interpolates data from ocean dataset from both Eulerian and Lagrangian perspective. "
Expand Down
8 changes: 4 additions & 4 deletions seaduck/eulerian.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -162,13 +162,13 @@ def from_latlon(self, x=None, y=None, z=None, t=None, data=None):
raise ValueError("Shapes of input coordinates are not compatible")

if isinstance(x, float):
x = np.array([1.0]) * x
x = np.ones(self.N, float) * x
if isinstance(y, float):
y = np.array([1.0]) * y
y = np.ones(self.N, float) * y
if isinstance(z, float):
z = np.array([1.0]) * z
z = np.ones(self.N, float) * z
if isinstance(t, float):
t = np.array([1.0]) * t
t = np.ones(self.N, float) * t

for thing in [x, y, z, t]:
if thing is None:
Expand Down
296 changes: 159 additions & 137 deletions seaduck/lagrangian.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -535,7 +535,7 @@ def _out_of_bound(self): # pragma: no cover
z_out = False
return np.logical_or(np.logical_or(x_out, y_out), z_out)

def trim(self, tol=1e-6):
def trim(self, tol=1e-12):
"""Move the particles from outside the cell into the cell.

At the same time change the velocity accordingly.
Expand Down Expand Up @@ -654,125 +654,15 @@ def _contract(self): # pragma: no cover

self.t[out] += contract_time

def update_after_cell_change(self):
"""Update properties after particles cross wall.

A wall event is triggered when particle reached the wall.
This method handle the coords translation as a particle cross
a wall.
"""
if self.face is not None:
self.face = self.face.astype(int)
self.iy = self.iy.astype(int)
self.ix = self.ix.astype(int)
if self.iz is not None:
self.iz = self.iz.astype(int)
if self.izl_lin is not None:
self.izl_lin = self.izl_lin.astype(int)

if self.ocedata.readiness["Z"]:
self.rel.update(self.ocedata.find_rel_v(self.dep))
if self.ocedata.readiness["h"] == "local_cartesian":
if self.face is not None:
self.bx, self.by = (
self.ocedata.XC[self.face, self.iy, self.ix],
self.ocedata.YC[self.face, self.iy, self.ix],
)
self.cs, self.sn = (
self.ocedata.CS[self.face, self.iy, self.ix],
self.ocedata.SN[self.face, self.iy, self.ix],
)
self.dx, self.dy = (
self.ocedata.dX[self.face, self.iy, self.ix],
self.ocedata.dY[self.face, self.iy, self.ix],
)

else:
self.bx, self.by = (
self.ocedata.XC[self.iy, self.ix],
self.ocedata.YC[self.iy, self.ix],
)
self.cs, self.sn = (
self.ocedata.CS[self.iy, self.ix],
self.ocedata.SN[self.iy, self.ix],
)
self.dx, self.dy = (
self.ocedata.dX[self.iy, self.ix],
self.ocedata.dY[self.iy, self.ix],
)
elif self.ocedata.readiness["h"] == "oceanparcel":
if self.face is not None:
self.bx, self.by = (
self.ocedata.XC[self.face, self.iy, self.ix],
self.ocedata.YC[self.face, self.iy, self.ix],
)

else: # pragema: no cover
self.bx, self.by = (
self.ocedata.XC[self.iy, self.ix],
self.ocedata.YC[self.iy, self.ix],
)

elif self.ocedata.readiness["h"] == "rectilinear":
self.bx = self.ocedata.lon[self.ix]
self.by = self.ocedata.lat[self.iy]
self.cs = np.ones_like(self.bx)
self.sn = np.zeros_like(self.bx)
self.dx = self.ocedata.dlon[self.ix] * np.cos(self.by * np.pi / 180)
self.dy = self.ocedata.dlat[self.iy]

if self.izl_lin is not None:
self.bzl_lin = self.ocedata.Zl[self.izl_lin]
self.dzl_lin = self.ocedata.dZl[self.izl_lin - 1]
if self.dz is not None:
self.dz = self.ocedata.dZ[self.iz]
try:
self.px, self.py = self.get_px_py()
self.rx, self.ry = find_rx_ry_oceanparcel(
self.lon, self.lat, self.px, self.py
)
except AttributeError:
# if True:
dlon = to_180(self.lon - self.bx)
dlat = to_180(self.lat - self.by)
self.rx = (
(dlon * np.cos(self.by * np.pi / 180) * self.cs + dlat * self.sn)
* DEG2M
/ self.dx
)
self.ry = (
(dlat * self.cs - dlon * self.sn * np.cos(self.by * np.pi / 180))
* DEG2M
/ self.dy
)
if self.rzl_lin is not None:
self.rzl_lin = (self.dep - self.bzl_lin) / self.dzl_lin

def analytical_step(self, tf, which=None):
"""Integrate the particle with velocity.

The core method.
A set of particles trying to integrate for time tf
(could be negative).
at the end of the call, every particle are either:
1. ended up somewhere within the cell after time tf.
2. ended up on a cell wall before
(if tf is negative, then "after") tf.

Parameters
----------
tf: float, numpy.ndarray
The longest duration of the simulation for each particle.
which: numpy.ndarray, optional
Boolean or int array that specify the subset of points to
do the operation.
"""
def _subset_velocity_position(self, tf, which):
"""See analytical_step."""
if which is None:
which = np.ones(self.N).astype(bool)
if isinstance(tf, float):
tf = np.array([tf for i in range(self.N)])

tf = tf[which]
t_now = self.t[which]

if self.rzl_lin is not None:
xs = [self.rx[which], self.ry[which], self.rzl_lin[which] - 1 / 2]
Expand All @@ -781,46 +671,32 @@ def analytical_step(self, tf, which=None):
xs = [x_temp, self.ry[which], np.zeros_like(x_temp)]
us = [self.u[which], self.v[which], self.w[which]]
dus = [self.du[which], self.dv[which], self.dw[which]]
return which, tf, t_now, us, dus, xs

ts = time2wall(xs, us, dus)

tend, the_t = which_early(tf, ts)
self.t[which] += the_t

def _move_within_cell(self, the_t, t_now, us, dus, xs):
t_now += the_t
new_x = []
new_u = []
for i in range(3):
x_move = stationary(the_t, us[i], dus[i], 0)
new_u.append(us[i] + dus[i] * x_move)
new_x.append(x_move + xs[i])

# for rr in new_x:
# if np.logical_or(rr>0.6,rr<-0.6).any():
# where = np.where(np.logical_or(rr>0.6,rr<-0.6))[0][0]
# raise ValueError(
# f'Particle way out of bound.tend = {tend[where]},'
# f' the_t = {the_t[where]},'
# f' rx = {new_x[0][where]},ry = {new_x[1][where]},rz = {new_x[2][where]}'
# f'start with u = {self.u[where]}, du = {self.du[where]}, x={self.rx[where]}'
# f'start with v = {self.v[where]}, dv = {self.dv[where]}, y={self.ry[where]}'
# f'start with w = {self.w[where]}, dw = {self.dv[where]}, z={self.rzl_lin[where]}'
# )
for rr in new_x:
if np.logical_or(rr > 0.6, rr < -0.6).any():
where = np.where(np.logical_or(rr > 0.6, rr < -0.6))[0][0]
raise ValueError(
f"Particle way out of bound.tend = {tend[where]},"
f"Particle way out of bound."
# f"tend = {tend[where]},"
f" the_t = {the_t[where]},"
f" rx = {new_x[0][where]},ry = {new_x[1][where]},rz = {new_x[2][where]}"
f"start with u = {self.u[where]}, du = {self.du[where]}, x={self.rx[where]}"
f"start with v = {self.v[where]}, dv = {self.dv[where]}, y={self.ry[where]}"
f"start with w = {self.w[where]}, dw = {self.dv[where]}, z={self.rzl_lin[where]}"
)
self.rx[which], self.ry[which], temp = new_x
if self.rzl_lin is not None:
self.rzl_lin[which] = temp + 1 / 2
return t_now, new_x, new_u

self.u[which], self.v[which], self.w[which] = new_u
def _sync_latlondep_before_cross(self):
try:
px, py = self.px, self.py
w = self.get_f_node_weight()
Expand Down Expand Up @@ -850,10 +726,53 @@ def analytical_step(self, tf, which=None):
self.by,
bzl_lin,
)

def analytical_step(self, tf, which=None):
"""Integrate the particle with velocity.

The core method.
A set of particles trying to integrate for time tf
(could be negative).
at the end of the call, every particle are either:
1. ended up somewhere within the cell after time tf.
2. ended up on a cell wall before
(if tf is negative, then "after") tf.

Parameters
----------
tf: float, numpy.ndarray
The longest duration of the simulation for each particle.
which: numpy.ndarray, optional
Boolean or int array that specify the subset of points to
do the operation.
"""
which, tf, t_now, us, dus, xs = self._subset_velocity_position(tf, which)

ts = time2wall(xs, us, dus)

tend, the_t = which_early(tf, ts)

t_now, new_x, new_u = self._move_within_cell(the_t, t_now, us, dus, xs)

# Could potentially move this block all the way back
self.t[which] += t_now
self.rx[which], self.ry[which], temp = new_x
if self.rzl_lin is not None:
self.rzl_lin[which] = temp + 1 / 2

self.u[which], self.v[which], self.w[which] = new_u

self._sync_latlondep_before_cross()

if self.save_raw:
# record the moment just before crossing the wall
# or the moment reaching destination.
self.note_taking(which)
return tend

def _cross_cell_wall_index(self, tend, which):
if which is None:
which = np.ones(self.N).astype(bool)
type1 = tend <= 3
translate = {0: 2, 1: 3, 2: 1, 3: 0}
# left # right # down # up
Expand Down Expand Up @@ -906,6 +825,109 @@ def analytical_step(self, tf, which=None):
if self.izl_lin is not None:
self.izl_lin[which] = tiz

def _cross_cell_wall_read(self):
"""Update properties after particles cross wall.

A wall event is triggered when particle reached the wall.
This method handle the coords translation as a particle cross
a wall.
"""
if self.face is not None:
self.face = self.face.astype(int)
self.iy = self.iy.astype(int)
self.ix = self.ix.astype(int)
if self.iz is not None:
self.iz = self.iz.astype(int)
if self.izl_lin is not None:
self.izl_lin = self.izl_lin.astype(int)

if self.ocedata.readiness["Z"]:
self.rel.update(self.ocedata.find_rel_v(self.dep))
if self.ocedata.readiness["h"] == "local_cartesian":
if self.face is not None:
self.bx, self.by = (
self.ocedata.XC[self.face, self.iy, self.ix],
self.ocedata.YC[self.face, self.iy, self.ix],
)
self.cs, self.sn = (
self.ocedata.CS[self.face, self.iy, self.ix],
self.ocedata.SN[self.face, self.iy, self.ix],
)
self.dx, self.dy = (
self.ocedata.dX[self.face, self.iy, self.ix],
self.ocedata.dY[self.face, self.iy, self.ix],
)

else:
self.bx, self.by = (
self.ocedata.XC[self.iy, self.ix],
self.ocedata.YC[self.iy, self.ix],
)
self.cs, self.sn = (
self.ocedata.CS[self.iy, self.ix],
self.ocedata.SN[self.iy, self.ix],
)
self.dx, self.dy = (
self.ocedata.dX[self.iy, self.ix],
self.ocedata.dY[self.iy, self.ix],
)
elif self.ocedata.readiness["h"] == "oceanparcel":
if self.face is not None:
self.bx, self.by = (
self.ocedata.XC[self.face, self.iy, self.ix],
self.ocedata.YC[self.face, self.iy, self.ix],
)

else: # pragema: no cover
self.bx, self.by = (
self.ocedata.XC[self.iy, self.ix],
self.ocedata.YC[self.iy, self.ix],
)

elif self.ocedata.readiness["h"] == "rectilinear":
self.bx = self.ocedata.lon[self.ix]
self.by = self.ocedata.lat[self.iy]
self.cs = np.ones_like(self.bx)
self.sn = np.zeros_like(self.bx)
self.dx = self.ocedata.dlon[self.ix] * np.cos(self.by * np.pi / 180)
self.dy = self.ocedata.dlat[self.iy]

if self.izl_lin is not None:
self.bzl_lin = self.ocedata.Zl[self.izl_lin]
self.dzl_lin = self.ocedata.dZl[self.izl_lin - 1]
if self.dz is not None:
self.dz = self.ocedata.dZ[self.iz]

def _cross_cell_wall_rel(self):
try:
self.px, self.py = self.get_px_py()
self.rx, self.ry = find_rx_ry_oceanparcel(
self.lon, self.lat, self.px, self.py
)
except AttributeError:
# if True:
dlon = to_180(self.lon - self.bx)
dlat = to_180(self.lat - self.by)
self.rx = (
(dlon * np.cos(self.by * np.pi / 180) * self.cs + dlat * self.sn)
* DEG2M
/ self.dx
)
self.ry = (
(dlat * self.cs - dlon * self.sn * np.cos(self.by * np.pi / 180))
* DEG2M
/ self.dy
)
if self.rzl_lin is not None:
self.rzl_lin = (self.dep - self.bzl_lin) / self.dzl_lin

def cross_cell_wall(self, tend, which=None):
if which is None:
which = np.ones(self.N).astype(bool)
self._cross_cell_wall_index(tend, which)
self._cross_cell_wall_read()
self._cross_cell_wall_rel()

def deepcopy(self):
"""Return a clone of the object.

Expand Down Expand Up @@ -955,8 +977,8 @@ def to_next_stop(self, t1):
trim_tol = 1e-10
self.trim(tol=trim_tol)
logging.debug(sum(todo), "left")
self.analytical_step(tf, todo)
self.update_after_cell_change()
tend = self.analytical_step(tf, todo)
self.cross_cell_wall(tend, which=todo)
if self.transport:
self.get_vol()
self.get_u_du(todo)
Expand Down
Loading