diff --git a/docs/notebook/AVISO.ipynb b/docs/notebook/AVISO.ipynb index 59b2d9b8..2eb8ac3e 100644 --- a/docs/notebook/AVISO.ipynb +++ b/docs/notebook/AVISO.ipynb @@ -38,7 +38,13 @@ "attachments": {}, "cell_type": "markdown", "id": "256affe3", - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "source": [ "## Loading data\n", "\n", @@ -311,6 +317,10 @@ "execution_count": null, "id": "124718ad", "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, "tags": [] }, "outputs": [], @@ -347,6 +357,10 @@ "execution_count": null, "id": "6d8ee9cc", "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, "tags": [] }, "outputs": [], @@ -371,6 +385,10 @@ "execution_count": null, "id": "25f862a7", "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, "tags": [] }, "outputs": [], @@ -393,6 +411,10 @@ "execution_count": null, "id": "0e0fc0c8", "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, "tags": [ "hide-input" ] @@ -412,7 +434,13 @@ "attachments": {}, "cell_type": "markdown", "id": "c7046590", - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "source": [ "Voila!" ] @@ -422,6 +450,10 @@ "execution_count": null, "id": "5338b2f2", "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, "tags": [ "hide-input" ] @@ -444,6 +476,232 @@ "source": [ "**Fig.2** The trajectories of particles advected by AVISO-derived surface velocity field." ] + }, + { + "cell_type": "markdown", + "id": "072f7015-6e4e-473e-85e3-ba15388963b9", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "## Advanced use of `Particle`\n", + "\n", + "In this subsection, we are going to demonstrate how to access the analytical trajectories of particles. We are also going to demonstrate the flexibility of Particle release.\n", + "\n", + "Most of the notebooks release particles at the same time. We are going to do it differently this time. They are going to be released at 64W and 55S and 64S. Crucially, all of the particles are released at a different time." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b1a404b2-c59c-4d1a-87be-a49230639ebe", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [ + "hide-input" + ] + }, + "outputs": [], + "source": [ + "N = 777\n", + "t_earliest = sd.utils.convert_time(\"1970-01-01\")\n", + "t_final = sd.utils.convert_time(\"1970-02-01\")\n", + "t_latest = sd.utils.convert_time(\"1970-03-03\") # damn you, February.\n", + "t = np.linspace(t_earliest, t_latest, N)\n", + "\n", + "number_of_loop = 6.18\n", + "# y = -65.0+np.abs(10-np.linspace(0,number_of_loop*20,N)%20)\n", + "y = -59.5 + 4.5 * np.sin(number_of_loop * np.linspace(-3.14, 3.14, N))\n", + "x = np.ones(N) * (-64.0)\n", + "z = -np.ones(N)" + ] + }, + { + "cell_type": "markdown", + "id": "19c51858-ecf6-4229-adf9-0c6051d00405", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "To better demonstrate, here is a plot showing the release pattern of the particle. The time is referenced against the final time of the simulation. Note that every single particle will be released at a different time and there are going to be both forward and backward particles in the same simulation.\n", + "\n", + "This is as if there is a ship is commuting in drake passage while releasing particles." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1f5c63f8-6074-41d7-b0f6-513ee484e89a", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [ + "hide-input" + ] + }, + "outputs": [], + "source": [ + "plt.plot(t, y, \"o\", markersize=1)\n", + "plt.xticks(\n", + " [t_earliest, t_final - 15 * 86400, t_final, t_final + 15 * 86400, t_latest],\n", + " [\"-31 days\", \"-15 days\", \"0\", \"+15 days\", \"+31 days\"],\n", + ")\n", + "plt.ylabel(\"Latitude\")\n", + "plt.title(\"Latitude of particles released at 64W\")\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "c4c7c3de-7c11-428e-b342-c854cc38e4ed", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "**Fig.3** Pattern of particle release. Time is relative to the final time." + ] + }, + { + "cell_type": "markdown", + "id": "2b5e335f-198a-4ba1-b3ff-85f5cc39d31b", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "The only difference in preparation between this one and the previous simulation is that we have `save_raw = True`. This means the particles will record all the necessary informations to reconstruct the analytical trajectory." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e329dca7-f916-4b08-89da-e99264f2aa76", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "p = sd.Particle(\n", + " x=x,\n", + " y=y,\n", + " z=z,\n", + " t=t,\n", + " data=bathtub,\n", + " uname=\"u\",\n", + " vname=\"v\",\n", + " wname=None,\n", + " callback=interested_in,\n", + " save_raw=True,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "1ef2b1d9-eceb-4418-9911-683116a8beb1", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "No tricks need to be played while execution." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "eedf2bfd-5e4c-4387-b1d1-c7d81ed209f8", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "%%time\n", + "p.to_next_stop(t_final)" + ] + }, + { + "cell_type": "markdown", + "id": "850b27c8-5bcd-4d5e-bba0-5b88a380cf2b", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "When `save_raw = True` is selected, the `Particle` object records location and velocity information everytime velocity is updated. The following plot is plotted from longitude and latitudes of particles crossing cell walls. Each trajectory is colored based on the time of release, Purple is the earliest, red is the latest (furthest into the future)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0d6a4387-b10e-489c-a527-47d54b3f101f", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [ + "hide-input" + ] + }, + "outputs": [], + "source": [ + "plt.figure(figsize=(9, 16))\n", + "ax = plt.axes(projection=ccrs.PlateCarree())\n", + "rainbow = plt.get_cmap(\"rainbow\")\n", + "for i in range(0, N):\n", + " color = rainbow(t[i] / 2 / t_final)\n", + " ax.plot([x[i]] + p.xxlist[i], [y[i]] + p.yylist[i], color=color, lw=0.2)\n", + "ax.coastlines()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "80a5ee7e-4791-46cb-9635-2084c262449d", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "**Fig.4** Trajectory of particles released at different time. Warm color are particles released after the final time, and cold colors are those released before the final time." + ] } ], "metadata": { @@ -463,7 +721,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.4" + "version": "3.11.3" } }, "nbformat": 4, diff --git a/seaduck/eulerian.py b/seaduck/eulerian.py index be7cdef6..41137107 100644 --- a/seaduck/eulerian.py +++ b/seaduck/eulerian.py @@ -501,51 +501,6 @@ def get_f_node_weight(self): """Find weight for the corner points interpolation.""" return weight_f_node(self.rx, self.ry) - def _get_lon_lat(self): # pragma: no cover - """Return the lat-lon value based on relative coordinate. - - This method only work if the dataset has readiness['h'] == 'oceanparcel'. - """ - px, py = self.get_px_py() - w = self.get_f_node_weight() - lon = np.einsum("nj,nj->n", w, px.T) - lat = np.einsum("nj,nj->n", w, py.T) - return lon, lat - - def _get_needed( - self, var_name, knw, required=None, prefetched=None, **kwarg - ): # pragma: no cover - if required is None: - required = self.ocedata._ds[var_name].dims - ind = self.fatten(knw, required=required, **kwarg) - if len(ind) != len(self.ocedata._ds[var_name].dims): - raise IndexError( - "Dimension mismatch. " - "Please check if the Position objects " - "have all the dimensions needed" - ) - if prefetched is None: - return smart_read(self.ocedata[var_name], ind) - else: - return prefetched[ind] - - def _get_masked(self, knw, cuvwg="C", **kwarg): # pragma: no cover - ind = self.fatten(knw, four_d=True, **kwarg) - if self.it is not None: - ind = ind[1:] - if len(ind) != len(self.ocedata._ds["maskC"].dims): - raise IndexError( - "Dimension mismatch. " - "Please check if the Position objects " - "have all the dimensions needed" - ) - return get_masked(self.ocedata, ind, cuvwg=cuvwg) - - def _find_pk4d(self, knw, cuvwg="C"): # pragma: no cover - masked = self._get_masked(knw, cuvwg=cuvwg) - pk4d = find_pk_4d(masked, inheritance=knw.inheritance) - return pk4d - def _register_interpolation_input( self, var_name, knw, prefetched=None, prefetch_prefix=None ): diff --git a/seaduck/kernel_weight.py b/seaduck/kernel_weight.py index f3dde34d..3eb9f6a2 100644 --- a/seaduck/kernel_weight.py +++ b/seaduck/kernel_weight.py @@ -8,24 +8,23 @@ from seaduck.runtime_conf import compileable # default kernel for interpolation. -default_kernel = np.array( +DEFAULT_KERNEL = np.array( [[0, 0], [0, 1], [0, 2], [0, -1], [0, -2], [-1, 0], [-2, 0], [1, 0], [2, 0]] ) -default_inheritance = [ +DEFAULT_INHERITANCE = [ [0, 1, 2, 3, 4, 5, 6, 7, 8], [0, 1, 2, 3, 5, 7, 8], [0, 1, 3, 5, 7], [0], ] -weight_func: dict[int, str] = {} -default_kernels = [ - np.array([default_kernel[i] for i in doll]) for doll in default_inheritance +DEFAULT_KERNELS = [ + np.array([DEFAULT_KERNEL[i] for i in doll]) for doll in DEFAULT_INHERITANCE ] # It just tell you what the kernels look like -def show_kernels(kernels=default_kernels): +def show_kernels(kernels=DEFAULT_KERNELS): """Plot a small scatter plot of the shape of a list of kernel. Parameters @@ -476,7 +475,7 @@ def kernel_weight(kernel, kernel_type="interp", order=0): raise NotImplementedError("The shape of the kernel is neither cross or square") -default_interp_funcs = [kernel_weight_x(a_kernel) for a_kernel in default_kernels] +default_interp_funcs = [kernel_weight_x(a_kernel) for a_kernel in DEFAULT_KERNELS] def find_which_points_for_each_kernel(masked, inheritance="default"): @@ -510,7 +509,7 @@ def find_which_points_for_each_kernel(masked, inheritance="default"): return """ if inheritance == "default": - inheritance = default_inheritance + inheritance = DEFAULT_INHERITANCE already_wet = [] for i, doll in enumerate(inheritance): wet_1d = masked[:, np.array(doll)].all(axis=1) @@ -527,8 +526,8 @@ def get_weight_cascade( rx, ry, pk, - kernel_large=default_kernel, - inheritance=default_inheritance, + kernel_large=DEFAULT_KERNEL, + inheritance=DEFAULT_INHERITANCE, funcs=default_interp_funcs, ): """Compute the weight. @@ -569,7 +568,7 @@ def get_weight_cascade( return weight -def find_pk_4d(masked, inheritance=default_inheritance): +def find_pk_4d(masked, inheritance=DEFAULT_INHERITANCE): """Find the masking condition for 4D space time. See find_which_points_for_each_kernel @@ -675,7 +674,7 @@ class KnW: def __init__( self, - kernel=default_kernel, + kernel=DEFAULT_KERNEL, inheritance="auto", # None, or list of lists hkernel="interp", # dx,dy vkernel="nearest", # linear,dz diff --git a/seaduck/lagrangian.py b/seaduck/lagrangian.py index 4b1eb656..45886a61 100644 --- a/seaduck/lagrangian.py +++ b/seaduck/lagrangian.py @@ -102,10 +102,11 @@ def time2wall(pos_list, u_list, du_list, tf): for i in range(3): tl, tr = stationary_time(u_list[i], du_list[i], pos_list[i]) ul, ur = uleftright_from_udu(u_list[i], du_list[i], pos_list[i]) - cannot_left = ul * tf >= 0 - tl[cannot_left] = -np.sign(tf[cannot_left]) - cannot_right = ur * tf <= 0 - tr[cannot_right] = -np.sign(tf[cannot_right]) + sign = np.sign(tf) + cannot_left = -ul * sign <= 1e-12 # aroung 30000 years + tl[cannot_left] = -sign[cannot_left] + cannot_right = ur * tf <= 1e-12 + tr[cannot_right] = -sign[cannot_right] ts.append(tl) ts.append(tr) return ts @@ -415,7 +416,7 @@ def fillna(self): np.nan_to_num(self.dv, copy=False) np.nan_to_num(self.dw, copy=False) - def note_taking(self, subset_index=None): + def note_taking(self, subset_index=None, stamp=-1): """Record raw data into list of lists. This method is only called in save_raw = True particles. @@ -437,6 +438,11 @@ def note_taking(self, subset_index=None): raise AttributeError("This is not a particle_rawlist object") from exc if subset_index is None: subset_index = np.arange(self.N) + else: + if len(subset_index) != self.N: + raise IndexError( + "The subset used for notetaking is not generated from the int indexes" + ) for isub, ifull in enumerate(subset_index): if self.face is not None: self.fclist[ifull].append(self.face[isub]) @@ -459,6 +465,10 @@ def note_taking(self, subset_index=None): self.xxlist[ifull].append(self.lon[isub]) self.yylist[ifull].append(self.lat[isub]) self.zzlist[ifull].append(self.dep[isub]) + if isinstance(stamp, int): + self.vslist[ifull].append(stamp) + else: + self.vslist[ifull].append(stamp[isub]) def empty_lists(self): """Empty/Create the lists. @@ -489,6 +499,7 @@ def empty_lists(self): self.xxlist = [[] for i in range(self.N)] self.yylist = [[] for i in range(self.N)] self.zzlist = [[] for i in range(self.N)] + self.vslist = [[] for i in range(self.N)] def _out_of_bound(self): # pragma: no cover """Return particles that are out of the cell bound. @@ -856,22 +867,20 @@ def to_next_stop(self, t_stop): if self.callback is not None: bool_todo = np.logical_and(bool_todo, self.callback(self)) int_todo = np.where(bool_todo)[0] + sub = self.subset(int_todo) if len(int_todo) == 0: logging.info("Nothing left to simulate") return tf_used = tf[int_todo] - trim_tol = 1e-12 + trim_tol = 1e-14 for i in range(self.max_iteration): - if i > self.max_iteration * 0.95: - trim_tol = 1e-3 logging.debug(len(int_todo), "left") - sub = self.subset(int_todo) sub.trim(tol=trim_tol) tend = sub.analytical_step(tf_used) if self.save_raw: # record the moment just before crossing the wall # or the moment reaching destination. - self.note_taking(int_todo) + sub.note_taking(int_todo, stamp=tend) sub.cross_cell_wall(tend) if self.transport: @@ -883,12 +892,13 @@ def to_next_stop(self, t_stop): if self.callback is not None: bool_todo = np.logical_and(bool_todo, self.callback(sub)) int_todo = int_todo[bool_todo] - tf_used = tf_used[bool_todo] if len(int_todo) == 0: break + sub = self.subset(int_todo) + tf_used = tf_used[bool_todo] if self.save_raw: # record those who cross the wall - self.note_taking(int_todo) + sub.note_taking(int_todo, stamp=7) if i == self.max_iteration - 1: warnings.warn("maximum iteration count reached") @@ -971,7 +981,7 @@ def to_list_of_time( logging.info(np.datetime64(round(tl), "s")) if self.save_raw: # save the very start of everything. - self.note_taking() + self.note_taking(stamp=15) self.to_next_stop(tl) if update[i]: if self.too_large: # pragma: no cover diff --git a/seaduck/topology.py b/seaduck/topology.py index bb37b737..b6a1fe6d 100644 --- a/seaduck/topology.py +++ b/seaduck/topology.py @@ -5,9 +5,9 @@ from seaduck.runtime_conf import compileable -legal_tends = [0, 1, 2, 3] # up,down,left,right #list(llc_face_connect.columns) +# legal_tends = [0, 1, 2, 3] # up,down,left,right #list(LLC_FACE_CONNECT.columns) -llc_face_connect = np.array( +LLC_FACE_CONNECT = np.array( [ [1, 42, 12, 3], [2, 0, 11, 4], @@ -25,7 +25,7 @@ ] ) -directions = np.array([np.pi / 2, -np.pi / 2, np.pi, 0]) +DIRECTIONS = np.array([np.pi / 2, -np.pi / 2, np.pi, 0]) @compileable @@ -35,8 +35,8 @@ def llc_mutual_direction(face, neighbor_face, transitive=False): The compileable version of mutual direction for llc grid. See Topology.mutual direction for more detail. """ - connected_to_face = np.where(llc_face_connect[face] == neighbor_face) - connected_to_neigh = np.where(llc_face_connect[neighbor_face] == face) + connected_to_face = np.where(LLC_FACE_CONNECT[face] == neighbor_face) + connected_to_neigh = np.where(LLC_FACE_CONNECT[neighbor_face] == face) if len(connected_to_face[0]) == 0: found = False else: @@ -47,8 +47,8 @@ def llc_mutual_direction(face, neighbor_face, transitive=False): return connected_to_face[0][0], connected_to_neigh[0][0] elif transitive: common = -1 - for i in llc_face_connect[face]: - if i in llc_face_connect[neighbor_face]: + for i in LLC_FACE_CONNECT[face]: + if i in LLC_FACE_CONNECT[neighbor_face]: common = i break if common < 0: @@ -56,10 +56,10 @@ def llc_mutual_direction(face, neighbor_face, transitive=False): "The two faces does not share common face, transitive did not help" ) else: - edge_1 = np.where(llc_face_connect[face] == common)[0][0] - new_edge_1 = np.where(llc_face_connect[common] == face)[0][0] - edge_2 = np.where(llc_face_connect[common] == neighbor_face)[0][0] - new_edge_2 = np.where(llc_face_connect[neighbor_face] == common)[0][0] + edge_1 = np.where(LLC_FACE_CONNECT[face] == common)[0][0] + new_edge_1 = np.where(LLC_FACE_CONNECT[common] == face)[0][0] + edge_2 = np.where(LLC_FACE_CONNECT[common] == neighbor_face)[0][0] + new_edge_2 = np.where(LLC_FACE_CONNECT[neighbor_face] == common)[0][0] if (edge_1 in [0, 1] and new_edge_1 in [0, 1]) or ( edge_1 in [2, 3] and new_edge_1 in [2, 3] ): @@ -83,7 +83,7 @@ def llc_get_the_other_edge(face, edge): The compileable version of get_the_other_edge for llc grid. See Topology.get_the_other_edge for more detail. """ - face_connect = llc_face_connect + face_connect = LLC_FACE_CONNECT neighbor_face = face_connect[face, edge] if neighbor_face == 42: raise IndexError( @@ -234,7 +234,7 @@ def llc_get_uv_mask_from_face(faces): edge, new_edge = llc_mutual_direction( faces[0], faces[i], transitive=True ) - rot = np.pi - directions[edge] + directions[new_edge] + rot = np.pi - DIRECTIONS[edge] + DIRECTIONS[new_edge] # you can think of this as a rotation matrix UfromUvel[i] = np.cos(rot) UfromVvel[i] = np.sin(rot) @@ -442,7 +442,7 @@ def ind_moves(self, ind, moves, **kwarg): edge, new_edge = self.mutual_direction( face, ind[0], transitive=True ) - rot = (np.pi - directions[edge] + directions[new_edge]) % ( + rot = (np.pi - DIRECTIONS[edge] + DIRECTIONS[new_edge]) % ( np.pi * 2 ) if np.isclose(rot, 0):