Add support for MultiField on GPU

CMakeLists.txt

@@ -18,6 +18,9 @@ find_package( ecbuild 3.4 REQUIRED HINTS ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CUR
 
 project( atlas LANGUAGES CXX )
 
+set(CMAKE_CXX_STANDARD 17)
+set(CMAKE_CXX_STANDARD_REQUIRED ON)
+
 set( ATLAS_BITS_GLOBAL 64 CACHE STRING "bits used to define a global index (atlas::gidx_t)" )
 set( ATLAS_BITS_LOCAL  32 CACHE STRING "bits used to define a local index (atlas::idx_t)" )
 

atlas_io/src/atlas_io/CMakeLists.txt

@@ -103,4 +103,4 @@ ecbuild_add_library( TARGET atlas_io
         ${CMAKE_CURRENT_BINARY_DIR}/detail/defines.h
 )
 
-target_compile_features( atlas_io PUBLIC cxx_std_17 )
+#target_compile_features( atlas_io PUBLIC cxx_std_17 )

cmake/features/ACC.cmake

@@ -31,7 +31,15 @@ if( atlas_HAVE_ACC )
     set( ACC_Fortran_FLAGS ${OpenACC_Fortran_FLAGS} )
     set( ACC_C_FLAGS       ${OpenACC_C_FLAGS} )
   endif()
-endif()
+
+  if (ATLAS_FIND_OPENACC)
+    enable_language(C)  # So that find_package(OpenACC) looks C components
+    find_package(OpenACC)
+  endif()
+
+  ecbuild_add_option( FEATURE ACC
+                      DESCRIPTION  "OpenACC capable data structures"
+                      CONDITION HAVE_CUDA AND OpenACC_C_FOUND )
 
 else()
   set( HAVE_ACC 0 )

cmake/project_summary.cmake

@@ -42,9 +42,9 @@ endif()
 if( atlas_HAVE_ACC )
 
     ecbuild_info( "ACC" )
-    ecbuild_info( "    ACC_C_COMPILER     : [${ACC_C_COMPILER}]" )
-    ecbuild_info( "    ACC_C_FLAGS        : [${ACC_C_FLAGS}]" )
-    ecbuild_info( "    ACC_Fortran_FLAGS  : [${ACC_Fortran_FLAGS}]" )
+    ecbuild_info( "    OpenACC_C_FLAGS       : [${OpenACC_C_FLAGS}]" )
+    ecbuild_info( "    OpenACC_CXX_FLAGS     : [${OpenACC_CXX_FLAGS}]" )
+    ecbuild_info( "    OpenACC_Fortran_FLAGS : [${OpenACC_Fortran_FLAGS}]" )
 
 endif()
 

src/atlas/CMakeLists.txt

@@ -472,12 +472,24 @@ field/FieldSet.cc
 field/FieldSet.h
 field/MissingValue.cc
 field/MissingValue.h
+field/MultiField.cc
+field/MultiField.h
+field/MultiFieldCreator.cc
+field/MultiFieldCreator.h
+field/MultiFieldCreatorIFS.cc
+field/MultiFieldCreatorIFS.h
+field/MultiFieldCreatorArray.cc
+field/MultiFieldCreatorArray.h
 field/State.cc
 field/State.h
 field/detail/FieldImpl.cc
 field/detail/FieldImpl.h
 field/detail/FieldInterface.cc
 field/detail/FieldInterface.h
+field/detail/MultiFieldImpl.cc
+field/detail/MultiFieldImpl.h
+field/detail/MultiFieldInterface.cc
+field/detail/MultiFieldInterface.h
 field/detail/MissingValue.cc
 field/detail/MissingValue.h
 )
@@ -811,6 +823,8 @@ util/Topology.h
 util/Allocate.h
 util/Allocate.cc
 util/GPUClonable.h
+util/RegisterPointerInfo.cc
+util/RegisterPointerInfo.h
 
 util/Constants.h
 util/ConvexSphericalPolygon.cc

src/atlas/array/Array.h

@@ -13,14 +13,21 @@
 #include <memory>
 #include <vector>
 
+#include "eckit/config/Parametrisation.h"
+
 #include "atlas/library/config.h"
+#include "atlas/util/Config.h"
 
 #include "atlas/util/Object.h"
 
 #include "atlas/array/ArrayDataStore.h"
 #include "atlas/array/DataType.h"
 #include "atlas/array_fwd.h"
 
+namespace eckit {
+class Parametrisation;
+}
+
 namespace atlas {
 namespace array {
 
@@ -85,12 +92,30 @@ class Array : public util::Object {
 
     const ArrayStrides& strides() const { return spec_.strides(); }
 
+    const ArrayStrides& device_strides() const { 
+        if( mapped_ ) {
+            return spec_.strides();
+        }
+        else {
+            return spec_.device_strides();
+        }
+    }
+
     const ArrayShape& shape() const { return spec_.shape(); }
 
     const std::vector<int>& shapef() const { return spec_.shapef(); }
 
     const std::vector<int>& stridesf() const { return spec_.stridesf(); }
 
+    const std::vector<int>& device_stridesf() const {
+        if( mapped_ ) {
+            return spec_.stridesf();
+        }
+        else {
+            return spec_.device_stridesf();
+        }
+    }
+
     bool contiguous() const { return spec_.contiguous(); }
 
     bool hasDefaultLayout() const { return spec_.hasDefaultLayout(); }
@@ -190,6 +215,7 @@ class Array : public util::Object {
     Array(ArraySpec&& spec): spec_(std::move(spec)) {}
     ArraySpec spec_;
     std::unique_ptr<ArrayDataStore> data_store_;
+    bool mapped_{false};
 
     void replace(Array& array) {
         data_store_.swap(array.data_store_);

src/atlas/array/ArrayDataStore.cc

@@ -16,7 +16,7 @@
 namespace atlas {
 namespace array {
 
-void throw_OutOfRange(const std::string& class_name, char idx_str, int idx, int max) {
+void throw_OutOfRange(std::string_view class_name, char idx_str, int idx, int max) {
     std::ostringstream msg;
     msg << class_name << " index " << idx << " out of bounds: " << idx << " >= " << max;
     throw_Exception(msg.str(), Here());

src/atlas/array/ArrayDataStore.h

@@ -10,7 +10,9 @@
 
 #pragma once
 
-#include <string>
+#include <string_view>
+
+#include <eckit/config/Parametrisation.h>
 
 #include "atlas/array/ArrayIdx.h"
 #include "atlas/array/ArrayLayout.h"
@@ -81,7 +83,7 @@ static constexpr char array_dim() {
     return Dim == 0 ? 'i' : (Dim == 1 ? 'j' : (Dim == 2 ? 'k' : (Dim == 3 ? 'l' : (Dim == 4 ? 'm' : ('*')))));
 }
 
-void throw_OutOfRange(const std::string& class_name, char idx_str, int idx, int max);
+void throw_OutOfRange(std::string_view class_name, char idx_str, int idx, int max);
 #endif
 
 //------------------------------------------------------------------------------------------------------

src/atlas/array/ArraySpec.cc

@@ -44,15 +44,19 @@ ArraySpec::ArraySpec(const ArrayShape& shape, const ArrayAlignment& alignment):
     rank_           = static_cast<int>(shape.size());
     size_           = 1;
     allocated_size_ = 1;
+    allocated_device_size_ = 1;
     shape_.resize(rank_);
     strides_.resize(rank_);
+    device_strides_.resize(rank_);
     layout_.resize(rank_);
     for (int j = rank_ - 1; j >= 0; --j) {
         shape_[j]   = shape[j];
         strides_[j] = allocated_size_;
+        device_strides_[j] = allocated_device_size_;
         layout_[j]  = j;
         size_ *= size_t(shape_[j]);
         allocated_size_ *= size_t(aligned_shape[j]);
+        allocated_device_size_ *= size_t(shape[j]);
     }
     ATLAS_ASSERT(allocated_size_ == compute_aligned_size(size_t(shape_[0]) * size_t(strides_[0]), size_t(alignment)));
     contiguous_     = (size_ == allocated_size_);
@@ -80,10 +84,12 @@ ArraySpec::ArraySpec(const ArrayShape& shape, const ArrayStrides& strides, const
     size_ = 1;
     shape_.resize(rank_);
     strides_.resize(rank_);
+    device_strides_.resize(rank_);
     layout_.resize(rank_);
     for (int j = rank_ - 1; j >= 0; --j) {
         shape_[j]   = shape[j];
         strides_[j] = strides[j];
+        device_strides_[j] = size_;
         layout_[j]  = j;
         size_ *= size_t(shape_[j]);
     }
@@ -120,11 +126,13 @@ ArraySpec::ArraySpec(const ArrayShape& shape, const ArrayStrides& strides, const
     size_ = 1;
     shape_.resize(rank_);
     strides_.resize(rank_);
+    device_strides_.resize(rank_);
     layout_.resize(rank_);
     default_layout_ = true;
     for (int j = rank_ - 1; j >= 0; --j) {
         shape_[j]   = shape[j];
         strides_[j] = strides[j];
+        device_strides_[j] = size_;
         layout_[j]  = layout[j];
         size_ *= size_t(shape_[j]);
         if (layout_[j] != idx_t(j)) {
@@ -144,6 +152,7 @@ ArraySpec::ArraySpec(DataType datatype, const ArrayShape& shape, const ArrayStri
     ArraySpec(shape, strides, layout, alignment) {
     datatype_ = datatype;
 }
+
 const std::vector<int>& ArraySpec::shapef() const {
     return shapef_;
 }
@@ -152,12 +161,19 @@ const std::vector<int>& ArraySpec::stridesf() const {
     return stridesf_;
 }
 
+const std::vector<int>& ArraySpec::device_stridesf() const {
+    return device_stridesf_;
+}
+
 void ArraySpec::allocate_fortran_specs() {
     shapef_.resize(rank_);
     stridesf_.resize(rank_);
+    device_stridesf_.resize(rank_);
+    device_stridesf_[rank_ - 1] = stridesf_[rank_ - 1];
     for (idx_t j = 0; j < rank_; ++j) {
         shapef_[j] = shape_[rank_ - 1 - layout_[j]];
         stridesf_[j] = strides_[rank_ -1 - layout_[j]];
+        device_stridesf_[rank_ - j - 1] = device_stridesf_[rank_ - j] * shapef_[rank_ - j - 1];
     }
 }
 

src/atlas/array/ArraySpec.h

@@ -30,13 +30,16 @@ class ArraySpec {
     size_t size_;
     idx_t rank_;
     size_t allocated_size_;
+    size_t allocated_device_size_;
     DataType datatype_;
     ArrayShape shape_;
     ArrayStrides strides_;
+    ArrayStrides device_strides_;
     ArrayLayout layout_;
     ArrayAlignment alignment_;
     std::vector<int> shapef_;
     std::vector<int> stridesf_;
+    std::vector<int> device_stridesf_;
     bool contiguous_;
     bool default_layout_;
 
@@ -61,9 +64,11 @@ class ArraySpec {
     const ArrayShape& shape() const { return shape_; }
     const ArrayAlignment& alignment() const { return alignment_; }
     const ArrayStrides& strides() const { return strides_; }
+    const ArrayStrides& device_strides() const { return device_strides_; }
     const ArrayLayout& layout() const { return layout_; }
     const std::vector<int>& shapef() const;
     const std::vector<int>& stridesf() const;
+    const std::vector<int>& device_stridesf() const;
     bool contiguous() const { return contiguous_; }
     bool hasDefaultLayout() const { return default_layout_; }
 

src/atlas/array/ArrayStrides.h

@@ -28,6 +28,7 @@ class ArrayStrides : public std::vector<idx_t> {
     ArrayStrides() {}
     ArrayStrides(std::initializer_list<idx_t> list): Base(list) {}
     ArrayStrides(Base&& base): Base(std::forward<Base>(base)) {}
+    ArrayStrides(const std::vector<idx_t>& list): Base(list.begin(), list.end()) {}
 };
 
 namespace detail {

src/atlas/array/native/NativeArray.cc

@@ -10,12 +10,15 @@
 
 #include <iostream>
 
+//#include "eckit/config/Parametrisation.h"
+
 #include "atlas/array.h"
 #include "atlas/array/ArrayDataStore.h"
 #include "atlas/array/MakeView.h"
 #include "atlas/array/helpers/ArrayInitializer.h"
 #include "atlas/array/helpers/ArrayWriter.h"
 #include "atlas/array/native/NativeDataStore.h"
+#include "atlas/util/Config.h"
 #include "atlas/runtime/Exception.h"
 
 using namespace atlas::array::helpers;