No mesh 2

cpputil/usid/.clang-format

@@ -0,0 +1,60 @@
+---
+Language:        Cpp
+# BasedOnStyle:  LLVM
+AccessModifierOffset: -2
+ConstructorInitializerIndentWidth: 4
+#AlignEscapedNewlinesLeft: false
+AlignTrailingComments: true
+AllowAllParametersOfDeclarationOnNextLine: true
+AllowShortBlocksOnASingleLine: false
+AllowShortIfStatementsOnASingleLine: false
+AllowShortLoopsOnASingleLine: false
+AllowShortFunctionsOnASingleLine: All
+AlwaysBreakTemplateDeclarations: true
+AlwaysBreakBeforeMultilineStrings: false
+BreakBeforeBinaryOperators: false
+BreakBeforeTernaryOperators: true
+BreakConstructorInitializersBeforeComma: false
+BinPackParameters: false
+BinPackArguments: false
+ColumnLimit:     120
+ConstructorInitializerAllOnOneLineOrOnePerLine: false
+DerivePointerAlignment: false
+ExperimentalAutoDetectBinPacking: false
+IndentCaseLabels: false
+IndentWrappedFunctionNames: false
+IndentFunctionDeclarationAfterType: false
+MaxEmptyLinesToKeep: 1
+KeepEmptyLinesAtTheStartOfBlocks: true
+NamespaceIndentation: All
+ObjCSpaceAfterProperty: false
+ObjCSpaceBeforeProtocolList: true
+PenaltyBreakBeforeFirstCallParameter: 19
+PenaltyBreakComment: 300
+PenaltyBreakString: 1000
+PenaltyBreakFirstLessLess: 120
+PenaltyExcessCharacter: 1000000
+PenaltyReturnTypeOnItsOwnLine: 60
+PointerAlignment: Right
+SpacesBeforeTrailingComments: 1
+Cpp11BracedListStyle: true
+Standard:        Cpp11
+IndentWidth:     4
+TabWidth:        8
+UseTab:          Never
+BreakBeforeBraces: Attach
+SpacesInParentheses: false
+SpacesInAngles:  false
+SpaceInEmptyParentheses: false
+SpacesInCStyleCastParentheses: false
+SpacesInContainerLiterals: true
+SpaceBeforeAssignmentOperators: true
+ContinuationIndentWidth: 4
+CommentPragmas:  '*'
+ForEachMacros:   [ foreach, Q_FOREACH, BOOST_FOREACH ]
+SpaceBeforeParens: ControlStatements
+DisableFormat:   false
+AlignAfterOpenBracket: false
+AlignEscapedNewlinesLeft: true
+...
+

cpputil/usid/CMakeLists.txt

@@ -0,0 +1,56 @@
+cmake_minimum_required(VERSION 3.14.5)
+project(usid LANGUAGES CXX)
+
+set(CMAKE_CXX_STANDARD 17)
+
+add_compile_options(-Wall)
+
+include(FetchContent)
+FetchContent_Declare(GridTools
+    GIT_REPOSITORY https://github.com/GridTools/gridtools.git
+    GIT_TAG        master
+)
+FetchContent_MakeAvailable(GridTools)
+
+add_library(usid_naive_helper INTERFACE)
+add_library(usid::usid_naive_helper ALIAS usid_naive_helper)
+target_include_directories(usid_naive_helper INTERFACE ${PROJECT_SOURCE_DIR}/include)
+target_link_libraries(usid_naive_helper INTERFACE GridTools::gridtools)
+
+include(CheckLanguage)
+check_language(CUDA)
+if(CMAKE_CUDA_COMPILER)
+    enable_language(CUDA)
+
+    add_library(usid_cuda_helper INTERFACE)
+    add_library(usid::usid_cuda_helper ALIAS usid_cuda_helper)
+    target_include_directories(usid_cuda_helper INTERFACE ${PROJECT_SOURCE_DIR}/include)
+    target_link_libraries(usid_cuda_helper INTERFACE GridTools::gridtools GridTools::stencil_gpu)
+endif()
+
+set(BUILD_SHARED_LIBS OFF)
+
+FetchContent_Declare(GoogleTest
+    GIT_REPOSITORY https://github.com/google/googletest.git
+    GIT_TAG        release-1.10.0
+)
+FetchContent_MakeAvailable(GoogleTest)
+
+find_package(eckit REQUIRED)
+find_package(Atlas REQUIRED)
+
+add_library(fvm_nabla_driver INTERFACE)
+add_library(usid::fvm_nabla_driver ALIAS fvm_nabla_driver)
+target_include_directories(fvm_nabla_driver INTERFACE ${PROJECT_SOURCE_DIR}/include)
+target_link_libraries(fvm_nabla_driver INTERFACE atlas eckit GridTools::gridtools gtest)
+if(CMAKE_CUDA_COMPILER)
+    target_link_libraries(fvm_nabla_driver INTERFACE GridTools::stencil_gpu)
+endif()
+
+
+if(CMAKE_PROJECT_NAME STREQUAL PROJECT_NAME)
+    include(CTest)
+    if(BUILD_TESTING)
+        add_subdirectory(tests)
+    endif()
+endif()

cpputil/usid/include/gridtools/usid/atlas.hpp

@@ -0,0 +1,19 @@
+#pragma once
+
+#include <atlas/mesh/Connectivity.h>
+
+#include <gridtools/storage/builder.hpp>
+
+namespace atlas::mesh {
+template <class Connectivity, class MaxNeighbors>
+auto make_storage_producer(MaxNeighbors max_neighbors,
+                           Connectivity const &src) {
+  return [&src, max_neighbors](auto traits) {
+    return gridtools::storage::builder<decltype(traits)>
+            .template type<int>()
+            .dimensions(src.rows(), max_neighbors)
+            .initializer([&src](auto row, auto col) { return col < src.cols(row) ? src(row, col) : -1; })
+            .build();
+  };
+}
+} // namespace atlas::mesh

cpputil/usid/include/gridtools/usid/cuda_helpers.hpp

@@ -0,0 +1,52 @@
+#pragma once
+
+#ifndef __CUDACC__
+#error Tried to compile CUDA code with a regular C++ compiler.
+#endif
+
+#include <gridtools/common/cuda_util.hpp>
+#include <gridtools/common/hymap.hpp>
+#include <gridtools/common/pair.hpp>
+#include <gridtools/common/tuple_util.hpp>
+#include <gridtools/sid/allocator.hpp>
+#include <gridtools/sid/concept.hpp>
+#include <gridtools/storage/gpu.hpp>
+
+#include "dim.hpp"
+#include "helpers.hpp"
+
+namespace gridtools::usid::cuda {
+using traits_t = storage::gpu;
+
+inline auto make_allocator() {
+  return sid::device::make_cached_allocator(&cuda_util::cuda_malloc<char[]>);
+}
+
+template <class Kernel, class Ptr, class Strides, class... Neighbors>
+__global__ void kernel(int h_size, Ptr ptr_holder, Strides strides,
+                       Neighbors... neighbors) {
+  auto h = blockIdx.x * blockDim.x + threadIdx.x;
+  if (h >= h_size)
+    return;
+  Kernel()()(
+      sid::shifted(ptr_holder(), device::at_key<dim::h>(strides), h), strides,
+      tuple_util::device::make<pair>(neighbors.first(), neighbors.second)...);
+}
+
+template <class Kernel, class Size, class Sid, class... Sids>
+void call_kernel(Size size, Sid &&fields, Sids &&...neighbor_fields) {
+  int threads_per_block = 32;
+  int blocks = (size + threads_per_block - 1) / threads_per_block;
+  kernel<Kernel><<<blocks, threads_per_block>>>(
+      size, sid::get_origin(fields), sid::get_strides(fields),
+      tuple_util::make<pair>(
+          sid::get_origin(neighbor_fields),
+          at_key<dim::h>(sid::get_strides(neighbor_fields)))...);
+  GT_CUDA_CHECK(cudaGetLastError());
+}
+
+template <class Tag, class Ptr>
+__device__ decltype(auto) field(Ptr const &ptr) {
+  return *device::at_key<Tag>(ptr);
+}
+} // namespace gridtools::usid::cuda

cpputil/usid/include/gridtools/usid/dim.hpp

@@ -0,0 +1,16 @@
+#pragma once
+
+#include <gridtools/common/defs.hpp>
+#include <gridtools/common/integral_constant.hpp>
+
+namespace gridtools::usid::dim {
+using horizontal = integral_constant<int_t, 0>;
+using vertical = integral_constant<int_t, 1>;
+using neighbor = integral_constant<int_t, 1>;
+using sparse = integral_constant<int_t, 2>;
+
+using h = horizontal;
+using k = vertical;
+using n = neighbor;
+using s = sparse;
+} // namespace gridtools::usid::dim

cpputil/usid/include/gridtools/usid/domain.hpp

@@ -0,0 +1,10 @@
+#pragma once
+
+namespace gridtools::usid {
+struct domain {
+  int vertex;
+  int edge;
+  int cell;
+  int k;
+};
+} // namespace gridtools::usid

cpputil/usid/include/gridtools/usid/helpers.hpp

@@ -0,0 +1,115 @@
+#pragma once
+
+#include <limits>
+#include <type_traits>
+#include <utility>
+
+#include <gridtools/common/defs.hpp>
+#include <gridtools/common/host_device.hpp>
+#include <gridtools/common/hymap.hpp>
+#include <gridtools/common/integral_constant.hpp>
+#include <gridtools/common/utility.hpp>
+#include <gridtools/meta/id.hpp>
+#include <gridtools/sid/composite.hpp>
+#include <gridtools/sid/concept.hpp>
+#include <gridtools/sid/contiguous.hpp>
+#include <gridtools/sid/loop.hpp>
+#include <gridtools/sid/rename_dimensions.hpp>
+
+#include "dim.hpp"
+#include "domain.hpp"
+
+namespace gridtools::usid {
+    template <class T, class Alloc, class HSize, class KSize>
+    auto make_simple_tmp_storage(HSize h_size, KSize k_size, Alloc &alloc) {
+        return sid::make_contiguous<T>(alloc, hymap::keys<dim::h, dim::k>::values<HSize, KSize>(h_size, k_size));
+    }
+
+    template <int_t N, bool HasSkipValues>
+    struct connectivity {
+        using max_neighbors_t = integral_constant<int_t, N>;
+        using has_skip_values_t = std::bool_constant<HasSkipValues>;
+    };
+
+    template <class Connectivity>
+    struct sparse_field {
+        using connectivity_t = Connectivity;
+    };
+
+    template <class, class = void>
+    struct is_sparse_field : std::false_type {};
+
+    template <class T>
+    struct is_sparse_field<T, std::enable_if_t<std::is_base_of_v<sparse_field<typename T::connectivity_t>, T>>>
+        : std::true_type {};
+
+    template <class T, class Conncectivity, class F, class Init, class G, class Ptr, class Strides, class Neighbors>
+    GT_FUNCTION T fold_neighbors(F f, Init init, G g, Ptr &&ptr, Strides &&strides, Neighbors &&neighbors) {
+        T acc = init(meta::lazy::id<T>());
+        sid::make_loop<Conncectivity>(typename Conncectivity::max_neighbors_t())([&](auto const &ptr, auto &&) {
+            auto i = *host_device::at_key<Conncectivity>(ptr);
+            if constexpr (Conncectivity::has_skip_values_t::value)
+                if (i < 0)
+                    return;
+            acc = f(acc, g(ptr, sid::shifted(neighbors.first, neighbors.second, i)));
+        })(wstd::forward<decltype(ptr)>(ptr), strides);
+        return acc;
+    }
+
+    template <class T, class Conncectivity, class F, class... Args>
+    GT_FUNCTION T sum_neighbors(F f, Args &&...args) {
+        return fold_neighbors<T, Conncectivity>([](auto x, auto y) { return x + y; },
+            [](auto z) -> typename decltype(z)::type { return 0; },
+            f,
+            wstd::forward<Args>(args)...);
+    }
+
+    template <class T, class Conncectivity, class F, class... Args>
+    GT_FUNCTION T product_neighbors(F f, Args &&...args) {
+        return fold_neighbors<T, Conncectivity>([](auto x, auto y) { return x * y; },
+            [](auto z) -> typename decltype(z)::type { return 1; },
+            f,
+            wstd::forward<Args>(args)...);
+    }
+
+    template <class T, class Conncectivity, class F, class... Args>
+    GT_FUNCTION T min_neighbors(F f, Args &&...args) {
+        return fold_neighbors<T, Conncectivity>([](auto x, auto y) { return x < y ? x : y; },
+            [](auto z) -> typename decltype(z)::type {
+                constexpr auto res = std::numeric_limits<typename decltype(z)::type>::max();
+                return res;
+            },
+            f,
+            wstd::forward<Args>(args)...);
+    }
+
+    template <class T, class Conncectivity, class F, class... Args>
+    GT_FUNCTION T max_neighbors(F f, Args &&...args) {
+        return fold_neighbors<T, Conncectivity>([](auto x, auto y) { return x > y ? x : y; },
+            [](auto z) -> typename decltype(z)::type {
+                constexpr auto res = std::numeric_limits<typename decltype(z)::type>::min();
+                return res;
+            },
+            f,
+            wstd::forward<Args>(args)...);
+    }
+
+    template <class Tag, class Val>
+    decltype(auto) composite_item(Tag, Val &&val) {
+        if constexpr (is_sparse_field<Tag>::value)
+            return sid::rename_dimensions<dim::s, typename Tag::connectivity_t>(std::forward<Val>(val));
+        else
+            return std::forward<Val>(val);
+    }
+
+    template <class... Tags>
+    struct make_composite_f {
+        template <class... Vals>
+        auto operator()(Vals &&...vals) const {
+            return sid::composite::make<Tags...>(composite_item(Tags(), std::forward<Vals>(vals))...);
+        }
+    };
+
+    template <class... Tags>
+    constexpr make_composite_f<Tags...> make_composite = {};
+} // namespace gridtools::usid

cpputil/usid/include/gridtools/usid/naive_helpers.hpp

@@ -0,0 +1,38 @@
+#pragma once
+
+#include <memory>
+#include <tuple>
+#include <utility>
+
+#include <gridtools/common/hymap.hpp>
+#include <gridtools/sid/allocator.hpp>
+#include <gridtools/sid/concept.hpp>
+#include <gridtools/sid/loop.hpp>
+#include <gridtools/storage/cpu_ifirst.hpp>
+
+#include "dim.hpp"
+#include "helpers.hpp"
+
+namespace gridtools::usid::naive {
+using traits_t = storage::cpu_ifirst;
+
+inline auto make_allocator() {
+  return sid::make_cached_allocator(&std::make_unique<char[]>);
+}
+
+template <class Kernel, class Size, class Sid, class... Sids>
+void call_kernel(Size size, Sid &&fields, Sids &&...neighbor_fields) {
+  sid::make_loop<dim::h>(size)(
+      [params = std::make_tuple(std::make_pair(
+           sid::get_origin(neighbor_fields)(),
+           sid::get_stride<dim::h>(sid::get_strides(neighbor_fields)))...)](
+          auto &ptr, auto const &strides) {
+        std::apply(Kernel()(),
+                   std::tuple_cat(std::forward_as_tuple(ptr, strides), params));
+      })(sid::get_origin(fields)(), sid::get_strides(fields));
+}
+
+template <class Tag, class Ptr> decltype(auto) field(Ptr const &ptr) {
+  return *at_key<Tag>(ptr);
+}
+} // namespace gridtools::usid::naive