From f91b8978e4d423bed70529bf36da851f9610cca5 Mon Sep 17 00:00:00 2001
From: Pradipta Samanta <samanta@dkrz.de>
Date: Thu, 20 Mar 2025 18:46:08 +0100
Subject: [PATCH] added the first test for mo_lib_gradients
---
test/c/CMakeLists.txt | 1 +
test/c/test_horizontal_gradients.cpp | 211 +++++++++++++++++++++++++++
2 files changed, 212 insertions(+)
create mode 100644 test/c/test_horizontal_gradients.cpp
diff --git a/test/c/CMakeLists.txt b/test/c/CMakeLists.txt
index 90ab1e3..43c9ac4 100644
--- a/test/c/CMakeLists.txt
+++ b/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
diff --git a/test/c/test_horizontal_gradients.cpp b/test/c/test_horizontal_gradients.cpp
new file mode 100644
index 0000000..d5a3697
--- /dev/null
+++ b/test/c/test_horizontal_gradients.cpp
@@ -0,0 +1,211 @@
+// 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;
+ }
+ }
+}
--
GitLab