From e10372a13a9c598bee0435e787414759fe1a5536 Mon Sep 17 00:00:00 2001 From: Wenrui Jiang Date: Sun, 2 Jul 2023 15:29:31 -0400 Subject: [PATCH 01/11] visa and stamp system for notetaking --- seaduck/lagrangian.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/seaduck/lagrangian.py b/seaduck/lagrangian.py index 4b1eb656..0fefa4dd 100644 --- a/seaduck/lagrangian.py +++ b/seaduck/lagrangian.py @@ -415,7 +415,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=0): """Record raw data into list of lists. This method is only called in save_raw = True particles. @@ -459,6 +459,7 @@ 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]) + self.vslist[ifull].append(stamp) def empty_lists(self): """Empty/Create the lists. @@ -489,6 +490,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. @@ -871,7 +873,7 @@ def to_next_stop(self, t_stop): if self.save_raw: # record the moment just before crossing the wall # or the moment reaching destination. - self.note_taking(int_todo) + self.note_taking(int_todo, stamp=0) sub.cross_cell_wall(tend) if self.transport: @@ -888,7 +890,7 @@ def to_next_stop(self, t_stop): break if self.save_raw: # record those who cross the wall - self.note_taking(int_todo) + self.note_taking(int_todo, stamp=1) if i == self.max_iteration - 1: warnings.warn("maximum iteration count reached") @@ -971,7 +973,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=2) self.to_next_stop(tl) if update[i]: if self.too_large: # pragma: no cover From f1fe6742679381016da5a78b4cf5db8ec51154ab Mon Sep 17 00:00:00 2001 From: Wenrui Jiang Date: Sun, 2 Jul 2023 15:31:48 -0400 Subject: [PATCH 02/11] It seems like some of the pragma no cover funcs are actually not useful even for post process, thus removed --- seaduck/eulerian.py | 45 --------------------------------------------- 1 file changed, 45 deletions(-) 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 ): From 36af6aee8a0ff9b92a0fbb1c8a797117e01d4fd5 Mon Sep 17 00:00:00 2001 From: Wenrui Jiang Date: Thu, 6 Jul 2023 11:44:31 -0400 Subject: [PATCH 03/11] should be subset doing notetaking --- seaduck/lagrangian.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/seaduck/lagrangian.py b/seaduck/lagrangian.py index 0fefa4dd..e2644f99 100644 --- a/seaduck/lagrangian.py +++ b/seaduck/lagrangian.py @@ -873,7 +873,7 @@ def to_next_stop(self, t_stop): if self.save_raw: # record the moment just before crossing the wall # or the moment reaching destination. - self.note_taking(int_todo, stamp=0) + sub.note_taking(int_todo, stamp=0) sub.cross_cell_wall(tend) if self.transport: @@ -890,7 +890,7 @@ def to_next_stop(self, t_stop): break if self.save_raw: # record those who cross the wall - self.note_taking(int_todo, stamp=1) + sub.note_taking(int_todo, stamp=1) if i == self.max_iteration - 1: warnings.warn("maximum iteration count reached") From b58f2f8b0f946e2aac319d954211ebad58a1b11a Mon Sep 17 00:00:00 2001 From: Wenrui Jiang Date: Thu, 6 Jul 2023 19:14:52 -0400 Subject: [PATCH 04/11] add a section in AVISO and debug the note_taking after cross_cell_Wall --- docs/notebook/AVISO.ipynb | 262 +++++++++++++++++++++++++++++++++++++- seaduck/lagrangian.py | 12 +- 2 files changed, 268 insertions(+), 6 deletions(-) diff --git a/docs/notebook/AVISO.ipynb b/docs/notebook/AVISO.ipynb index 59b2d9b8..26de177f 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,230 @@ "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 a special feature of this package: users can release particles once at a time without hurting performance.\n", + "\n", + "Here is how we are going to release the particles. We are going to release them all at 64W and going back and forth between 55S and 64S. This is as if there is a ship is commuting in drake passage while releasing particles." + ] + }, + { + "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." + ] + }, + { + "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 +719,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.4" + "version": "3.11.3" } }, "nbformat": 4, diff --git a/seaduck/lagrangian.py b/seaduck/lagrangian.py index e2644f99..c86cdf23 100644 --- a/seaduck/lagrangian.py +++ b/seaduck/lagrangian.py @@ -415,7 +415,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, stamp=0): + 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 +437,11 @@ def note_taking(self, subset_index=None, stamp=0): 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]) @@ -858,6 +863,7 @@ 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 @@ -867,7 +873,6 @@ def to_next_stop(self, t_stop): 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: @@ -885,9 +890,10 @@ 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 sub.note_taking(int_todo, stamp=1) From 179b115a93e9fc34ef36da506dde11334add5e38 Mon Sep 17 00:00:00 2001 From: Wenrui Jiang Date: Thu, 6 Jul 2023 22:37:16 -0400 Subject: [PATCH 05/11] improve stamp and deal with weird underflow prob --- seaduck/lagrangian.py | 24 +++++++++++++----------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/seaduck/lagrangian.py b/seaduck/lagrangian.py index c86cdf23..101d09aa 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-14 + tl[cannot_left] = -sign[cannot_left] + cannot_right = ur * tf <= 1e-14 + tr[cannot_right] = -sign[cannot_right] ts.append(tl) ts.append(tr) return ts @@ -464,7 +465,10 @@ def note_taking(self, subset_index=None, stamp=-1): self.xxlist[ifull].append(self.lon[isub]) self.yylist[ifull].append(self.lat[isub]) self.zzlist[ifull].append(self.dep[isub]) - self.vslist[ifull].append(stamp) + if isinstance(stamp,int): + self.vslist[ifull].append(stamp) + else: + self.vslist[ifull].append(stamp[isub]) def empty_lists(self): """Empty/Create the lists. @@ -868,17 +872,15 @@ def to_next_stop(self, t_stop): 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.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. - sub.note_taking(int_todo, stamp=0) + sub.note_taking(int_todo, stamp=tend) sub.cross_cell_wall(tend) if self.transport: @@ -896,7 +898,7 @@ def to_next_stop(self, t_stop): tf_used = tf_used[bool_todo] if self.save_raw: # record those who cross the wall - sub.note_taking(int_todo, stamp=1) + sub.note_taking(int_todo, stamp=7) if i == self.max_iteration - 1: warnings.warn("maximum iteration count reached") @@ -979,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(stamp=2) + self.note_taking(stamp=15) self.to_next_stop(tl) if update[i]: if self.too_large: # pragma: no cover From 654ebb0f4b31d68736faf3d5d2a2ce370b953a57 Mon Sep 17 00:00:00 2001 From: Wenrui Jiang Date: Thu, 6 Jul 2023 22:39:35 -0400 Subject: [PATCH 06/11] style fix, confirming that this result in particles unstuck --- seaduck/lagrangian.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/seaduck/lagrangian.py b/seaduck/lagrangian.py index 101d09aa..1740a2b1 100644 --- a/seaduck/lagrangian.py +++ b/seaduck/lagrangian.py @@ -465,7 +465,7 @@ def note_taking(self, subset_index=None, stamp=-1): 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): + if isinstance(stamp, int): self.vslist[ifull].append(stamp) else: self.vslist[ifull].append(stamp[isub]) From a43064e986a2df2d44fcbc55d59c51e258b4e0db Mon Sep 17 00:00:00 2001 From: Wenrui Jiang Date: Sat, 8 Jul 2023 17:32:31 -0400 Subject: [PATCH 07/11] change the words in AVISO --- docs/notebook/AVISO.ipynb | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/docs/notebook/AVISO.ipynb b/docs/notebook/AVISO.ipynb index 26de177f..2eb8ac3e 100644 --- a/docs/notebook/AVISO.ipynb +++ b/docs/notebook/AVISO.ipynb @@ -490,9 +490,9 @@ "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 a special feature of this package: users can release particles once at a time without hurting performance.\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", - "Here is how we are going to release the particles. We are going to release them all at 64W and going back and forth between 55S and 64S. This is as if there is a ship is commuting in drake passage while releasing particles." + "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." ] }, { @@ -534,7 +534,9 @@ "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." + "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." ] }, { From 78b62349c17ecbae0b0512919b51d48a06d9c27b Mon Sep 17 00:00:00 2001 From: Wenrui Jiang Date: Sun, 9 Jul 2023 17:40:58 -0400 Subject: [PATCH 08/11] seems like the previous fix was too lenient on underflow --- seaduck/lagrangian.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/seaduck/lagrangian.py b/seaduck/lagrangian.py index 1740a2b1..2c80817a 100644 --- a/seaduck/lagrangian.py +++ b/seaduck/lagrangian.py @@ -103,9 +103,9 @@ def time2wall(pos_list, u_list, du_list, tf): 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]) sign = np.sign(tf) - cannot_left = -ul * sign <= 1e-14 + cannot_left = -ul * sign <= 1e-11 # aroung 3000 years tl[cannot_left] = -sign[cannot_left] - cannot_right = ur * tf <= 1e-14 + cannot_right = ur * tf <= 1e-11 tr[cannot_right] = -sign[cannot_right] ts.append(tl) ts.append(tr) From 550a000e72c5e17726378b0b09bde539204ce5a7 Mon Sep 17 00:00:00 2001 From: Wenrui Jiang Date: Tue, 11 Jul 2023 21:50:07 -0400 Subject: [PATCH 09/11] change name of global variable --- seaduck/kernel_weight.py | 23 +++++++++++------------ seaduck/lagrangian.py | 4 ++-- seaduck/topology.py | 32 ++++++++++++++++---------------- 3 files changed, 29 insertions(+), 30 deletions(-) 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 2c80817a..3a5e0e54 100644 --- a/seaduck/lagrangian.py +++ b/seaduck/lagrangian.py @@ -103,9 +103,9 @@ def time2wall(pos_list, u_list, du_list, tf): 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]) sign = np.sign(tf) - cannot_left = -ul * sign <= 1e-11 # aroung 3000 years + cannot_left = -ul * sign <= 1e-12 # aroung 30000 years tl[cannot_left] = -sign[cannot_left] - cannot_right = ur * tf <= 1e-11 + cannot_right = ur * tf <= 1e-12 tr[cannot_right] = -sign[cannot_right] ts.append(tl) ts.append(tr) diff --git a/seaduck/topology.py b/seaduck/topology.py index bb37b737..ee233786 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,20 +35,20 @@ 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: - found = connected_to_face[0][0] in [0, 1, 2, 3] and connected_to_neigh[0][ + found = connected_to_face[0][0] in LEGAL_TENDS and connected_to_neigh[0][ 0 - ] in [0, 1, 2, 3] + ] in LEGAL_TENDS if found: 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): From 5e5b46d389d11e11a31b60839bc6fc7aac2f1585 Mon Sep 17 00:00:00 2001 From: Wenrui Jiang Date: Tue, 11 Jul 2023 22:00:02 -0400 Subject: [PATCH 10/11] topology somehow breaks stuff --- seaduck/lagrangian.py | 2 +- seaduck/topology.py | 32 ++++++++++++++++---------------- 2 files changed, 17 insertions(+), 17 deletions(-) diff --git a/seaduck/lagrangian.py b/seaduck/lagrangian.py index 3a5e0e54..45886a61 100644 --- a/seaduck/lagrangian.py +++ b/seaduck/lagrangian.py @@ -103,7 +103,7 @@ def time2wall(pos_list, u_list, du_list, tf): 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]) sign = np.sign(tf) - cannot_left = -ul * sign <= 1e-12 # aroung 30000 years + 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] diff --git a/seaduck/topology.py b/seaduck/topology.py index ee233786..bb37b737 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,20 +35,20 @@ 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: - found = connected_to_face[0][0] in LEGAL_TENDS and connected_to_neigh[0][ + found = connected_to_face[0][0] in [0, 1, 2, 3] and connected_to_neigh[0][ 0 - ] in LEGAL_TENDS + ] in [0, 1, 2, 3] if found: 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): From 68dd6c92138e52edee9f52963485107dfcbe42d3 Mon Sep 17 00:00:00 2001 From: Wenrui Jiang Date: Tue, 11 Jul 2023 22:08:08 -0400 Subject: [PATCH 11/11] ready to push after fixing topology --- seaduck/topology.py | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) 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):