implement a procedure to carry out a regridding from high to low reso…

…lution (binning) (#191)

* implemented a procedure to evaluate the areas of the grid cell (grid type: CS-LFR)

* implemented a procedure to carry out a regridding from high to low resolution (binning)

* added test for the regridding procedure (high to low resolution)

* updated build procedure

* update 03.05.24-01

* updated build procedure

* added check to identify and ignore an empty row

* updated setup procedure for the method 'binning'

* update 10.05.24-01

* added flag to control the adjoint operation

* added tests to battery of tests

* update 10.05.24-02

* added halo exchange operation before generating Gmsh data

* removed variable grid_type_

* removed unnecessary namespaces

* update 07.06.24-01

* update 07.06.24-02

* update 07.06.24-03

* update 07.06.24-04

* Consistency with SphericalVector: rename 'ancillary_scheme' to 'scheme'

* Fix intent of halo_exchange and adjoint in Binning method

* Remove SphericalVector::adjoint_ and use base class

* No casting to CubedSphereProjection necessary

* Reduce dimensionality of area field

* redefined the procedure to evaluate the areas of the grid cells

* updated build procedure

* update 14.06.24-01

---------

Co-authored-by: Willem Deconinck <willem.deconinck@ecmwf.int>

src/atlas/CMakeLists.txt

@@ -376,8 +376,6 @@ mesh/detail/AccumulateFacets.cc
 util/Unique.h
 util/Unique.cc
 
-mesh/actions/ExtendNodesGlobal.h
-mesh/actions/ExtendNodesGlobal.cc
 mesh/actions/BuildDualMesh.h
 mesh/actions/BuildCellCentres.cc
 mesh/actions/BuildCellCentres.h
@@ -401,6 +399,10 @@ mesh/actions/BuildStatistics.cc
 mesh/actions/BuildStatistics.h
 mesh/actions/BuildXYZField.cc
 mesh/actions/BuildXYZField.h
+mesh/actions/ExtendNodesGlobal.h
+mesh/actions/ExtendNodesGlobal.cc
+mesh/actions/GetCubedSphereNodalArea.cc
+mesh/actions/GetCubedSphereNodalArea.h
 mesh/actions/WriteLoadBalanceReport.cc
 mesh/actions/BuildTorusXYZField.h
 mesh/actions/BuildTorusXYZField.cc
@@ -618,6 +620,8 @@ interpolation/method/PointIndex3.cc
 interpolation/method/PointIndex3.h
 interpolation/method/PointIndex2.cc
 interpolation/method/PointIndex2.h
+interpolation/method/binning/Binning.cc
+interpolation/method/binning/Binning.h
 interpolation/method/cubedsphere/CubedSphereBilinear.cc
 interpolation/method/cubedsphere/CubedSphereBilinear.h
 interpolation/method/knn/GridBox.cc

src/atlas/grid/detail/grid/CubedSphere.h

@@ -241,7 +241,7 @@ class CubedSphere : public Grid {
     virtual std::string name() const override;
     virtual std::string type() const override;
 
-    // Return number of faces on cube
+    // Return N_, where (N_ * N_) is the number of cells on a tile
     inline idx_t N() const { return N_; }
 
     // Access to the tile class
@@ -436,7 +436,7 @@ class CubedSphere : public Grid {
     void xyt2xy(const double xyt[], double xy[]) const;
 
 protected:
-    // Number of faces on tile
+    // (N_ * N_) = number of cells on a tile
     idx_t N_;
 
     // Number of tiles

src/atlas/interpolation/method/Method.h

@@ -160,10 +160,10 @@ class Method : public util::Object {
     interpolation::MatrixCache matrix_cache_;
     NonLinear nonLinear_;
     std::string linalg_backend_;
-    bool adjoint_{false};
     Matrix matrix_transpose_;
 
 protected:
+    bool adjoint_{false};
     bool allow_halo_exchange_{true};
     std::vector<idx_t> missing_;
 };

src/atlas/interpolation/method/MethodFactory.cc

@@ -11,6 +11,7 @@
 #include "MethodFactory.h"
 
 // for static linking
+#include "binning/Binning.h"
 #include "cubedsphere/CubedSphereBilinear.h"
 #include "knn/GridBoxAverage.h"
 #include "knn/GridBoxMaximum.h"
@@ -49,6 +50,7 @@ void force_link() {
             MethodBuilder<method::GridBoxMaximum>();
             MethodBuilder<method::CubedSphereBilinear>();
             MethodBuilder<method::SphericalVector>();
+            MethodBuilder<method::Binning>();
         }
     } link;
 }

src/atlas/interpolation/method/binning/Binning.cc

@@ -0,0 +1,186 @@
+/*
+ * (C) Crown Copyright 2024 Met Office
+ *
+ * This software is licensed under the terms of the Apache Licence Version 2.0
+ * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+ */
+
+
+#include "atlas/functionspace/NodeColumns.h"
+#include "atlas/functionspace/StructuredColumns.h"
+#include "atlas/grid.h"
+#include "atlas/interpolation/Interpolation.h"
+#include "atlas/interpolation/method/binning/Binning.h"
+#include "atlas/interpolation/method/MethodFactory.h"
+#include "atlas/mesh.h"
+#include "atlas/mesh/actions/GetCubedSphereNodalArea.h"
+#include "atlas/runtime/Trace.h"
+
+#include "eckit/config/LocalConfiguration.h"
+#include "eckit/linalg/SparseMatrix.h"
+#include "eckit/linalg/Triplet.h"
+#include "eckit/mpi/Comm.h"
+
+
+namespace atlas {
+namespace interpolation {
+namespace method {
+
+
+namespace {
+
+MethodBuilder<Binning> __builder("binning");
+
+}
+
+
+Binning::Binning(const Config& config) : Method(config) {
+  const auto* conf = dynamic_cast<const eckit::LocalConfiguration*>(&config);
+  ATLAS_ASSERT(conf, "config must be derived from eckit::LocalConfiguration");
+  interpAncillaryScheme_ = conf->getSubConfiguration("scheme");
+  // enabling or disabling the adjoint operation
+  adjoint_ = conf->getBool("adjoint", false);
+  // enabling or disabling the halo exchange
+  allow_halo_exchange_ = conf->getBool("halo_exchange", true);
+}
+
+
+void Binning::do_setup(const Grid& source,
+                       const Grid& target,
+                       const Cache&) {
+  ATLAS_NOTIMPLEMENTED;
+}
+
+
+void Binning::do_setup(const FunctionSpace& source,
+                       const FunctionSpace& target) {
+  ATLAS_TRACE("atlas::interpolation::method::Binning::do_setup()");
+
+  using Index = eckit::linalg::Index;
+  using Triplet = eckit::linalg::Triplet;
+  using SMatrix = eckit::linalg::SparseMatrix;
+
+  source_ = source;
+  target_ = target;
+
+  if (target_.size() == 0) {
+    return;
+  }
+
+  // note that the 'source' grid for the low-to-high regridding (interpolation)
+  // is the 'target' grid for high-to-low regridding (binning) and
+  // the 'target' grid for the low-to-high regridding (interpolation) is the
+  // 'source' grid for the for high-to-low regridding (binning)
+  const auto& fs_source_interp = target_;
+  const auto& fs_target_interp = source_;
+
+  const auto interp = Interpolation(
+    interpAncillaryScheme_, fs_source_interp, fs_target_interp);
+  auto smx_interp_cache = interpolation::MatrixCache(interp);
+
+  auto smx_interp = smx_interp_cache.matrix();
+
+  auto smx_interp_tr = smx_interp.transpose();
+
+  const auto rows_tamx = smx_interp_tr.rows();
+  const auto cols_tamx = smx_interp_tr.cols();
+
+  const double* ptr_tamx_data = smx_interp_tr.data();
+  const Index* ptr_tamx_idxs_col = smx_interp_tr.inner();
+  const Index* ptr_tamx_o = smx_interp_tr.outer();
+
+  // diagonal of 'area weights matrix', W
+  auto ds_aweights = getAreaWeights(source_);
+
+  auto smx_binning_els = std::vector<Triplet>{};
+  size_t idx_row_next = 0;
+
+  for (size_t idx_row = 0; idx_row < rows_tamx; ++idx_row) {
+    idx_row_next = (idx_row+1);
+    // start of the indexes associated with the row 'i'
+    size_t lbound = ptr_tamx_o[idx_row];
+    // start of the indexes associated with the row 'i+1'
+    size_t ubound = ptr_tamx_o[idx_row_next];
+
+    if (lbound == ubound)
+      continue;
+
+    double sum_row = 0;
+    for (size_t i = lbound; i < ubound; ++i) {
+      sum_row += (ptr_tamx_data[i] * ds_aweights.at(ptr_tamx_idxs_col[i]));
+    }
+
+    // normalization factor
+    double nfactor = 1/sum_row;
+
+    for (size_t i = lbound; i < ubound; ++i) {
+      // evaluating the non-zero elements of the binning matrix
+      smx_binning_els.emplace_back(
+        idx_row, ptr_tamx_idxs_col[i],
+        (nfactor * (ptr_tamx_data[i] * ds_aweights.at(ptr_tamx_idxs_col[i]))));
+    }
+  }
+
+  // 'binning matrix' (sparse matrix), B = N A^T W
+  SMatrix smx_binning{rows_tamx, cols_tamx, smx_binning_els};
+  setMatrix(smx_binning);
+}
+
+
+void Binning::print(std::ostream&) const {
+  ATLAS_NOTIMPLEMENTED;
+}
+
+
+std::vector<double> Binning::getAreaWeights(const FunctionSpace& fspace) const {
+  // diagonal of 'area weights matrix', W
+  std::vector<double> ds_aweights;
+
+  bool is_cubed_sphere {false};
+  if (auto csfs = functionspace::NodeColumns(fspace)) {
+    if (CubedSphereGrid(csfs.mesh().grid())) {
+      is_cubed_sphere = true;
+    }
+  }
+
+  if (is_cubed_sphere) { 
+
+    const auto csfs = functionspace::NodeColumns(fspace); 
+    auto csmesh = csfs.mesh();
+
+    // areas of the cells (geographic coord. system)
+    auto gcell_areas = mesh::actions::GetCubedSphereNodalArea()(csmesh);
+    auto gcell_areas_view = array::make_view<double, 1>(gcell_areas);
+
+    auto is_ghost = array::make_view<int, 1>(csfs.ghost());
+
+    double total_area {0.};
+    for (idx_t i = 0; i < csfs.size(); i++) {
+      if (!is_ghost[i]) {
+        total_area += gcell_areas_view(i);
+      }
+    }
+    eckit::mpi::comm().allReduceInPlace(total_area, eckit::mpi::Operation::SUM);
+
+    double aweight_temp {0.};
+    for (idx_t i = 0; i < csfs.size(); i++) {
+      if (!is_ghost[i]) {
+        aweight_temp = gcell_areas_view(i)/total_area;
+        ds_aweights.emplace_back(aweight_temp);
+      }
+    }
+
+  } else {
+
+    // area weights (default)
+    ds_aweights.assign(fspace.size(), 1.);
+
+  }
+
+  return ds_aweights;
+}
+
+
+}  // namespace method
+}  // namespace interpolation
+}  // namespace atlas

src/atlas/interpolation/method/binning/Binning.h

@@ -0,0 +1,74 @@
+/*
+ * (C) Crown Copyright 2024 Met Office
+ *
+ * This software is licensed under the terms of the Apache Licence Version 2.0
+ * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+ */
+
+#pragma once
+
+#include <string>
+#include <vector>
+
+#include "atlas/functionspace/FunctionSpace.h"
+#include "atlas/interpolation/Cache.h"
+#include "atlas/interpolation/method/Method.h"
+#include "atlas/grid/Grid.h"
+
+#include "eckit/config/Configuration.h"
+
+
+namespace atlas {
+namespace interpolation {
+namespace method {
+
+
+class Binning : public Method {
+ public:
+  /// @brief   Binning procedure to carry out a regridding from high
+  ///          to low resolution
+  ///
+  /// @details This method is part of the family of the Interpolation operations;
+  ///          it relies on the evaluation of the transpose of an interpolation matrix.
+  ///          The rigridding from high to low resolution is carried out by using
+  ///          a 'binning matrix':
+  ///                         binning matrix: B = N A^T W
+  ///                    area weights matrix: W
+  ///                   interpolation matrix: A
+  ///            normalization factor matrix: N
+  ///
+  ///          Setup, configuration variables:
+  ///                      <type>: method used to evaluate the 'B' matrix;
+  ///                       value: 'binning'
+  ///                    <scheme>: method used to evaluate the 'A' matrix;
+  ///                       value: 'cubedsphere-bilinear', 'structured-bilinear', ...
+  ///             <halo_exchange>: flag to control the halo exchange procedure
+  ///                       value: 'true', 'false'
+  ///                   <adjoint>: flag to control the adjoint operation
+  ///                       value: 'true', 'false'
+  ///
+  Binning(const Config& config);
+  ~Binning() override {}
+
+  void print(std::ostream&) const override;
+  const FunctionSpace& source() const override { return source_; }
+  const FunctionSpace& target() const override { return target_; }
+
+ private:
+  using Method::do_setup;
+  void do_setup(const FunctionSpace& source,
+                const FunctionSpace& target) override;
+  void do_setup(const Grid& source, const Grid& target, const Cache&) override;
+
+  std::vector<double> getAreaWeights(const FunctionSpace& source) const;
+
+  eckit::LocalConfiguration interpAncillaryScheme_{};
+
+  FunctionSpace source_{};
+  FunctionSpace target_{};
+};
+
+
+}  // namespace method
+}  // namespace interpolation
+}  // namespace atlas

src/atlas/interpolation/method/sphericalvector/SphericalVector.h

@@ -72,7 +72,6 @@ class SphericalVector : public Method {
   void do_setup(const Grid& source, const Grid& target, const Cache&) override;
 
   eckit::LocalConfiguration interpolationScheme_{};
-  bool adjoint_{};
 
   FunctionSpace source_{};
   FunctionSpace target_{};