Tested the c++ codes run on GPUs and fixed unit tests associated with them

What is the bug

The c++ codes were not running in parallel using Kokkos and were not tested properly

How do you fix it

How urgent is the bugfix

• I need it as soon as possible
• I can wait for a couple of days
• None of my current codes is directly affected

Mandatory steps before review

• Gitlab CI passes (Hint: use make format for linting)
• Bugfix is covered by additional unit tests
• Mark the merge request as ready by removing Draft:

Mandatory steps before merge

• Reviewed by a maintainer
• Incorporate review suggestions
• Remember to edit the commit message and select the proper changelog category (feature/bugfix/other)
• Prior to merging, please remove any boilerplate from the MR description, retaining only the What is the bug and How do you fix it section to maintain

You are not supposed to merge this request by yourself, the maintainers of libiconmath take care of this action!

test/c/test_horizontal_div.cpp

@@ -68,7 +68,7 @@ protected:
         f4dout("f4dout", dim_combine(nproma, nlev, nblks_c, dim4d)),
         cell_neighbor_idx("cell_neighbor_idx", dim_combine(nproma, nblks_c, 3)),
         cell_neighbor_blk("cell_neighbor_blk", dim_combine(nproma, nblks_c, 3)),
-        avg_coeff("avg_coeff", dim_combine(nproma, nlev, nblks_c)),
+        avg_coeff("avg_coeff", dim_combine(nproma, 4, nblks_c)),
         opt_in2("opt_in2", dim_combine(nproma, nlev, nblks_e)),
         opt_out2("opt_out2", dim_combine(nproma, nlev, nblks_c))
   {
@@ -628,7 +628,6 @@ TYPED_TEST(HorizontalDivTest, TestDiv4DRandom) {
   }
 }
 
-/*
 TYPED_TEST(HorizontalDivTest, TestDivAvgSpecific) {
   constexpr int nproma = this->nproma;
   constexpr int nlev = this->nlev;
@@ -640,6 +639,8 @@ TYPED_TEST(HorizontalDivTest, TestDivAvgSpecific) {
   const auto &cell_edge_at = at<nproma, nblks_c, 3>;
   const auto &geofac_div_at = at<nproma, 3, nblks_c>;
   const auto &div_vec_c_at = at<nproma, nlev, nblks_c>;
+  const auto &cell_neighbor_at = at<nproma, nblks_c, 3>;
+  const auto &avg_coeff_at = at<nproma, 4, nblks_c>;
 
   // Vectors for additional parameters
   // Vectors for block and index ranges
@@ -653,51 +654,71 @@ TYPED_TEST(HorizontalDivTest, TestDivAvgSpecific) {
   bool l_limited_area = true;
   bool l2fields = true;
 
-  const auto &cell_neighbor_at = at<nproma, nblks_c, 3>;
-  const auto &avg_coeff_at = at<nproma, 4, nblks_c>;
+  // Create mirror views to store data on the host
+  auto vec_e_h = Kokkos::create_mirror_view(this->vec_e);
+  auto cell_neighbor_idx_h = Kokkos::create_mirror_view(this->cell_neighbor_idx);
+  auto cell_neighbor_blk_h = Kokkos::create_mirror_view(this->cell_neighbor_blk);
+  auto cell_edge_idx_h = Kokkos::create_mirror_view(this->cell_edge_idx);
+  auto cell_edge_blk_h = Kokkos::create_mirror_view(this->cell_edge_blk);
+  auto geofac_div_h = Kokkos::create_mirror_view(this->geofac_div);
+  auto avg_coeff_h = Kokkos::create_mirror_view(this->avg_coeff);
+  auto div_vec_c_h = Kokkos::create_mirror_view(this->div_vec_c);
+  auto opt_in2_h = Kokkos::create_mirror_view(this->opt_in2);
+  auto opt_out2_h = Kokkos::create_mirror_view(this->opt_out2);
 
   // Initialize the vectors with specific values
   for (int i = 0; i < nproma; ++i) {
     for (int k = 0; k < nlev; ++k) {
-      this->vec_e[vec_e_at(i, k, 0)] = (i + 1) * (k + 1); // Simple pattern
-      this->opt_in2[vec_e_at(i, k, 0)] =
-          (i + 1) * (k + 1) * 0.5; // Half of vec_e
+      vec_e_h[vec_e_at(i, k, 0)] = static_cast<TypeParam>((i + 1) * (k + 1)); // Simple pattern
+      opt_in2_h[vec_e_at(i, k, 0)] = static_cast<TypeParam>((i + 1) * (k + 1) * 0.5); // Half of vec_e
     }
 
     // Set edge indices to point to specific cells
-    this->cell_edge_idx[cell_edge_at(i, 0, 0)] = i;
-    this->cell_edge_idx[cell_edge_at(i, 0, 1)] = (i + 1) % nproma;
-    this->cell_edge_idx[cell_edge_at(i, 0, 2)] = (i + 2) % nproma;
+    cell_edge_idx_h[cell_edge_at(i, 0, 0)] = i;
+    cell_edge_idx_h[cell_edge_at(i, 0, 1)] = (i + 1) % nproma;
+    cell_edge_idx_h[cell_edge_at(i, 0, 2)] = (i + 2) % nproma;
 
     // Set neighbor indices similarly
-    this->cell_neighbor_idx[cell_neighbor_at(i, 0, 0)] = i;
-    this->cell_neighbor_idx[cell_neighbor_at(i, 0, 1)] = (i + 1) % nproma;
-    this->cell_neighbor_idx[cell_neighbor_at(i, 0, 2)] = (i + 2) % nproma;
+    cell_neighbor_idx_h[cell_neighbor_at(i, 0, 0)] = i;
+    cell_neighbor_idx_h[cell_neighbor_at(i, 0, 1)] = (i + 1) % nproma;
+    cell_neighbor_idx_h[cell_neighbor_at(i, 0, 2)] = (i + 2) % nproma;
 
     // All edges and neighbors are in the same block for this test
     for (int j = 0; j < 3; ++j) {
-      this->cell_edge_blk[cell_edge_at(i, 0, j)] = 0;
-      this->cell_neighbor_blk[cell_neighbor_at(i, 0, j)] = 0;
+      cell_edge_blk_h[cell_edge_at(i, 0, j)] = 0;
+      cell_neighbor_blk_h[cell_neighbor_at(i, 0, j)] = 0;
     }
 
     // Geometric factors
-    this->geofac_div[geofac_div_at(i, 0, 0)] = 0.5;
-    this->geofac_div[geofac_div_at(i, 1, 0)] = 0.3;
-    this->geofac_div[geofac_div_at(i, 2, 0)] = 0.2;
+    geofac_div_h[geofac_div_at(i, 0, 0)] = static_cast<TypeParam>(0.5);
+    geofac_div_h[geofac_div_at(i, 1, 0)] = static_cast<TypeParam>(0.3);
+    geofac_div_h[geofac_div_at(i, 2, 0)] = static_cast<TypeParam>(0.2);
 
     // Average coefficients
-    this->avg_coeff[avg_coeff_at(i, 0, 0)] = 0.4; // Self
-    this->avg_coeff[avg_coeff_at(i, 1, 0)] = 0.2; // First neighbor
-    this->avg_coeff[avg_coeff_at(i, 2, 0)] = 0.2; // Second neighbor
-    this->avg_coeff[avg_coeff_at(i, 3, 0)] = 0.2; // Third neighbor
+    avg_coeff_h[avg_coeff_at(i, 0, 0)] = static_cast<TypeParam>(0.4); // Self
+    avg_coeff_h[avg_coeff_at(i, 1, 0)] = static_cast<TypeParam>(0.2); // First neighbor
+    avg_coeff_h[avg_coeff_at(i, 2, 0)] = static_cast<TypeParam>(0.2); // Second neighbor
+    avg_coeff_h[avg_coeff_at(i, 3, 0)] = static_cast<TypeParam>(0.2); // Third neighbor
 
     // Initialize div_vec_c and opt_out2 to zero
     for (int k = 0; k < nlev; ++k) {
-      this->div_vec_c[div_vec_c_at(i, k, 0)] = 0.0;
-      this->opt_out2[div_vec_c_at(i, k, 0)] = 0.0;
+      div_vec_c_h[div_vec_c_at(i, k, 0)] = static_cast<TypeParam>(0.0);
+      opt_out2_h[div_vec_c_at(i, k, 0)] = static_cast<TypeParam>(0.0);
     }
   }
 
+  // Copy initialized data to device
+  Kokkos::deep_copy(this->vec_e, vec_e_h);
+  Kokkos::deep_copy(this->cell_neighbor_idx, cell_neighbor_idx_h);
+  Kokkos::deep_copy(this->cell_neighbor_blk, cell_neighbor_blk_h);
+  Kokkos::deep_copy(this->cell_edge_idx, cell_edge_idx_h);
+  Kokkos::deep_copy(this->cell_edge_blk, cell_edge_blk_h);
+  Kokkos::deep_copy(this->geofac_div, geofac_div_h);
+  Kokkos::deep_copy(this->avg_coeff, avg_coeff_h);
+  Kokkos::deep_copy(this->div_vec_c, div_vec_c_h);
+  Kokkos::deep_copy(this->opt_in2, opt_in2_h);
+  Kokkos::deep_copy(this->opt_out2, opt_out2_h);
+
   // Call the div_avg function
   div_avg<TypeParam>(
       this->vec_e.data(), this->cell_neighbor_idx.data(),
@@ -709,19 +730,25 @@ TYPED_TEST(HorizontalDivTest, TestDivAvgSpecific) {
       this->nproma, patch_id, l_limited_area, l2fields, this->lacc, this->nlev,
       this->nblks_c, this->nblks_e);
 
-  EXPECT_NEAR(this->div_vec_c[div_vec_c_at(0, 0, 0)], 1.88, 1e-6);
-  EXPECT_NEAR(this->div_vec_c[div_vec_c_at(0, 1, 0)], 3.76, 1e-6);
-  EXPECT_NEAR(this->div_vec_c[div_vec_c_at(1, 0, 0)], 2.04, 1e-6);
-  EXPECT_NEAR(this->div_vec_c[div_vec_c_at(1, 1, 0)], 4.08, 1e-6);
-  EXPECT_NEAR(this->div_vec_c[div_vec_c_at(2, 0, 0)], 2.08, 1e-6);
-  EXPECT_NEAR(this->div_vec_c[div_vec_c_at(2, 1, 0)], 4.16, 1e-6);
-
-  EXPECT_NEAR(this->opt_out2[div_vec_c_at(0, 0, 0)], 0.94, 1e-6);
-  EXPECT_NEAR(this->opt_out2[div_vec_c_at(0, 1, 0)], 1.88, 1e-6);
-  EXPECT_NEAR(this->opt_out2[div_vec_c_at(1, 0, 0)], 1.02, 1e-6);
-  EXPECT_NEAR(this->opt_out2[div_vec_c_at(1, 1, 0)], 2.04, 1e-6);
-  EXPECT_NEAR(this->opt_out2[div_vec_c_at(2, 0, 0)], 1.04, 1e-6);
-  EXPECT_NEAR(this->opt_out2[div_vec_c_at(2, 1, 0)], 2.08, 1e-6);
+  // Copy results back to host for verification
+  Kokkos::deep_copy(div_vec_c_h, this->div_vec_c);
+  Kokkos::deep_copy(opt_out2_h, this->opt_out2);
+
+  // Verify first field results
+  EXPECT_NEAR(div_vec_c_h[div_vec_c_at(0, 0, 0)], static_cast<TypeParam>(1.88), 1e-6);
+  EXPECT_NEAR(div_vec_c_h[div_vec_c_at(0, 1, 0)], static_cast<TypeParam>(3.76), 1e-6);
+  EXPECT_NEAR(div_vec_c_h[div_vec_c_at(1, 0, 0)], static_cast<TypeParam>(2.04), 1e-6);
+  EXPECT_NEAR(div_vec_c_h[div_vec_c_at(1, 1, 0)], static_cast<TypeParam>(4.08), 1e-6);
+  EXPECT_NEAR(div_vec_c_h[div_vec_c_at(2, 0, 0)], static_cast<TypeParam>(2.08), 1e-6);
+  EXPECT_NEAR(div_vec_c_h[div_vec_c_at(2, 1, 0)], static_cast<TypeParam>(4.16), 1e-6);
+
+  // Verify second field results
+  EXPECT_NEAR(opt_out2_h[div_vec_c_at(0, 0, 0)], static_cast<TypeParam>(0.94), 1e-6);
+  EXPECT_NEAR(opt_out2_h[div_vec_c_at(0, 1, 0)], static_cast<TypeParam>(1.88), 1e-6);
+  EXPECT_NEAR(opt_out2_h[div_vec_c_at(1, 0, 0)], static_cast<TypeParam>(1.02), 1e-6);
+  EXPECT_NEAR(opt_out2_h[div_vec_c_at(1, 1, 0)], static_cast<TypeParam>(2.04), 1e-6);
+  EXPECT_NEAR(opt_out2_h[div_vec_c_at(2, 0, 0)], static_cast<TypeParam>(1.04), 1e-6);
+  EXPECT_NEAR(opt_out2_h[div_vec_c_at(2, 1, 0)], static_cast<TypeParam>(2.08), 1e-6);
 }
 
 TYPED_TEST(HorizontalDivTest, TestDivAvgRandom) {
@@ -749,47 +776,69 @@ TYPED_TEST(HorizontalDivTest, TestDivAvgRandom) {
   const auto &cell_neighbor_at = at<nproma, nblks_c, 3>;
   const auto &avg_coeff_at = at<nproma, 4, nblks_c>;
 
+  // Create mirror views to store data on the host
+  auto vec_e_h = Kokkos::create_mirror_view(this->vec_e);
+  auto cell_neighbor_idx_h = Kokkos::create_mirror_view(this->cell_neighbor_idx);
+  auto cell_neighbor_blk_h = Kokkos::create_mirror_view(this->cell_neighbor_blk);
+  auto cell_edge_idx_h = Kokkos::create_mirror_view(this->cell_edge_idx);
+  auto cell_edge_blk_h = Kokkos::create_mirror_view(this->cell_edge_blk);
+  auto geofac_div_h = Kokkos::create_mirror_view(this->geofac_div);
+  auto avg_coeff_h = Kokkos::create_mirror_view(this->avg_coeff);
+  auto div_vec_c_h = Kokkos::create_mirror_view(this->div_vec_c);
+  auto opt_in2_h = Kokkos::create_mirror_view(this->opt_in2);
+  auto opt_out2_h = Kokkos::create_mirror_view(this->opt_out2);
+
   // 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);
+  std::uniform_real_distribution<TypeParam> real_distrib(-1.0, 1.0);
 
   // Initialize with random values
   for (int i = 0; i < nproma; ++i) {
     for (int k = 0; k < nlev; ++k) {
-      this->vec_e[vec_e_at(i, k, 0)] = real_distrib(gen);
-      this->opt_in2[vec_e_at(i, k, 0)] = real_distrib(gen);
+      vec_e_h[vec_e_at(i, k, 0)] = real_distrib(gen);
+      opt_in2_h[vec_e_at(i, k, 0)] = real_distrib(gen);
     }
 
     // Set random edge indices
     for (int j = 0; j < 3; ++j) {
-      this->cell_edge_idx[cell_edge_at(i, 0, j)] = int_distrib(gen);
-      this->cell_edge_blk[cell_edge_at(i, 0, j)] =
-          0; // Keep in same block for simplicity
+      cell_edge_idx_h[cell_edge_at(i, 0, j)] = int_distrib(gen);
+      cell_edge_blk_h[cell_edge_at(i, 0, j)] = 0; // Keep in same block for simplicity
 
-      this->cell_neighbor_idx[cell_neighbor_at(i, 0, j)] = int_distrib(gen);
-      this->cell_neighbor_blk[cell_neighbor_at(i, 0, j)] =
-          0; // Keep in same block for simplicity
+      cell_neighbor_idx_h[cell_neighbor_at(i, 0, j)] = int_distrib(gen);
+      cell_neighbor_blk_h[cell_neighbor_at(i, 0, j)] = 0; // Keep in same block for simplicity
     }
 
     // Random geometric factors
     for (int j = 0; j < 3; ++j) {
-      this->geofac_div[geofac_div_at(i, j, 0)] = real_distrib(gen);
+      geofac_div_h[geofac_div_at(i, j, 0)] = real_distrib(gen);
     }
 
     // Random average coefficients
     for (int j = 0; j < 4; ++j) {
-      this->avg_coeff[avg_coeff_at(i, j, 0)] = real_distrib(gen);
+      avg_coeff_h[avg_coeff_at(i, j, 0)] = real_distrib(gen);
     }
 
     // Random initial values for div_vec_c and opt_out2
     for (int k = 0; k < nlev; ++k) {
-      this->div_vec_c[div_vec_c_at(i, k, 0)] = real_distrib(gen);
-      this->opt_out2[div_vec_c_at(i, k, 0)] = real_distrib(gen);
+      div_vec_c_h[div_vec_c_at(i, k, 0)] = static_cast<TypeParam>(0.0);
+      opt_out2_h[div_vec_c_at(i, k, 0)] = static_cast<TypeParam>(0.0);
     }
   }
 
+  // Copy initialized data to device
+  Kokkos::deep_copy(this->vec_e, vec_e_h);
+  Kokkos::deep_copy(this->cell_neighbor_idx, cell_neighbor_idx_h);
+  Kokkos::deep_copy(this->cell_neighbor_blk, cell_neighbor_blk_h);
+  Kokkos::deep_copy(this->cell_edge_idx, cell_edge_idx_h);
+  Kokkos::deep_copy(this->cell_edge_blk, cell_edge_blk_h);
+  Kokkos::deep_copy(this->geofac_div, geofac_div_h);
+  Kokkos::deep_copy(this->avg_coeff, avg_coeff_h);
+  Kokkos::deep_copy(this->div_vec_c, div_vec_c_h);
+  Kokkos::deep_copy(this->opt_in2, opt_in2_h);
+  Kokkos::deep_copy(this->opt_out2, opt_out2_h);
+
   // Call the div_avg function
   div_avg<TypeParam>(
       this->vec_e.data(), this->cell_neighbor_idx.data(),
@@ -801,6 +850,10 @@ TYPED_TEST(HorizontalDivTest, TestDivAvgRandom) {
       this->nproma, patch_id, l_limited_area, l2fields, this->lacc, this->nlev,
       this->nblks_c, this->nblks_e);
 
+  // Copy results back to host for verification
+  Kokkos::deep_copy(div_vec_c_h, this->div_vec_c);
+  Kokkos::deep_copy(opt_out2_h, this->opt_out2);
+
   // Calculate reference values manually
   std::vector<TypeParam> aux_c(dim_combine(nproma, nlev, nblks_c));
   std::vector<TypeParam> aux_c2(dim_combine(nproma, nlev, nblks_c));
@@ -816,32 +869,32 @@ TYPED_TEST(HorizontalDivTest, TestDivAvgRandom) {
     for (int jk = this->slev[0]; jk < this->elev[0]; ++jk) {
       for (int jc = i_startidx; jc < i_endidx; ++jc) {
         aux_c[div_vec_c_at(jc, jk, jb)] =
-            this->vec_e[vec_e_at(
-                this->cell_edge_idx[cell_edge_at(jc, jb, 0)], jk,
-                this->cell_edge_blk[cell_edge_at(jc, jb, 0)])] *
-                this->geofac_div[geofac_div_at(jc, 0, jb)] +
-            this->vec_e[vec_e_at(
-                this->cell_edge_idx[cell_edge_at(jc, jb, 1)], jk,
-                this->cell_edge_blk[cell_edge_at(jc, jb, 1)])] *
-                this->geofac_div[geofac_div_at(jc, 1, jb)] +
-            this->vec_e[vec_e_at(
-                this->cell_edge_idx[cell_edge_at(jc, jb, 2)], jk,
-                this->cell_edge_blk[cell_edge_at(jc, jb, 2)])] *
-                this->geofac_div[geofac_div_at(jc, 2, jb)];
+            vec_e_h[vec_e_at(
+                cell_edge_idx_h[cell_edge_at(jc, jb, 0)], jk,
+                cell_edge_blk_h[cell_edge_at(jc, jb, 0)])] *
+                geofac_div_h[geofac_div_at(jc, 0, jb)] +
+            vec_e_h[vec_e_at(
+                cell_edge_idx_h[cell_edge_at(jc, jb, 1)], jk,
+                cell_edge_blk_h[cell_edge_at(jc, jb, 1)])] *
+                geofac_div_h[geofac_div_at(jc, 1, jb)] +
+            vec_e_h[vec_e_at(
+                cell_edge_idx_h[cell_edge_at(jc, jb, 2)], jk,
+                cell_edge_blk_h[cell_edge_at(jc, jb, 2)])] *
+                geofac_div_h[geofac_div_at(jc, 2, jb)];
 
         aux_c2[div_vec_c_at(jc, jk, jb)] =
-            this->opt_in2[vec_e_at(
-                this->cell_edge_idx[cell_edge_at(jc, jb, 0)], jk,
-                this->cell_edge_blk[cell_edge_at(jc, jb, 0)])] *
-                this->geofac_div[geofac_div_at(jc, 0, jb)] +
-            this->opt_in2[vec_e_at(
-                this->cell_edge_idx[cell_edge_at(jc, jb, 1)], jk,
-                this->cell_edge_blk[cell_edge_at(jc, jb, 1)])] *
-                this->geofac_div[geofac_div_at(jc, 1, jb)] +
-            this->opt_in2[vec_e_at(
-                this->cell_edge_idx[cell_edge_at(jc, jb, 2)], jk,
-                this->cell_edge_blk[cell_edge_at(jc, jb, 2)])] *
-                this->geofac_div[geofac_div_at(jc, 2, jb)];
+            opt_in2_h[vec_e_at(
+                cell_edge_idx_h[cell_edge_at(jc, jb, 0)], jk,
+                cell_edge_blk_h[cell_edge_at(jc, jb, 0)])] *
+                geofac_div_h[geofac_div_at(jc, 0, jb)] +
+            opt_in2_h[vec_e_at(
+                cell_edge_idx_h[cell_edge_at(jc, jb, 1)], jk,
+                cell_edge_blk_h[cell_edge_at(jc, jb, 1)])] *
+                geofac_div_h[geofac_div_at(jc, 1, jb)] +
+            opt_in2_h[vec_e_at(
+                cell_edge_idx_h[cell_edge_at(jc, jb, 2)], jk,
+                cell_edge_blk_h[cell_edge_at(jc, jb, 2)])] *
+                geofac_div_h[geofac_div_at(jc, 2, jb)];
       }
     }
   }
@@ -872,35 +925,35 @@ TYPED_TEST(HorizontalDivTest, TestDivAvgRandom) {
       for (int jc = i_startidx; jc < i_endidx; ++jc) {
         ref_div_vec_c[div_vec_c_at(jc, jk, jb)] =
             aux_c[div_vec_c_at(jc, jk, jb)] *
-                this->avg_coeff[avg_coeff_at(jc, 0, jb)] +
+                avg_coeff_h[avg_coeff_at(jc, 0, jb)] +
             aux_c[div_vec_c_at(
-                this->cell_neighbor_idx[cell_neighbor_at(jc, jb, 0)], jk,
-                this->cell_neighbor_blk[cell_neighbor_at(jc, jb, 0)])] *
-                this->avg_coeff[avg_coeff_at(jc, 1, jb)] +
+                cell_neighbor_idx_h[cell_neighbor_at(jc, jb, 0)], jk,
+                cell_neighbor_blk_h[cell_neighbor_at(jc, jb, 0)])] *
+                avg_coeff_h[avg_coeff_at(jc, 1, jb)] +
             aux_c[div_vec_c_at(
-                this->cell_neighbor_idx[cell_neighbor_at(jc, jb, 1)], jk,
-                this->cell_neighbor_blk[cell_neighbor_at(jc, jb, 1)])] *
-                this->avg_coeff[avg_coeff_at(jc, 2, jb)] +
+                cell_neighbor_idx_h[cell_neighbor_at(jc, jb, 1)], jk,
+                cell_neighbor_blk_h[cell_neighbor_at(jc, jb, 1)])] *
+                avg_coeff_h[avg_coeff_at(jc, 2, jb)] +
             aux_c[div_vec_c_at(
-                this->cell_neighbor_idx[cell_neighbor_at(jc, jb, 2)], jk,
-                this->cell_neighbor_blk[cell_neighbor_at(jc, jb, 2)])] *
-                this->avg_coeff[avg_coeff_at(jc, 3, jb)];
+                cell_neighbor_idx_h[cell_neighbor_at(jc, jb, 2)], jk,
+                cell_neighbor_blk_h[cell_neighbor_at(jc, jb, 2)])] *
+                avg_coeff_h[avg_coeff_at(jc, 3, jb)];
 
         ref_opt_out2[div_vec_c_at(jc, jk, jb)] =
             aux_c2[div_vec_c_at(jc, jk, jb)] *
-                this->avg_coeff[avg_coeff_at(jc, 0, jb)] +
+                avg_coeff_h[avg_coeff_at(jc, 0, jb)] +
             aux_c2[div_vec_c_at(
-                this->cell_neighbor_idx[cell_neighbor_at(jc, jb, 0)], jk,
-                this->cell_neighbor_blk[cell_neighbor_at(jc, jb, 0)])] *
-                this->avg_coeff[avg_coeff_at(jc, 1, jb)] +
+                cell_neighbor_idx_h[cell_neighbor_at(jc, jb, 0)], jk,
+                cell_neighbor_blk_h[cell_neighbor_at(jc, jb, 0)])] *
+                avg_coeff_h[avg_coeff_at(jc, 1, jb)] +
             aux_c2[div_vec_c_at(
-                this->cell_neighbor_idx[cell_neighbor_at(jc, jb, 1)], jk,
-                this->cell_neighbor_blk[cell_neighbor_at(jc, jb, 1)])] *
-                this->avg_coeff[avg_coeff_at(jc, 2, jb)] +
+                cell_neighbor_idx_h[cell_neighbor_at(jc, jb, 1)], jk,
+                cell_neighbor_blk_h[cell_neighbor_at(jc, jb, 1)])] *
+                avg_coeff_h[avg_coeff_at(jc, 2, jb)] +
             aux_c2[div_vec_c_at(
-                this->cell_neighbor_idx[cell_neighbor_at(jc, jb, 2)], jk,
-                this->cell_neighbor_blk[cell_neighbor_at(jc, jb, 2)])] *
-                this->avg_coeff[avg_coeff_at(jc, 3, jb)];
+                cell_neighbor_idx_h[cell_neighbor_at(jc, jb, 2)], jk,
+                cell_neighbor_blk_h[cell_neighbor_at(jc, jb, 2)])] *
+                avg_coeff_h[avg_coeff_at(jc, 3, jb)];
       }
     }
   }
@@ -908,11 +961,11 @@ TYPED_TEST(HorizontalDivTest, TestDivAvgRandom) {
   // Verify results
   for (int i = 0; i < nproma; ++i) {
     for (int k = 0; k < nlev; ++k) {
-      EXPECT_NEAR(this->div_vec_c[div_vec_c_at(i, k, 0)],
+      EXPECT_NEAR(div_vec_c_h[div_vec_c_at(i, k, 0)],
                   ref_div_vec_c[div_vec_c_at(i, k, 0)], 1e-5)
           << "div_vec_c results differ at i=" << i << ", k=" << k;
 
-      EXPECT_NEAR(this->opt_out2[div_vec_c_at(i, k, 0)],
+      EXPECT_NEAR(opt_out2_h[div_vec_c_at(i, k, 0)],
                   ref_opt_out2[div_vec_c_at(i, k, 0)], 1e-5)
           << "opt_out2 results differ at i=" << i << ", k=" << k;
     }
@@ -945,48 +998,71 @@ TYPED_TEST(HorizontalDivTest, TestDivAvgSpecificNoL2fields) {
   const auto &cell_neighbor_at = at<nproma, nblks_c, 3>;
   const auto &avg_coeff_at = at<nproma, 4, nblks_c>;
 
+  // Create mirror views to store data on the host
+  auto vec_e_h = Kokkos::create_mirror_view(this->vec_e);
+  auto cell_neighbor_idx_h = Kokkos::create_mirror_view(this->cell_neighbor_idx);
+  auto cell_neighbor_blk_h = Kokkos::create_mirror_view(this->cell_neighbor_blk);
+  auto cell_edge_idx_h = Kokkos::create_mirror_view(this->cell_edge_idx);
+  auto cell_edge_blk_h = Kokkos::create_mirror_view(this->cell_edge_blk);
+  auto geofac_div_h = Kokkos::create_mirror_view(this->geofac_div);
+  auto avg_coeff_h = Kokkos::create_mirror_view(this->avg_coeff);
+  auto div_vec_c_h = Kokkos::create_mirror_view(this->div_vec_c);
+  auto opt_in2_h = Kokkos::create_mirror_view(this->opt_in2);
+  auto opt_out2_h = Kokkos::create_mirror_view(this->opt_out2);
+
   // Initialize the vectors with specific values
   for (int i = 0; i < nproma; ++i) {
     for (int k = 0; k < nlev; ++k) {
-      this->vec_e[vec_e_at(i, k, 0)] = (i + 1) * (k + 1); // Simple pattern
-      this->opt_in2[vec_e_at(i, k, 0)] =
-          (i + 1) * (k + 1) * 0.5; // Half of vec_e
+      vec_e_h[vec_e_at(i, k, 0)] = static_cast<TypeParam>((i + 1) * (k + 1)); // Simple pattern
+      opt_in2_h[vec_e_at(i, k, 0)] = static_cast<TypeParam>((i + 1) * (k + 1) * 0.5); // Half of vec_e
     }
 
     // Set edge indices to point to specific cells
-    this->cell_edge_idx[cell_edge_at(i, 0, 0)] = i;
-    this->cell_edge_idx[cell_edge_at(i, 0, 1)] = (i + 1) % nproma;
-    this->cell_edge_idx[cell_edge_at(i, 0, 2)] = (i + 2) % nproma;
+    cell_edge_idx_h[cell_edge_at(i, 0, 0)] = i;
+    cell_edge_idx_h[cell_edge_at(i, 0, 1)] = (i + 1) % nproma;
+    cell_edge_idx_h[cell_edge_at(i, 0, 2)] = (i + 2) % nproma;
 
     // Set neighbor indices similarly
-    this->cell_neighbor_idx[cell_neighbor_at(i, 0, 0)] = i;
-    this->cell_neighbor_idx[cell_neighbor_at(i, 0, 1)] = (i + 1) % nproma;
-    this->cell_neighbor_idx[cell_neighbor_at(i, 0, 2)] = (i + 2) % nproma;
+    cell_neighbor_idx_h[cell_neighbor_at(i, 0, 0)] = i;
+    cell_neighbor_idx_h[cell_neighbor_at(i, 0, 1)] = (i + 1) % nproma;
+    cell_neighbor_idx_h[cell_neighbor_at(i, 0, 2)] = (i + 2) % nproma;
 
     // All edges and neighbors are in the same block for this test
     for (int j = 0; j < 3; ++j) {
-      this->cell_edge_blk[cell_edge_at(i, 0, j)] = 0;
-      this->cell_neighbor_blk[cell_neighbor_at(i, 0, j)] = 0;
+      cell_edge_blk_h[cell_edge_at(i, 0, j)] = 0;
+      cell_neighbor_blk_h[cell_neighbor_at(i, 0, j)] = 0;
     }
 
     // Geometric factors
-    this->geofac_div[geofac_div_at(i, 0, 0)] = 0.5;
-    this->geofac_div[geofac_div_at(i, 1, 0)] = 0.3;
-    this->geofac_div[geofac_div_at(i, 2, 0)] = 0.2;
+    geofac_div_h[geofac_div_at(i, 0, 0)] = static_cast<TypeParam>(0.5);
+    geofac_div_h[geofac_div_at(i, 1, 0)] = static_cast<TypeParam>(0.3);
+    geofac_div_h[geofac_div_at(i, 2, 0)] = static_cast<TypeParam>(0.2);
 
     // Average coefficients
-    this->avg_coeff[avg_coeff_at(i, 0, 0)] = 0.4; // Self
-    this->avg_coeff[avg_coeff_at(i, 1, 0)] = 0.2; // First neighbor
-    this->avg_coeff[avg_coeff_at(i, 2, 0)] = 0.2; // Second neighbor
-    this->avg_coeff[avg_coeff_at(i, 3, 0)] = 0.2; // Third neighbor
+    avg_coeff_h[avg_coeff_at(i, 0, 0)] = static_cast<TypeParam>(0.4); // Self
+    avg_coeff_h[avg_coeff_at(i, 1, 0)] = static_cast<TypeParam>(0.2); // First neighbor
+    avg_coeff_h[avg_coeff_at(i, 2, 0)] = static_cast<TypeParam>(0.2); // Second neighbor
+    avg_coeff_h[avg_coeff_at(i, 3, 0)] = static_cast<TypeParam>(0.2); // Third neighbor
 
     // Initialize div_vec_c and opt_out2 to zero
     for (int k = 0; k < nlev; ++k) {
-      this->div_vec_c[div_vec_c_at(i, k, 0)] = 0.0;
-      this->opt_out2[div_vec_c_at(i, k, 0)] = 0.0;
+      div_vec_c_h[div_vec_c_at(i, k, 0)] = static_cast<TypeParam>(0.0);
+      opt_out2_h[div_vec_c_at(i, k, 0)] = static_cast<TypeParam>(0.0);
     }
   }
 
+  // Copy initialized data to device
+  Kokkos::deep_copy(this->vec_e, vec_e_h);
+  Kokkos::deep_copy(this->cell_neighbor_idx, cell_neighbor_idx_h);
+  Kokkos::deep_copy(this->cell_neighbor_blk, cell_neighbor_blk_h);
+  Kokkos::deep_copy(this->cell_edge_idx, cell_edge_idx_h);
+  Kokkos::deep_copy(this->cell_edge_blk, cell_edge_blk_h);
+  Kokkos::deep_copy(this->geofac_div, geofac_div_h);
+  Kokkos::deep_copy(this->avg_coeff, avg_coeff_h);
+  Kokkos::deep_copy(this->div_vec_c, div_vec_c_h);
+  Kokkos::deep_copy(this->opt_in2, opt_in2_h);
+  Kokkos::deep_copy(this->opt_out2, opt_out2_h);
+
   // Call the div_avg function
   div_avg<TypeParam>(
       this->vec_e.data(), this->cell_neighbor_idx.data(),
@@ -998,19 +1074,25 @@ TYPED_TEST(HorizontalDivTest, TestDivAvgSpecificNoL2fields) {
       this->nproma, patch_id, l_limited_area, l2fields, this->lacc, this->nlev,
       this->nblks_c, this->nblks_e);
 
-  EXPECT_NEAR(this->div_vec_c[div_vec_c_at(0, 0, 0)], 1.88, 1e-6);
-  EXPECT_NEAR(this->div_vec_c[div_vec_c_at(0, 1, 0)], 3.76, 1e-6);
-  EXPECT_NEAR(this->div_vec_c[div_vec_c_at(1, 0, 0)], 2.04, 1e-6);
-  EXPECT_NEAR(this->div_vec_c[div_vec_c_at(1, 1, 0)], 4.08, 1e-6);
-  EXPECT_NEAR(this->div_vec_c[div_vec_c_at(2, 0, 0)], 2.08, 1e-6);
-  EXPECT_NEAR(this->div_vec_c[div_vec_c_at(2, 1, 0)], 4.16, 1e-6);
-
-  EXPECT_NEAR(this->opt_out2[div_vec_c_at(0, 0, 0)], 0.0, 1e-6);
-  EXPECT_NEAR(this->opt_out2[div_vec_c_at(0, 1, 0)], 0.0, 1e-6);
-  EXPECT_NEAR(this->opt_out2[div_vec_c_at(1, 0, 0)], 0.0, 1e-6);
-  EXPECT_NEAR(this->opt_out2[div_vec_c_at(1, 1, 0)], 0.0, 1e-6);
-  EXPECT_NEAR(this->opt_out2[div_vec_c_at(2, 0, 0)], 0.0, 1e-6);
-  EXPECT_NEAR(this->opt_out2[div_vec_c_at(2, 1, 0)], 0.0, 1e-6);
+  // Copy results back to host for verification
+  Kokkos::deep_copy(div_vec_c_h, this->div_vec_c);
+  Kokkos::deep_copy(opt_out2_h, this->opt_out2);
+
+  // Verify first field results
+  EXPECT_NEAR(div_vec_c_h[div_vec_c_at(0, 0, 0)], static_cast<TypeParam>(1.88), 1e-6);
+  EXPECT_NEAR(div_vec_c_h[div_vec_c_at(0, 1, 0)], static_cast<TypeParam>(3.76), 1e-6);
+  EXPECT_NEAR(div_vec_c_h[div_vec_c_at(1, 0, 0)], static_cast<TypeParam>(2.04), 1e-6);
+  EXPECT_NEAR(div_vec_c_h[div_vec_c_at(1, 1, 0)], static_cast<TypeParam>(4.08), 1e-6);
+  EXPECT_NEAR(div_vec_c_h[div_vec_c_at(2, 0, 0)], static_cast<TypeParam>(2.08), 1e-6);
+  EXPECT_NEAR(div_vec_c_h[div_vec_c_at(2, 1, 0)], static_cast<TypeParam>(4.16), 1e-6);
+
+  // Since l2fields=false, opt_out2 should not be modified
+  EXPECT_NEAR(opt_out2_h[div_vec_c_at(0, 0, 0)], static_cast<TypeParam>(0.0), 1e-6);
+  EXPECT_NEAR(opt_out2_h[div_vec_c_at(0, 1, 0)], static_cast<TypeParam>(0.0), 1e-6);
+  EXPECT_NEAR(opt_out2_h[div_vec_c_at(1, 0, 0)], static_cast<TypeParam>(0.0), 1e-6);
+  EXPECT_NEAR(opt_out2_h[div_vec_c_at(1, 1, 0)], static_cast<TypeParam>(0.0), 1e-6);
+  EXPECT_NEAR(opt_out2_h[div_vec_c_at(2, 0, 0)], static_cast<TypeParam>(0.0), 1e-6);
+  EXPECT_NEAR(opt_out2_h[div_vec_c_at(2, 1, 0)], static_cast<TypeParam>(0.0), 1e-6);
 }
 
 TYPED_TEST(HorizontalDivTest, TestDivAvgRandomNoL2fields) {
@@ -1038,49 +1120,69 @@ TYPED_TEST(HorizontalDivTest, TestDivAvgRandomNoL2fields) {
   const auto &cell_neighbor_at = at<nproma, nblks_c, 3>;
   const auto &avg_coeff_at = at<nproma, 4, nblks_c>;
 
+  // Create mirror views to store data on the host
+  auto vec_e_h = Kokkos::create_mirror_view(this->vec_e);
+  auto cell_neighbor_idx_h = Kokkos::create_mirror_view(this->cell_neighbor_idx);
+  auto cell_neighbor_blk_h = Kokkos::create_mirror_view(this->cell_neighbor_blk);
+  auto cell_edge_idx_h = Kokkos::create_mirror_view(this->cell_edge_idx);
+  auto cell_edge_blk_h = Kokkos::create_mirror_view(this->cell_edge_blk);
+  auto geofac_div_h = Kokkos::create_mirror_view(this->geofac_div);
+  auto avg_coeff_h = Kokkos::create_mirror_view(this->avg_coeff);
+  auto div_vec_c_h = Kokkos::create_mirror_view(this->div_vec_c);
+  auto opt_in2_h = Kokkos::create_mirror_view(this->opt_in2);
+  auto opt_out2_h = Kokkos::create_mirror_view(this->opt_out2);
+
   // 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);
+  std::uniform_real_distribution<TypeParam> real_distrib(-1.0, 1.0);
 
   // Initialize with random values
   for (int i = 0; i < nproma; ++i) {
     for (int k = 0; k < nlev; ++k) {
-      this->vec_e[vec_e_at(i, k, 0)] = real_distrib(gen);
-      this->opt_in2[vec_e_at(i, k, 0)] =
-          real_distrib(gen); // Not used but initialize anyway
+      vec_e_h[vec_e_at(i, k, 0)] = real_distrib(gen);
+      opt_in2_h[vec_e_at(i, k, 0)] = real_distrib(gen); // Not used but initialize anyway
     }
 
     // Set random edge indices
     for (int j = 0; j < 3; ++j) {
-      this->cell_edge_idx[cell_edge_at(i, 0, j)] = int_distrib(gen);
-      this->cell_edge_blk[cell_edge_at(i, 0, j)] =
-          0; // Keep in same block for simplicity
+      cell_edge_idx_h[cell_edge_at(i, 0, j)] = int_distrib(gen);
+      cell_edge_blk_h[cell_edge_at(i, 0, j)] = 0; // Keep in same block for simplicity
 
-      this->cell_neighbor_idx[cell_neighbor_at(i, 0, j)] = int_distrib(gen);
-      this->cell_neighbor_blk[cell_neighbor_at(i, 0, j)] =
-          0; // Keep in same block for simplicity
+      cell_neighbor_idx_h[cell_neighbor_at(i, 0, j)] = int_distrib(gen);
+      cell_neighbor_blk_h[cell_neighbor_at(i, 0, j)] = 0; // Keep in same block for simplicity
     }
 
     // Random geometric factors
     for (int j = 0; j < 3; ++j) {
-      this->geofac_div[geofac_div_at(i, j, 0)] = real_distrib(gen);
+      geofac_div_h[geofac_div_at(i, j, 0)] = real_distrib(gen);
     }
 
     // Random average coefficients
     for (int j = 0; j < 4; ++j) {
-      this->avg_coeff[avg_coeff_at(i, j, 0)] = real_distrib(gen);
+      avg_coeff_h[avg_coeff_at(i, j, 0)] = real_distrib(gen);
     }
 
     // Random initial values for div_vec_c and opt_out2
     for (int k = 0; k < nlev; ++k) {
-      this->div_vec_c[div_vec_c_at(i, k, 0)] = real_distrib(gen);
-      this->opt_out2[div_vec_c_at(i, k, 0)] =
-          real_distrib(gen); // Not used but initialize anyway
+      div_vec_c_h[div_vec_c_at(i, k, 0)] = static_cast<TypeParam>(0.0);
+      opt_out2_h[div_vec_c_at(i, k, 0)] = static_cast<TypeParam>(0.0); // Not used but initialize anyway
     }
   }
 
+  // Copy initialized data to device
+  Kokkos::deep_copy(this->vec_e, vec_e_h);
+  Kokkos::deep_copy(this->cell_neighbor_idx, cell_neighbor_idx_h);
+  Kokkos::deep_copy(this->cell_neighbor_blk, cell_neighbor_blk_h);
+  Kokkos::deep_copy(this->cell_edge_idx, cell_edge_idx_h);
+  Kokkos::deep_copy(this->cell_edge_blk, cell_edge_blk_h);
+  Kokkos::deep_copy(this->geofac_div, geofac_div_h);
+  Kokkos::deep_copy(this->avg_coeff, avg_coeff_h);
+  Kokkos::deep_copy(this->div_vec_c, div_vec_c_h);
+  Kokkos::deep_copy(this->opt_in2, opt_in2_h);
+  Kokkos::deep_copy(this->opt_out2, opt_out2_h);
+
   // Call the div_avg function with l2fields=false
   div_avg<TypeParam>(
       this->vec_e.data(), this->cell_neighbor_idx.data(),
@@ -1092,6 +1194,10 @@ TYPED_TEST(HorizontalDivTest, TestDivAvgRandomNoL2fields) {
       this->nproma, patch_id, l_limited_area, l2fields, this->lacc, this->nlev,
       this->nblks_c, this->nblks_e);
 
+  // Copy results back to host for verification
+  Kokkos::deep_copy(div_vec_c_h, this->div_vec_c);
+  Kokkos::deep_copy(opt_out2_h, this->opt_out2);
+
   // Calculate reference values manually
   std::vector<TypeParam> aux_c(dim_combine(nproma, nlev, nblks_c));
   std::vector<TypeParam> ref_div_vec_c(dim_combine(nproma, nlev, nblks_c));
@@ -1105,18 +1211,18 @@ TYPED_TEST(HorizontalDivTest, TestDivAvgRandomNoL2fields) {
     for (int jk = this->slev[0]; jk < this->elev[0]; ++jk) {
       for (int jc = i_startidx; jc < i_endidx; ++jc) {
         aux_c[div_vec_c_at(jc, jk, jb)] =
-            this->vec_e[vec_e_at(
-                this->cell_edge_idx[cell_edge_at(jc, jb, 0)], jk,
-                this->cell_edge_blk[cell_edge_at(jc, jb, 0)])] *
-                this->geofac_div[geofac_div_at(jc, 0, jb)] +
-            this->vec_e[vec_e_at(
-                this->cell_edge_idx[cell_edge_at(jc, jb, 1)], jk,
-                this->cell_edge_blk[cell_edge_at(jc, jb, 1)])] *
-                this->geofac_div[geofac_div_at(jc, 1, jb)] +
-            this->vec_e[vec_e_at(
-                this->cell_edge_idx[cell_edge_at(jc, jb, 2)], jk,
-                this->cell_edge_blk[cell_edge_at(jc, jb, 2)])] *
-                this->geofac_div[geofac_div_at(jc, 2, jb)];
+            vec_e_h[vec_e_at(
+                cell_edge_idx_h[cell_edge_at(jc, jb, 0)], jk,
+                cell_edge_blk_h[cell_edge_at(jc, jb, 0)])] *
+                geofac_div_h[geofac_div_at(jc, 0, jb)] +
+            vec_e_h[vec_e_at(
+                cell_edge_idx_h[cell_edge_at(jc, jb, 1)], jk,
+                cell_edge_blk_h[cell_edge_at(jc, jb, 1)])] *
+                geofac_div_h[geofac_div_at(jc, 1, jb)] +
+            vec_e_h[vec_e_at(
+                cell_edge_idx_h[cell_edge_at(jc, jb, 2)], jk,
+                cell_edge_blk_h[cell_edge_at(jc, jb, 2)])] *
+                geofac_div_h[geofac_div_at(jc, 2, jb)];
       }
     }
   }
@@ -1147,19 +1253,19 @@ TYPED_TEST(HorizontalDivTest, TestDivAvgRandomNoL2fields) {
       for (int jc = i_startidx; jc < i_endidx; ++jc) {
         ref_div_vec_c[div_vec_c_at(jc, jk, jb)] =
             aux_c[div_vec_c_at(jc, jk, jb)] *
-                this->avg_coeff[avg_coeff_at(jc, 0, jb)] +
+                avg_coeff_h[avg_coeff_at(jc, 0, jb)] +
             aux_c[div_vec_c_at(
-                this->cell_neighbor_idx[cell_neighbor_at(jc, jb, 0)], jk,
-                this->cell_neighbor_blk[cell_neighbor_at(jc, jb, 0)])] *
-                this->avg_coeff[avg_coeff_at(jc, 1, jb)] +
+                cell_neighbor_idx_h[cell_neighbor_at(jc, jb, 0)], jk,
+                cell_neighbor_blk_h[cell_neighbor_at(jc, jb, 0)])] *
+                avg_coeff_h[avg_coeff_at(jc, 1, jb)] +
             aux_c[div_vec_c_at(
-                this->cell_neighbor_idx[cell_neighbor_at(jc, jb, 1)], jk,
-                this->cell_neighbor_blk[cell_neighbor_at(jc, jb, 1)])] *
-                this->avg_coeff[avg_coeff_at(jc, 2, jb)] +
+                cell_neighbor_idx_h[cell_neighbor_at(jc, jb, 1)], jk,
+                cell_neighbor_blk_h[cell_neighbor_at(jc, jb, 1)])] *
+                avg_coeff_h[avg_coeff_at(jc, 2, jb)] +
             aux_c[div_vec_c_at(
-                this->cell_neighbor_idx[cell_neighbor_at(jc, jb, 2)], jk,
-                this->cell_neighbor_blk[cell_neighbor_at(jc, jb, 2)])] *
-                this->avg_coeff[avg_coeff_at(jc, 3, jb)];
+                cell_neighbor_idx_h[cell_neighbor_at(jc, jb, 2)], jk,
+                cell_neighbor_blk_h[cell_neighbor_at(jc, jb, 2)])] *
+                avg_coeff_h[avg_coeff_at(jc, 3, jb)];
       }
     }
   }
@@ -1168,10 +1274,9 @@ TYPED_TEST(HorizontalDivTest, TestDivAvgRandomNoL2fields) {
   // isn't updated
   for (int i = 0; i < nproma; ++i) {
     for (int k = 0; k < nlev; ++k) {
-      EXPECT_NEAR(this->div_vec_c[div_vec_c_at(i, k, 0)],
+      EXPECT_NEAR(div_vec_c_h[div_vec_c_at(i, k, 0)],
                   ref_div_vec_c[div_vec_c_at(i, k, 0)], 1e-5)
           << "div_vec_c results differ at i=" << i << ", k=" << k;
     }
   }
 }
-*/