Skip to content

Commit c5f59a0

Browse files
Separating interpolators into x and ux files
1 parent d056ed5 commit c5f59a0

3 files changed

Lines changed: 111 additions & 81 deletions

File tree

Lines changed: 37 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,37 @@
1+
from .uxinterpolators import (
2+
Ux_Velocity,
3+
UxPiecewiseConstantFace,
4+
UxPiecewiseLinearNode,
5+
)
6+
from .xinterpolators import (
7+
CGrid_Tracer,
8+
CGrid_Velocity,
9+
XConstantField,
10+
XFreeslip,
11+
XLinear,
12+
XLinear_Velocity,
13+
XLinearInvdistLandTracer,
14+
XNearest,
15+
XPartialslip,
16+
ZeroInterpolator,
17+
ZeroInterpolator_Vector,
18+
)
19+
20+
__all__ = [ # noqa: RUF022
21+
# xinterpolators
22+
"CGrid_Tracer",
23+
"CGrid_Velocity",
24+
"XConstantField",
25+
"XFreeslip",
26+
"XLinear",
27+
"XLinearInvdistLandTracer",
28+
"XLinear_Velocity",
29+
"XNearest",
30+
"XPartialslip",
31+
"ZeroInterpolator",
32+
"ZeroInterpolator_Vector",
33+
# uxinterpolators
34+
"UxPiecewiseConstantFace",
35+
"UxPiecewiseLinearNode",
36+
"Ux_Velocity",
37+
]
Lines changed: 73 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,73 @@
1+
"""Collection of pre-built interpolation kernels for unstructured grids."""
2+
3+
from __future__ import annotations
4+
5+
from typing import TYPE_CHECKING
6+
7+
import numpy as np
8+
9+
if TYPE_CHECKING:
10+
from parcels._core.field import Field, VectorField
11+
from parcels._core.uxgrid import _UXGRID_AXES
12+
13+
14+
def UxPiecewiseConstantFace(
15+
particle_positions: dict[str, float | np.ndarray],
16+
grid_positions: dict[_UXGRID_AXES, dict[str, int | float | np.ndarray]],
17+
field: Field,
18+
):
19+
"""
20+
Piecewise constant interpolation kernel for face registered data.
21+
This interpolation method is appropriate for fields that are
22+
face registered, such as u,v in FESOM.
23+
"""
24+
return field.data.values[
25+
grid_positions["T"]["index"], grid_positions["Z"]["index"], grid_positions["FACE"]["index"]
26+
]
27+
28+
29+
def UxPiecewiseLinearNode(
30+
particle_positions: dict[str, float | np.ndarray],
31+
grid_positions: dict[_UXGRID_AXES, dict[str, int | float | np.ndarray]],
32+
field: Field,
33+
):
34+
"""
35+
Piecewise linear interpolation kernel for node registered data located at vertical interface levels.
36+
This interpolation method is appropriate for fields that are node registered such as the vertical
37+
velocity W in FESOM2. Effectively, it applies barycentric interpolation in the lateral direction
38+
and piecewise linear interpolation in the vertical direction.
39+
"""
40+
ti = grid_positions["T"]["index"]
41+
zi, fi = grid_positions["Z"]["index"], grid_positions["FACE"]["index"]
42+
z = particle_positions["z"]
43+
bcoords = grid_positions["FACE"]["bcoord"]
44+
node_ids = field.grid.uxgrid.face_node_connectivity[fi, :].values
45+
# The zi refers to the vertical layer index. The field in this routine are assumed to be defined at the vertical interface levels.
46+
# For interface zi, the interface indices are [zi, zi+1], so we need to use the values at zi and zi+1.
47+
# First, do barycentric interpolation in the lateral direction for each interface level
48+
fzk = np.sum(field.data.values[ti[:, None], zi[:, None], node_ids] * bcoords, axis=-1)
49+
fzkp1 = np.sum(field.data.values[ti[:, None], zi[:, None] + 1, node_ids] * bcoords, axis=-1)
50+
51+
# Then, do piecewise linear interpolation in the vertical direction
52+
zk = field.grid.z.values[zi]
53+
zkp1 = field.grid.z.values[zi + 1]
54+
return (fzk * (zkp1 - z) + fzkp1 * (z - zk)) / (zkp1 - zk) # Linear interpolation in the vertical direction
55+
56+
57+
def Ux_Velocity(
58+
particle_positions: dict[str, float | np.ndarray],
59+
grid_positions: dict[_UXGRID_AXES, dict[str, int | float | np.ndarray]],
60+
vectorfield: VectorField,
61+
):
62+
"""Interpolation kernel for Vectorfields of velocity on a UxGrid."""
63+
u = vectorfield.U._interp_method(particle_positions, grid_positions, vectorfield.U)
64+
v = vectorfield.V._interp_method(particle_positions, grid_positions, vectorfield.V)
65+
if vectorfield.grid._mesh == "spherical":
66+
u /= 1852 * 60 * np.cos(np.deg2rad(particle_positions["lat"]))
67+
v /= 1852 * 60
68+
69+
if "3D" in vectorfield.vector_type:
70+
w = vectorfield.W._interp_method(particle_positions, grid_positions, vectorfield.W)
71+
else:
72+
w = 0.0
73+
return u, v, w

src/parcels/interpolators/xinterpolators.py

Lines changed: 1 addition & 81 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
"""Collection of pre-built interpolation kernels."""
1+
"""Collection of pre-built interpolation kernels for structured grids."""
22

33
from __future__ import annotations
44

@@ -12,26 +12,8 @@
1212

1313
if TYPE_CHECKING:
1414
from parcels._core.field import Field, VectorField
15-
from parcels._core.uxgrid import _UXGRID_AXES
1615
from parcels._core.xgrid import _XGRID_AXES
1716

18-
__all__ = [
19-
"CGrid_Tracer",
20-
"CGrid_Velocity",
21-
"UxPiecewiseConstantFace",
22-
"UxPiecewiseLinearNode",
23-
"Ux_Velocity",
24-
"XConstantField",
25-
"XFreeslip",
26-
"XLinear",
27-
"XLinearInvdistLandTracer",
28-
"XLinear_Velocity",
29-
"XNearest",
30-
"XPartialslip",
31-
"ZeroInterpolator",
32-
"ZeroInterpolator_Vector",
33-
]
34-
3517

3618
def ZeroInterpolator(
3719
particle_positions: dict[str, float | np.ndarray],
@@ -666,65 +648,3 @@ def XLinearInvdistLandTracer(
666648
values[exact_particles] = exact_vals[exact_particles]
667649

668650
return values.compute() if is_dask_collection(values) else values
669-
670-
671-
def UxPiecewiseConstantFace(
672-
particle_positions: dict[str, float | np.ndarray],
673-
grid_positions: dict[_UXGRID_AXES, dict[str, int | float | np.ndarray]],
674-
field: Field,
675-
):
676-
"""
677-
Piecewise constant interpolation kernel for face registered data.
678-
This interpolation method is appropriate for fields that are
679-
face registered, such as u,v in FESOM.
680-
"""
681-
return field.data.values[
682-
grid_positions["T"]["index"], grid_positions["Z"]["index"], grid_positions["FACE"]["index"]
683-
]
684-
685-
686-
def UxPiecewiseLinearNode(
687-
particle_positions: dict[str, float | np.ndarray],
688-
grid_positions: dict[_UXGRID_AXES, dict[str, int | float | np.ndarray]],
689-
field: Field,
690-
):
691-
"""
692-
Piecewise linear interpolation kernel for node registered data located at vertical interface levels.
693-
This interpolation method is appropriate for fields that are node registered such as the vertical
694-
velocity W in FESOM2. Effectively, it applies barycentric interpolation in the lateral direction
695-
and piecewise linear interpolation in the vertical direction.
696-
"""
697-
ti = grid_positions["T"]["index"]
698-
zi, fi = grid_positions["Z"]["index"], grid_positions["FACE"]["index"]
699-
z = particle_positions["z"]
700-
bcoords = grid_positions["FACE"]["bcoord"]
701-
node_ids = field.grid.uxgrid.face_node_connectivity[fi, :].values
702-
# The zi refers to the vertical layer index. The field in this routine are assumed to be defined at the vertical interface levels.
703-
# For interface zi, the interface indices are [zi, zi+1], so we need to use the values at zi and zi+1.
704-
# First, do barycentric interpolation in the lateral direction for each interface level
705-
fzk = np.sum(field.data.values[ti[:, None], zi[:, None], node_ids] * bcoords, axis=-1)
706-
fzkp1 = np.sum(field.data.values[ti[:, None], zi[:, None] + 1, node_ids] * bcoords, axis=-1)
707-
708-
# Then, do piecewise linear interpolation in the vertical direction
709-
zk = field.grid.z.values[zi]
710-
zkp1 = field.grid.z.values[zi + 1]
711-
return (fzk * (zkp1 - z) + fzkp1 * (z - zk)) / (zkp1 - zk) # Linear interpolation in the vertical direction
712-
713-
714-
def Ux_Velocity(
715-
particle_positions: dict[str, float | np.ndarray],
716-
grid_positions: dict[_UXGRID_AXES, dict[str, int | float | np.ndarray]],
717-
vectorfield: VectorField,
718-
):
719-
"""Interpolation kernel for Vectorfields of velocity on a UxGrid."""
720-
u = vectorfield.U._interp_method(particle_positions, grid_positions, vectorfield.U)
721-
v = vectorfield.V._interp_method(particle_positions, grid_positions, vectorfield.V)
722-
if vectorfield.grid._mesh == "spherical":
723-
u /= 1852 * 60 * np.cos(np.deg2rad(particle_positions["lat"]))
724-
v /= 1852 * 60
725-
726-
if "3D" in vectorfield.vector_type:
727-
w = vectorfield.W._interp_method(particle_positions, grid_positions, vectorfield.W)
728-
else:
729-
w = 0.0
730-
return u, v, w

0 commit comments

Comments
 (0)