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

@@ -213,3 +213,80 @@ void grad_fe_cell_adv_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).
+// 2D version, i.e. for a single vertical level
+// LITERATURE:
+// Fish. J and T. Belytschko, 2007: A first course in finite elements,
+//                                  John Wiley and Sons
+template <typename T>
+void grad_fe_cell_adv_2d_lib(const T* p_cc,
+                             const int* cell_neighbor_idx, const int* cell_neighbor_blk,
+                             const T* gradc_bmat,
+                             T* p_grad,
+                             // Integer parameters
+                             int i_startblk, int i_endblk,
+                             int i_startidx_in, int i_endidx_in,
+                             int nproma, int patch_id,
+                             // Extra array dimension
+                             int nblks_c,
+                             // OpenACC flag
+                             bool lacc)
+{
+    // Host or device check
+
+    // Wrap raw pointers in unmanaged Kokkos views
+    typedef Kokkos::View<const T**, Kokkos::LayoutLeft, Kokkos::MemoryUnmanaged> UnmanagedConstT2D;
+    typedef Kokkos::View<const int***, Kokkos::LayoutLeft, Kokkos::MemoryUnmanaged> UnmanagedConstInt3D;
+    typedef Kokkos::View<const T****, Kokkos::LayoutLeft, Kokkos::MemoryUnmanaged> UnmanagedConstT4D;
+    typedef Kokkos::View<T***, Kokkos::LayoutLeft, Kokkos::MemoryUnmanaged> UnmanagedT3D;
+
+    UnmanagedConstT2D p_cc_view(p_cc, nproma, 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);
+    UnmanagedT3D p_grad_view(p_grad, 2, nproma, 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::RangePolicy innerPolicy(i_startidx, i_endidx);
+        Kokkos::parallel_for("grad_fe_cell_adv_2d_lib_inner", innerPolicy,
+                             KOKKOS_LAMBDA(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 precomputed grid geometry factor.
+
+            // zonal(u)-component of gradient
+            p_grad_view(0, jc, jb) =
+                gradc_bmat_view(jc, 0, 0, jb) *
+                p_cc_view(cell_neighbor_idx_view(jc, jb, 0), cell_neighbor_blk_view(jc, jb, 0)) +
+                gradc_bmat_view(jc, 0, 1, jb) *
+                p_cc_view(cell_neighbor_idx_view(jc, jb, 1), cell_neighbor_blk_view(jc, jb, 1)) +
+                gradc_bmat_view(jc, 0, 2, jb) *
+                p_cc_view(cell_neighbor_idx_view(jc, jb, 2), cell_neighbor_blk_view(jc, jb, 2));
+
+            // meridional(v)-component of gradient
+            p_grad_view(1, jc, jb) =
+                gradc_bmat_view(jc, 1, 0, jb) *
+                p_cc_view(cell_neighbor_idx_view(jc, jb, 0), cell_neighbor_blk_view(jc, jb, 0)) +
+                gradc_bmat_view(jc, 1, 1, jb) *
+                p_cc_view(cell_neighbor_idx_view(jc, jb, 1), cell_neighbor_blk_view(jc, jb, 1)) +
+                gradc_bmat_view(jc, 1, 2, jb) *
+                p_cc_view(cell_neighbor_idx_view(jc, jb, 2), cell_neighbor_blk_view(jc, jb, 2));
+        });
+        Kokkos::fence();
+    }
+}