Rbf interpolation stencils improvement

model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/solve_nonhydro.py

@@ -404,7 +404,10 @@ def _validate(self):
             raise NotImplementedError("divdamp_order can only be 24")
 
         if self.divdamp_type == DivergenceDampingType.TWO_DIMENSIONAL:
-            raise NotImplementedError("`DivergenceDampingType.TWO_DIMENSIONAL` (2) is not yet implemented")
+            raise NotImplementedError(
+                "`DivergenceDampingType.TWO_DIMENSIONAL` (2) is not yet implemented"
+            )
+
 
 class NonHydrostaticParams:
     """Calculates derived quantities depending on the NonHydrostaticConfig."""

model/atmosphere/dycore/tests/dycore_stencil_tests/test_compute_advection_in_horizontal_momentum_equation.py

@@ -64,7 +64,9 @@ def reference(
         normal_wind_advective_tendency_cp = normal_wind_advective_tendency.copy()
         k = np.arange(nlev)
 
-        upward_vorticity_at_vertices = mo_math_divrot_rot_vertex_ri_dsl_numpy(connectivities, vn, geofac_rot)
+        upward_vorticity_at_vertices = mo_math_divrot_rot_vertex_ri_dsl_numpy(
+            connectivities, vn, geofac_rot
+        )
 
         normal_wind_advective_tendency = compute_advective_normal_wind_tendency_numpy(
             connectivities,
@@ -79,7 +81,7 @@ def reference(
             vn_on_half_levels,
             ddqz_z_full_e,
         )
-        
+
         condition = (np.maximum(3, nrdmax - 2) - 1 <= k) & (k < nlev - 4)
         normal_wind_advective_tendency_extra_diffu = (
             add_extra_diffusion_for_normal_wind_tendency_approaching_cfl_numpy(

model/common/src/icon4py/model/common/grid/geometry.py

@@ -464,6 +464,34 @@ def _register_computed_fields(self):
             pairs=(("u_cell_1", "u_cell_2"), ("v_cell_1", "v_cell_2")),
         )
         self.register_provider(tangent_cell_wrapper)
+        cartesian_edge_centers = factory.FieldOperatorProvider(
+            func=math_helpers.geographical_to_cartesian_on_edges.with_backend(self.backend),
+            domain=(dims.EdgeDim,),
+            fields={
+                attrs.EDGE_CENTER_X: attrs.EDGE_CENTER_X,
+                attrs.EDGE_CENTER_Y: attrs.EDGE_CENTER_Y,
+                attrs.EDGE_CENTER_Z: attrs.EDGE_CENTER_Z,
+            },
+            deps={
+                "lat": attrs.EDGE_LAT,
+                "lon": attrs.EDGE_LON,
+            },
+        )
+        self.register_provider(cartesian_edge_centers)
+        cartesian_cell_centers = factory.FieldOperatorProvider(
+            func=math_helpers.geographical_to_cartesian_on_cells.with_backend(self.backend),
+            domain=(dims.CellDim,),
+            fields={
+                attrs.CELL_CENTER_X: attrs.CELL_CENTER_X,
+                attrs.CELL_CENTER_Y: attrs.CELL_CENTER_Y,
+                attrs.CELL_CENTER_Z: attrs.CELL_CENTER_Z,
+            },
+            deps={
+                "lat": attrs.CELL_LAT,
+                "lon": attrs.CELL_LON,
+            },
+        )
+        self.register_provider(cartesian_cell_centers)
 
     def _inverse_field_provider(self, field_name: str):
         meta = attrs.metadata_for_inverse(attrs.attrs[field_name])

model/common/src/icon4py/model/common/grid/geometry_attributes.py

@@ -27,6 +27,11 @@
 
 CELL_LON: Final[str] = "grid_longitude_of_cell_center"
 CELL_LAT: Final[str] = "grid_latitude_of_cell_center"
+CELL_CENTER_X: Final[str] = "x_component_of_cell_center"
+CELL_CENTER_Y: Final[str] = "y_component_of_cell_center"
+
+CELL_CENTER_Z: Final[str] = "z_component_of_cell_center"
+
 CELL_AREA: Final[str] = "cell_area"
 EDGE_AREA: Final[str] = "edge_area"
 DUAL_AREA: Final[str] = "dual_area"
@@ -62,6 +67,10 @@
 EDGE_NORMAL_VERTEX_V: Final[str] = "northward_component_of_edge_normal_on_vertex"
 EDGE_NORMAL_CELL_U: Final[str] = "eastward_component_of_edge_normal_on_cell"
 EDGE_NORMAL_CELL_V: Final[str] = "northward_component_of_edge_normal_on_cell"
+EDGE_CENTER_X: Final[str] = "x_coordinate_of_edge_center"
+EDGE_CENTER_Y: Final[str] = "y_coordinate_of_edge_center"
+EDGE_CENTER_Z: Final[str] = "z_coordinate_of_edge_center"
+
 
 attrs: dict[str, model.FieldMetaData] = {
     CELL_LAT: dict(
@@ -169,6 +178,30 @@
         icon_var_name="t_grid_cells%area",
         dtype=ta.wpfloat,
     ),
+    CELL_CENTER_X: dict(
+        standard_name=CELL_CENTER_X,
+        long_name="x component of cell center",
+        units="",
+        dims=(dims.CellDim,),
+        icon_var_name="t_grid_cells%%cartesian_center%x(1)",
+        dtype=ta.wpfloat,
+    ),
+    CELL_CENTER_Y: dict(
+        standard_name=CELL_CENTER_Y,
+        long_name="y-component of cell center",
+        units="",
+        dims=(dims.CellDim,),
+        icon_var_name="t_grid_cells%%cartesian_center%x(2)",
+        dtype=ta.wpfloat,
+    ),
+    CELL_CENTER_Z: dict(
+        standard_name=CELL_CENTER_Z,
+        long_name="z-component of cell center",
+        units="",
+        dims=(dims.CellDim,),
+        icon_var_name="t_grid_cells%%cartesian_center%x(3)",
+        dtype=ta.wpfloat,
+    ),
     DUAL_AREA: dict(
         standard_name=DUAL_AREA,
         long_name="area of the dual grid cell (hexagon cell)",
@@ -345,6 +378,30 @@
         icon_var_name="ptr_patch%edges%dual_normal%v2",
         dtype=ta.wpfloat,
     ),
+    EDGE_CENTER_X: dict(
+        standard_name=EDGE_CENTER_X,
+        long_name="x component of cartesian coordinates of edge centers",
+        units="1",
+        dims=(dims.EdgeDim,),
+        icon_var_name="t_grid_edges%cartesian_center%x(1)",
+        dtype=ta.wpfloat,
+    ),
+    EDGE_CENTER_Y: dict(
+        standard_name=EDGE_CENTER_Y,
+        long_name="x component of cartesian coordinates of edge centers",
+        units="1",
+        dims=(dims.EdgeDim,),
+        icon_var_name="t_grid_edges%cartesian_center%x(2)",
+        dtype=ta.wpfloat,
+    ),
+    EDGE_CENTER_Z: dict(
+        standard_name=EDGE_CENTER_Z,
+        long_name="z component of cartesian coordinates of edge centers",
+        units="1",
+        dims=(dims.EdgeDim,),
+        icon_var_name="t_grid_edges%cartesian_center%x(3)",
+        dtype=ta.wpfloat,
+    ),
 }
 
 

model/common/src/icon4py/model/common/grid/geometry_stencils.py

@@ -15,7 +15,7 @@
 from icon4py.model.common import dimension as dims, field_type_aliases as fa, type_alias as ta
 from icon4py.model.common.dimension import E2C, E2C2V, E2V, EdgeDim
 from icon4py.model.common.math.helpers import (
-    arc_length,
+    arc_length_on_edges,
     cross_product_on_edges,
     geographical_to_cartesian_on_edges,
     geographical_to_cartesian_on_vertex,
@@ -362,7 +362,7 @@ def cell_center_arc_distance(
     x0, y0, z0 = geographical_to_cartesian_on_edges(lat_neighbor_0, lon_neighbor_0)
     x1, y1, z1 = geographical_to_cartesian_on_edges(lat_neighbor_1, lon_neighbor_1)
     # (xi, yi, zi) are normalized by construction
-    arc = arc_length(x0, x1, y0, y1, z0, z1, radius)
+    arc = arc_length_on_edges(x0, x1, y0, y1, z0, z1, radius)
     return arc
 
 
@@ -406,7 +406,7 @@ def arc_distance_of_far_edges_in_diamond(
     z3 = z(E2C2V[3])
     # (xi, yi, zi) are normalized by construction
 
-    far_vertex_vertex_length = arc_length(x2, x3, y2, y3, z2, z3, radius)
+    far_vertex_vertex_length = arc_length_on_edges(x2, x3, y2, y3, z2, z3, radius)
     return far_vertex_vertex_length
 
 
@@ -440,7 +440,7 @@ def edge_length(
     z1 = z(E2V[1])
     # (xi, yi, zi) are normalized by construction
 
-    length = arc_length(x0, x1, y0, y1, z0, z1, radius)
+    length = arc_length_on_edges(x0, x1, y0, y1, z0, z1, radius)
     return length