diff --git a/src/mir/method/ProxyWeightedMethod.cc b/src/mir/method/ProxyWeightedMethod.cc index 680ee5e6d..8583910a7 100644 --- a/src/mir/method/ProxyWeightedMethod.cc +++ b/src/mir/method/ProxyWeightedMethod.cc @@ -13,7 +13,6 @@ #include "mir/method/ProxyWeightedMethod.h" #include -#include #include "eckit/log/JSON.h" #include "eckit/utils/MD5.h" @@ -33,7 +32,7 @@ namespace { struct StructuredBicubic final : public ProxyWeightedMethod { explicit StructuredBicubic(const param::MIRParametrisation& param) : - ProxyWeightedMethod(param, "structured-bicubic", "serial-halo to serial") {} + ProxyWeightedMethod(param, "structured-bicubic") {} }; @@ -45,12 +44,35 @@ static const MethodFactory* METHODS[]{ }; -ProxyWeightedMethod::ProxyWeightedMethod(const param::MIRParametrisation& param, const std::string& interpolation_type, - const std::string& renumber_type) : +ProxyWeightedMethod::ProxyWeightedMethod(const param::MIRParametrisation& param, + const std::string& interpolation_type) : MethodWeighted(param), type_(interpolation_type) { interpol_.set("matrix_free", false); interpol_.set("type", interpolation_type); - renumber_.set("type", renumber_type); +} + + +void ProxyWeightedMethod::foldSourceHalo(const atlas::Interpolation& interpol, size_t Nr, size_t Nc, + WeightMatrix& W) const { + ASSERT(Nr == W.rows()); + ASSERT(Nc < W.cols()); + + const auto& fs = interpol.source(); + const auto global_index{atlas::array::make_view(fs.global_index())}; + ASSERT(global_index.size() == W.cols()); + + auto* a = const_cast(W.data()); + for (auto c = Nc; c < W.cols(); ++c) { + ASSERT(1 <= global_index[c] && global_index[c] < Nc + 1); // (global indexing is 1-based) + a[global_index[c] - 1] += a[c]; + a[c] = 0.; + } + + W.prune(); + WeightMatrix M(new util::allocator::InPlaceAllocator{ + Nr, Nc, W.nonZeros(), const_cast(W.outer()), + const_cast(W.inner()), const_cast(W.data())}); + W.swap(M); } @@ -60,7 +82,7 @@ const char* ProxyWeightedMethod::name() const { int ProxyWeightedMethod::version() const { - return -1; + return 1; } @@ -79,12 +101,12 @@ bool ProxyWeightedMethod::sameAs(const Method& other) const { const auto* o = dynamic_cast(&other); return (o != nullptr) && name() == std::string{o->name()} && digest(interpol_) == digest(o->interpol_) && - digest(renumber_) == digest(o->renumber_) && MethodWeighted::sameAs(*o); + MethodWeighted::sameAs(*o); } void ProxyWeightedMethod::print(std::ostream& out) const { - out << "ProxyWeightedMethod[interpolation=" << interpol_ << ",renumber=" << renumber_ << ","; + out << "ProxyWeightedMethod[interpolation=" << interpol_ << ","; MethodWeighted::print(out); out << "]"; } @@ -92,7 +114,7 @@ void ProxyWeightedMethod::print(std::ostream& out) const { void ProxyWeightedMethod::json(eckit::JSON& j) const { j.startObject(); - j << "options" << interpol_.json(eckit::JSON::Formatting::compact()); + j << "interpolation" << interpol_.json(eckit::JSON::Formatting::compact()); MethodWeighted::json(j); j.endObject(); } @@ -100,43 +122,13 @@ void ProxyWeightedMethod::json(eckit::JSON& j) const { void ProxyWeightedMethod::assemble(util::MIRStatistics&, WeightMatrix& W, const repres::Representation& in, const repres::Representation& out) const { - // interpolation matrix build and assign (with a halo on the source grid) + // build matrix (with a halo on the source grid), move out of cache atlas::Interpolation interpol{interpol_, in.atlasGrid(), out.atlasGrid()}; - { - atlas::interpolation::MatrixCache cache{interpol}; - const auto& M = cache.matrix(); - W.swap(const_cast(M)); - } - - // fold serial+halo into serial - const auto Nr = out.numberOfPoints(); - const auto Nc = in.numberOfPoints(); - ASSERT(Nr == W.rows()); - ASSERT(Nc < W.cols()); - - { - const auto& fs = interpol.source(); - const auto gidx{atlas::array::make_view(fs.global_index())}; - - auto* a = const_cast(W.data()); - for (auto c = Nc; c < W.cols(); ++c) { - ASSERT(1 <= gidx[c] && gidx[c] < Nc + 1); // (global indexing is 1-based) - a[gidx[c] - 1] += a[c]; - a[c] = 0.; - } - W.prune(); - } - - // reshape matrix - WeightMatrix M(new util::allocator::InPlaceAllocator{ - Nr, Nc, W.nonZeros(), const_cast(W.outer()), - const_cast(W.inner()), const_cast(W.data())}); - W.swap(M); -} - + atlas::interpolation::MatrixCache cache{interpol}; + W.swap(const_cast(cache.matrix())); -bool ProxyWeightedMethod::validateMatrixWeights() const { - return true; + // fold source grid halo (from serial + halo into serial) + foldSourceHalo(interpol, out.numberOfPoints(), in.numberOfPoints(), W); } diff --git a/src/mir/method/ProxyWeightedMethod.h b/src/mir/method/ProxyWeightedMethod.h index 19b32cd72..410c5d272 100644 --- a/src/mir/method/ProxyWeightedMethod.h +++ b/src/mir/method/ProxyWeightedMethod.h @@ -24,19 +24,22 @@ class ProxyWeightedMethod : public MethodWeighted { protected: // -- Constructor - ProxyWeightedMethod(const param::MIRParametrisation&, const std::string& interpolation_type, - const std::string& renumber_type = ""); + ProxyWeightedMethod(const param::MIRParametrisation&, const std::string& interpolation_type); // -- Destructor ~ProxyWeightedMethod() override = default; +protected: + // -- Methods + + void foldSourceHalo(const atlas::Interpolation&, size_t Nr, size_t Nc, WeightMatrix&) const; + private: // -- Members const std::string type_; atlas::util::Config interpol_; - atlas::util::Config renumber_; // -- Overridden methods @@ -51,8 +54,6 @@ class ProxyWeightedMethod : public MethodWeighted { void assemble(util::MIRStatistics&, WeightMatrix&, const repres::Representation& in, const repres::Representation& out) const override; - - bool validateMatrixWeights() const override; };