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

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

@@ -35,7 +35,9 @@ MethodBuilder<Binning> __builder("binning");
 Binning::Binning(const Config& config) : Method(config) {
   const auto* ptr_config = dynamic_cast<const eckit::LocalConfiguration*>(&config);
   interpAncillaryScheme_ = ptr_config->getSubConfiguration("ancillary_scheme");
+  grid_type_ = ptr_config->getString("grid_type", "ATLAS");
   halo_exchange_ = ptr_config->getBool("halo_exchange", true);
+  adjoint_ = ptr_config->getBool("adjoint", false);
 }
 
 
@@ -61,8 +63,10 @@ void Binning::do_setup(const FunctionSpace& source,
     return;
   }
 
-  // enabling or disabling halo exchange
+  // enabling or disabling the halo exchange
   this->allow_halo_exchange_ = halo_exchange_;
+  // enabling or disabling the adjoint operation
+  this->adjoint_ = adjoint_;
 
   // note that the 'source' grid for the low-to-high regridding (interpolation)
   // is the 'target' grid for high-to-low regridding (binning) and
@@ -99,6 +103,9 @@ void Binning::do_setup(const FunctionSpace& source,
     // 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]));
@@ -130,7 +137,8 @@ std::vector<double> Binning::getAreaWeights(const FunctionSpace& fspace) const {
   // diagonal of 'area weights matrix', W
   std::vector<double> ds_aweights;
 
-  if (fspace.type().compare(functionspace::detail::NodeColumns::static_type()) == 0) {
+  if (grid_type_ == "CS-LFR") {
+
     const auto csfs = atlas::functionspace::NodeColumns(fspace);
 
     const auto csgrid = atlas::CubedSphereGrid(csfs.mesh().grid());
@@ -155,17 +163,6 @@ std::vector<double> Binning::getAreaWeights(const FunctionSpace& fspace) const {
         ds_aweights.emplace_back(aweight_temp);
       }
     }
-
-  } else if (fspace.type().compare(functionspace::detail::StructuredColumns::static_type()) == 0) {
-
-    const auto scfs = atlas::functionspace::StructuredColumns(fspace);
-
-    const auto scgrid = scfs.grid();
-    // area weight (appoximated)
-    //   area weight = area_cell/area_total
-    //               = area_cell/(no_cells*area_cell) = 1/no_cells
-    const double aweight = 1/static_cast<double>(scgrid.size());
-    ds_aweights.assign(fspace.size(), aweight);
 
   } else {
 

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

@@ -7,6 +7,9 @@
 
 #pragma once
 
+#include <string>
+#include <vector>
+
 #include "atlas/functionspace/FunctionSpace.h"
 #include "atlas/interpolation/Cache.h"
 #include "atlas/interpolation/method/Method.h"
@@ -34,6 +37,18 @@ class Binning : public Method {
   ///                   interpolation matrix: A
   ///            normalization factor matrix: N
   ///
+  ///          Setup, configuration variables:
+  ///                      <type>: method used to evaluate the 'B' matrix;
+  ///                       value: 'binning'
+  ///          <ancillary_scheme>: method used to evaluate the 'A' matrix;
+  ///                       value: 'cubedsphere-bilinear', 'structured-bilinear', ...
+  ///                 <grid_type>: type of grid
+  ///                       value: 'CS-LFR', 'O', ...
+  ///             <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 {}
 
@@ -54,7 +69,9 @@ class Binning : public Method {
   FunctionSpace source_{};
   FunctionSpace target_{};
 
+  std::string grid_type_{"ATLAS"};
   bool halo_exchange_{true};
+  bool adjoint_{false};
 };
 
 

src/tests/interpolation/CMakeLists.txt

@@ -132,6 +132,7 @@ ecbuild_add_test( TARGET atlas_test_interpolation_binning
   SOURCES  test_interpolation_binning.cc
   LIBS     atlas
   MPI      6
+  CONDITION eckit_HAVE_MPI AND MPI_SLOTS GREATER_EQUAL 6
   ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT}
 )
 

src/tests/interpolation/test_interpolation_binning.cc

@@ -13,12 +13,14 @@
 #include "atlas/field/FieldSet.h"
 #include "atlas/functionspace/CubedSphereColumns.h"
 #include "atlas/functionspace/NodeColumns.h"
+#include "atlas/functionspace/StructuredColumns.h"
 #include "atlas/grid.h"
 #include "atlas/interpolation.h"
 #include "atlas/mesh.h"
 #include "atlas/meshgenerator/MeshGenerator.h"
 #include "atlas/option/Options.h"
 #include "atlas/output/Gmsh.h"
+#include "atlas/runtime/Log.h"
 #include "atlas/util/Config.h"
 #include "atlas/util/CoordinateEnums.h"
 #include "atlas/util/function/VortexRollup.h"
@@ -54,9 +56,9 @@ Field getVortexField(
 
 /// function to write a field set in a Gmsh file
 void makeGmshOutput(const std::string& file_name,
+                    const Mesh& mesh,
                     const FieldSet& fields) {
   const auto& fs = fields[0].functionspace();
-  const auto& mesh = functionspace::NodeColumns(fs).mesh();
 
   const auto config_gmsh =
     util::Config("coordinates", "xyz") |
@@ -68,6 +70,27 @@ void makeGmshOutput(const std::string& file_name,
 }
 
 
+/// function to carry out a dot product
+double dotProd(const atlas::Field& field01, const atlas::Field& field02) {
+  double dprod{};
+
+  const auto field01_view = array::make_view<double, 2>(field01);
+  const auto field02_view = array::make_view<double, 2>(field02);
+
+  for (atlas::idx_t l=0; l < field01_view.shape(1); ++l) {
+    for (atlas::idx_t i=0; i < field01_view.shape(0); ++i) {
+      dprod += field01_view(i, l) * field02_view(i, l);
+    }
+  }
+  eckit::mpi::comm().allReduceInPlace(dprod, eckit::mpi::Operation::SUM);
+
+  return dprod;
+}
+
+
+//--
+
+
 /// test to carry out the rigridding from 'high' to 'low' resolution,
 /// for a given type of grid (CS-LFR)
 ///
@@ -101,7 +124,8 @@ CASE("rigridding from high to low resolution; grid type: CS-LFR") {
 
 
   const auto scheme = util::Config("type", "binning") |
-    util::Config("ancillary_scheme", option::type("cubedsphere-bilinear"));
+                      util::Config("ancillary_scheme", option::type("cubedsphere-bilinear")) |
+                      util::Config("grid_type", "CS-LFR");
 
   Interpolation regrid_high2low(scheme, csfs_s, csfs_t);
 
@@ -111,13 +135,147 @@ CASE("rigridding from high to low resolution; grid type: CS-LFR") {
 
   //--
 
-  const std::string fname_s("regridding_h2l_field_s.msh");
+  const std::string fname_s("regridding_h2l_cs_s.msh");
+
+  makeGmshOutput(fname_s, csmesh_s, fs_s["field_01_s"]);
+
+  const std::string fname_t("regridding_h2l_cs_t.msh");
+
+  makeGmshOutput(fname_t, csmesh_t, fs_t["field_01_t"]);
+}
+
+
+/// test to carry out the rigridding from 'high' to 'low' resolution,
+/// for a given type of grid (O)
+///
+/// after the regridding, the field set that has been generated is stored
+/// in a Gmsh file (data visualization);
+///
+CASE("rigridding from high to low resolution; grid type: O") {
+
+  // source grid (high res.)
+  auto ncgrid_s = Grid("O32");
+  auto ncmesh_s = MeshGenerator("structured").generate(ncgrid_s);
+  auto ncfs_s = functionspace::StructuredColumns(ncgrid_s, option::halo(3));
+
+  // target grid (low res.)
+  auto ncgrid_t = Grid("O16");
+  auto ncmesh_t = MeshGenerator("structured").generate(ncgrid_t);
+  auto ncfs_t = functionspace::StructuredColumns(ncgrid_t, option::halo(3));
+
+  size_t nb_levels = 1;
+
+  auto field_01_s = getVortexField(ncfs_s, "field_01_s", nb_levels);
+
+  atlas::FieldSet fs_s;
+  fs_s.add(field_01_s);
+
+  auto field_01_t = ncfs_t.createField<double>(option::name("field_01_t") |
+                                               option::levels(nb_levels));
+
+  FieldSet fs_t;
+  fs_t.add(field_01_t);
+
+
+  const auto scheme = util::Config("type", "binning") |
+                      util::Config("ancillary_scheme", option::type("structured-bilinear")) |
+                      util::Config("grid_type", "O");
+
+  Interpolation regrid_high2low(scheme, ncfs_s, ncfs_t);
+
+  // performing the regridding from high to low resolution
+  regrid_high2low.execute(fs_s, fs_t);
+
+
+  //--
+
+  const std::string fname_s("regridding_h2l_nc_s.msh");
+
+  makeGmshOutput(fname_s, ncmesh_s, fs_s["field_01_s"]);
+
+  const std::string fname_t("regridding_h2l_nc_t.msh");
+
+  makeGmshOutput(fname_t, ncmesh_t, fs_t["field_01_t"]);
+}
+
+
+/// test to carry out the 'dot-product' test for the rigridding from
+/// 'high' to 'low' resolution, for a given type of grid (CS-LFR)
+///
+CASE("dot-product test for the rigridding from high to low resolution; grid type: CS-LFR") {
+
+  // source grid (high res.)
+  auto csgrid_s = Grid("CS-LFR-100");
+  auto csmesh_s = MeshGenerator("cubedsphere_dual").generate(csgrid_s);
+  auto csfs_s = functionspace::NodeColumns(csmesh_s);
+
+  // target grid (low res.)
+  auto csgrid_t = Grid("CS-LFR-50");
+  auto csmesh_t = MeshGenerator("cubedsphere_dual").generate(csgrid_t);
+  auto csfs_t = functionspace::NodeColumns(csmesh_t);
+
+  size_t nb_levels = 1;
+
+  // source field
+  auto field_01_s = getVortexField(csfs_s, "field_01_s", nb_levels);
+
+  atlas::FieldSet fs_s;
+  fs_s.add(field_01_s);
+
+  // target field
+  auto field_01_t = csfs_t.createField<double>(option::name("field_01_t") |
+                                               option::levels(nb_levels));
+
+  atlas::FieldSet fs_t;
+  fs_t.add(field_01_t);
+
+  const auto scheme = util::Config("type", "binning") |
+                      util::Config("ancillary_scheme", option::type("cubedsphere-bilinear")) |
+                      util::Config("grid_type", "CS-LFR") |
+                      util::Config("adjoint", true);
+
+  Interpolation regrid_high2low(scheme, csfs_s, csfs_t);
+
+  // performing the regridding from high to low resolution
+  regrid_high2low.execute(fs_s, fs_t);
+
+
+  fs_t["field_01_t"].haloExchange();
+
+  // target field (adjoint)
+  auto field_01_ad_t = csfs_t.createField<double>(option::name("field_01_ad_t") |
+                                                  option::levels(nb_levels));
+  atlas::array::make_view<double, 2>(field_01_ad_t).assign(
+    atlas::array::make_view<double, 2>(field_01_t));
+  field_01_ad_t.adjointHaloExchange();
+
+  atlas::FieldSet fs_ad_t;
+  fs_ad_t.add(field_01_ad_t);
+
+  // source field (adjoint)
+  auto field_01_ad_s = csfs_s.createField<double>(option::name("field_01_ad_s") |
+                                                  option::levels(nb_levels));
+  atlas::array::make_view<double, 2>(field_01_ad_s).assign(0.);
+
+  atlas::FieldSet fs_ad_s;
+  fs_ad_s.add(field_01_ad_s);
+
+  // performing adjoint operation
+  regrid_high2low.execute_adjoint(fs_ad_s, fs_ad_t);
+
+
+  const auto t_dot_t = dotProd(fs_t["field_01_t"], fs_t["field_01_t"]);
+  const auto s_dot_ad_s = dotProd(fs_s["field_01_s"], fs_ad_s["field_01_ad_s"]);
 
-  makeGmshOutput(fname_s, fs_s["field_01_s"]);
+  double scaled_diff = std::abs(t_dot_t - s_dot_ad_s)/std::abs(t_dot_t);
 
-  const std::string fname_t("regridding_h2l_field_t.msh");
+  // carrrying out a dot-product test ...
+  Log::info() << "\n- dot-product test:\n"
+              << "(Ax) . (Ax) = " << t_dot_t << "; "
+              << "x . (A^t A x) = " << s_dot_ad_s << "; "
+              << "scaled difference = " << scaled_diff << "\n" << std::endl;
 
-  makeGmshOutput(fname_t, fs_t["field_01_t"]);
+  EXPECT(scaled_diff < 1e-12);
 }