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

src/atlas/CMakeLists.txt

@@ -618,6 +618,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/CubedSphereGrid.h

@@ -9,6 +9,8 @@
 
 #include <string>
 
+#include "atlas/field.h"
+#include "atlas/functionspace/FunctionSpace.h"
 #include "atlas/grid/Grid.h"
 #include "atlas/grid/Tiles.h"
 #include "atlas/grid/detail/grid/CubedSphere.h"
@@ -297,6 +299,9 @@ class CubedSphereGrid : public Grid {
     // Return the size of the cubed sphere grid, where N is the number of grid boxes along the edge of a tile
     inline int N() const { return grid_->N(); }
 
+    // Return the areas of the cells
+    inline Field gridCellArea(const FunctionSpace& fspace) const { return grid_->gridCellArea(fspace); };
+
     /// @brief return tiles object.
     inline CubedSphereTiles tiles() const { return grid_->tiles(); }
 

src/atlas/grid/detail/grid/CubedSphere.cc

@@ -7,6 +7,8 @@
 
 #include "CubedSphere.h"
 
+#include <math.h>
+
 #include <algorithm>
 #include <iomanip>
 #include <limits>
@@ -16,7 +18,9 @@
 #include "eckit/utils/Hash.h"
 #include "eckit/utils/Translator.h"
 
+#include "atlas/array.h"
 #include "atlas/domain/Domain.h"
+#include "atlas/functionspace/NodeColumns.h"
 #include "atlas/grid/CubedSphereGrid.h"
 #include "atlas/grid/Tiles.h"
 #include "atlas/grid/detail/grid/GridBuilder.h"
@@ -303,6 +307,41 @@ void CubedSphere::xyt2xy(const double xyt[], double xy[]) const {
     xy[YY]             = normalisedY * 90. + tiles_offsets_ab2xy_[LAT][t];
 }
 
+// Provide the areas of the cells
+Field CubedSphere::gridCellArea(const FunctionSpace& fspace) const {
+
+  constexpr auto degrees2rads = M_PI / 180.;
+  const auto ncfs = functionspace::NodeColumns(fspace);
+
+  const auto csgrid = CubedSphereGrid(ncfs.mesh().grid());
+
+  auto lonlat = array::make_view<double, 2>(fspace.lonlat());
+
+  // (grid_res * grid_res) = no. of cells on a tile
+  auto grid_res = csgrid.N();
+
+  const auto& proj = csgrid.cubedSphereProjection();
+
+  // area of a grid cell (cubed-sphere coord. system)
+  double gcell_area_cs = M_PI/(2*grid_res) * M_PI/(2*grid_res);
+
+  auto gcell_area_field = ncfs.createField<double>(
+    atlas::option::name("grid_cell_areas") | atlas::option::levels(1));
+
+  auto gcell_area_fview = array::make_view<double, 2>(gcell_area_field);
+
+  for (size_t i = 0; i < gcell_area_fview.size(); i++) {
+    PointLonLat loc = PointLonLat(lonlat(i, atlas::LON), lonlat(i, atlas::LAT));
+    double cos_lat = std::cos(degrees2rads * loc.lat());
+    double grid_jac_det = 1/proj.jacobian(loc).determinant();
+    // area of a grid cell (geographic coord. system)
+    gcell_area_fview(i, 0) = grid_jac_det * gcell_area_cs * cos_lat;
+  }
+
+  return gcell_area_field;
+}
+
+
 // ------------------------------------------
 
 namespace {

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

@@ -16,6 +16,8 @@
 
 #include "eckit/types/Types.h"
 
+#include "atlas/field.h"
+#include "atlas/functionspace/FunctionSpace.h"
 #include "atlas/grid/Spacing.h"
 #include "atlas/grid/Tiles.h"
 #include "atlas/grid/detail/grid/Grid.h"
@@ -241,9 +243,12 @@ 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_; }
 
+    // Return the areas of the faces/cells
+    Field gridCellArea(const FunctionSpace&) const;
+
     // Access to the tile class
     inline atlas::grid::CubedSphereTiles tiles() const { return tiles_; }
 
@@ -436,7 +441,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,184 @@
+/*
+ * (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/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);
+
+    const auto csgrid = CubedSphereGrid(csfs.mesh().grid());
+    // areas of the cells (geographic coord. system)
+    const auto gcell_areas = csgrid.gridCellArea(fspace);
+
+    auto gcell_areas_view = array::make_view<double, 2>(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, 0);
+      }
+    }
+    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, 0)/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