added the first test for mo_lib_gradients

test/c/CMakeLists.txt

@@ -35,6 +35,7 @@ set(SOURCES
   test_horizontal_div.cpp
   test_horizontal_recon.cpp
   test_horizontal_rot.cpp
+  test_horizontal_gradients.cpp
   test_tdma_solver.cpp
   test_interpolation_vector.cpp
   test_intp_rbf.cpp

test/c/test_horizontal_gradients.cpp

+// ICON
+//
+// ---------------------------------------------------------------
+// Copyright (C) 2004-2025, DWD, MPI-M, DKRZ, KIT, ETH, MeteoSwiss
+// Contact information: icon-model.org
+//
+// See AUTHORS.TXT for a list of authors
+// See LICENSES/ for license information
+// SPDX-License-Identifier: BSD-3-Clause
+// ---------------------------------------------------------------
+
+#include <iostream>
+#include <random>
+#include <vector>
+
+#include <Kokkos_Core.hpp>
+#include <gtest/gtest.h>
+#include <dim_helper.hpp>
+#include <horizontal/mo_lib_gradients.hpp>
+#include <support/mo_lib_loopindices.hpp>
+
+#include <gtest/gtest.h>
+#include <random>
+#include <vector>
+
+template <typename ValueType>
+class GradientNormalTest : public ::testing::Test {
+protected:
+  static constexpr int nproma = 3;   // inner loop length
+  static constexpr int nlev = 2;     // number of vertical levels
+  static constexpr int nblks_c = 1;  // number of cell blocks
+  static constexpr int nblks_e = 1;  // number of edge blocks
+
+  int i_startblk = 0;
+  int i_endblk = nblks_e;      // Test blocks [0 .. nblks_e-1]
+  int i_startidx_in = 0;
+  int i_endidx_in = nproma;    // Full range: 0 .. nproma-1
+  std::vector<int> slev;
+  std::vector<int> elev;
+  bool lacc = false;           // Not using ACC-specific behavior.
+
+  std::vector<ValueType> psi_c;
+  std::vector<int> edge_cell_idx;
+  std::vector<int> edge_cell_blk;
+  std::vector<ValueType> inv_dual_edge_length;
+  std::vector<ValueType> grad_norm_psi_e;
+
+  GradientNormalTest() {
+    slev.resize(1, 0);
+    elev.resize(1, nlev);      // Full vertical range (0 .. nlev-1)
+
+    psi_c.resize(dim_combine(nproma, nlev, nblks_c));
+    edge_cell_idx.resize(dim_combine(nproma, nblks_e, 2));
+    edge_cell_blk.resize(dim_combine(nproma, nblks_e, 2));
+    inv_dual_edge_length.resize(dim_combine(nproma, nblks_e));
+    grad_norm_psi_e.resize(dim_combine(nproma, nlev, nblks_e));
+  }
+};
+
+template <typename ValueType>
+class GradientNormal3DTest : public GradientNormalTest<ValueType> {};
+
+// Define value types to test (e.g., float and double)
+using ValueTypes = ::testing::Types<float, double>;
+TYPED_TEST_SUITE(GradientNormal3DTest, ValueTypes);
+
+TYPED_TEST(GradientNormal3DTest, TestSpecific) {
+  constexpr int nproma = this->nproma;
+  constexpr int nlev = this->nlev;
+  constexpr int nblks_c = this->nblks_c;
+  constexpr int nblks_e = this->nblks_e;
+
+  const auto &psi_c_at = at<nproma, nlev, nblks_c>;
+  const auto &edge_cell_at = at<nproma, nblks_e, 2>;
+  const auto &inv_edge_at = at<nproma, nblks_e>;
+  const auto &grad_psi_at = at<nproma, nlev, nblks_e>;
+
+  // Initialization with specific values
+  for (int i = 0; i < nproma; ++i) {
+    for (int k = 0; k < nlev; ++k) {
+      this->psi_c[psi_c_at(i, k, 0)] = (i + 1) * (k + 1) * 2.0; // Simple pattern
+    }
+
+    // Set edge cell indices to point to specific cells
+    this->edge_cell_idx[edge_cell_at(i, 0, 0)] = i;
+    this->edge_cell_idx[edge_cell_at(i, 0, 1)] = (i + 1) % nproma;
+
+    // All cells are in the same block for this test
+    this->edge_cell_blk[edge_cell_at(i, 0, 0)] = 0;
+    this->edge_cell_blk[edge_cell_at(i, 0, 1)] = 0;
+
+    // Set inverse dual edge length
+    this->inv_dual_edge_length[inv_edge_at(i, 0)] = 0.5;
+
+    // Initialize grad_norm_psi_e to zero
+    for (int k = 0; k < nlev; ++k) {
+      this->grad_norm_psi_e[grad_psi_at(i, k, 0)] = 0.0;
+    }
+  }
+
+  // Call the gradient function
+  grad_fd_norm_lib<TypeParam>(
+      this->psi_c.data(), this->edge_cell_idx.data(),
+      this->edge_cell_blk.data(), this->inv_dual_edge_length.data(),
+      this->grad_norm_psi_e.data(), this->i_startblk, this->i_endblk,
+      this->i_startidx_in, this->i_endidx_in, this->slev[0],
+      this->elev[0], this->nproma, this->nlev, this->nblks_c,
+      this->nblks_e, this->lacc);
+
+  // Expected values based on the initialization
+  // For edge 0: (cell 1 value - cell 0 value) * inv_dual_edge_length
+  EXPECT_NEAR(this->grad_norm_psi_e[grad_psi_at(0, 0, 0)],
+              (this->psi_c[psi_c_at(1, 0, 0)] - this->psi_c[psi_c_at(0, 0, 0)]) * 0.5,
+              1e-6);
+  EXPECT_NEAR(this->grad_norm_psi_e[grad_psi_at(0, 1, 0)],
+              (this->psi_c[psi_c_at(1, 1, 0)] - this->psi_c[psi_c_at(0, 1, 0)]) * 0.5,
+              1e-6);
+
+  // For edge 1: (cell 2 value - cell 1 value) * inv_dual_edge_length
+  EXPECT_NEAR(this->grad_norm_psi_e[grad_psi_at(1, 0, 0)],
+              (this->psi_c[psi_c_at(2, 0, 0)] - this->psi_c[psi_c_at(1, 0, 0)]) * 0.5,
+              1e-6);
+  EXPECT_NEAR(this->grad_norm_psi_e[grad_psi_at(1, 1, 0)],
+              (this->psi_c[psi_c_at(2, 1, 0)] - this->psi_c[psi_c_at(1, 1, 0)]) * 0.5,
+              1e-6);
+
+  // For edge 2: (cell 0 value - cell 2 value) * inv_dual_edge_length
+  EXPECT_NEAR(this->grad_norm_psi_e[grad_psi_at(2, 0, 0)],
+              (this->psi_c[psi_c_at(0, 0, 0)] - this->psi_c[psi_c_at(2, 0, 0)]) * 0.5,
+              1e-6);
+  EXPECT_NEAR(this->grad_norm_psi_e[grad_psi_at(2, 1, 0)],
+              (this->psi_c[psi_c_at(0, 1, 0)] - this->psi_c[psi_c_at(2, 1, 0)]) * 0.5,
+              1e-6);
+}
+
+TYPED_TEST(GradientNormal3DTest, TestRandom) {
+  constexpr int nproma = this->nproma;
+  constexpr int nlev = this->nlev;
+  constexpr int nblks_c = this->nblks_c;
+  constexpr int nblks_e = this->nblks_e;
+
+  const auto &psi_c_at = at<nproma, nlev, nblks_c>;
+  const auto &edge_cell_at = at<nproma, nblks_e, 2>;
+  const auto &inv_edge_at = at<nproma, nblks_e>;
+  const auto &grad_psi_at = at<nproma, nlev, nblks_e>;
+
+  // Set up random number generators
+  std::random_device rd;
+  std::mt19937 gen(rd());
+  std::uniform_int_distribution<int> int_distrib(0, nproma - 1);
+  std::uniform_real_distribution<TypeParam> real_distrib(-10.0, 10.0);
+
+  // Initialization with random values
+  for (int i = 0; i < nproma; ++i) {
+    for (int k = 0; k < nlev; ++k) {
+      this->psi_c[psi_c_at(i, k, 0)] = real_distrib(gen);
+    }
+
+    // Set random edge cell indices
+    for (int j = 0; j < 2; ++j) {
+      this->edge_cell_idx[edge_cell_at(i, 0, j)] = int_distrib(gen);
+      this->edge_cell_blk[edge_cell_at(i, 0, j)] = 0; // Keep in same block for simplicity
+    }
+
+    // Random inverse dual edge length
+    this->inv_dual_edge_length[inv_edge_at(i, 0)] = real_distrib(gen);
+
+    // Initialize grad_norm_psi_e to random values
+    for (int k = 0; k < nlev; ++k) {
+      this->grad_norm_psi_e[grad_psi_at(i, k, 0)] = real_distrib(gen);
+    }
+  }
+
+  // Call the gradient function
+  grad_fd_norm_lib<TypeParam>(
+      this->psi_c.data(), this->edge_cell_idx.data(),
+      this->edge_cell_blk.data(), this->inv_dual_edge_length.data(),
+      this->grad_norm_psi_e.data(), this->i_startblk, this->i_endblk,
+      this->i_startidx_in, this->i_endidx_in, this->slev[0],
+      this->elev[0], this->nproma, this->nlev, this->nblks_c,
+      this->nblks_e, this->lacc);
+
+  // Calculate reference values separately and verify results
+  std::vector<TypeParam> ref_grad_psi_e(nproma * nlev * nblks_e, 0.0);
+
+  for (int jb = this->i_startblk; jb < this->i_endblk; ++jb) {
+    int i_startidx, i_endidx;
+    get_indices_e_lib(this->i_startidx_in, this->i_endidx_in, nproma, jb,
+                     this->i_startblk, this->i_endblk, i_startidx, i_endidx);
+
+    for (int jk = this->slev[0]; jk < this->elev[0]; ++jk) {
+      for (int je = i_startidx; je < i_endidx; ++je) {
+        ref_grad_psi_e[grad_psi_at(je, jk, jb)] =
+            (this->psi_c[psi_c_at(this->edge_cell_idx[edge_cell_at(je, jb, 1)], jk,
+                               this->edge_cell_blk[edge_cell_at(je, jb, 1)])] -
+             this->psi_c[psi_c_at(this->edge_cell_idx[edge_cell_at(je, jb, 0)], jk,
+                               this->edge_cell_blk[edge_cell_at(je, jb, 0)])]) *
+            this->inv_dual_edge_length[inv_edge_at(je, jb)];
+      }
+    }
+  }
+
+  // Verify results
+  for (int i = 0; i < nproma; ++i) {
+    for (int k = 0; k < nlev; ++k) {
+      EXPECT_NEAR(this->grad_norm_psi_e[grad_psi_at(i, k, 0)],
+                  ref_grad_psi_e[grad_psi_at(i, k, 0)], 1e-5)
+          << "Results differ at i=" << i << ", k=" << k;
+    }
+  }
+}