dawn-ico · Stagno · Mar 17, 2020 · Mar 17, 2020
diff --git a/.gitignore b/.gitignore
@@ -0,0 +1 @@
+build*
diff --git a/README.md b/README.md
@@ -6,7 +6,7 @@ This repoistory contains some utilities to transform Atlas meshes as well as the
 
 Requirements:
 
-* libnetcdc-c++4
+* libnetcdf-c++4-dev
 * Atlas
 * eckit
 
@@ -60,8 +60,8 @@ python3 convergence_plot.py ./(mylib|atlas)IconLaplaceDriver
 
 Various error norms are collected for divergence, curl and (normal) vector Laplacian. 
 
-Various Octave scripts are provided to assist in visualizing debug and output data. For example, to get a convergence plot:
+Various Octave scripts (inside `scripts` directory) are provided to assist in visualizing debug and output data. For example, to get a convergence plot:
 
 ```
-octave:1> convergence_plot_dec('<path/to/conv.csv') 
+octave:1> convergence_plot_dec('path/to/conv.csv') 
 ```
diff --git a/scripts/convergence_plot_dec.m b/scripts/convergence_plot_dec.m
@@ -12,7 +12,7 @@ function convergence_plot_dec(fname)
   hcb = legend(['L_{\infty}';'L_1';'L_2';'lin. conv';'quad. conv']);
   xlabel('log(\Delta x)[°]')
   ylabel('log(L) [-]')
-  set(gca,'FontSize',50)
-  set(hcb,'FontSize',50)
+  set(gca,'FontSize',20)
+  set(hcb,'FontSize',20)
   ylim([-6,0]);
  endfunction
diff --git a/stencils/CMakeLists.txt b/stencils/CMakeLists.txt
@@ -1,8 +1,11 @@
 add_executable(atlasIconLaplaceDriver atlasIconLaplaceDriver.cpp)
 target_link_libraries(atlasIconLaplaceDriver atlas eckit atlasUtilsLib)
+target_include_directories(atlasIconLaplaceDriver PRIVATE ${CMAKE_CURRENT_SOURCE_DIR})
 
 add_executable(atlasShallowWater shallowWater.cpp)
 target_link_libraries(atlasShallowWater atlas eckit atlasUtilsLib)
+target_include_directories(atlasShallowWater PRIVATE ${CMAKE_CURRENT_SOURCE_DIR})
 
 add_executable(mylibIconLaplaceDriver mylibIconLaplaceDriver.cpp)
-target_link_libraries(mylibIconLaplaceDriver atlasUtilsLib myLib)
+target_link_libraries(mylibIconLaplaceDriver atlasUtilsLib myLib)
+target_include_directories(mylibIconLaplaceDriver PRIVATE ${CMAKE_CURRENT_SOURCE_DIR})
diff --git a/stencils/atlasIconLaplaceDriver.cpp b/stencils/atlasIconLaplaceDriver.cpp
@@ -115,7 +115,7 @@ int main(int argc, char const* argv[]) {
 
   // dump a whole bunch of debug output (meant to be visualized using Octave, but gnuplot and the
   // like will certainly work too)
-  const bool dbg_out = false;
+  const bool dbg_out = true;
   const bool readMeshFromDisk = false;
 
   atlas::Mesh mesh;
@@ -490,8 +490,8 @@ int main(int argc, char const* argv[]) {
       primal_edge_length, dual_edge_length, tangent_orientation, geofac_rot, geofac_div)
       .run();
   clock_t end = clock();
-  std::cout << "run time Laplacian at resolution " << w << " "
-            << (end - start) / (double CLOCKS_PER_SEC) << "\n";
+  // std::cout << "run time Laplacian at resolution " << w << " "
+  //           << (end - start) / (double CLOCKS_PER_SEC) << "\n";
 
   if(dbg_out) {
     dumpEdgeField("laplICONatlas_nabla2t1.txt", mesh, wrapper, nabla2t1_vec, level,

diff --git a/stencils/interfaces/defs.hpp → stencils/driver-includes/defs.hpp b/stencils/interfaces/defs.hpp → stencils/driver-includes/defs.hpp
diff --git a/stencils/driver-includes/extent.hpp b/stencils/driver-includes/extent.hpp
@@ -0,0 +1,58 @@
+//===--------------------------------------------------------------------------------*- C++ -*-===//
+//                          _
+//                         | |
+//                       __| | __ ___      ___ ___
+//                      / _` |/ _` \ \ /\ / / '_  |
+//                     | (_| | (_| |\ V  V /| | | |
+//                      \__,_|\__,_| \_/\_/ |_| |_| - Compiler Toolchain
+//
+//
+//  This file is distributed under the MIT License (MIT).
+//  See LICENSE.txt for details.
+//
+//===------------------------------------------------------------------------------------------===//
+
+#pragma once
+
+#include <array>
+
+namespace dawn {
+namespace driver {
+
+// Access dimensions with `auto dim_extent = get<D>();`, where D = 0,1,2 (=i,j,k)
+// and minus/plus component with `get<C>(dim_extent)`, where C = 0 (minus), 1 (plus)
+// (Note, need `using std::get` for ADL to kick-in if std < c++20)
+using cartesian_extent = std::array<std::array<int, 2>, 3>;
+
+// Access dimensions with get<D>(), where D = 0,1 (=horizontal, vertical)
+// for D == 0 it is either true (have extent) or false (no extent)
+// for D == 1, see Cartesian case.
+struct unstructured_extent {
+  bool horizontal;
+  std::array<int, 2> vertical;
+};
+
+namespace detail {
+template <int I>
+struct get_impl;
+
+template <>
+struct get_impl<0> {
+  using type = bool;
+  constexpr type operator()(unstructured_extent const& extent) { return extent.horizontal; }
+};
+template <>
+struct get_impl<1> {
+  using type = std::array<int, 2>;
+  constexpr type operator()(unstructured_extent const& extent) { return extent.vertical; }
+};
+
+} // namespace detail
+
+template <int I>
+constexpr auto get(unstructured_extent const& extent) -> decltype(detail::get_impl<I>{}(extent)) {
+  return detail::get_impl<I>{}(extent);
+}
+
+} // namespace driver
+} // namespace dawn
diff --git a/...ils/interfaces/unstructured_interface.hpp → ...river-includes/unstructured_interface.hpp b/...ils/interfaces/unstructured_interface.hpp → ...river-includes/unstructured_interface.hpp
@@ -15,6 +15,7 @@
 #pragma once
 
 #include "defs.hpp"
+#include "extent.hpp"
 
 namespace dawn {
 

diff --git a/stencils/generated_iconLaplace.hpp b/stencils/generated_iconLaplace.hpp
@@ -1,12 +1,18 @@
+//---- Preprocessor defines ----
 #define DAWN_GENERATED 1
+#undef DAWN_BACKEND_T
 #define DAWN_BACKEND_T CXXNAIVEICO
-#include "interfaces/unstructured_interface.hpp"
+#include <driver-includes/unstructured_interface.hpp>
+
+//---- Globals ----
+
+//---- Stencils ----
 namespace dawn_generated {
 namespace cxxnaiveico {
 template <typename LibTag>
 class icon {
 private:
-  struct stencil_399 {
+  struct stencil_64 {
     dawn::mesh_t<LibTag> const& m_mesh;
     int m_k_size;
     dawn::edge_field_t<LibTag, double>& m_vec;
@@ -22,91 +28,109 @@ class icon {
     dawn::sparse_cell_field_t<LibTag, double>& m_geofac_div;
 
   public:
-    stencil_399(dawn::mesh_t<LibTag> const& mesh, int k_size,
-                dawn::edge_field_t<LibTag, double>& vec,
-                dawn::cell_field_t<LibTag, double>& div_vec,
-                dawn::vertex_field_t<LibTag, double>& rot_vec,
-                dawn::edge_field_t<LibTag, double>& nabla2t1_vec,
-                dawn::edge_field_t<LibTag, double>& nabla2t2_vec,
-                dawn::edge_field_t<LibTag, double>& nabla2_vec,
-                dawn::edge_field_t<LibTag, double>& primal_edge_length,
-                dawn::edge_field_t<LibTag, double>& dual_edge_length,
-                dawn::edge_field_t<LibTag, double>& tangent_orientation,
-                dawn::sparse_vertex_field_t<LibTag, double>& geofac_rot,
-                dawn::sparse_cell_field_t<LibTag, double>& geofac_div)
+    stencil_64(dawn::mesh_t<LibTag> const& mesh, int k_size,
+               dawn::edge_field_t<LibTag, double>& vec, dawn::cell_field_t<LibTag, double>& div_vec,
+               dawn::vertex_field_t<LibTag, double>& rot_vec,
+               dawn::edge_field_t<LibTag, double>& nabla2t1_vec,
+               dawn::edge_field_t<LibTag, double>& nabla2t2_vec,
+               dawn::edge_field_t<LibTag, double>& nabla2_vec,
+               dawn::edge_field_t<LibTag, double>& primal_edge_length,
+               dawn::edge_field_t<LibTag, double>& dual_edge_length,
+               dawn::edge_field_t<LibTag, double>& tangent_orientation,
+               dawn::sparse_vertex_field_t<LibTag, double>& geofac_rot,
+               dawn::sparse_cell_field_t<LibTag, double>& geofac_div)
         : m_mesh(mesh), m_k_size(k_size), m_vec(vec), m_div_vec(div_vec), m_rot_vec(rot_vec),
           m_nabla2t1_vec(nabla2t1_vec), m_nabla2t2_vec(nabla2t2_vec), m_nabla2_vec(nabla2_vec),
           m_primal_edge_length(primal_edge_length), m_dual_edge_length(dual_edge_length),
           m_tangent_orientation(tangent_orientation), m_geofac_rot(geofac_rot),
           m_geofac_div(geofac_div) {}
 
-    ~stencil_399() {}
+    ~stencil_64() {}
 
     void sync_storages() {}
+    static constexpr dawn::driver::unstructured_extent vec_extent = {false, 0, 0};
+    static constexpr dawn::driver::unstructured_extent div_vec_extent = {false, 0, 0};
+    static constexpr dawn::driver::unstructured_extent rot_vec_extent = {false, 0, 0};
+    static constexpr dawn::driver::unstructured_extent nabla2t1_vec_extent = {false, 0, 0};
+    static constexpr dawn::driver::unstructured_extent nabla2t2_vec_extent = {false, 0, 0};
+    static constexpr dawn::driver::unstructured_extent nabla2_vec_extent = {false, 0, 0};
+    static constexpr dawn::driver::unstructured_extent primal_edge_length_extent = {false, 0, 0};
+    static constexpr dawn::driver::unstructured_extent dual_edge_length_extent = {false, 0, 0};
+    static constexpr dawn::driver::unstructured_extent tangent_orientation_extent = {false, 0, 0};
+    static constexpr dawn::driver::unstructured_extent geofac_rot_extent = {false, 0, 0};
+    static constexpr dawn::driver::unstructured_extent geofac_div_extent = {false, 0, 0};
 
     void run() {
       using dawn::deref;
-      ;
       {
         for(int k = 0 + 0; k <= (m_k_size == 0 ? 0 : (m_k_size - 1)) + 0 + 0; ++k) {
-          for(auto const& loc : getVertices(LibTag{}, m_mesh)) {
+          for(auto const& loc : getCells(LibTag{}, m_mesh)) {
             int m_sparse_dimension_idx = 0;
-            m_rot_vec(deref(LibTag{}, loc), k + 0) = reduceEdgeToVertex(
-                LibTag{}, m_mesh, loc, (::dawn::float_type)0.000000,
+            m_div_vec(deref(LibTag{}, loc), k + 0) = reduceEdgeToCell(
+                LibTag{}, m_mesh, loc, (::dawn::float_type)0.0,
                 [&](auto& lhs, auto const& red_loc) {
                   lhs += (m_vec(deref(LibTag{}, red_loc), k + 0) *
-                          m_geofac_rot(deref(LibTag{}, loc), m_sparse_dimension_idx, k + 0));
+                          m_geofac_div(deref(LibTag{}, loc), m_sparse_dimension_idx, k + 0));
                   m_sparse_dimension_idx++;
                   return lhs;
                 });
           }
-          for(auto const& loc : getCells(LibTag{}, m_mesh)) {
+          for(auto const& loc : getEdges(LibTag{}, m_mesh)) {
             int m_sparse_dimension_idx = 0;
-            m_div_vec(deref(LibTag{}, loc), k + 0) = reduceEdgeToCell(
-                LibTag{}, m_mesh, loc, (::dawn::float_type)0.000000,
+            m_nabla2t2_vec(deref(LibTag{}, loc), k + 0) =
+                reduceCellToEdge(LibTag{}, m_mesh, loc, (::dawn::float_type)0.0,
+                                 [&](auto& lhs, auto const& red_loc, auto const& weight) {
+                                   lhs += weight * m_div_vec(deref(LibTag{}, red_loc), k + 0);
+                                   m_sparse_dimension_idx++;
+                                   return lhs;
+                                 },
+                                 std::vector<double>({-1.000000, 1.000000}));
+          }
+          for(auto const& loc : getEdges(LibTag{}, m_mesh)) {
+            m_nabla2t2_vec(deref(LibTag{}, loc), k + 0) =
+                (m_nabla2t2_vec(deref(LibTag{}, loc), k + 0) /
+                 m_dual_edge_length(deref(LibTag{}, loc), k + 0));
+          }
+          for(auto const& loc : getVertices(LibTag{}, m_mesh)) {
+            int m_sparse_dimension_idx = 0;
+            m_rot_vec(deref(LibTag{}, loc), k + 0) = reduceEdgeToVertex(
+                LibTag{}, m_mesh, loc, (::dawn::float_type)0.0,
                 [&](auto& lhs, auto const& red_loc) {
                   lhs += (m_vec(deref(LibTag{}, red_loc), k + 0) *
-                          m_geofac_div(deref(LibTag{}, loc), m_sparse_dimension_idx, k + 0));
+                          m_geofac_rot(deref(LibTag{}, loc), m_sparse_dimension_idx, k + 0));
                   m_sparse_dimension_idx++;
                   return lhs;
                 });
           }
           for(auto const& loc : getEdges(LibTag{}, m_mesh)) {
             int m_sparse_dimension_idx = 0;
-            m_nabla2t1_vec(deref(LibTag{}, loc), k + 0) = reduceVertexToEdge(
-                LibTag{}, m_mesh, loc, (::dawn::float_type)0.000000,
-                [&](auto& lhs, auto const& red_loc, auto const& weight) {
-                  lhs += weight * m_rot_vec(deref(LibTag{}, red_loc), k + 0);
-                  m_sparse_dimension_idx++;
-                  return lhs;
-                },
-                std::vector<double>({-1.000000, 1.000000}));
+            m_nabla2t1_vec(deref(LibTag{}, loc), k + 0) =
+                reduceVertexToEdge(LibTag{}, m_mesh, loc, (::dawn::float_type)0.0,
+                                   [&](auto& lhs, auto const& red_loc, auto const& weight) {
+                                     lhs += weight * m_rot_vec(deref(LibTag{}, red_loc), k + 0);
+                                     m_sparse_dimension_idx++;
+                                     return lhs;
+                                   },
+                                   std::vector<double>({-1.000000, 1.000000}));
+          }
+          for(auto const& loc : getEdges(LibTag{}, m_mesh)) {
             m_nabla2t1_vec(deref(LibTag{}, loc), k + 0) =
                 ((m_tangent_orientation(deref(LibTag{}, loc), k + 0) *
                   m_nabla2t1_vec(deref(LibTag{}, loc), k + 0)) /
                  m_primal_edge_length(deref(LibTag{}, loc), k + 0));
-            m_nabla2t2_vec(deref(LibTag{}, loc), k + 0) = reduceCellToEdge(
-                LibTag{}, m_mesh, loc, (::dawn::float_type)0.000000,
-                [&](auto& lhs, auto const& red_loc, auto const& weight) {
-                  lhs += weight * m_div_vec(deref(LibTag{}, red_loc), k + 0);
-                  m_sparse_dimension_idx++;
-                  return lhs;
-                },
-                std::vector<double>({-1.000000, 1.000000}));
-            m_nabla2t2_vec(deref(LibTag{}, loc), k + 0) =
-                (m_nabla2t2_vec(deref(LibTag{}, loc), k + 0) /
-                 m_dual_edge_length(deref(LibTag{}, loc), k + 0));
+          }
+          for(auto const& loc : getEdges(LibTag{}, m_mesh)) {
             m_nabla2_vec(deref(LibTag{}, loc), k + 0) =
-                (-m_nabla2t1_vec(deref(LibTag{}, loc), k + 0) +
-                 m_nabla2t2_vec(deref(LibTag{}, loc), k + 0));
+                (m_nabla2t2_vec(deref(LibTag{}, loc), k + 0) -
+                 m_nabla2t1_vec(deref(LibTag{}, loc), k + 0));
           }
         }
       }
       sync_storages();
     }
   };
   static constexpr const char* s_name = "icon";
-  stencil_399 m_stencil_399;
+  stencil_64 m_stencil_64;
 
 public:
   icon(const icon&) = delete;
@@ -123,12 +147,12 @@ class icon {
        dawn::edge_field_t<LibTag, double>& tangent_orientation,
        dawn::sparse_vertex_field_t<LibTag, double>& geofac_rot,
        dawn::sparse_cell_field_t<LibTag, double>& geofac_div)
-      : m_stencil_399(mesh, k_size, vec, div_vec, rot_vec, nabla2t1_vec, nabla2t2_vec, nabla2_vec,
-                      primal_edge_length, dual_edge_length, tangent_orientation, geofac_rot,
-                      geofac_div) {}
+      : m_stencil_64(mesh, k_size, vec, div_vec, rot_vec, nabla2t1_vec, nabla2t2_vec, nabla2_vec,
+                     primal_edge_length, dual_edge_length, tangent_orientation, geofac_rot,
+                     geofac_div) {}
 
   void run() {
-    m_stencil_399.run();
+    m_stencil_64.run();
     ;
   }
 };

diff --git a/stencils/interfaces/mylib_interface.hpp b/stencils/interfaces/mylib_interface.hpp
@@ -1,8 +1,8 @@
 #ifndef DAWN_PROTOTYPE_MYLIB_INTERFACE_H_
 #define DAWN_PROTOTYPE_MYLIB_INTERFACE_H_
 
+#include "driver-includes/unstructured_interface.hpp"
 #include "mylib.hpp"
-#include "unstructured_interface.hpp"
 #include <functional>
 
 namespace mylibInterface {