Fix velocity advection stencil test and references

ci/benchmark.yml

@@ -4,8 +4,7 @@ include:
 .benchmark_model_stencils:
   stage: benchmark
   script:
-    # force execution of tests where validation is expected to fail, because the reason for failure is wrong numpy reference
-    - nox -s benchmark_model-3.10  -- --backend=$BACKEND --grid=$GRID --runxfail
+    - nox -s benchmark_model-3.10  -- --backend=$BACKEND --grid=$GRID
   parallel:
     matrix:
       - BACKEND: [gtfn_cpu, gtfn_gpu]

model/atmosphere/dycore/src/icon4py/model/atmosphere/dycore/stencils/fused_velocity_advection_stencil_8_to_13.py

@@ -57,7 +57,7 @@ def _fused_velocity_advection_stencil_8_to_13_predictor(
 
     w_concorr_c = where(
         nflatlev + 1 <= k < nlev,
-        _interpolate_to_half_levels_vp(interpolant=z_w_concorr_mc, wgtfac_c=wgtfac_c),
+        _interpolate_to_half_levels_vp(wgtfac_c=wgtfac_c, interpolant=z_w_concorr_mc),
         w_concorr_c,
     )
 

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

@@ -13,9 +13,21 @@
     compute_horizontal_advection_term_for_vertical_velocity,
 )
 from icon4py.model.common import dimension as dims
+from icon4py.model.common.grid import horizontal as h_grid, icon as icon_grid
+from icon4py.model.common.type_alias import vpfloat, wpfloat
 from icon4py.model.common.utils.data_allocation import random_field, zero_field
 from icon4py.model.testing.helpers import StencilTest
-from icon4py.model.common.type_alias import vpfloat, wpfloat
+
+
+def _horizontal_range(grid):
+    if isinstance(grid, icon_grid.IconGrid):
+        # For the ICON grid we use the proper domain bounds (otherwise we will run into non-protected skip values)
+        edge_domain = h_grid.domain(dims.EdgeDim)
+        return grid.start_index(edge_domain(h_grid.Zone.LATERAL_BOUNDARY_LEVEL_7)), grid.end_index(
+            edge_domain(h_grid.Zone.HALO)
+        )
+    else:
+        return 0, gtx.int32(grid.num_edges)
 
 
 def compute_horizontal_advection_term_for_vertical_velocity_numpy(
@@ -59,27 +71,26 @@ def reference(
         inv_primal_edge_length: np.array,
         tangent_orientation: np.array,
         z_w_v: np.array,
+        z_v_grad_w: np.array,
         **kwargs,
     ) -> dict:
-        z_v_grad_w = compute_horizontal_advection_term_for_vertical_velocity_numpy(
-            grid,
-            vn_ie,
-            inv_dual_edge_length,
-            w,
-            z_vt_ie,
-            inv_primal_edge_length,
-            tangent_orientation,
-            z_w_v,
-        )
+        horizontal_start, horizontal_end = _horizontal_range(grid)
+        z_v_grad_w[horizontal_start:horizontal_end, :] = (
+            compute_horizontal_advection_term_for_vertical_velocity_numpy(
+                grid,
+                vn_ie,
+                inv_dual_edge_length,
+                w,
+                z_vt_ie,
+                inv_primal_edge_length,
+                tangent_orientation,
+                z_w_v,
+            )
+        )[horizontal_start:horizontal_end, :]
         return dict(z_v_grad_w=z_v_grad_w)
 
     @pytest.fixture
     def input_data(self, grid):
-        if np.any(grid.connectivities[dims.E2CDim] == -1) or np.any(
-            grid.connectivities[dims.E2VDim] == -1
-        ):
-            pytest.xfail("Stencil does not support missing neighbors.")
-
         vn_ie = random_field(grid, dims.EdgeDim, dims.KDim, dtype=vpfloat)
         inv_dual_edge_length = random_field(grid, dims.EdgeDim, dtype=wpfloat)
         w = random_field(grid, dims.CellDim, dims.KDim, dtype=wpfloat)
@@ -89,6 +100,8 @@ def input_data(self, grid):
         z_w_v = random_field(grid, dims.VertexDim, dims.KDim, dtype=vpfloat)
         z_v_grad_w = zero_field(grid, dims.EdgeDim, dims.KDim, dtype=vpfloat)
 
+        horizontal_start, horizontal_end = _horizontal_range(grid)
+
         return dict(
             vn_ie=vn_ie,
             inv_dual_edge_length=inv_dual_edge_length,
@@ -98,8 +111,8 @@ def input_data(self, grid):
             tangent_orientation=tangent_orientation,
             z_w_v=z_w_v,
             z_v_grad_w=z_v_grad_w,
-            horizontal_start=0,
-            horizontal_end=gtx.int32(grid.num_edges),
+            horizontal_start=horizontal_start,
+            horizontal_end=horizontal_end,
             vertical_start=0,
             vertical_end=gtx.int32(grid.num_levels),
         )

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

@@ -19,7 +19,7 @@
 
 
 def copy_cell_kdim_field_to_vp_numpy(field: np.array) -> np.array:
-    field_copy = field
+    field_copy = field.copy()
     return field_copy
 
 

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

@@ -19,7 +19,7 @@
 
 
 def correct_contravariant_vertical_velocity_numpy(
-    w_concorr_c: np.array, z_w_con_c: np.array
+    z_w_con_c: np.array, w_concorr_c: np.array
 ) -> np.array:
     z_w_con_c = z_w_con_c - w_concorr_c
     return z_w_con_c
@@ -31,7 +31,7 @@ class TestCorrectContravariantVerticalVelocity(StencilTest):
 
     @staticmethod
     def reference(grid, w_concorr_c: np.array, z_w_con_c: np.array, **kwargs) -> dict:
-        z_w_con_c = correct_contravariant_vertical_velocity_numpy(w_concorr_c, z_w_con_c)
+        z_w_con_c = correct_contravariant_vertical_velocity_numpy(z_w_con_c, w_concorr_c)
         return dict(z_w_con_c=z_w_con_c)
 
     @pytest.fixture

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

@@ -13,8 +13,9 @@
     fused_velocity_advection_stencil_1_to_7,
 )
 from icon4py.model.common import dimension as dims
-from icon4py.model.testing.helpers import StencilTest
+from icon4py.model.common.grid import base as base_grid, horizontal as h_grid, icon as icon_grid
 from icon4py.model.common.utils import data_allocation as data_alloc
+from icon4py.model.testing.helpers import StencilTest
 
 from .test_compute_contravariant_correction import compute_contravariant_correction_numpy
 from .test_compute_horizontal_advection_term_for_vertical_velocity import (
@@ -106,7 +107,7 @@ def _fused_velocity_advection_stencil_1_to_6_numpy(
             z_w_concorr_me,
         )
 
-        return vt, vn_ie, z_kin_hor_e, z_w_concorr_me
+        return vt, vn_ie, z_vt_ie, z_kin_hor_e, z_w_concorr_me
 
     @classmethod
     def reference(
@@ -145,6 +146,7 @@ def reference(
             (
                 vt,
                 vn_ie,
+                z_vt_ie,
                 z_kin_hor_e,
                 z_w_concorr_me,
             ) = cls._fused_velocity_advection_stencil_1_to_6_numpy(
@@ -170,7 +172,7 @@ def reference(
 
         condition_mask = (lateral_boundary_7 <= edge) & (edge < halo_1) & (k_nlev < nlev)
 
-        z_v_w = mo_icon_interpolation_scalar_cells2verts_scalar_ri_dsl_numpy(grid, w, c_intp)
+        z_w_v = mo_icon_interpolation_scalar_cells2verts_scalar_ri_dsl_numpy(grid, w, c_intp)
 
         if not lvn_only:
             z_v_grad_w = np.where(
@@ -183,7 +185,7 @@ def reference(
                     z_vt_ie,
                     inv_primal_edge_length,
                     tangent_orientation,
-                    z_v_w,
+                    z_w_v,
                 ),
                 z_v_grad_w,
             )
@@ -197,11 +199,7 @@ def reference(
         )
 
     @pytest.fixture
-    def input_data(self, grid):
-        pytest.xfail(
-            "Verification of z_v_grad_w currently not working, because numpy version incorrect."
-        )
-
+    def input_data(self, grid: base_grid.BaseGrid):
         c_intp = data_alloc.random_field(grid, dims.VertexDim, dims.V2CDim)
         vn = data_alloc.random_field(grid, dims.EdgeDim, dims.KDim)
         rbf_vec_coeff_e = data_alloc.random_field(grid, dims.EdgeDim, dims.E2C2EDim)
@@ -232,8 +230,14 @@ def input_data(self, grid):
         istep = 1
         lvn_only = False
 
-        lateral_boundary_7 = 0
-        halo_1 = grid.num_edges
+        if isinstance(grid, icon_grid.IconGrid):
+            # For the ICON grid we use the proper domain bounds (otherwise we will run into non-protected skip values)
+            edge_domain = h_grid.domain(dims.EdgeDim)
+            lateral_boundary_7 = grid.start_index(edge_domain(h_grid.Zone.LATERAL_BOUNDARY_LEVEL_7))
+            halo_1 = grid.end_index(edge_domain(h_grid.Zone.HALO))
+        else:
+            lateral_boundary_7 = 0
+            halo_1 = grid.num_edges
 
         horizontal_start = 0
         horizontal_end = grid.num_edges

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

@@ -14,10 +14,7 @@
 )
 from icon4py.model.common import dimension as dims
 from icon4py.model.testing.helpers import StencilTest
-from icon4py.model.common.utils.data_allocation import (
-    random_field,
-    zero_field
-)
+from icon4py.model.common.utils.data_allocation import random_field, zero_field
 from icon4py.model.common.utils import data_allocation as data_alloc
 
 from .test_copy_cell_kdim_field_to_vp import copy_cell_kdim_field_to_vp_numpy
@@ -64,15 +61,13 @@ def reference(
         )
 
         if istep == 1:
-            z_w_concorr_mc = np.where(
-                (nflatlev <= k_nlev) & (k_nlev < nlev),
-                interpolate_to_cell_center_numpy(grid, z_w_concorr_me, e_bln_c_s),
-                z_w_concorr_mc,
-            )
+            z_w_concorr_mc = interpolate_to_cell_center_numpy(grid, z_w_concorr_me, e_bln_c_s)
 
             w_concorr_c = np.where(
                 (nflatlev + 1 <= k_nlev) & (k_nlev < nlev),
-                interpolate_to_half_levels_vp_numpy(grid, z_w_concorr_mc, wgtfac_c),
+                interpolate_to_half_levels_vp_numpy(
+                    grid, wgtfac_c=wgtfac_c, interpolant=z_w_concorr_mc
+                ),
                 w_concorr_c,
             )
 
@@ -96,9 +91,6 @@ def reference(
 
     @pytest.fixture
     def input_data(self, grid):
-        pytest.xfail(
-            "Verification of w_concorr_c currently not working, because numpy version is incorrect."
-        )
         z_kin_hor_e = random_field(grid, dims.EdgeDim, dims.KDim)
         e_bln_c_s = random_field(grid, dims.CellDim, dims.C2EDim)
         z_ekinh = zero_field(grid, dims.CellDim, dims.KDim)

model/testing/src/icon4py/model/testing/helpers.py

@@ -21,7 +21,6 @@
 from icon4py.model.common.utils import data_allocation as data_alloc
 
 
-
 try:
     import pytest_benchmark
 except ModuleNotFoundError:
@@ -96,11 +95,12 @@ def _test_validation(self, grid, backend, input_data):
             else (out, (slice(None),), (slice(None),))
         )
 
-        assert np.allclose(
+        np.testing.assert_allclose(
             input_data[name].asnumpy()[gtslice],
             reference_outputs[name][refslice],
             equal_nan=True,
-        ), f"Validation failed for '{name}'"
+            err_msg=f"Validation failed for '{name}'",
+        )
 
 
 if pytest_benchmark: