Draft: C++ port of horizontal/mo_lib_gradients.F90

Port the functions of mo_lib_gradients.F90 to C++, with C bindings. Untested, but compiles.

src/horizontal/mo_lib_gradients.cpp

@@ -290,3 +290,100 @@ void grad_fe_cell_adv_2d_lib(const T* p_cc,
         Kokkos::fence();
     }
 }
+
+
+// Computes the cell centered gradient in geographical coordinates.
+// The gradient is computed by taking the derivative of the shape functions
+// for a three-node triangular element (Finite Element thinking).
+// Special dycore version, which handles two fields at a time.
+// LITERATURE:
+// Fish. J and T. Belytschko, 2007: A first course in finite elements,
+//                                  John Wiley and Sons
+template <typename T, typename U>
+void grad_fe_cell_dycore_lib(const U* p_ccpr,
+                             const int* cell_neighbor_idx, const int* cell_neighbor_blk,
+                             const T* gradc_bmat,
+                             U* p_grad,
+                             // Integer parameters
+                             int i_startblk, int i_endblk,
+                             int i_startidx_in, int i_endidx_in,
+                             int slev, int elev, int nproma, int patch_id,
+                             // Extra array dimension
+                             int nlev, int nblks_c,
+                             // OpenACC flag
+                             bool lacc)
+{
+    // Host or device check
+
+    // Wrap raw pointers in unmanaged Kokkos views
+    typedef Kokkos::View<const U****, Kokkos::LayoutLeft, Kokkos::MemoryUnmanaged> UnmanagedConstU4D;
+    typedef Kokkos::View<const int***, Kokkos::LayoutLeft, Kokkos::MemoryUnmanaged> UnmanagedConstInt3D;
+    typedef Kokkos::View<const T****, Kokkos::LayoutLeft, Kokkos::MemoryUnmanaged> UnmanagedConstT4D;
+    typedef Kokkos::View<U****, Kokkos::LayoutLeft, Kokkos::MemoryUnmanaged> UnmanagedU4D;
+
+    UnmanagedConstU4D p_ccpr_view(p_ccpr, nproma, 2, nlev, nblks_c);
+    UnmanagedConstInt3D cell_neighbor_idx_view(cell_neighbor_idx, nproma, nblks_c, 3);
+    UnmanagedConstInt3D cell_neighbor_blk_view(cell_neighbor_blk, nproma, nblks_c, 3);
+    UnmanagedConstT4D gradc_bmat_view(gradc_bmat, nproma, 2, 3, nblks_c);
+    UnmanagedU4D p_grad_view(p_grad, 2, nproma, nlev, nblks_c);
+
+    if (patch_id > 1)
+    {
+        // Fill nest boundaries with zeros to avoid trouble with MPI synchronization
+    }
+
+    for (int jb = i_startblk; jb < i_endblk; ++jb)
+    {
+        int i_startidx, i_endidx;
+        get_indices_c_lib(i_startidx_in, i_endidx_in, nproma, jb,
+                          i_startblk, i_endblk, i_startidx, i_endidx);
+
+        Kokkos::MDRangePolicy<Kokkos::Rank<2>> innerPolicy({slev, i_startidx},
+                                                           {elev, i_endidx});
+        Kokkos::parallel_for("grad_fe_cell_dycore_lib_inner", innerPolicy,
+                             KOKKOS_LAMBDA(const int jk, const int jc)
+        {
+            // We do not make use of the intrinsic function DOT_PRODUCT on purpose,
+            // since it is extremely slow on the SX9, when combined with indirect
+            // addressing.
+            // multiply cell-based input values with shape function derivatives
+
+            // zonal(u)-component of gradient, field 0
+            p_grad_view(0, jc, jk, jb) =
+                gradc_bmat_view(jc, 0, 0, jb) *
+                p_ccpr_view(0, cell_neighbor_idx_view(jc, jb, 0), jk, cell_neighbor_blk_view(jc, jb, 0)) +
+                gradc_bmat_view(jc, 0, 1, jb) *
+                p_ccpr_view(0, cell_neighbor_idx_view(jc, jb, 1), jk, cell_neighbor_blk_view(jc, jb, 1)) +
+                gradc_bmat_view(jc, 0, 2, jb) *
+                p_ccpr_view(0, cell_neighbor_idx_view(jc, jb, 2), jk, cell_neighbor_blk_view(jc, jb, 2));
+
+            // meridional(v)-component of gradient, field 0
+            p_grad_view(1, jc, jk, jb) =
+                gradc_bmat_view(jc, 1, 0, jb) *
+                p_ccpr_view(0, cell_neighbor_idx_view(jc, jb, 0), jk, cell_neighbor_blk_view(jc, jb, 0)) +
+                gradc_bmat_view(jc, 1, 1, jb) *
+                p_ccpr_view(0, cell_neighbor_idx_view(jc, jb, 1), jk, cell_neighbor_blk_view(jc, jb, 1)) +
+                gradc_bmat_view(jc, 1, 2, jb) *
+                p_ccpr_view(0, cell_neighbor_idx_view(jc, jb, 2), jk, cell_neighbor_blk_view(jc, jb, 2));
+
+            // zonal(u)-component of gradient, field 1
+            p_grad_view(2, jc, jk, jb) =
+                gradc_bmat_view(jc, 0, 0, jb) *
+                p_ccpr_view(1, cell_neighbor_idx_view(jc, jb, 0), jk, cell_neighbor_blk_view(jc, jb, 0)) +
+                gradc_bmat_view(jc, 0, 1, jb) *
+                p_ccpr_view(1, cell_neighbor_idx_view(jc, jb, 1), jk, cell_neighbor_blk_view(jc, jb, 1)) +
+                gradc_bmat_view(jc, 0, 2, jb) *
+                p_ccpr_view(1, cell_neighbor_idx_view(jc, jb, 2), jk, cell_neighbor_blk_view(jc, jb, 2));
+
+            // meridional(v)-component of gradient, field 1
+            p_grad_view(3, jc, jk, jb) =
+                gradc_bmat_view(jc, 1, 0, jb) *
+                p_ccpr_view(1, cell_neighbor_idx_view(jc, jb, 0), jk, cell_neighbor_blk_view(jc, jb, 0)) +
+                gradc_bmat_view(jc, 1, 1, jb) *
+                p_ccpr_view(1, cell_neighbor_idx_view(jc, jb, 1), jk, cell_neighbor_blk_view(jc, jb, 1)) +
+                gradc_bmat_view(jc, 1, 2, jb) *
+                p_ccpr_view(1, cell_neighbor_idx_view(jc, jb, 2), jk, cell_neighbor_blk_view(jc, jb, 2));
+        });
+        Kokkos::fence();
+    }
+}