Add SphericalVector to MethodFactory

ecmwf · Dec 7, 2023 · a96057e · a96057e
1 parent ee4c3ad
commit a96057e
Show file tree

Hide file tree

Showing 3 changed files with 23 additions and 6 deletions.
diff --git a/src/atlas/interpolation/method/MethodFactory.cc b/src/atlas/interpolation/method/MethodFactory.cc
@@ -16,6 +16,7 @@
 #include "knn/GridBoxMaximum.h"
 #include "knn/KNearestNeighbours.h"
 #include "knn/NearestNeighbour.h"
+#include "sphericalvector/SphericalVector.h"
 #include "structured/Cubic2D.h"
 #include "structured/Cubic3D.h"
 #include "structured/Linear2D.h"
@@ -46,6 +47,7 @@ void force_link() {
             MethodBuilder<method::GridBoxAverage>();
             MethodBuilder<method::GridBoxMaximum>();
             MethodBuilder<method::CubedSphereBilinear>();
+            MethodBuilder<method::SphericalVector>();
         }
     } link;
 }

diff --git a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc
@@ -6,7 +6,6 @@
  */
 
 #include "atlas/library/defines.h"
-#if ATLAS_HAVE_EIGEN
 
 #include <cmath>
 #include <tuple>
@@ -37,17 +36,21 @@ namespace method {
 
 using Complex = SphericalVector::Complex;
 
+#if ATLAS_HAVE_EIGEN
 template <typename Value>
 using SparseMatrix = SphericalVector::SparseMatrix<Value>;
 using RealMatrixMap = Eigen::Map<const SparseMatrix<double>>;
 using ComplexTriplets = std::vector<Eigen::Triplet<Complex>>;
 using RealTriplets = std::vector<Eigen::Triplet<double>>;
+#endif
+
 using EckitMatrix = eckit::linalg::SparseMatrix;
 
 namespace {
 
 MethodBuilder<SphericalVector> __builder("spherical-vector");
 
+#if ATLAS_HAVE_EIGEN
 RealMatrixMap makeMatrixMap(const EckitMatrix& baseMatrix) {
   return RealMatrixMap(baseMatrix.rows(), baseMatrix.cols(),
                        baseMatrix.nonZeros(), baseMatrix.outer(),
@@ -110,6 +113,7 @@ void matrixMultiply(const SourceView& sourceView, TargetView& targetView,
 
   sparseMatrixForEach(multiplyColumn, matrices...);
 }
+#endif
 
 }  // namespace
 
@@ -128,6 +132,8 @@ void SphericalVector::do_setup(const FunctionSpace& source,
     return;
   }
 
+#if ATLAS_HAVE_EIGEN
+
   setMatrix(Interpolation(interpolationScheme_, source_, target_));
 
   // Get matrix data.
@@ -169,6 +175,10 @@ void SphericalVector::do_setup(const FunctionSpace& source,
   realWeights_->setFromTriplets(realTriplets.begin(), realTriplets.end());
 
   ATLAS_ASSERT(complexWeights_->nonZeros() == matrix().nonZeros());
+
+#else
+  ATLAS_THROW_EXCEPTION("atlas has been compiled without Eigen");
+#endif
 }
 
 void SphericalVector::print(std::ostream&) const { ATLAS_NOTIMPLEMENTED; }
@@ -241,6 +251,7 @@ void SphericalVector::interpolate_vector_field(const Field& sourceField,
   auto targetView = array::make_view<Value, Rank>(targetField);
   targetView.assign(0.);
 
+#if ATLAS_HAVE_EIGEN
   const auto horizontalComponent = [](const auto& sourceVars, auto& targetVars,
                                       const auto& complexWeight) {
     const auto sourceVector = Complex(sourceVars(0), sourceVars(1));
@@ -267,12 +278,13 @@ void SphericalVector::interpolate_vector_field(const Field& sourceField,
 
     return;
   }
+#else
+  ATLAS_THROW_EXCEPTION("atlas has been compiled without Eigen");
+#endif
 
   ATLAS_NOTIMPLEMENTED;
 }
 
 }  // namespace method
 }  // namespace interpolation
 }  // namespace atlas
-
-#endif
diff --git a/src/atlas/interpolation/method/sphericalvector/SphericalVector.h b/src/atlas/interpolation/method/sphericalvector/SphericalVector.h
@@ -6,14 +6,15 @@
  */
 
 #include "atlas/library/defines.h"
-#if ATLAS_HAVE_EIGEN
 
 #pragma once
 
 #include <complex>
 #include <memory>
 
+#if ATLAS_HAVE_EIGEN
 #include <Eigen/Sparse>
+#endif
 
 #include "atlas/functionspace/FunctionSpace.h"
 #include "atlas/interpolation/method/Method.h"
@@ -28,10 +29,12 @@ class SphericalVector : public Method {
  public:
   using Complex = std::complex<double>;
 
+#if ATLAS_HAVE_EIGEN
   template <typename Value>
   using SparseMatrix = Eigen::SparseMatrix<Value, Eigen::RowMajor>;
   using ComplexMatrix = SparseMatrix<Complex>;
   using RealMatrix = SparseMatrix<double>;
+#endif
 
   /// @brief   Interpolation post-processor for vector field interpolation
   ///
@@ -90,12 +93,12 @@ class SphericalVector : public Method {
   FunctionSpace source_;
   FunctionSpace target_;
 
+#if ATLAS_HAVE_EIGEN
   std::shared_ptr<ComplexMatrix> complexWeights_;
   std::shared_ptr<RealMatrix> realWeights_;
+#endif
 };
 
 }  // namespace method
 }  // namespace interpolation
 }  // namespace atlas
-
-#endif