Merge branch 'develop' into feature/cubed_sphere_interp_omp

Author: wdeconinck

.github/ci-hpc-config.yml

@@ -7,21 +7,25 @@ mpi_on:
     modules:
       - ninja
     modules_package:
-      - atlas:fftw,eigen,openmpi
+      - atlas:fftw,eigen,qhull,python3/3.12,openmpi
+      - fckit:python3/3.12,openmpi
       - eckit:openmpi
     dependencies:
       - ecmwf/ecbuild@develop
       - ecmwf/eckit@develop
       - ecmwf/fckit@develop
     parallel: 64
     ntasks: 16
+    env:
+      - ATLAS_FPE=0
 
 mpi_off:
   build:
     modules:
       - ninja
     modules_package:
-      - atlas:fftw,eigen
+      - atlas:fftw,eigen,qhull,python3/3.12
+      - fckit:python3/3.12
     dependencies:
       - ecmwf/ecbuild@develop
       - ecmwf/eckit@develop

.github/workflows/build-hpc.yml

@@ -93,6 +93,7 @@ jobs:
                   - -DENABLE_ECKIT_SQL=OFF
                   - -DENABLE_ECKIT_CMD=OFF
                   - -DENABLE_AIO=OFF
+                  - -DENABLE_ECKIT_GEO=OFF
                   - -DENABLE_WARNINGS=OFF
                   - ${{ matrix.cmake_options || '' }}
               ecmwf/fckit:

pluto/src/pluto_f/pluto_module.F90

@@ -7,6 +7,10 @@ module pluto_module
 public :: pluto, pluto_memory_resource, pluto_allocator
 
 interface
+    function c_pluto_devices() result(devices) bind(c)
+        use iso_c_binding, only: c_int
+        integer(c_int) :: devices
+    end function
     subroutine c_pluto_host_set_default_resource_name(name, name_size) bind(c)
         use iso_c_binding, only: c_ptr, c_int
         type(c_ptr), value, intent(in) :: name
@@ -249,6 +253,7 @@ function c_pluto_trace_enabled() result(enabled) bind(c)
     type(pluto_scope_t)  :: scope
     type(pluto_trace_t)  :: trace
 contains
+    procedure, nopass :: devices => pluto_devices
     procedure, nopass :: get_registered_resource => pluto_get_registered_resource
     procedure, nopass :: new_delete_resource => pluto_new_delete_resource
     procedure, nopass :: null_memory_resource => pluto_null_memory_resource
@@ -328,6 +333,11 @@ function pluto_trace_enabled()
     endif
 end function
 
+function pluto_devices()
+    integer(c_int) :: pluto_devices
+    pluto_devices = c_pluto_devices()
+end function
+
 function pluto_get_registered_resource(name) result(memory_resource)
     type(pluto_memory_resource) :: memory_resource
     character(len=*), target, intent(in) :: name

pluto/src/pluto_f/pluto_module.cc

@@ -109,4 +109,8 @@ pluto::memory_pool_resource* c_pluto_device_pool_resource() {
 pluto::memory_pool_resource* c_pluto_managed_pool_resource() {
     return pluto::managed_pool_resource();
 }
+int c_pluto_devices() {
+    return pluto::devices();
+}
+
 }

src/atlas/array/native/NativeMakeView.cc

@@ -10,6 +10,8 @@
 
 #include <sstream>
 
+#include "pluto/pluto.h"
+
 #include "atlas/array.h"
 #include "atlas/array/ArrayView.h"
 #include "atlas/array/IndexView.h"
@@ -50,7 +52,9 @@ ArrayView<const Value, Rank> make_host_view(const Array& array) {
 template <typename Value, int Rank>
 ArrayView<Value, Rank> make_device_view(Array& array) {
 #if ATLAS_HAVE_GPU
-    ATLAS_ASSERT(array.deviceAllocated(),"make_device_view: Array not allocated on device");
+    if (pluto::devices()) {
+        ATLAS_ASSERT(array.deviceAllocated() || pluto::devices(), "make_device_view: Array not allocated on device");
+    }
     return ArrayView<Value, Rank>((array.device_data<Value>()), array.shape(), array.device_strides());
 #else
     return make_host_view<Value, Rank>(array);
@@ -60,7 +64,9 @@ ArrayView<Value, Rank> make_device_view(Array& array) {
 template <typename Value, int Rank>
 ArrayView<const Value, Rank> make_device_view(const Array& array) {
 #if ATLAS_HAVE_GPU
-    ATLAS_ASSERT(array.deviceAllocated(),"make_device_view: Array not allocated on device");
+    if (pluto::devices()) {
+        ATLAS_ASSERT(array.deviceAllocated(), "make_device_view: Array not allocated on device");
+    }
     return ArrayView<const Value, Rank>(array.device_data<const Value>(), array.shape(), array.device_strides());
 #else
     return make_host_view<const Value, Rank>(array);

src/atlas/field/detail/FieldInterface.cc

@@ -46,7 +46,12 @@ void atlas__Field__device_data_specf(FieldImpl* This, Value*& data, int& rank, i
     }
     data     = This->array().device_data<Value>();
     shapef   = const_cast<int*>(This->shapef().data());
-    stridesf = const_cast<int*>(This->device_stridesf().data());
+    if (data == This->array().host_data<Value>()) {
+        stridesf = const_cast<int*>(This->stridesf().data());
+    }
+    else {
+        stridesf = const_cast<int*>(This->device_stridesf().data());
+    }
     rank     = This->shapef().size();
 }
 

src/atlas/functionspace/NodeColumns.cc

@@ -292,6 +292,10 @@ void NodeColumns::set_field_metadata(const eckit::Configuration& config, Field&
     if (config.has("type")) {
       field.metadata().set("type", config.getString("type"));
     }
+
+    if (config.has("vector_component")) {
+        field.metadata().set("vector_component", config.getSubConfiguration("vector_component"));
+    }
 }
 
 array::DataType NodeColumns::config_datatype(const eckit::Configuration& config) const {

src/atlas/functionspace/PointCloud.cc

@@ -81,6 +81,7 @@ static std::string get_mpi_comm(const eckit::Configuration& config) {
 
 template <>
 PointCloud::PointCloud(const std::vector<PointXY>& points, const eckit::Configuration& config) {
+    ATLAS_TRACE("PointCloud(std::vector<PointXY>, eckit::Configuration)");
     mpi_comm_ = get_mpi_comm(config);
     lonlat_     = Field("lonlat", array::make_datatype<double>(), array::make_shape(points.size(), 2));
     auto lonlat = array::make_view<double, 2>(lonlat_);
@@ -92,6 +93,7 @@ PointCloud::PointCloud(const std::vector<PointXY>& points, const eckit::Configur
 
 template <>
 PointCloud::PointCloud(const std::vector<PointXYZ>& points, const eckit::Configuration& config) {
+    ATLAS_TRACE("PointCloud(std::vector<PointXYZ>, eckit::Configuration)");
     mpi_comm_ = get_mpi_comm(config);
     lonlat_       = Field("lonlat", array::make_datatype<double>(), array::make_shape(points.size(), 2));
     vertical_     = Field("vertical", array::make_datatype<double>(), array::make_shape(points.size()));
@@ -105,16 +107,18 @@ PointCloud::PointCloud(const std::vector<PointXYZ>& points, const eckit::Configu
 }
 
 PointCloud::PointCloud(const Field& lonlat, const eckit::Configuration& config): lonlat_(lonlat) {
+        ATLAS_TRACE("PointCloud(Field lonlat, eckit::Configuration)");
         mpi_comm_ = get_mpi_comm(config);
 }
 
 PointCloud::PointCloud(const Field& lonlat, const Field& ghost, const eckit::Configuration& config): lonlat_(lonlat), ghost_(ghost) {
+    ATLAS_TRACE("PointCloud(Field lonlat, Field ghost, eckit::Configuration)");
     mpi_comm_ = get_mpi_comm(config);
-    setupHaloExchange();
-    setupGatherScatter();
+    setupParallel();
 }
 
 PointCloud::PointCloud(const FieldSet& flds, const eckit::Configuration& config): lonlat_(flds["lonlat"]) {
+    ATLAS_TRACE("PointCloud(Fieldset, eckit::Configuration)");
     mpi_comm_ = get_mpi_comm(config);
     if (flds.has("ghost")) {
         ghost_ = flds["ghost"];
@@ -129,8 +133,7 @@ PointCloud::PointCloud(const FieldSet& flds, const eckit::Configuration& config)
         global_index_ = flds["global_index"];
     }
     if( ghost_ && remote_index_ && partition_ ) {
-        setupHaloExchange();
-        setupGatherScatter();
+        setupParallel();
     }
 }
 
@@ -165,6 +168,8 @@ PointCloud::PointCloud(const Grid& grid, const grid::Distribution& distribution,
     auto size_owned = distribution.nb_pts()[part_];
     size_owned_ = size_owned;
 
+    size_global_ = grid.size();
+
     if (halo_radius == 0. || nb_partitions_ == 1) {
         idx_t size_halo = size_owned;
         ATLAS_ASSERT(size_owned > 0);
@@ -180,10 +185,13 @@ PointCloud::PointCloud(const Grid& grid, const grid::Distribution& distribution,
         array::make_view<int,1>(ghost_).assign(0);
         array::make_view<int,1>(partition_).assign(part_);
 
+        ATLAS_ASSERT(grid.size() == distribution.size());
+
         idx_t j{0};
         gidx_t g{0};
         for (auto p : grid.lonlat()) {
             if( distribution.partition(g) == part_ ) {
+                ATLAS_ASSERT(j < size_halo);
                 gidx(j) = g+1;
                 ridx(j) = j;
                 lonlat(j, 0) = p.lon();
@@ -262,13 +270,16 @@ PointCloud::PointCloud(const Grid& grid, const grid::Distribution& distribution,
         }
 
     }
-
-    setupHaloExchange();
-    setupGatherScatter();
+    setupParallel();
 }
 
 PointCloud::PointCloud(const Grid& grid, const grid::Partitioner& _partitioner, const eckit::Configuration& config):
-PointCloud(grid, ((_partitioner) ? _partitioner : grid::Partitioner("equal_regions", util::Config("mpi_comm",get_mpi_comm(config)))).partition(grid), config) {
+PointCloud(
+    grid,
+    grid::Distribution{grid, (_partitioner) ? _partitioner :
+        grid::Partitioner{grid.partitioner() | util::Config("mpi_comm",get_mpi_comm(config))}
+    },
+    config) {
     ATLAS_TRACE("PointCloud(grid,partitioner,config)");
 }
 
@@ -286,18 +297,18 @@ const Grid& PointCloud::grid() const {
     }
 
     std::vector<PointXY> points;
-    points.reserve(size_global_);
+    points.reserve(size_global());
     if (nb_partitions_ == 1) {
         for (const auto& point : iterate().xy()) {
             points.push_back(point);
         }
     }
     else {
         std::vector<int> gidx;
-        gidx.reserve(size_global_);
+        gidx.reserve(size_global());
         std::vector<double> x, y;
-        x.reserve(size_global_);
-        y.reserve(size_global_);
+        x.reserve(size_global());
+        y.reserve(size_global());
         const auto gidxView = array::make_view<gidx_t, 1>(global_index_);
         const auto ghostView = array::make_view<int, 1>(ghost_);
         int i = 0;
@@ -333,7 +344,7 @@ array::ArrayShape PointCloud::config_shape(const eckit::Configuration& config) c
             idx_t owner(0);
             config.get("owner", owner);
             idx_t rank = mpi::comm(mpi_comm()).rank();
-            _size = (rank == owner ? size_global_ : 0);
+            _size = (rank == owner ? size_global() : 0);
         }
     }
 
@@ -381,6 +392,11 @@ std::string PointCloud::config_name(const eckit::Configuration& config) const {
 }
 
 const parallel::HaloExchange& PointCloud::halo_exchange() const {
+    if (halo_exchange_) {
+        return *halo_exchange_;
+    }
+    const_cast<PointCloud&>(*this).setupHaloExchange();
+    ATLAS_ASSERT(halo_exchange_);
     return *halo_exchange_;
 }
 
@@ -428,10 +444,18 @@ void PointCloud::gather(const Field& local, Field& global) const {
     gather(local_fields, global_fields);
 }
 const parallel::GatherScatter& PointCloud::gather() const {
+    if (gather_scatter_) {
+        return *gather_scatter_;
+    }
+    const_cast<PointCloud&>(*this).setupGatherScatter();
     ATLAS_ASSERT(gather_scatter_);
     return *gather_scatter_;
 }
 const parallel::GatherScatter& PointCloud::scatter() const {
+    if (gather_scatter_) {
+        return *gather_scatter_;
+    }
+    const_cast<PointCloud&>(*this).setupGatherScatter();
     ATLAS_ASSERT(gather_scatter_);
     return *gather_scatter_;
 }
@@ -512,6 +536,10 @@ void PointCloud::set_field_metadata(const eckit::Configuration& config, Field& f
     if (config.has("type")) {
         field.metadata().set("type", config.getString("type"));
     }
+
+    if (config.has("vector_component")) {
+        field.metadata().set("vector_component", config.getSubConfiguration("vector_component"));
+    }
 }
 
 Field PointCloud::createField(const eckit::Configuration& options) const {
@@ -618,7 +646,7 @@ void dispatch_adjointHaloExchange(Field& field, const parallel::HaloExchange& ha
 }  // namespace
 
 void PointCloud::haloExchange(const FieldSet& fieldset, bool on_device) const {
-    if (halo_exchange_) {
+    if (parallel_) {
         for (idx_t f = 0; f < fieldset.size(); ++f) {
             Field& field = const_cast<FieldSet&>(fieldset)[f];
             switch (field.rank()) {
@@ -891,7 +919,7 @@ void PointCloud::create_remote_index() const {
     }
 }
 
-void PointCloud::setupHaloExchange() {
+void PointCloud::setupParallel() {
     ATLAS_TRACE();
     if (ghost_ and partition_ and global_index_ and not remote_index_) {
         create_remote_index();
@@ -1025,16 +1053,23 @@ void PointCloud::setupHaloExchange() {
     ATLAS_ASSERT(remote_index_);
     ATLAS_ASSERT(ghost_.size() == remote_index_.size());
     ATLAS_ASSERT(ghost_.size() == partition_.size());
- 
-    halo_exchange_.reset(new parallel::HaloExchange());
-    halo_exchange_->setup(mpi_comm_,
-                          array::make_view<int, 1>(partition_).data(),
-                          array::make_view<idx_t, 1>(remote_index_).data(),
-                          REMOTE_IDX_BASE,
-                          ghost_.size());
+    parallel_ = true;
+}
+
+void PointCloud::setupHaloExchange() {
+    ATLAS_TRACE();
+    if (ghost_ and partition_ and remote_index_) {
+        halo_exchange_.reset(new parallel::HaloExchange());
+        halo_exchange_->setup(mpi_comm_,
+                            array::make_view<int, 1>(partition_).data(),
+                            array::make_view<idx_t, 1>(remote_index_).data(),
+                            REMOTE_IDX_BASE,
+                            ghost_.size());
+    }
 }
 
 void PointCloud::setupGatherScatter() {
+    ATLAS_TRACE();
     if (ghost_ and partition_ and remote_index_ and global_index_) {
         gather_scatter_.reset(new parallel::GatherScatter());
         gather_scatter_->setup(mpi_comm_,
@@ -1044,13 +1079,29 @@ void PointCloud::setupGatherScatter() {
                             array::make_view<gidx_t, 1>(global_index_).data(),
                             array::make_view<int, 1>(ghost_).data(),
                             ghost_.size());
-        size_global_ = gather_scatter_->glb_dof();
     }
 }
 
+idx_t PointCloud::size_global() const {
+    if (size_global_ == -1) {
+        if (not parallel_) {
+            size_global_ = lonlat().size();
+        }
+        else {
+            if( !gather_scatter_ ) {
+                const_cast<PointCloud&>(*this).setupGatherScatter();
+            }
+            if (gather_scatter_) {
+                size_global_ = gather_scatter_->glb_dof();
+            }
+        }
+    }
+    return size_global_;
+}
+
 
 void PointCloud::adjointHaloExchange(const FieldSet& fieldset, bool on_device) const {
-    if (halo_exchange_) {
+    if (parallel_) {
         for (idx_t f = 0; f < fieldset.size(); ++f) {
             Field& field = const_cast<FieldSet&>(fieldset)[f];
             switch (field.rank()) {