diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index 0bddfc529..b37ee6934 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -98,6 +98,10 @@ jobs: - name: Checkout Repository uses: actions/checkout@v2 + - name: Set up Homebrew + id: set-up-homebrew + uses: Homebrew/actions/setup-homebrew@master + - name: Environment run: | echo "DEPS_DIR=${{ runner.temp }}/deps" >> $GITHUB_ENV diff --git a/CHANGELOG.md b/CHANGELOG.md index 6fc68a51c..abc131be9 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,25 @@ This project adheres to [Semantic Versioning](http://semver.org/spec/v2.0.0.html ## [Unreleased] +## [0.31.0] - 2022-11-10 +### Added +- Extend PointCloud functionspace to do halo exchanges, including Fortran API +- Add FunctionSpace::gather and FunctionSpace::scatter abstraction + +### Changed +- Improve performance of MatchingMeshPartitionerLonLatPolygon using OpenMP +- Improve performance of BuildHalo using OpenMP and unordered_map +- Improve performance of StructuredMeshGenerator using OpenMP +- Improve performance of ATLAS_TRACE +- Reduce memory peak in GatherScatter setup +- Global element numbering of RegularLonLat grids is now following rows independent of partitioning + +### Fixed +- Running CI with Github Actions +- Fix output of atlas-grids y-range for shifted grids +- Fix building with ATLAS_BITS_LOCAL=64 + + ## [0.30.0] - 2022-08-22 ### Added - Fortran API for Interpolation::execute_adjoint() @@ -386,6 +405,7 @@ This project adheres to [Semantic Versioning](http://semver.org/spec/v2.0.0.html ## 0.13.0 - 2018-02-16 [Unreleased]: https://github.com/ecmwf/atlas/compare/master...develop +[0.31.0]: https://github.com/ecmwf/atlas/compare/0.30.0...0.31.0 [0.30.0]: https://github.com/ecmwf/atlas/compare/0.29.0...0.30.0 [0.29.0]: https://github.com/ecmwf/atlas/compare/0.28.1...0.29.0 [0.28.1]: https://github.com/ecmwf/atlas/compare/0.28.0...0.28.1 diff --git a/VERSION b/VERSION index c25c8e5b7..26bea73e8 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -0.30.0 +0.31.0 diff --git a/atlas_io/src/atlas_io/Data.h b/atlas_io/src/atlas_io/Data.h index f57181d94..7431fd7f4 100644 --- a/atlas_io/src/atlas_io/Data.h +++ b/atlas_io/src/atlas_io/Data.h @@ -23,7 +23,7 @@ class Data { public: Data() = default; Data(void*, size_t); - Data(Data&&) = default; + Data(Data&&) = default; Data& operator=(Data&&) = default; operator const void*() const { return data(); } diff --git a/atlas_io/src/atlas_io/Metadata.h b/atlas_io/src/atlas_io/Metadata.h index dce380105..b6f2ab123 100644 --- a/atlas_io/src/atlas_io/Metadata.h +++ b/atlas_io/src/atlas_io/Metadata.h @@ -83,6 +83,14 @@ class Metadata : public eckit::LocalConfiguration { root.remove(name); return *this; } + + + std::vector keys() const { + // Preserves order of keys + std::vector result; + eckit::fromValue(result, get().keys()); + return result; + } }; //--------------------------------------------------------------------------------------------------------------------- diff --git a/atlas_io/src/atlas_io/Record.cc b/atlas_io/src/atlas_io/Record.cc index 9eff80b01..12ca10c7e 100644 --- a/atlas_io/src/atlas_io/Record.cc +++ b/atlas_io/src/atlas_io/Record.cc @@ -132,6 +132,20 @@ Record::operator const ParsedRecord&() const { //--------------------------------------------------------------------------------------------------------------------- +static void parse_record(ParsedRecord& record, const std::string& key, const Metadata& metadata) { + if (metadata.type() || metadata.link()) { + record.items.emplace(key, metadata); + record.keys.emplace_back(key); + } + else { + for (auto& next_key : metadata.keys()) { + parse_record(record, key + "." + next_key, metadata.getSubConfiguration(next_key)); + } + } +} + +//--------------------------------------------------------------------------------------------------------------------- + Record& Record::read(Stream& in, bool read_to_end) { if (not empty()) { return *this; @@ -206,9 +220,9 @@ Record& Record::read(Stream& in, bool read_to_end) { ATLAS_IO_ASSERT(r.metadata_format == "yaml"); Metadata metadata = eckit::YAMLConfiguration(metadata_str); - record_->keys = metadata.keys(); - for (const auto& key : record_->keys) { - record_->items.emplace(key, metadata.getSubConfiguration(key)); + + for (auto& key : metadata.keys()) { + parse_record(*record_, key, metadata.getSubConfiguration(key)); } diff --git a/atlas_io/src/atlas_io/RecordWriter.h b/atlas_io/src/atlas_io/RecordWriter.h index ef30883bd..14d6ae075 100644 --- a/atlas_io/src/atlas_io/RecordWriter.h +++ b/atlas_io/src/atlas_io/RecordWriter.h @@ -24,6 +24,9 @@ #include "atlas_io/detail/TypeTraits.h" #include "atlas_io/types/array/ArrayReference.h" +#include "atlas_io/types/scalar.h" +#include "atlas_io/types/string.h" + #include "atlas_io/detail/Defaults.h" @@ -88,6 +91,14 @@ class RecordWriter { set(key, Encoder{value}, config); } + void set(const Key& key, const char* value, const eckit::Configuration& config = NoConfig()) { + set(key, Encoder{std::string(value)}, config); + } + + void set(const Key& key, const std::string& value, const eckit::Configuration& config = NoConfig()) { + set(key, Encoder{std::string(value)}, config); + } + /// @brief Write new record to path size_t write(const eckit::PathName&, Mode = Mode::write) const; diff --git a/atlas_io/src/atlas_io/Trace.cc b/atlas_io/src/atlas_io/Trace.cc index 468626cae..5e52b7dc9 100644 --- a/atlas_io/src/atlas_io/Trace.cc +++ b/atlas_io/src/atlas_io/Trace.cc @@ -16,20 +16,26 @@ namespace atlas { namespace io { atlas::io::Trace::Trace(const eckit::CodeLocation& loc) { - for (auto& hook : TraceHookRegistry::instance().hooks) { - hooks_.emplace_back(hook(loc, loc.func())); + for (size_t id = 0; id < TraceHookRegistry::size(); ++id) { + if (TraceHookRegistry::enabled(id)) { + hooks_.emplace_back(TraceHookRegistry::hook(id)(loc, loc.func())); + } } } Trace::Trace(const eckit::CodeLocation& loc, const std::string& title) { - for (auto& hook : TraceHookRegistry::instance().hooks) { - hooks_.emplace_back(hook(loc, title)); + for (size_t id = 0; id < TraceHookRegistry::size(); ++id) { + if (TraceHookRegistry::enabled(id)) { + hooks_.emplace_back(TraceHookRegistry::hook(id)(loc, title)); + } } } Trace::Trace(const eckit::CodeLocation& loc, const std::string& title, const Labels& labels) { - for (auto& hook : TraceHookRegistry::instance().hooks) { - hooks_.emplace_back(hook(loc, title)); + for (size_t id = 0; id < TraceHookRegistry::size(); ++id) { + if (TraceHookRegistry::enabled(id)) { + hooks_.emplace_back(TraceHookRegistry::hook(id)(loc, title)); + } } } diff --git a/atlas_io/src/atlas_io/Trace.h b/atlas_io/src/atlas_io/Trace.h index bf622edc0..f08e1d2bc 100644 --- a/atlas_io/src/atlas_io/Trace.h +++ b/atlas_io/src/atlas_io/Trace.h @@ -31,12 +31,26 @@ struct TraceHook { struct TraceHookRegistry { using TraceHookBuilder = std::function(const eckit::CodeLocation&, const std::string&)>; std::vector hooks; + std::vector enabled_; static TraceHookRegistry& instance() { static TraceHookRegistry instance; return instance; } - static void add(TraceHookBuilder&& hook) { instance().hooks.emplace_back(hook); } - static void add(const TraceHookBuilder& hook) { instance().hooks.emplace_back(hook); } + static size_t add(TraceHookBuilder&& hook) { + instance().hooks.emplace_back(hook); + instance().enabled_.emplace_back(true); + return instance().hooks.size() - 1; + } + static size_t add(const TraceHookBuilder& hook) { + instance().hooks.emplace_back(hook); + instance().enabled_.emplace_back(true); + return instance().hooks.size() - 1; + } + static void enable(size_t id) { instance().enabled_[id] = true; } + static void disable(size_t id) { instance().enabled_[id] = false; } + static bool enabled(size_t id) { return instance().enabled_[id]; } + static size_t size() { return instance().hooks.size(); } + static TraceHookBuilder& hook(size_t id) { return instance().hooks[id]; } private: TraceHookRegistry() = default; diff --git a/atlas_io/src/atlas_io/detail/Encoder.h b/atlas_io/src/atlas_io/detail/Encoder.h index d6a7ae05f..5e203c5d5 100644 --- a/atlas_io/src/atlas_io/detail/Encoder.h +++ b/atlas_io/src/atlas_io/detail/Encoder.h @@ -38,6 +38,10 @@ class Encoder { Encoder(Encoder&& other): self_(std::move(other.self_)) {} + template = 0> + explicit Encoder(const T& x): self_(new EncodableValue(x)) {} + + Encoder& operator=(Encoder&& rhs) { self_ = std::move(rhs.self_); return *this; @@ -61,6 +65,11 @@ class Encoder { struct EncodableValue : Encodable { EncodableValue(Value&& v): value_{std::move(v)} { sfinae::encode_metadata(value_, metadata_, data_size_); } + template = 0> + EncodableValue(const Value& v): value_{v} { + sfinae::encode_metadata(value_, metadata_, data_size_); + } + size_t encode_metadata_(atlas::io::Metadata& metadata) const override { metadata.set(metadata_); return data_size_; diff --git a/atlas_io/src/atlas_io/detail/TypeTraits.h b/atlas_io/src/atlas_io/detail/TypeTraits.h index 93150a247..2ecf759d3 100644 --- a/atlas_io/src/atlas_io/detail/TypeTraits.h +++ b/atlas_io/src/atlas_io/detail/TypeTraits.h @@ -12,6 +12,7 @@ #include +#include "atlas_io/detail/DataType.h" namespace atlas { namespace io { @@ -134,6 +135,20 @@ template using enable_if_move_constructible_decodable_rvalue_t = enable_if_t() && std::is_rvalue_reference() && std::is_move_constructible()>; +template +using enable_if_scalar_t = enable_if_t::value && EnableBool>; + + +template +constexpr bool is_array_datatype() { + return std::is_same::value || std::is_same::value || + std::is_same::value || std::is_same::value || + std::is_same::value || std::is_same::value; +} + +template +using enable_if_array_datatype = typename std::enable_if(), int>::type; + } // namespace io } // namespace atlas diff --git a/atlas_io/src/atlas_io/print/JSONFormat.cc b/atlas_io/src/atlas_io/print/JSONFormat.cc index 51446a160..caafb9a6b 100644 --- a/atlas_io/src/atlas_io/print/JSONFormat.cc +++ b/atlas_io/src/atlas_io/print/JSONFormat.cc @@ -52,11 +52,13 @@ void JSONFormat::print(std::ostream& out) const { } if (print_details_) { - m.set("data.compression.type", item.data.compression()); - m.set("data.compression.size", item.data.compressed_size()); - m.set("data.size", item.data.size()); - m.set("data.byte_order", (item.data.endian() == Endian::little) ? "little endian" : "big endian"); - m.set("data.checksum", item.data.checksum().str()); + if (item.data.size()) { + m.set("data.compression.type", item.data.compression()); + m.set("data.compression.size", item.data.compressed_size()); + m.set("data.size", item.data.size()); + m.set("data.byte_order", (item.data.endian() == Endian::little) ? "little endian" : "big endian"); + m.set("data.checksum", item.data.checksum().str()); + } m.set("version", item.record.version().str()); m.set("created", item.record.created().str()); } diff --git a/atlas_io/src/atlas_io/print/TableFormat.cc b/atlas_io/src/atlas_io/print/TableFormat.cc index 968257341..44d418e84 100644 --- a/atlas_io/src/atlas_io/print/TableFormat.cc +++ b/atlas_io/src/atlas_io/print/TableFormat.cc @@ -135,7 +135,6 @@ class ScalarMetadataPrettyPrint : public MetadataPrettyPrintBase { std::string type = metadata_.getString("type"); ATLAS_IO_ASSERT(type == "scalar"); std::string datatype = metadata_.getString("datatype"); - std::string base64 = metadata_.getString("base64"); out << std::setw(7) << std::left << datatype << ": "; if (datatype == DataType::str()) { out << decode(); diff --git a/atlas_io/src/atlas_io/types/array/ArrayReference.cc b/atlas_io/src/atlas_io/types/array/ArrayReference.cc index f880a7f83..05d173963 100644 --- a/atlas_io/src/atlas_io/types/array/ArrayReference.cc +++ b/atlas_io/src/atlas_io/types/array/ArrayReference.cc @@ -75,8 +75,8 @@ ArrayReference::ArrayReference(ArrayReference&& other): ArrayMetadata(std::move( ArrayReference& ArrayReference::operator=(ArrayReference&& rhs) { ArrayMetadata::operator=(std::move(rhs)); - data_ = rhs.data_; - rhs.data_ = nullptr; + data_ = rhs.data_; + rhs.data_ = nullptr; return *this; } diff --git a/atlas_io/src/atlas_io/types/array/adaptors/StdVectorAdaptor.h b/atlas_io/src/atlas_io/types/array/adaptors/StdVectorAdaptor.h index 18aaeb47f..b2be76a32 100644 --- a/atlas_io/src/atlas_io/types/array/adaptors/StdVectorAdaptor.h +++ b/atlas_io/src/atlas_io/types/array/adaptors/StdVectorAdaptor.h @@ -15,6 +15,7 @@ #include "atlas_io/Data.h" #include "atlas_io/Exceptions.h" #include "atlas_io/Metadata.h" +#include "atlas_io/detail/TypeTraits.h" #include "atlas_io/types/array/ArrayMetadata.h" #include "atlas_io/types/array/ArrayReference.h" @@ -22,7 +23,7 @@ namespace std { //--------------------------------------------------------------------------------------------------------------------- -template +template = 0> void interprete(const std::vector& vector, atlas::io::ArrayReference& out) { using atlas::io::ArrayReference; out = ArrayReference{vector.data(), {int(vector.size())}}; @@ -30,7 +31,7 @@ void interprete(const std::vector& vector, atlas::io::ArrayReference& out) { //--------------------------------------------------------------------------------------------------------------------- -template +template = 0> void decode(const atlas::io::Metadata& m, const atlas::io::Data& encoded, std::vector& out) { atlas::io::ArrayMetadata array(m); if (array.datatype().kind() != atlas::io::ArrayMetadata::DataType::kind()) { diff --git a/atlas_io/src/atlas_io/types/scalar.cc b/atlas_io/src/atlas_io/types/scalar.cc index 0e4ea1ea4..61448c82c 100644 --- a/atlas_io/src/atlas_io/types/scalar.cc +++ b/atlas_io/src/atlas_io/types/scalar.cc @@ -135,16 +135,16 @@ void encode_data(const double&, atlas::io::Data&) {} //--------------------------------------------------------------------------------------------------------------------- void decode(const atlas::io::Metadata& metadata, const atlas::io::Data&, int& value) { - decode_scalar_b64(metadata, value); + decode_scalar(metadata, value); } void decode(const atlas::io::Metadata& metadata, const atlas::io::Data&, long& value) { - decode_scalar_b64(metadata, value); + decode_scalar(metadata, value); } void decode(const atlas::io::Metadata& metadata, const atlas::io::Data&, long long& value) { - decode_scalar_b64(metadata, value); + decode_scalar(metadata, value); } void decode(const atlas::io::Metadata& metadata, const atlas::io::Data&, unsigned long& value) { - decode_scalar_b64(metadata, value); + decode_scalar(metadata, value); } void decode(const atlas::io::Metadata& metadata, const atlas::io::Data&, unsigned long long& value) { decode_scalar_b64(metadata, value); diff --git a/atlas_io/tests/test_io_encoding.cc b/atlas_io/tests/test_io_encoding.cc index 9a248d129..9f327dc77 100644 --- a/atlas_io/tests/test_io_encoding.cc +++ b/atlas_io/tests/test_io_encoding.cc @@ -635,11 +635,21 @@ void test_encode_decode_scalar() { CASE("Encode/Decode scalar") { // bit identical encoding via Base64 string within the metadata! - SECTION("int32") { test_encode_decode_scalar(); } - SECTION("int64") { test_encode_decode_scalar(); } - SECTION("real32") { test_encode_decode_scalar(); } - SECTION("real64") { test_encode_decode_scalar(); } - SECTION("uint64") { test_encode_decode_scalar(); } + SECTION("int32") { + test_encode_decode_scalar(); + } + SECTION("int64") { + test_encode_decode_scalar(); + } + SECTION("real32") { + test_encode_decode_scalar(); + } + SECTION("real64") { + test_encode_decode_scalar(); + } + SECTION("uint64") { + test_encode_decode_scalar(); + } } // ------------------------------------------------------------------------------------------------------- diff --git a/atlas_io/tests/test_io_record.cc b/atlas_io/tests/test_io_record.cc index a1e057cfa..f6d56fa6c 100644 --- a/atlas_io/tests/test_io_record.cc +++ b/atlas_io/tests/test_io_record.cc @@ -137,9 +137,15 @@ CASE("Write records, each in separate file (offset=0)") { write_length(length, path + ".length"); }; - SECTION("record1.atlas" + suffix()) { write_record(globals::record1.data, "record1.atlas" + suffix()); } - SECTION("record2.atlas" + suffix()) { write_record(globals::record2.data, "record2.atlas" + suffix()); } - SECTION("record3.atlas" + suffix()) { write_record(globals::record3.data, "record3.atlas" + suffix()); } + SECTION("record1.atlas" + suffix()) { + write_record(globals::record1.data, "record1.atlas" + suffix()); + } + SECTION("record2.atlas" + suffix()) { + write_record(globals::record2.data, "record2.atlas" + suffix()); + } + SECTION("record3.atlas" + suffix()) { + write_record(globals::record3.data, "record3.atlas" + suffix()); + } } //----------------------------------------------------------------------------- diff --git a/cmake/features/TRANS.cmake b/cmake/features/TRANS.cmake index 71e19487c..4136d0f26 100644 --- a/cmake/features/TRANS.cmake +++ b/cmake/features/TRANS.cmake @@ -2,7 +2,7 @@ set( atlas_HAVE_ECTRANS 0 ) if( ENABLE_TRANS OR NOT DEFINED ENABLE_TRANS ) - find_package( ectrans 1.0 COMPONENTS transi double QUIET ) + find_package( ectrans 1.1 COMPONENTS transi double QUIET ) if( TARGET transi_dp ) set( transi_FOUND TRUE ) if( NOT TARGET transi ) diff --git a/src/apps/atlas-grids.cc b/src/apps/atlas-grids.cc index 47b798ea1..953795985 100644 --- a/src/apps/atlas-grids.cc +++ b/src/apps/atlas-grids.cc @@ -282,9 +282,10 @@ int AtlasGrids::execute(const Args& args) { Log::info() << " x : [ " << std::setw(10) << std::fixed << structuredgrid.xspace().min() << " , " << std::setw(10) << std::fixed << structuredgrid.xspace().max() << " ] deg" << std::endl; - Log::info() << " y : [ " << std::setw(10) << std::fixed << structuredgrid.yspace().min() << " , " - << std::setw(10) << std::fixed << structuredgrid.yspace().max() << " ] deg" - << std::endl; + double ymin = std::min(structuredgrid.yspace().front(), structuredgrid.yspace().back()); + double ymax = std::max(structuredgrid.yspace().front(), structuredgrid.yspace().back()); + Log::info() << " y : [ " << std::setw(10) << std::fixed << ymin << " , " << std::setw(10) + << std::fixed << ymax << " ] deg" << std::endl; } auto it = grid.lonlat().begin(); Log::info() << " lonlat(first) : " << *it << std::endl; diff --git a/src/apps/atlas-meshgen.cc b/src/apps/atlas-meshgen.cc index f5d6d635c..ecc478e5e 100644 --- a/src/apps/atlas-meshgen.cc +++ b/src/apps/atlas-meshgen.cc @@ -88,6 +88,9 @@ MeshGenerator make_meshgenerator(const Grid& grid, const AtlasTool::Args& args) config.set("halo", args.getInt("halo")); } + if (args.has("angle")) { + config.set("angle", args.getDouble("angle")); + } return MeshGenerator{config}; } diff --git a/src/atlas/array/helpers/ArrayAssigner.h b/src/atlas/array/helpers/ArrayAssigner.h index 46190fc96..bd20d13c7 100644 --- a/src/atlas/array/helpers/ArrayAssigner.h +++ b/src/atlas/array/helpers/ArrayAssigner.h @@ -83,6 +83,14 @@ struct array_assigner { ATLAS_ASSERT(it = iterable.end()); } + template + static void apply(IndexView& arr, const Iterable& iterable) { + typename Iterable::const_iterator it = iterable.begin(); + array_assigner_impl::apply(arr, it); + ATLAS_ASSERT(it = iterable.end()); + } + + static void apply(LocalView& arr, Value value) { array_assigner_impl::apply(arr, value); // Note: no need to apply variadic pack (idxs...) diff --git a/src/atlas/array/native/NativeIndexView.cc b/src/atlas/array/native/NativeIndexView.cc index 87f6cd99a..af7b3d809 100644 --- a/src/atlas/array/native/NativeIndexView.cc +++ b/src/atlas/array/native/NativeIndexView.cc @@ -11,12 +11,17 @@ #include #include "atlas/array/native/NativeIndexView.h" +#include "atlas/array/helpers/ArrayAssigner.h" //------------------------------------------------------------------------------------------------------ namespace atlas { namespace array { +#undef ENABLE_IF_NON_CONST +#define ENABLE_IF_NON_CONST \ + template ::value && EnableBool), int>::type*> + //------------------------------------------------------------------------------------------------------ template @@ -42,12 +47,20 @@ void IndexView::dump(std::ostream& os) const { os << "]" << std::endl; } +template +ENABLE_IF_NON_CONST void IndexView::assign(const std::initializer_list& list) { + helpers::array_assigner::apply(*this, list); +} + + //------------------------------------------------------------------------------------------------------ // Explicit template instatiation #define EXPLICIT_TEMPLATE_INSTANTIATION_TYPE_RANK(TYPE, RANK) \ template class IndexView; \ - template class IndexView; + template class IndexView; \ + template void IndexView::assign(std::initializer_list const&); + #define EXPLICIT_TEMPLATE_INSTANTIATION(RANK) \ EXPLICIT_TEMPLATE_INSTANTIATION_TYPE_RANK(int, RANK) \ diff --git a/src/atlas/array/native/NativeIndexView.h b/src/atlas/array/native/NativeIndexView.h index 7b2b44abf..37b3b2cfe 100644 --- a/src/atlas/array/native/NativeIndexView.h +++ b/src/atlas/array/native/NativeIndexView.h @@ -39,6 +39,8 @@ #pragma once #include +#include +#include #include "atlas/array/ArrayUtil.h" #include "atlas/library/config.h" @@ -108,6 +110,11 @@ class FortranIndex { #define TO_FORTRAN #endif + +#define ENABLE_IF_NON_CONST \ + template ::value && EnableBool), int>::type* = nullptr> + //------------------------------------------------------------------------------------------------------ ///@brief Multidimensional access to Array or Field objects containing index fields @@ -152,6 +159,14 @@ class IndexView { return data_[index(idx...)] FROM_FORTRAN; } + ENABLE_IF_NON_CONST + void assign(const std::initializer_list& list); + + template + idx_t shape(Int idx) const { + return shape_[idx]; + } + private: // -- Private methods @@ -223,7 +238,7 @@ class LocalIndexView : public IndexView { #undef INDEX_REF #undef FROM_FORTRAN #undef TO_FORTRAN - +#undef ENABLE_IF_NON_CONST //------------------------------------------------------------------------------------------------------ diff --git a/src/atlas/functionspace/CellColumns.cc b/src/atlas/functionspace/CellColumns.cc index 7da6a2080..85bb5c062 100644 --- a/src/atlas/functionspace/CellColumns.cc +++ b/src/atlas/functionspace/CellColumns.cc @@ -823,30 +823,6 @@ const parallel::HaloExchange& CellColumns::halo_exchange() const { return functionspace_->halo_exchange(); } -void CellColumns::gather(const FieldSet& local, FieldSet& global) const { - functionspace_->gather(local, global); -} - -void CellColumns::gather(const Field& local, Field& global) const { - functionspace_->gather(local, global); -} - -const parallel::GatherScatter& CellColumns::gather() const { - return functionspace_->gather(); -} - -void CellColumns::scatter(const FieldSet& global, FieldSet& local) const { - functionspace_->scatter(global, local); -} - -void CellColumns::scatter(const Field& global, Field& local) const { - functionspace_->scatter(global, local); -} - -const parallel::GatherScatter& CellColumns::scatter() const { - return functionspace_->scatter(); -} - std::string CellColumns::checksum(const FieldSet& fieldset) const { return functionspace_->checksum(fieldset); } diff --git a/src/atlas/functionspace/CellColumns.h b/src/atlas/functionspace/CellColumns.h index 1f84c6ea5..ff620d9c9 100644 --- a/src/atlas/functionspace/CellColumns.h +++ b/src/atlas/functionspace/CellColumns.h @@ -43,12 +43,12 @@ class CellColumns : public functionspace::FunctionSpaceImpl { public: CellColumns(const Mesh&, const eckit::Configuration& = util::NoConfig()); - virtual ~CellColumns() override; + ~CellColumns() override; static std::string static_type() { return "CellColumns"; } - virtual std::string type() const override { return static_type(); } + std::string type() const override { return static_type(); } - virtual std::string distribution() const override; + std::string distribution() const override; idx_t nb_cells() const; idx_t nb_cells_global() const; // Only on MPI rank 0, will this be different from 0 @@ -63,29 +63,29 @@ class CellColumns : public functionspace::FunctionSpaceImpl { // -- Field creation methods - virtual Field createField(const eckit::Configuration&) const override; + Field createField(const eckit::Configuration&) const override; - virtual Field createField(const Field&, const eckit::Configuration&) const override; + Field createField(const Field&, const eckit::Configuration&) const override; // -- Parallelisation aware methods - virtual void haloExchange(const FieldSet&, bool on_device = false) const override; - virtual void haloExchange(const Field&, bool on_device = false) const override; + void haloExchange(const FieldSet&, bool on_device = false) const override; + void haloExchange(const Field&, bool on_device = false) const override; const parallel::HaloExchange& halo_exchange() const; - void gather(const FieldSet&, FieldSet&) const; - void gather(const Field&, Field&) const; - const parallel::GatherScatter& gather() const; + void gather(const FieldSet&, FieldSet&) const override; + void gather(const Field&, Field&) const override; + const parallel::GatherScatter& gather() const override; - void scatter(const FieldSet&, FieldSet&) const; - void scatter(const Field&, Field&) const; - const parallel::GatherScatter& scatter() const; + void scatter(const FieldSet&, FieldSet&) const override; + void scatter(const Field&, Field&) const override; + const parallel::GatherScatter& scatter() const override; std::string checksum(const FieldSet&) const; std::string checksum(const Field&) const; const parallel::Checksum& checksum() const; - virtual idx_t size() const override { return nb_cells_; } + idx_t size() const override { return nb_cells_; } Field lonlat() const override; @@ -102,7 +102,7 @@ class CellColumns : public functionspace::FunctionSpaceImpl { idx_t config_levels(const eckit::Configuration&) const; array::ArrayShape config_shape(const eckit::Configuration&) const; void set_field_metadata(const eckit::Configuration&, Field&) const; - virtual size_t footprint() const override; + size_t footprint() const override; private: // data Mesh mesh_; // non-const because functionspace may modify mesh @@ -183,14 +183,6 @@ class CellColumns : public FunctionSpace { // -- Parallelisation aware methods const parallel::HaloExchange& halo_exchange() const; - void gather(const FieldSet&, FieldSet&) const; - void gather(const Field&, Field&) const; - const parallel::GatherScatter& gather() const; - - void scatter(const FieldSet&, FieldSet&) const; - void scatter(const Field&, Field&) const; - const parallel::GatherScatter& scatter() const; - std::string checksum(const FieldSet&) const; std::string checksum(const Field&) const; const parallel::Checksum& checksum() const; diff --git a/src/atlas/functionspace/EdgeColumns.cc b/src/atlas/functionspace/EdgeColumns.cc index 281e866de..122f98efb 100644 --- a/src/atlas/functionspace/EdgeColumns.cc +++ b/src/atlas/functionspace/EdgeColumns.cc @@ -806,30 +806,6 @@ const parallel::HaloExchange& EdgeColumns::halo_exchange() const { return functionspace_->halo_exchange(); } -void EdgeColumns::gather(const FieldSet& local, FieldSet& global) const { - functionspace_->gather(local, global); -} - -void EdgeColumns::gather(const Field& local, Field& global) const { - functionspace_->gather(local, global); -} - -const parallel::GatherScatter& EdgeColumns::gather() const { - return functionspace_->gather(); -} - -void EdgeColumns::scatter(const FieldSet& global, FieldSet& local) const { - functionspace_->scatter(global, local); -} - -void EdgeColumns::scatter(const Field& global, Field& local) const { - functionspace_->scatter(global, local); -} - -const parallel::GatherScatter& EdgeColumns::scatter() const { - return functionspace_->scatter(); -} - std::string EdgeColumns::checksum(const FieldSet& fieldset) const { return functionspace_->checksum(fieldset); } diff --git a/src/atlas/functionspace/EdgeColumns.h b/src/atlas/functionspace/EdgeColumns.h index 9a5ae18d1..c25fd085c 100644 --- a/src/atlas/functionspace/EdgeColumns.h +++ b/src/atlas/functionspace/EdgeColumns.h @@ -43,11 +43,11 @@ class EdgeColumns : public functionspace::FunctionSpaceImpl { public: EdgeColumns(const Mesh&, const eckit::Configuration& = util::NoConfig()); - virtual ~EdgeColumns() override; + ~EdgeColumns() override; - virtual std::string type() const override { return "Edges"; } + std::string type() const override { return "Edges"; } - virtual std::string distribution() const override; + std::string distribution() const override; idx_t nb_edges() const; idx_t nb_edges_global() const; // Only on MPI rank 0, will this be different from 0 @@ -61,29 +61,29 @@ class EdgeColumns : public functionspace::FunctionSpaceImpl { // -- Field creation methods - virtual Field createField(const eckit::Configuration&) const override; + Field createField(const eckit::Configuration&) const override; - virtual Field createField(const Field&, const eckit::Configuration&) const override; + Field createField(const Field&, const eckit::Configuration&) const override; // -- Parallelisation aware methods - virtual void haloExchange(const FieldSet&, bool on_device = false) const override; - virtual void haloExchange(const Field&, bool on_device = false) const override; + void haloExchange(const FieldSet&, bool on_device = false) const override; + void haloExchange(const Field&, bool on_device = false) const override; const parallel::HaloExchange& halo_exchange() const; - void gather(const FieldSet&, FieldSet&) const; - void gather(const Field&, Field&) const; - const parallel::GatherScatter& gather() const; + void gather(const FieldSet&, FieldSet&) const override; + void gather(const Field&, Field&) const override; + const parallel::GatherScatter& gather() const override; - void scatter(const FieldSet&, FieldSet&) const; - void scatter(const Field&, Field&) const; - const parallel::GatherScatter& scatter() const; + void scatter(const FieldSet&, FieldSet&) const override; + void scatter(const Field&, Field&) const override; + const parallel::GatherScatter& scatter() const override; std::string checksum(const FieldSet&) const; std::string checksum(const Field&) const; const parallel::Checksum& checksum() const; - virtual idx_t size() const override { return nb_edges_; } + idx_t size() const override { return nb_edges_; } Field lonlat() const override; @@ -98,7 +98,7 @@ class EdgeColumns : public functionspace::FunctionSpaceImpl { idx_t config_levels(const eckit::Configuration&) const; array::ArrayShape config_shape(const eckit::Configuration&) const; void set_field_metadata(const eckit::Configuration&, Field&) const; - virtual size_t footprint() const override; + size_t footprint() const override; private: // data Mesh mesh_; // non-const because functionspace may modify mesh @@ -173,14 +173,6 @@ class EdgeColumns : public FunctionSpace { // -- Parallelisation aware methods const parallel::HaloExchange& halo_exchange() const; - void gather(const FieldSet&, FieldSet&) const; - void gather(const Field&, Field&) const; - const parallel::GatherScatter& gather() const; - - void scatter(const FieldSet&, FieldSet&) const; - void scatter(const Field&, Field&) const; - const parallel::GatherScatter& scatter() const; - std::string checksum(const FieldSet&) const; std::string checksum(const Field&) const; const parallel::Checksum& checksum() const; diff --git a/src/atlas/functionspace/FunctionSpace.cc b/src/atlas/functionspace/FunctionSpace.cc index 4a85c01af..d56a30dc9 100644 --- a/src/atlas/functionspace/FunctionSpace.cc +++ b/src/atlas/functionspace/FunctionSpace.cc @@ -86,6 +86,30 @@ void FunctionSpace::adjointHaloExchange(const FieldSet& fields, bool on_device) get()->adjointHaloExchange(fields, on_device); } +void FunctionSpace::gather(const FieldSet& local, FieldSet& global) const { + get()->gather(local, global); +} + +void FunctionSpace::gather(const Field& local, Field& global) const { + get()->gather(local, global); +} + +void FunctionSpace::scatter(const FieldSet& global, FieldSet& local) const { + get()->scatter(global, local); +} + +void FunctionSpace::scatter(const Field& global, Field& local) const { + get()->scatter(global, local); +} + +const parallel::GatherScatter& FunctionSpace::gather() const { + return get()->gather(); +} + +const parallel::GatherScatter& FunctionSpace::scatter() const { + return get()->scatter(); +} + const util::PartitionPolygon& FunctionSpace::polygon(idx_t halo) const { return get()->polygon(halo); } diff --git a/src/atlas/functionspace/FunctionSpace.h b/src/atlas/functionspace/FunctionSpace.h index 63388fd33..220adf4b1 100644 --- a/src/atlas/functionspace/FunctionSpace.h +++ b/src/atlas/functionspace/FunctionSpace.h @@ -26,6 +26,10 @@ class Projection; namespace functionspace { class FunctionSpaceImpl; } +namespace parallel { +class GatherScatter; +} // namespace parallel + namespace util { class PartitionPolygon; class PartitionPolygons; @@ -63,6 +67,12 @@ class FunctionSpace : DOXYGEN_HIDE(public util::ObjectHandlehalo_exchange(); } -void NodeColumns::gather(const FieldSet& local, FieldSet& global) const { - functionspace_->gather(local, global); -} - -void NodeColumns::gather(const Field& local, Field& global) const { - functionspace_->gather(local, global); -} - -const parallel::GatherScatter& NodeColumns::gather() const { - return functionspace_->gather(); -} - -void NodeColumns::scatter(const FieldSet& global, FieldSet& local) const { - functionspace_->scatter(global, local); -} - -void NodeColumns::scatter(const Field& global, Field& local) const { - functionspace_->scatter(global, local); -} - -const parallel::GatherScatter& NodeColumns::scatter() const { - return functionspace_->scatter(); -} - std::string NodeColumns::checksum(const FieldSet& fieldset) const { return functionspace_->checksum(fieldset); } diff --git a/src/atlas/functionspace/NodeColumns.h b/src/atlas/functionspace/NodeColumns.h index 2f8217fe7..559905805 100644 --- a/src/atlas/functionspace/NodeColumns.h +++ b/src/atlas/functionspace/NodeColumns.h @@ -52,12 +52,12 @@ class NodeColumns : public functionspace::FunctionSpaceImpl { NodeColumns(Mesh mesh, const eckit::Configuration&); NodeColumns(Mesh mesh); - virtual ~NodeColumns() override; + ~NodeColumns() override; static std::string static_type() { return "NodeColumns"; } - virtual std::string type() const override { return static_type(); } + std::string type() const override { return static_type(); } - virtual std::string distribution() const override; + std::string distribution() const override; idx_t nb_nodes() const; idx_t nb_nodes_global() const; // All MPI ranks will have same output @@ -73,9 +73,9 @@ class NodeColumns : public functionspace::FunctionSpaceImpl { using FunctionSpaceImpl::createField; /// @brief Create a field - virtual Field createField(const eckit::Configuration&) const override; + Field createField(const eckit::Configuration&) const override; - virtual Field createField(const Field&, const eckit::Configuration&) const override; + Field createField(const Field&, const eckit::Configuration&) const override; // -- Parallelisation aware methods @@ -88,13 +88,13 @@ class NodeColumns : public functionspace::FunctionSpaceImpl { void adjointHaloExchange(const FieldSet&, bool on_device = false) const override; void adjointHaloExchange(const Field&, bool on_device = false) const override; - void gather(const FieldSet&, FieldSet&) const; - void gather(const Field&, Field&) const; - const parallel::GatherScatter& gather() const; + void gather(const FieldSet&, FieldSet&) const override; + void gather(const Field&, Field&) const override; + const parallel::GatherScatter& gather() const override; - void scatter(const FieldSet&, FieldSet&) const; - void scatter(const Field&, Field&) const; - const parallel::GatherScatter& scatter() const; + void scatter(const FieldSet&, FieldSet&) const override; + void scatter(const Field&, Field&) const override; + const parallel::GatherScatter& scatter() const override; std::string checksum(const FieldSet&) const; std::string checksum(const Field&) const; @@ -497,14 +497,6 @@ class NodeColumns : public FunctionSpace { void haloExchange(const Field&, bool on_device = false) const; const parallel::HaloExchange& halo_exchange() const; - void gather(const FieldSet&, FieldSet&) const; - void gather(const Field&, Field&) const; - const parallel::GatherScatter& gather() const; - - void scatter(const FieldSet&, FieldSet&) const; - void scatter(const Field&, Field&) const; - const parallel::GatherScatter& scatter() const; - std::string checksum(const FieldSet&) const; std::string checksum(const Field&) const; const parallel::Checksum& checksum() const; diff --git a/src/atlas/functionspace/PointCloud.cc b/src/atlas/functionspace/PointCloud.cc index c9cd23a95..6fd4fafef 100644 --- a/src/atlas/functionspace/PointCloud.cc +++ b/src/atlas/functionspace/PointCloud.cc @@ -12,9 +12,12 @@ #include "atlas/functionspace/PointCloud.h" #include "atlas/array.h" #include "atlas/field/Field.h" +#include "atlas/field/FieldSet.h" #include "atlas/grid/Grid.h" #include "atlas/grid/Iterator.h" #include "atlas/option/Options.h" +#include "atlas/parallel/HaloExchange.h" +#include "atlas/parallel/mpi/mpi.h" #include "atlas/runtime/Exception.h" #include "atlas/util/CoordinateEnums.h" #include "atlas/util/Metadata.h" @@ -57,6 +60,18 @@ PointCloud::PointCloud(const Field& lonlat): lonlat_(lonlat) {} PointCloud::PointCloud(const Field& lonlat, const Field& ghost): lonlat_(lonlat), ghost_(ghost) {} +PointCloud::PointCloud(const FieldSet & flds): lonlat_(flds["lonlat"]), + ghost_(flds["ghost"]), remote_index_(flds["remote_index"]), partition_(flds["partition"]) +{ + ATLAS_ASSERT(ghost_.size() == remote_index_.size()); + ATLAS_ASSERT(ghost_.size() == partition_.size()); + halo_exchange_.reset(new parallel::HaloExchange()); + halo_exchange_->setup(array::make_view( partition_).data(), + array::make_view(remote_index_).data(), + REMOTE_IDX_BASE, + ghost_.size()); +} + PointCloud::PointCloud(const Grid& grid) { lonlat_ = Field("lonlat", array::make_datatype(), array::make_shape(grid.size(), 2)); auto lonlat = array::make_view(lonlat_); @@ -121,6 +136,12 @@ std::string PointCloud::config_name(const eckit::Configuration& config) const { return name; } + +const parallel::HaloExchange& PointCloud::halo_exchange() const { + return *halo_exchange_; +} + + void PointCloud::set_field_metadata(const eckit::Configuration& config, Field& field) const { field.set_functionspace(this); @@ -161,7 +182,7 @@ Field PointCloud::createField(const Field& other, const eckit::Configuration& co } std::string PointCloud::distribution() const { - return std::string("serial"); + return (partition_ ? std::string("pointcloud") : std::string("serial")); } static const array::Array& get_dummy() { @@ -206,6 +227,113 @@ template class PointCloud::IteratorT; template class PointCloud::IteratorT; template class PointCloud::IteratorT; + +namespace { + +template +void dispatch_haloExchange(Field& field, const parallel::HaloExchange& halo_exchange, bool on_device) { + if (field.datatype() == array::DataType::kind()) { + halo_exchange.template execute(field.array(), on_device); + } + else if (field.datatype() == array::DataType::kind()) { + halo_exchange.template execute(field.array(), on_device); + } + else if (field.datatype() == array::DataType::kind()) { + halo_exchange.template execute(field.array(), on_device); + } + else if (field.datatype() == array::DataType::kind()) { + halo_exchange.template execute(field.array(), on_device); + } + else { + throw_Exception("datatype not supported", Here()); + } + field.set_dirty(false); +} + +template +void dispatch_adjointHaloExchange(Field& field, const parallel::HaloExchange& halo_exchange, bool on_device) { + if (field.datatype() == array::DataType::kind()) { + halo_exchange.template execute_adjoint(field.array(), on_device); + } + else if (field.datatype() == array::DataType::kind()) { + halo_exchange.template execute_adjoint(field.array(), on_device); + } + else if (field.datatype() == array::DataType::kind()) { + halo_exchange.template execute_adjoint(field.array(), on_device); + } + else if (field.datatype() == array::DataType::kind()) { + halo_exchange.template execute_adjoint(field.array(), on_device); + } + else { + throw_Exception("datatype not supported", Here()); + } + field.set_dirty(false); +} + +} // namespace + +void PointCloud::haloExchange(const FieldSet& fieldset, bool on_device) const { + if (halo_exchange_) { + for (idx_t f = 0; f < fieldset.size(); ++f) { + Field& field = const_cast(fieldset)[f]; + switch (field.rank()) { + case 1: + dispatch_haloExchange<1>(field, halo_exchange(), on_device); + break; + case 2: + dispatch_haloExchange<2>(field, halo_exchange(), on_device); + break; + case 3: + dispatch_haloExchange<3>(field, halo_exchange(), on_device); + break; + case 4: + dispatch_haloExchange<4>(field, halo_exchange(), on_device); + break; + default: + throw_Exception("Rank not supported", Here()); + } + field.set_dirty(false); + } + } +} + +void PointCloud::haloExchange(const Field& field, bool on_device) const { + FieldSet fieldset; + fieldset.add(field); + haloExchange(fieldset, on_device); +} + +void PointCloud::adjointHaloExchange(const FieldSet& fieldset, bool on_device) const { + if (halo_exchange_) { + for (idx_t f = 0; f < fieldset.size(); ++f) { + Field& field = const_cast(fieldset)[f]; + switch (field.rank()) { + case 1: + dispatch_adjointHaloExchange<1>(field, halo_exchange(), on_device); + break; + case 2: + dispatch_adjointHaloExchange<2>(field, halo_exchange(), on_device); + break; + case 3: + dispatch_adjointHaloExchange<3>(field, halo_exchange(), on_device); + break; + case 4: + dispatch_adjointHaloExchange<4>(field, halo_exchange(), on_device); + break; + default: + throw_Exception("Rank not supported", Here()); + } + } + } +} + +void PointCloud::adjointHaloExchange(const Field& field, bool) const { + FieldSet fieldset; + fieldset.add(field); + adjointHaloExchange(fieldset); +} + + } // namespace detail PointCloud::PointCloud(const FunctionSpace& functionspace): @@ -214,13 +342,15 @@ PointCloud::PointCloud(const FunctionSpace& functionspace): PointCloud::PointCloud(const Field& points): FunctionSpace(new detail::PointCloud(points)), functionspace_(dynamic_cast(get())) {} +PointCloud::PointCloud(const FieldSet& points): + FunctionSpace(new detail::PointCloud(points)), functionspace_(dynamic_cast(get())) {} + PointCloud::PointCloud(const std::vector& points): FunctionSpace(new detail::PointCloud(points)), functionspace_(dynamic_cast(get())) {} PointCloud::PointCloud(const std::vector& points): FunctionSpace(new detail::PointCloud(points)), functionspace_(dynamic_cast(get())) {} - PointCloud::PointCloud(const std::initializer_list>& points): FunctionSpace((points.begin()->size() == 2 ? new detail::PointCloud{std::vector(points.begin(), points.end())} diff --git a/src/atlas/functionspace/PointCloud.h b/src/atlas/functionspace/PointCloud.h index e70c2f3b1..03258f4bc 100644 --- a/src/atlas/functionspace/PointCloud.h +++ b/src/atlas/functionspace/PointCloud.h @@ -10,14 +10,24 @@ #pragma once +#include + #include "atlas/array/ArrayView.h" #include "atlas/field/Field.h" +#include "atlas/field/FieldSet.h" #include "atlas/functionspace/FunctionSpace.h" #include "atlas/functionspace/detail/FunctionSpaceImpl.h" +#include "atlas/parallel/HaloExchange.h" #include "atlas/runtime/Exception.h" #include "atlas/util/Config.h" #include "atlas/util/Point.h" +namespace atlas { +namespace parallel { +class HaloExchange; +} // namespace parallel +} // namespace atlas + namespace atlas { class Grid; @@ -33,23 +43,31 @@ class PointCloud : public functionspace::FunctionSpaceImpl { PointCloud(const std::vector&); PointCloud(const Field& lonlat); PointCloud(const Field& lonlat, const Field& ghost); + + PointCloud(const FieldSet&); // assuming lonlat ghost ridx and partition present. + PointCloud(const Grid&); - virtual ~PointCloud() override {} - virtual std::string type() const override { return "PointCloud"; } - virtual operator bool() const override { return true; } - virtual size_t footprint() const override { return sizeof(*this); } - virtual std::string distribution() const override; + ~PointCloud() override {} + std::string type() const override { return "PointCloud"; } + operator bool() const override { return true; } + size_t footprint() const override { return sizeof(*this); } + std::string distribution() const override; Field lonlat() const override { return lonlat_; } const Field& vertical() const { return vertical_; } Field ghost() const override; virtual idx_t size() const override { return lonlat_.shape(0); } using FunctionSpaceImpl::createField; - virtual Field createField(const eckit::Configuration&) const override; - virtual Field createField(const Field&, const eckit::Configuration&) const override; + Field createField(const eckit::Configuration&) const override; + Field createField(const Field&, const eckit::Configuration&) const override; + + void haloExchange(const FieldSet&, bool on_device = false) const override; + void haloExchange(const Field&, bool on_device = false) const override; + + void adjointHaloExchange(const FieldSet&, bool on_device = false) const override; + void adjointHaloExchange(const Field&, bool on_device = false) const override; - void haloExchange(const FieldSet&, bool /*on_device*/ = false) const override {} - void haloExchange(const Field&, bool /*on_device*/ = false) const override {} + const parallel::HaloExchange& halo_exchange() const; template class IteratorT { @@ -129,10 +147,15 @@ class PointCloud : public functionspace::FunctionSpaceImpl { void set_field_metadata(const eckit::Configuration& config, Field& field) const; + + private: Field lonlat_; Field vertical_; mutable Field ghost_; + Field remote_index_; + Field partition_; + std::unique_ptr halo_exchange_; idx_t levels_{0}; }; @@ -146,6 +169,7 @@ class PointCloud : public FunctionSpace { public: PointCloud(const FunctionSpace&); PointCloud(const Field& points); + PointCloud(const FieldSet& flds); PointCloud(const std::vector&); PointCloud(const std::vector&); PointCloud(const std::initializer_list>&); diff --git a/src/atlas/functionspace/Spectral.cc b/src/atlas/functionspace/Spectral.cc index 3564f542c..f5d401a1e 100644 --- a/src/atlas/functionspace/Spectral.cc +++ b/src/atlas/functionspace/Spectral.cc @@ -467,22 +467,6 @@ int Spectral::truncation() const { return functionspace_->truncation(); } -void Spectral::gather(const FieldSet& local_fieldset, FieldSet& global_fieldset) const { - functionspace_->gather(local_fieldset, global_fieldset); -} - -void Spectral::gather(const Field& local, Field& global) const { - functionspace_->gather(local, global); -} - -void Spectral::scatter(const FieldSet& global_fieldset, FieldSet& local_fieldset) const { - functionspace_->scatter(global_fieldset, local_fieldset); -} - -void Spectral::scatter(const Field& global, Field& local) const { - functionspace_->scatter(global, local); -} - std::string Spectral::checksum(const FieldSet& fieldset) const { return functionspace_->checksum(fieldset); } diff --git a/src/atlas/functionspace/Spectral.h b/src/atlas/functionspace/Spectral.h index 231bdaefa..0932c2261 100644 --- a/src/atlas/functionspace/Spectral.h +++ b/src/atlas/functionspace/Spectral.h @@ -65,21 +65,23 @@ class Spectral : public functionspace::FunctionSpaceImpl { Spectral(const trans::Trans&, const eckit::Configuration& = util::NoConfig()); - virtual ~Spectral() override; + ~Spectral() override; - virtual std::string type() const override { return "Spectral"; } + std::string type() const override { return "Spectral"; } - virtual std::string distribution() const override; + std::string distribution() const override; using FunctionSpaceImpl::createField; - virtual Field createField(const eckit::Configuration&) const override; - virtual Field createField(const Field&, const eckit::Configuration&) const override; + Field createField(const eckit::Configuration&) const override; + Field createField(const Field&, const eckit::Configuration&) const override; - void gather(const FieldSet&, FieldSet&) const; - void gather(const Field&, Field&) const; + using FunctionSpaceImpl::gather; + void gather(const FieldSet&, FieldSet&) const override; + void gather(const Field&, Field&) const override; - void scatter(const FieldSet&, FieldSet&) const; - void scatter(const Field&, Field&) const; + using FunctionSpaceImpl::scatter; + void scatter(const FieldSet&, FieldSet&) const override; + void scatter(const Field&, Field&) const override; std::string checksum(const FieldSet&) const; std::string checksum(const Field&) const; @@ -180,7 +182,7 @@ class Spectral : public functionspace::FunctionSpaceImpl { idx_t nb_spectral_coefficients_global() const; int truncation() const { return truncation_; } - virtual idx_t size() const override { return nb_spectral_coefficients(); } + idx_t size() const override { return nb_spectral_coefficients(); } private: // methods array::DataType config_datatype(const eckit::Configuration&) const; @@ -221,12 +223,6 @@ class Spectral : public FunctionSpace { operator bool() const { return valid(); } bool valid() const { return functionspace_; } - void gather(const FieldSet&, FieldSet&) const; - void gather(const Field&, Field&) const; - - void scatter(const FieldSet&, FieldSet&) const; - void scatter(const Field&, Field&) const; - std::string checksum(const FieldSet&) const; std::string checksum(const Field&) const; diff --git a/src/atlas/functionspace/StructuredColumns.cc b/src/atlas/functionspace/StructuredColumns.cc index db643cf72..a99dfdba9 100644 --- a/src/atlas/functionspace/StructuredColumns.cc +++ b/src/atlas/functionspace/StructuredColumns.cc @@ -48,22 +48,6 @@ StructuredColumns::StructuredColumns(const Grid& grid, const grid::Distribution& FunctionSpace(new detail::StructuredColumns(grid, distribution, vertical, config)), functionspace_(dynamic_cast(get())) {} -void StructuredColumns::gather(const FieldSet& local, FieldSet& global) const { - functionspace_->gather(local, global); -} - -void StructuredColumns::gather(const Field& local, Field& global) const { - functionspace_->gather(local, global); -} - -void StructuredColumns::scatter(const FieldSet& global, FieldSet& local) const { - functionspace_->scatter(global, local); -} - -void StructuredColumns::scatter(const Field& global, Field& local) const { - functionspace_->scatter(global, local); -} - std::string StructuredColumns::checksum(const FieldSet& fieldset) const { return functionspace_->checksum(fieldset); } diff --git a/src/atlas/functionspace/StructuredColumns.h b/src/atlas/functionspace/StructuredColumns.h index fabdb5aa1..ff458f6d5 100644 --- a/src/atlas/functionspace/StructuredColumns.h +++ b/src/atlas/functionspace/StructuredColumns.h @@ -50,12 +50,6 @@ class StructuredColumns : public FunctionSpace { const StructuredGrid& grid() const { return functionspace_->grid(); } - void gather(const FieldSet&, FieldSet&) const; - void gather(const Field&, Field&) const; - - void scatter(const FieldSet&, FieldSet&) const; - void scatter(const Field&, Field&) const; - std::string checksum(const FieldSet&) const; std::string checksum(const Field&) const; diff --git a/src/atlas/functionspace/detail/FunctionSpaceImpl.cc b/src/atlas/functionspace/detail/FunctionSpaceImpl.cc index 9d24a1453..0bcd33eaa 100644 --- a/src/atlas/functionspace/detail/FunctionSpaceImpl.cc +++ b/src/atlas/functionspace/detail/FunctionSpaceImpl.cc @@ -95,6 +95,30 @@ idx_t FunctionSpaceImpl::nb_partitions() const { ATLAS_NOTIMPLEMENTED; } +void FunctionSpaceImpl::gather(const FieldSet& local, FieldSet& global) const { + ATLAS_NOTIMPLEMENTED; +} + +void FunctionSpaceImpl::gather(const Field& local, Field& global) const { + ATLAS_NOTIMPLEMENTED; +} + +void FunctionSpaceImpl::scatter(const FieldSet& global, FieldSet& local) const { + ATLAS_NOTIMPLEMENTED; +} + +void FunctionSpaceImpl::scatter(const Field& global, Field& local) const { + ATLAS_NOTIMPLEMENTED; +} + +const parallel::GatherScatter& FunctionSpaceImpl::gather() const { + ATLAS_NOTIMPLEMENTED; +} + +const parallel::GatherScatter& FunctionSpaceImpl::scatter() const { + ATLAS_NOTIMPLEMENTED; +} + template Field FunctionSpaceImpl::createField() const; template Field FunctionSpaceImpl::createField() const; diff --git a/src/atlas/functionspace/detail/FunctionSpaceImpl.h b/src/atlas/functionspace/detail/FunctionSpaceImpl.h index 53db155b9..7b08e50d4 100644 --- a/src/atlas/functionspace/detail/FunctionSpaceImpl.h +++ b/src/atlas/functionspace/detail/FunctionSpaceImpl.h @@ -31,6 +31,10 @@ class Metadata; class PartitionPolygon; class PartitionPolygons; } // namespace util +namespace parallel { +class GatherScatter; +} // namespace parallel + } // namespace atlas namespace atlas { @@ -78,6 +82,15 @@ class FunctionSpaceImpl : public util::Object { virtual void adjointHaloExchange(const FieldSet&, bool /*on_device*/ = false) const; virtual void adjointHaloExchange(const Field&, bool /* on_device*/ = false) const; + virtual void gather(const FieldSet&, FieldSet&) const; + virtual void gather(const Field&, Field&) const; + + virtual void scatter(const FieldSet&, FieldSet&) const; + virtual void scatter(const Field&, Field&) const; + + virtual const parallel::GatherScatter& gather() const; + virtual const parallel::GatherScatter& scatter() const; + virtual idx_t size() const = 0; virtual idx_t nb_partitions() const; diff --git a/src/atlas/functionspace/detail/FunctionSpaceInterface.cc b/src/atlas/functionspace/detail/FunctionSpaceInterface.cc index d7c25510a..ed91bd780 100644 --- a/src/atlas/functionspace/detail/FunctionSpaceInterface.cc +++ b/src/atlas/functionspace/detail/FunctionSpaceInterface.cc @@ -84,7 +84,8 @@ void atlas__FunctionSpace__halo_exchange_field(const FunctionSpaceImpl* This, fi //------------------------------------------------------------------------------ -void atlas__FunctionSpace__halo_exchange_fieldset(const FunctionSpaceImpl* This, field::FieldSetImpl* fieldset) { +void atlas__FunctionSpace__halo_exchange_fieldset(const FunctionSpaceImpl* This, + field::FieldSetImpl* fieldset) { ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_FunctionSpace"); ATLAS_ASSERT(fieldset != nullptr, "Cannot access uninitialised atlas_FieldSet"); FieldSet f(fieldset); @@ -93,7 +94,8 @@ void atlas__FunctionSpace__halo_exchange_fieldset(const FunctionSpaceImpl* This, //------------------------------------------------------------------------------ -void atlas__FunctionSpace_adjoint_halo_exchange_field(const FunctionSpaceImpl* This, field::FieldImpl* field) { +void atlas__FunctionSpace__adjoint_halo_exchange_field(const FunctionSpaceImpl* This, + field::FieldImpl* field) { ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_FunctionSpace"); ATLAS_ASSERT(field != nullptr, "Cannot access uninitialised atlas_Field"); Field f(field); diff --git a/src/atlas/functionspace/detail/PointCloudInterface.cc b/src/atlas/functionspace/detail/PointCloudInterface.cc index fe7d53e69..9556bb7a0 100644 --- a/src/atlas/functionspace/detail/PointCloudInterface.cc +++ b/src/atlas/functionspace/detail/PointCloudInterface.cc @@ -35,6 +35,9 @@ const detail::PointCloud* atlas__functionspace__PointCloud__new__lonlat_ghost(co return new detail::PointCloud(Field(lonlat), Field(ghost)); } +const detail::PointCloud* atlas__functionspace__PointCloud__new__fieldset(const field::FieldSetImpl* fset) { + return new detail::PointCloud(FieldSet(fset)); +} const detail::PointCloud* atlas__functionspace__PointCloud__new__grid(const Grid::Implementation* grid) { return new detail::PointCloud(Grid(grid)); diff --git a/src/atlas/functionspace/detail/PointCloudInterface.h b/src/atlas/functionspace/detail/PointCloudInterface.h index 34f0cacf1..0d378796e 100644 --- a/src/atlas/functionspace/detail/PointCloudInterface.h +++ b/src/atlas/functionspace/detail/PointCloudInterface.h @@ -41,6 +41,8 @@ const detail::PointCloud* atlas__functionspace__PointCloud__new__lonlat(const fi const detail::PointCloud* atlas__functionspace__PointCloud__new__lonlat_ghost(const field::FieldImpl* lonlat, const field::FieldImpl* ghost); +const detail::PointCloud* atlas__functionspace__PointCloud__new__fieldset(const field::FieldSetImpl* fset); + const detail::PointCloud* atlas__functionspace__PointCloud__new__grid(const GridImpl* grid); void atlas__functionspace__PointCloud__delete(detail::PointCloud* This); diff --git a/src/atlas/functionspace/detail/StructuredColumns.h b/src/atlas/functionspace/detail/StructuredColumns.h index e66c4fa05..c575c318d 100644 --- a/src/atlas/functionspace/detail/StructuredColumns.h +++ b/src/atlas/functionspace/detail/StructuredColumns.h @@ -78,32 +78,32 @@ class StructuredColumns : public FunctionSpaceImpl { StructuredColumns(const Grid&, const Vertical&, const grid::Partitioner&, const eckit::Configuration& = util::NoConfig()); - virtual ~StructuredColumns() override; + ~StructuredColumns() override; static std::string static_type() { return "StructuredColumns"; } - virtual std::string type() const override { return static_type(); } - virtual std::string distribution() const override; + std::string type() const override { return static_type(); } + std::string distribution() const override; /// @brief Create a Structured field - virtual Field createField(const eckit::Configuration&) const override; + Field createField(const eckit::Configuration&) const override; - virtual Field createField(const Field&, const eckit::Configuration&) const override; + Field createField(const Field&, const eckit::Configuration&) const override; - void gather(const FieldSet&, FieldSet&) const; - void gather(const Field&, Field&) const; + void gather(const FieldSet&, FieldSet&) const override; + void gather(const Field&, Field&) const override; - void scatter(const FieldSet&, FieldSet&) const; - void scatter(const Field&, Field&) const; + void scatter(const FieldSet&, FieldSet&) const override; + void scatter(const Field&, Field&) const override; - virtual void haloExchange(const FieldSet&, bool on_device = false) const override; - virtual void haloExchange(const Field&, bool on_device = false) const override; + void haloExchange(const FieldSet&, bool on_device = false) const override; + void haloExchange(const Field&, bool on_device = false) const override; - virtual void adjointHaloExchange(const FieldSet&, bool on_device = false) const override; - virtual void adjointHaloExchange(const Field&, bool on_device = false) const override; + void adjointHaloExchange(const FieldSet&, bool on_device = false) const override; + void adjointHaloExchange(const Field&, bool on_device = false) const override; idx_t sizeOwned() const { return size_owned_; } idx_t sizeHalo() const { return size_halo_; } - virtual idx_t size() const override { return size_halo_; } + idx_t size() const override { return size_halo_; } idx_t levels() const { return nb_levels_; } @@ -192,8 +192,8 @@ class StructuredColumns : public FunctionSpaceImpl { } [[noreturn]] void throw_outofbounds(idx_t i, idx_t j) const; - const parallel::GatherScatter& gather() const; - const parallel::GatherScatter& scatter() const; + const parallel::GatherScatter& gather() const override; + const parallel::GatherScatter& scatter() const override; const parallel::Checksum& checksum() const; const parallel::HaloExchange& halo_exchange() const; diff --git a/src/atlas/functionspace/detail/StructuredColumns_setup.cc b/src/atlas/functionspace/detail/StructuredColumns_setup.cc index 7fd8b6087..cb348da60 100644 --- a/src/atlas/functionspace/detail/StructuredColumns_setup.cc +++ b/src/atlas/functionspace/detail/StructuredColumns_setup.cc @@ -229,7 +229,7 @@ void StructuredColumns::setup(const grid::Distribution& distribution, const ecki i_begin_halo_.resize(-halo, grid_->ny() - 1 + halo); i_end_halo_.resize(-halo, grid_->ny() - 1 + halo); if (regional) { - j_begin_halo_ = std::max(j_begin_halo_, 0); + j_begin_halo_ = std::max(j_begin_halo_, idx_t{0}); j_end_halo_ = std::min(j_end_halo_, grid_->ny()); } @@ -367,18 +367,18 @@ void StructuredColumns::setup(const grid::Distribution& distribution, const ecki double x = grid_->x(i, j); - double x_next = grid_->x(i + 1, j); - double x_prev = grid_->x(i - 1, j); + double x_next = grid_->x(i + idx_t{1}, j); + double x_prev = grid_->x(i - idx_t{1}, j); idx_t jj_min = j - halo; idx_t jj_max = j + halo; if (regional) { - jj_min = std::max(jj_min, 0); - jj_max = std::min(jj_max, grid_->nx(j) - 1); + jj_min = std::max(jj_min, idx_t{0}); + jj_max = std::min(jj_max, grid_->nx(j) - idx_t{1}); } for (idx_t jj = jj_min; jj <= jj_max; ++jj) { idx_t jjj = compute_j(jj); idx_t nx_jjj = grid_->nx(jjj); - idx_t last = grid_->nx(jjj) - 1; + idx_t last = grid_->nx(jjj) - idx_t{1}; if (i == grid_->nx(j)) { ++last; } @@ -426,14 +426,14 @@ void StructuredColumns::setup(const grid::Distribution& distribution, const ecki idx_t i_plus_halo = iii + halo; if (regional) { - i_minus_halo = std::max(i_minus_halo, 0); - i_plus_halo = std::min(i_plus_halo, grid_->nx(jj) - 1); + i_minus_halo = std::max(i_minus_halo, idx_t{0}); + i_plus_halo = std::min(i_plus_halo, grid_->nx(jj) - idx_t{1}); } imin = std::min(imin, i_minus_halo); imax = std::max(imax, i_plus_halo); i_begin_halo_(jj) = std::min(i_begin_halo_(jj), i_minus_halo); - i_end_halo_(jj) = std::max(i_end_halo_(jj), i_plus_halo + 1); + i_end_halo_(jj) = std::max(i_end_halo_(jj), i_plus_halo + idx_t{1}); } } } diff --git a/src/atlas/grid/detail/partitioner/MatchingMeshPartitionerLonLatPolygon.cc b/src/atlas/grid/detail/partitioner/MatchingMeshPartitionerLonLatPolygon.cc index 5dbc571fd..d06131320 100644 --- a/src/atlas/grid/detail/partitioner/MatchingMeshPartitionerLonLatPolygon.cc +++ b/src/atlas/grid/detail/partitioner/MatchingMeshPartitionerLonLatPolygon.cc @@ -22,6 +22,7 @@ #include "atlas/mesh/Nodes.h" #include "atlas/parallel/mpi/mpi.h" #include "atlas/parallel/omp/fill.h" +#include "atlas/parallel/omp/omp.h" #include "atlas/runtime/Exception.h" #include "atlas/runtime/Log.h" #include "atlas/util/CoordinateEnums.h" @@ -51,27 +52,65 @@ void MatchingMeshPartitionerLonLatPolygon::partition(const Grid& grid, int parti double west = poly.coordinatesMin().x(); double east = poly.coordinatesMax().x(); - comm.allReduceInPlace(west, eckit::mpi::Operation::MIN); - comm.allReduceInPlace(east, eckit::mpi::Operation::MAX); + ATLAS_TRACE_MPI(ALLREDUCE) { + comm.allReduceInPlace(west, eckit::mpi::Operation::MIN); + comm.allReduceInPlace(east, eckit::mpi::Operation::MAX); + } Projection projection = prePartitionedMesh_.projection(); omp::fill(partitioning, partitioning + grid.size(), -1); auto compute = [&](double west) { - size_t i = 0; - - for (PointLonLat P : grid.lonlat()) { - if (partitioning[i] < 0) { - projection.lonlat2xy(P); - P.normalise(west); - partitioning[i] = poly.contains(P) ? mpi_rank : -1; +#if !defined(__NVCOMPILER) + atlas_omp_parallel +#elif (__NVCOMPILER_MAJOR__ >= 21) && (__NVCOMPILER_MINOR__ > 9 ) + // Internal compiler error with nvhpc 21.9: + // + // NVC++-S-0000-Internal compiler error. BAD sptr in var_refsym 0 (MatchingMeshPartitionerLonLatPolygon.cc: 64) (=following line) + // NVC++/x86-64 Linux 21.9-0: compilation completed with severe errors + atlas_omp_parallel +#endif + { + const idx_t num_threads = atlas_omp_get_num_threads(); + const idx_t thread_num = atlas_omp_get_thread_num(); + const idx_t begin = static_cast(thread_num * size_t(grid.size()) / num_threads); + const idx_t end = + static_cast((thread_num + 1) * size_t(grid.size()) / num_threads); + size_t i = begin; + auto it = grid.lonlat().begin() + i; + for(; i < end; ++i, ++it) { + PointLonLat P = *it; + if (partitioning[i] < 0) { + projection.lonlat2xy(P); + P.normalise(west); + partitioning[i] = poly.contains(P) ? mpi_rank : -1; + } } - ++i; } // Synchronize partitioning - comm.allReduceInPlace(partitioning, grid.size(), eckit::mpi::Operation::MAX); + ATLAS_TRACE_MPI(ALLREDUCE) { + comm.allReduceInPlace(partitioning, grid.size(), eckit::mpi::Operation::MAX); + } - return *std::min_element(partitioning, partitioning + grid.size()); + std::vector thread_min(atlas_omp_get_max_threads(),std::numeric_limits::max()); +#if !defined(__NVCOMPILER) + atlas_omp_parallel +#elif (__NVCOMPILER_MAJOR__ >= 21) && (__NVCOMPILER_MINOR__ > 9 ) + // Internal compiler error with nvhpc 21.9: + // + // NVC++-S-0000-Internal compiler error. BAD sptr in var_refsym 0 (MatchingMeshPartitionerLonLatPolygon.cc: 64) (=following line) + // NVC++/x86-64 Linux 21.9-0: compilation completed with severe errors + atlas_omp_parallel +#endif + { + const idx_t num_threads = atlas_omp_get_num_threads(); + const idx_t thread_num = atlas_omp_get_thread_num(); + const idx_t begin = static_cast(thread_num * size_t(grid.size()) / num_threads); + const idx_t end = + static_cast((thread_num + 1) * size_t(grid.size()) / num_threads); + thread_min[thread_num] = *std::min_element(partitioning+begin,partitioning+end); + } + return *std::min_element(thread_min.begin(), thread_min.end()); }; int min = compute(east - 360.); diff --git a/src/atlas/interpolation/method/Method.cc b/src/atlas/interpolation/method/Method.cc index 5c402cb07..70aaeb847 100644 --- a/src/atlas/interpolation/method/Method.cc +++ b/src/atlas/interpolation/method/Method.cc @@ -242,6 +242,11 @@ void Method::check_compatibility(const Field& src, const Field& tgt, const Matri template void Method::interpolate_field(const Field& src, Field& tgt, const Matrix& W) const { + // do nothing if there are no observations to interpolate (W will be NULL + // and would fail the compatibility check) + if (tgt.shape(0) == 0) { + return; + } check_compatibility(src, tgt, W); if (src.rank() == 1) { @@ -260,6 +265,11 @@ void Method::interpolate_field(const Field& src, Field& tgt, const Matrix& W) co template void Method::adjoint_interpolate_field(Field& src, const Field& tgt, const Matrix& W) const { + // do nothing if there are no observations to interpolate (W will be NULL + // and would fail the compatibility check) + if (tgt.shape(0) == 0) { + return; + } check_compatibility(tgt, src, W); if (src.rank() == 1) { diff --git a/src/atlas/interpolation/method/unstructured/UnstructuredBilinearLonLat.cc b/src/atlas/interpolation/method/unstructured/UnstructuredBilinearLonLat.cc index d84bad6c2..9814560f2 100644 --- a/src/atlas/interpolation/method/unstructured/UnstructuredBilinearLonLat.cc +++ b/src/atlas/interpolation/method/unstructured/UnstructuredBilinearLonLat.cc @@ -225,6 +225,13 @@ void UnstructuredBilinearLonLat::setup(const FunctionSpace& source) { idx_t inp_npts = i_nodes.size(); idx_t out_npts = olonlat_->shape(0); + // return early if no output points on this partition reserve is called on + // the triplets but also during the sparseMatrix constructor. This won't + // work for empty matrices + if (out_npts == 0) { + return; + } + array::ArrayView out_ghosts = array::make_view(target_ghost_); idx_t Nelements = meshSource.cells().size(); diff --git a/src/atlas/library/Library.cc b/src/atlas/library/Library.cc index cfc70d41f..10c06063d 100644 --- a/src/atlas/library/Library.cc +++ b/src/atlas/library/Library.cc @@ -51,6 +51,7 @@ static bool feature_MKL() { #include "atlas/library/git_sha1.h" #include "atlas/library/version.h" #include "atlas/parallel/mpi/mpi.h" +#include "atlas/parallel/omp/omp.h" #include "atlas/runtime/Exception.h" #include "atlas/runtime/Log.h" #include "atlas/runtime/Trace.h" @@ -298,13 +299,14 @@ void Library::initialise(const eckit::Parametrisation& config) { init_data_paths(data_paths_); } - atlas::io::TraceHookRegistry::add([](const eckit::CodeLocation& loc, const std::string& title) { - struct Adaptor : public atlas::io::TraceHook { - Adaptor(const eckit::CodeLocation& loc, const std::string& title): trace{loc, title} {} - atlas::Trace trace; - }; - return std::unique_ptr(new Adaptor{loc, title}); - }); + atlas_io_trace_hook_ = + atlas::io::TraceHookRegistry::add([](const eckit::CodeLocation& loc, const std::string& title) { + struct Adaptor : public atlas::io::TraceHook { + Adaptor(const eckit::CodeLocation& loc, const std::string& title): trace{loc, title} {} + atlas::Trace trace; + }; + return std::unique_ptr(new Adaptor{loc, title}); + }); // Summary @@ -317,6 +319,8 @@ void Library::initialise(const eckit::Parametrisation& config) { out << " communicator [" << mpi::comm() << "] \n"; out << " size [" << mpi::size() << "] \n"; out << " rank [" << mpi::rank() << "] \n"; + out << " OMP\n"; + out << " max_threads [" << atlas_omp_get_max_threads() << "] \n"; out << " \n"; out << " log.info [" << str(info_) << "] \n"; out << " log.trace [" << str(trace()) << "] \n"; @@ -336,6 +340,8 @@ void Library::initialise() { } void Library::finalise() { + atlas::io::TraceHookRegistry::disable(atlas_io_trace_hook_); + if (ATLAS_HAVE_TRACE && trace_report_) { Log::info() << atlas::Trace::report() << std::endl; } diff --git a/src/atlas/library/Library.h b/src/atlas/library/Library.h index 4d3516099..e911ee3a7 100644 --- a/src/atlas/library/Library.h +++ b/src/atlas/library/Library.h @@ -105,6 +105,7 @@ class Library : public eckit::system::Library { private: std::vector plugins_; std::vector data_paths_; + size_t atlas_io_trace_hook_; }; using Atlas = Library; diff --git a/src/atlas/mesh/actions/BuildHalo.cc b/src/atlas/mesh/actions/BuildHalo.cc index 8d3b6c769..c83195fa9 100644 --- a/src/atlas/mesh/actions/BuildHalo.cc +++ b/src/atlas/mesh/actions/BuildHalo.cc @@ -14,6 +14,8 @@ #include #include #include +#include +#include #include "atlas/array.h" #include "atlas/array/IndexView.h" @@ -28,6 +30,8 @@ #include "atlas/mesh/detail/AccumulateFacets.h" #include "atlas/parallel/mpi/Buffer.h" #include "atlas/parallel/mpi/mpi.h" +#include "atlas/parallel/omp/omp.h" +#include "atlas/parallel/omp/sort.h" #include "atlas/runtime/Exception.h" #include "atlas/runtime/Log.h" #include "atlas/runtime/Trace.h" @@ -428,7 +432,9 @@ class Notification { std::vector notes; }; -using Uid2Node = std::map; +using Uid2Node = std::unordered_map; + + void build_lookup_uid2node(Mesh& mesh, Uid2Node& uid2node) { ATLAS_TRACE(); Notification notes; @@ -444,21 +450,21 @@ void build_lookup_uid2node(Mesh& mesh, Uid2Node& uid2node) { uid_t uid = compute_uid(jnode); bool inserted = uid2node.insert(std::make_pair(uid, jnode)).second; if (not inserted) { - int other = uid2node[uid]; - std::stringstream msg; - msg << std::setprecision(10) << std::fixed << "Node uid: " << uid << " " << glb_idx(jnode) << " xy(" - << xy(jnode, XX) << "," << xy(jnode, YY) << ")"; - if (nodes.has_field("ij")) { - auto ij = array::make_view(nodes.field("ij")); - msg << " ij(" << ij(jnode, XX) << "," << ij(jnode, YY) << ")"; - } - msg << " has already been added as node " << glb_idx(other) << " (" << xy(other, XX) << "," << xy(other, YY) - << ")"; - if (nodes.has_field("ij")) { - auto ij = array::make_view(nodes.field("ij")); - msg << " ij(" << ij(other, XX) << "," << ij(other, YY) << ")"; - } - notes.add_error(msg.str()); + int other = uid2node[uid]; + std::stringstream msg; + msg << std::setprecision(10) << std::fixed << "Node uid: " << uid << " " << glb_idx(jnode) << " xy(" + << xy(jnode, XX) << "," << xy(jnode, YY) << ")"; + if (nodes.has_field("ij")) { + auto ij = array::make_view(nodes.field("ij")); + msg << " ij(" << ij(jnode, XX) << "," << ij(jnode, YY) << ")"; + } + msg << " has already been added as node " << glb_idx(other) << " (" << xy(other, XX) << "," << xy(other, YY) + << ")"; + if (nodes.has_field("ij")) { + auto ij = array::make_view(nodes.field("ij")); + msg << " ij(" << ij(other, XX) << "," << ij(other, YY) << ")"; + } + notes.add_error(msg.str()); } } if (notes.error()) { @@ -466,10 +472,12 @@ void build_lookup_uid2node(Mesh& mesh, Uid2Node& uid2node) { } } +/// For given nodes, find new connected elements void accumulate_elements(const Mesh& mesh, const mpi::BufferView& request_node_uid, const Uid2Node& uid2node, const Node2Elem& node2elem, std::vector& found_elements, std::set& new_nodes_uid) { - // ATLAS_TRACE(); + + ATLAS_TRACE(); const mesh::HybridElements::Connectivity& elem_nodes = mesh.cells().node_connectivity(); const auto elem_part = array::make_view(mesh.cells().partition()); @@ -477,14 +485,15 @@ void accumulate_elements(const Mesh& mesh, const mpi::BufferView& request const idx_t nb_request_nodes = static_cast(request_node_uid.size()); const int mpi_rank = static_cast(mpi::rank()); - std::set found_elements_set; + std::unordered_set found_elements_set; + found_elements_set.reserve(nb_request_nodes*2); for (idx_t jnode = 0; jnode < nb_request_nodes; ++jnode) { uid_t uid = request_node_uid(jnode); idx_t inode = -1; // search and get node index for uid - Uid2Node::const_iterator found = uid2node.find(uid); + auto found = uid2node.find(uid); if (found != uid2node.end()) { inode = found->second; } @@ -682,7 +691,7 @@ class BuildHaloHelper { template void fill_sendbuffer(Buffers& buf, const NodeContainer& nodes_uid, const ElementContainer& elems, const int p) { - // ATLAS_TRACE(); + ATLAS_TRACE(); idx_t nb_nodes = static_cast(nodes_uid.size()); buf.node_glb_idx[p].resize(nb_nodes); @@ -692,11 +701,11 @@ class BuildHaloHelper { buf.node_xy[p].resize(2 * nb_nodes); idx_t jnode = 0; - typename NodeContainer::iterator it; + typename NodeContainer::const_iterator it; for (it = nodes_uid.begin(); it != nodes_uid.end(); ++it, ++jnode) { uid_t uid = *it; - Uid2Node::iterator found = uid2node.find(uid); + auto found = uid2node.find(uid); if (found != uid2node.end()) // Point exists inside domain { idx_t node = found->second; @@ -758,11 +767,11 @@ class BuildHaloHelper { buf.node_xy[p].resize(2 * nb_nodes); int jnode = 0; - typename NodeContainer::iterator it; + typename NodeContainer::const_iterator it; for (it = nodes_uid.begin(); it != nodes_uid.end(); ++it, ++jnode) { uid_t uid = *it; - Uid2Node::iterator found = uid2node.find(uid); + auto found = uid2node.find(uid); if (found != uid2node.end()) // Point exists inside domain { int node = found->second; @@ -908,7 +917,7 @@ class BuildHaloHelper { // make sure new node was not already there { uid_t uid = compute_uid(loc_idx); - Uid2Node::iterator found = uid2node.find(uid); + auto found = uid2node.find(uid); if (found != uid2node.end()) { int other = found->second; std::stringstream msg; @@ -938,11 +947,11 @@ class BuildHaloHelper { std::set new_elem_uid; { ATLAS_TRACE("compute elem_uid"); - for (int jelem = 0; jelem < nb_elems; ++jelem) { + atlas_omp_parallel_for (int jelem = 0; jelem < nb_elems; ++jelem) { elem_uid[jelem * 2 + 0] = -compute_uid(elem_nodes->row(jelem)); elem_uid[jelem * 2 + 1] = cell_gidx(jelem); } - std::sort(elem_uid.begin(), elem_uid.end()); + omp::sort(elem_uid.begin(), elem_uid.end()); } auto element_already_exists = [&elem_uid, &new_elem_uid](uid_t uid) -> bool { std::vector::iterator it = std::lower_bound(elem_uid.begin(), elem_uid.end(), uid); @@ -1202,15 +1211,17 @@ void increase_halo_interior(BuildHaloHelper& helper) { gather_bdry_nodes(helper, send_bdry_nodes_uid, recv_bdry_nodes_uid_from_parts); + { + runtime::trace::Barriers set_barriers(false); + runtime::trace::Logging set_logging(false); #ifndef ATLAS_103 /* deprecated */ - for (idx_t jpart = 0; jpart < mpi_size; ++jpart) + atlas_omp_parallel_for (idx_t jpart = 0; jpart < mpi_size; ++jpart) #else const Mesh::PartitionGraph::Neighbours neighbours = helper.mesh.nearestNeighbourPartitions(); for (idx_t jpart : neighbours) #endif { - // 3) Find elements and nodes completing these elements in // other tasks that have my nodes through its UID @@ -1225,6 +1236,7 @@ void increase_halo_interior(BuildHaloHelper& helper) { // 4) Fill node and element buffers to send back helper.fill_sendbuffer(sendmesh, found_bdry_nodes_uid, found_bdry_elems, jpart); } + } // 5) Now communicate all buffers helper.all_to_all(sendmesh, recvmesh); @@ -1312,9 +1324,12 @@ void increase_halo_periodic(BuildHaloHelper& helper, const PeriodicPoints& perio gather_bdry_nodes(helper, send_bdry_nodes_uid, recv_bdry_nodes_uid_from_parts, /* periodic = */ true); + { + runtime::trace::Barriers set_barriers(false); + runtime::trace::Logging set_logging(false); #ifndef ATLAS_103 /* deprecated */ - for (idx_t jpart = 0; jpart < mpi_size; ++jpart) + atlas_omp_parallel_for (idx_t jpart = 0; jpart < mpi_size; ++jpart) #else Mesh::PartitionGraph::Neighbours neighbours = helper.mesh.nearestNeighbourPartitions(); // add own rank to neighbours to allow periodicity with self (pole caps) @@ -1323,6 +1338,7 @@ void increase_halo_periodic(BuildHaloHelper& helper, const PeriodicPoints& perio for (idx_t jpart : neighbours) #endif { + // 3) Find elements and nodes completing these elements in // other tasks that have my nodes through its UID @@ -1337,6 +1353,7 @@ void increase_halo_periodic(BuildHaloHelper& helper, const PeriodicPoints& perio // 4) Fill node and element buffers to send back helper.fill_sendbuffer(sendmesh, found_bdry_nodes_uid, found_bdry_elems, transform, newflags, jpart); } + } // 5) Now communicate all buffers helper.all_to_all(sendmesh, recvmesh); diff --git a/src/atlas/meshgenerator/detail/HealpixMeshGenerator.cc b/src/atlas/meshgenerator/detail/HealpixMeshGenerator.cc index c95c0467e..1cbe8bda7 100644 --- a/src/atlas/meshgenerator/detail/HealpixMeshGenerator.cc +++ b/src/atlas/meshgenerator/detail/HealpixMeshGenerator.cc @@ -64,6 +64,15 @@ HealpixMeshGenerator::HealpixMeshGenerator(const eckit::Parametrisation& p) { options.set("3d", three_dimensional); } + std::string pole_elements{"quqds"}; + if (p.get("pole_elements", pole_elements)) { + if (pole_elements != "pentagons" and pole_elements != "quads") { + Log::warning() << "Atlas::HealpixMeshGenerator accepts \"pentagons\" or \"quads\" for \"pole_elements\"." + << "Defaulting to pole_elements = quads" << std::endl; + } + } + options.set("pole_elements", pole_elements); + std::string partitioner; if (p.get("partitioner", partitioner)) { if (not grid::Partitioner::exists(partitioner)) { @@ -82,9 +91,13 @@ void HealpixMeshGenerator::configure_defaults() { // This option sets the part that will be generated options.set("part", mpi::rank()); - // This option switches between original HEALPix (=1) or HEALPix with 8 points are the pole (=8) + // This option switches between original HEALPix with 1 point at the pole (3d -> true) + // or HEALPix with 8 or 4 points at the pole (3d -> false) options.set("3d", false); + // This options switches between pentagons and quads as the pole elements for (3d -> false) + options.set("pole_elements", "quads"); + // This options sets the default partitioner std::string partitioner; if (grid::Partitioner::exists("equal_regions") && mpi::size() > 1) { @@ -97,27 +110,32 @@ void HealpixMeshGenerator::configure_defaults() { } -// match glb_idx of node in (nb_pole_nodes==8)-mesh to glb_idx of nodes in (nb_pole_nodes==1)-mesh -gidx_t HealpixMeshGenerator::match_idx(gidx_t gidx, const int ns) const { +// match glb_idx of node in (nb_pole_nodes==8 and ==4)-meshes to glb_idx of nodes in (nb_pole_nodes==1)-mesh +gidx_t HealpixMeshGenerator::match_node_idx(const gidx_t& gidx, const int ns) const { const gidx_t nb_nodes_orig = 12 * ns * ns; - if (nb_pole_nodes_ == 1 or gidx >= nb_nodes_) { + if (gidx > nb_nodes_ - 1) { + // no change in index for the periodic nodes return gidx; } - if (gidx == 4) { - return 0; - } - if (gidx == nb_nodes_ - 4) { - return nb_nodes_orig + 1; - } - bool at_north_pole = (gidx < 8); - bool at_south_pole = (gidx <= nb_nodes_ and gidx >= nb_nodes_orig + nb_pole_nodes_); - if (at_north_pole) { - return nb_nodes_orig + 2 + gidx - (gidx > 4 ? 1 : 0); - } - if (at_south_pole) { - return gidx - nb_pole_nodes_ + 9 - (gidx > nb_nodes_orig + 12 ? 1 : 0); + if (nb_pole_nodes_ > 1) { + if (gidx == nb_pole_nodes_ / 2) { + return 0; + } + if (gidx == nb_nodes_ - nb_pole_nodes_ / 2) { + return nb_nodes_orig + 1; + } + bool at_north_pole = (gidx < nb_pole_nodes_); + bool at_south_pole = (gidx < nb_nodes_ and gidx >= nb_nodes_orig + nb_pole_nodes_); + if (at_north_pole) { + return nb_nodes_orig + 2 + gidx - (gidx > nb_pole_nodes_ / 2 ? 1 : 0); + } + if (at_south_pole) { + return gidx - nb_pole_nodes_ + nb_pole_nodes_ + 1 - (gidx > nb_nodes_orig + nb_pole_nodes_ * 3 / 2 ? 1 : 0); + } + return gidx - nb_pole_nodes_ + 1; } - return gidx - nb_pole_nodes_ + 1; + // no change for 3d healpix mesh with one node per pole, i.e. nb_pole_nodes = 1 + return gidx; } @@ -181,7 +199,7 @@ gidx_t HealpixMeshGenerator::up_idx(const int xidx, const int yidx, const int ns } else if (yidx == 1) { ATLAS_ASSERT(xidx < 4); - ret = (nb_pole_nodes_ == 8 ? 2 * xidx : 0); + ret = (nb_pole_nodes_ == 8 ? 2 * xidx : (nb_pole_nodes_ == 4 ? xidx : 0)); } else if (yidx < ns) { ATLAS_ASSERT(xidx < 4 * yidx); @@ -234,10 +252,12 @@ gidx_t HealpixMeshGenerator::up_idx(const int xidx, const int yidx, const int ns ATLAS_ASSERT(xidx < nb_pole_nodes_); if (ns == 1) { if (xidx != nb_pole_nodes_ - 1) { - ret = nb_nodes_orig + nb_pole_nodes_ - 4 + (xidx % 2 ? -4 + (xidx + 1) / 2 : xidx / 2); + ret = nb_nodes_orig + nb_pole_nodes_ - 4 + + (nb_pole_nodes_ != 4 ? (xidx % 2 ? -4 + (xidx + 1) / 2 : xidx / 2) : xidx); } else { - ret = (nb_pole_nodes_ == 8 ? ghostIdx(4 * ns - 2) : ghostIdx(4 * ns)); + ret = (nb_pole_nodes_ == 4 ? nb_nodes_orig + xidx + : (nb_pole_nodes_ == 8 ? ghostIdx(4 * ns - 2) : ghostIdx(4 * ns))); } } else { @@ -285,7 +305,7 @@ gidx_t HealpixMeshGenerator::down_idx(const int xidx, const int yidx, const int } else if (yidx == 2 * ns && ns == 1) { ATLAS_ASSERT(xidx < 4); - ret = (nb_pole_nodes_ == 8 ? 16 + xidx : 9 + xidx); + ret = (nb_pole_nodes_ == 8 ? 16 + xidx : (nb_pole_nodes_ == 4 ? 12 + xidx : 9 + xidx)); } else if (yidx < 3 * ns && ns > 1) { ATLAS_ASSERT(xidx < 4 * ns); @@ -303,6 +323,9 @@ gidx_t HealpixMeshGenerator::down_idx(const int xidx, const int yidx, const int if (nb_pole_nodes_ == 8) { ret = (xidx != 7 ? nb_nodes_orig + 4 + (xidx + 1) / 2 : ghostIdx(4 * ns - 1)); } + else if (nb_pole_nodes_ == 4) { + ret = (xidx != 7 ? nb_nodes_orig + (xidx + 1) / 2 : ghostIdx(4 * ns - 1)); + } else { ret = (xidx != 7 ? nb_nodes_orig - 3 + (xidx + 1) / 2 : ghostIdx(4 * ns - 1)); } @@ -312,6 +335,9 @@ gidx_t HealpixMeshGenerator::down_idx(const int xidx, const int yidx, const int if (nb_pole_nodes_ == 8) { ret = nb_nodes_orig + nb_pole_nodes_ + 2 * xidx; } + else if (nb_pole_nodes_ == 4) { + ret = nb_nodes_orig + nb_pole_nodes_ + xidx; + } else { ret = nb_nodes_orig + 1; } @@ -362,6 +388,9 @@ gidx_t HealpixMeshGenerator::right_idx(const int xidx, const int yidx, const int if (nb_pole_nodes_ == 8) { ret = 1 + 2 * xidx; } + else if (nb_pole_nodes_ == 4) { + ret = (xidx != 3 ? xidx + 1 : ghostIdx(0)); + } else { ret = (xidx != 3 ? xidx + 2 : ghostIdx(1)); } @@ -372,7 +401,15 @@ gidx_t HealpixMeshGenerator::right_idx(const int xidx, const int yidx, const int } else if (yidx == 3 && ns == 1) { ATLAS_ASSERT(xidx < 4); - ret = (nb_pole_nodes_ == 8 ? 21 + 2 * xidx : (xidx != 3 ? 10 + xidx : ghostIdx(yidx))); + if (nb_pole_nodes_ == 8) { + ret = 21 + 2 * xidx; + } + else if (nb_pole_nodes_ == 4) { + ret = (xidx != 3 ? 17 + xidx : ghostIdx(yidx + 1)); + } + else { + ret = (xidx != 3 ? 10 + xidx : ghostIdx(yidx)); + } } else if (yidx <= 3 * ns) { ATLAS_ASSERT(xidx < 4 * ns + 1); @@ -395,6 +432,9 @@ gidx_t HealpixMeshGenerator::right_idx(const int xidx, const int yidx, const int if (nb_pole_nodes_ == 8) { ret = nb_nodes_orig + nb_pole_nodes_ + 1 + 2 * xidx; } + else if (nb_pole_nodes_ == 4) { + ret = (xidx != 3 ? nb_nodes_orig + nb_pole_nodes_ + 1 + xidx : ghostIdx(yidx + 1)); + } else { ret = (xidx != 3 ? nb_nodes_orig - 2 + xidx : ghostIdx(yidx)); } @@ -411,6 +451,21 @@ gidx_t HealpixMeshGenerator::right_idx(const int xidx, const int yidx, const int return ret; } +// return global_id - 1 of the pentagon node "to the right of" (xidx,yidx) node +// pentagon points are only needed for yidx == 1 and yidx == 4 * ns - 1 +gidx_t HealpixMeshGenerator::pentagon_right_idx(const int xidx, const int yidx, const int ns) const { + auto ghostIdx = [ns, this](int latid) { return this->nb_nodes_ + latid; }; + if (yidx == 1) { + return (xidx != 3 ? nb_pole_nodes_ + xidx + 1 : ghostIdx(1)); + } + else if (yidx == 4 * ns - 1) { + return (xidx != 3 ? (12 * ns * ns + xidx + 1) : ghostIdx(yidx)); + } + else { + return -2; + } +} + void HealpixMeshGenerator::generate(const Grid& grid, Mesh& mesh) const { ATLAS_ASSERT(HealpixGrid(grid), "Grid could not be cast to a HealpixGrid"); ATLAS_ASSERT(!mesh.generated()); @@ -471,6 +526,7 @@ void HealpixMeshGenerator::generate_mesh(const StructuredGrid& grid, const grid: // array::ArrayView flags ( nodes.flags() ); // - define cells (only quadrilaterals for now) with // mesh.cells().add( new mesh::temporary::Quadrilateral(), nquads ); + // mesh.cells().add( new mesh::temporary::Pentagon(), npents ); // further define cells with // array::ArrayView cells_glb_idx( mesh.cells().global_index() // ); @@ -478,8 +534,8 @@ void HealpixMeshGenerator::generate_mesh(const StructuredGrid& grid, const grid: // - define connectivity with // mesh::HybridElements::Connectivity& node_connectivity = // mesh.cells().node_connectivity(); - // node_connectivity.set( jcell, quad_nodes ); - // where quad_nodes is a 4-element integer array containing the LOCAL + // node_connectivity.set( jcell, cell_nodes ); + // where cell_nodes is a 4-element or 5-element integer array containing the LOCAL // indices of the nodes // // The rule do determine if a cell belongs to a proc is the following with some @@ -490,13 +546,14 @@ void HealpixMeshGenerator::generate_mesh(const StructuredGrid& grid, const grid: ATLAS_ASSERT(HealpixGrid(grid)); - const int mypart = options.get("part"); - const int nparts = options.get("nb_parts"); - const bool three_dimensional = options.get("3d"); - const int nb_pole_nodes = three_dimensional ? 1 : 8; - const int ny = grid.ny() + 2; - const int ns = (ny - 1) / 4; - const int nvertices = 12 * ns * ns + 2 * nb_pole_nodes; + const int mypart = options.get("part"); + const int nparts = options.get("nb_parts"); + const bool three_dimensional = options.get("3d"); + const std::string pole_elements = options.get("pole_elements"); + const int nb_pole_nodes = (pole_elements == "pentagons") ? 4 : (three_dimensional ? 1 : 8); + const int ny = grid.ny() + 2; + const int ns = (ny - 1) / 4; + const int nvertices = 12 * ns * ns + 2 * nb_pole_nodes; nb_pole_nodes_ = nb_pole_nodes; @@ -516,6 +573,7 @@ void HealpixMeshGenerator::generate_mesh(const StructuredGrid& grid, const grid: std::vector local_idx(nvertices, -1); std::vector current_idx(nparts, 0); // index counter for each proc + // ANSATZ: requirement on the partitioner auto compute_part = [&](int iy, gidx_t ii_glb) -> int { // nodes at the pole belong to proc_0 (north) and proc_maxRank (south) // a node gets its proc rank from the element for which this node would be its west vertex @@ -526,9 +584,10 @@ void HealpixMeshGenerator::generate_mesh(const StructuredGrid& grid, const grid: for (iy = 0; iy < ny; iy++) { int nx = nb_lat_nodes(iy); for (ix = 0; ix < nx; ix++) { - Log::info() << "iy, ix, glb_idx, up_idx, down_idx, right_idx : " << iy << ", " << ix << ", " + Log::info() << "iy, ix, glb_idx, up_idx, down_idx, right_idx, pent_right_idx : " << iy << ", " << ix << ", " << idx_xy_to_x(ix, iy, ns) + 1 << ", " << up_idx(ix, iy, ns) + 1 << ", " - << down_idx(ix, iy, ns) + 1 << ", " << right_idx(ix, iy, ns) + 1 << std::endl; + << down_idx(ix, iy, ns) + 1 << ", " << right_idx(ix, iy, ns) + 1 << ", " + << pentagon_right_idx(ix, iy, ns) + 1 << std::endl; } Log::info() << std::endl; } @@ -652,65 +711,74 @@ void HealpixMeshGenerator::generate_mesh(const StructuredGrid& grid, const grid: // determine number of cells and number of nodes int nnodes = 0; - int ncells = 0; + int nquads = 0; + int npents = 0; ii = 0; - int iil; + + int glb2loc_ghost_offset = -nnodes_SB + iy_max + nb_nodes_ + 1; + auto get_local_id = [this, &parts_sidx, &glb2loc_ghost_offset](gidx_t gidx) { + return gidx - (gidx < nb_nodes_ ? parts_sidx : glb2loc_ghost_offset); + }; + for (iy = iy_min; iy <= iy_max; iy++) { int nx = nb_lat_nodes(iy); for (ix = 0; ix < nx; ix++) { - int not_duplicate_cell = (iy == 0 ? ix % 2 : 1) * (iy == ny - 1 ? ix % 2 : 1); - if ((iy == 0 or iy == ny - 1) and nb_pole_nodes == 1) { - not_duplicate_cell = false; - } - - if (not is_ghost_SB[ii] && not_duplicate_cell) { - // mark this node as being used - if (not is_node_SB[ii]) { - ++nnodes; - is_node_SB[ii] = true; - } - - ++ncells; - - const int glb2loc_ghost_offset = -nnodes_SB + iy_max + nb_nodes_ + 1; - // mark upper corner - iil = up_idx(ix, iy, ns); - iil -= (iil < nb_nodes_ ? parts_sidx : glb2loc_ghost_offset); - if (not is_node_SB[iil]) { - ++nnodes; - is_node_SB[iil] = true; + if (not is_ghost_SB[ii]) { + const bool at_pole = (iy == 0 or iy == ny - 1); + bool not_duplicate_cell = (at_pole ? (ix % 2) : 1); + if (at_pole and nb_pole_nodes < 8) { + // nodes at the poles do not own any pentagons + // if nb_pole_node=1, the node at the poles does not own any quads + not_duplicate_cell = false; } - // mark lower corner - iil = down_idx(ix, iy, ns); - iil -= (iil < nb_nodes_ ? parts_sidx : glb2loc_ghost_offset); - if (not is_node_SB[iil]) { - ++nnodes; - is_node_SB[iil] = true; - } - // mark right corner - iil = right_idx(ix, iy, ns); - iil -= (iil < nb_nodes_ ? parts_sidx : glb2loc_ghost_offset); - if (not is_node_SB[iil]) { - ++nnodes; - is_node_SB[iil] = true; + if (not at_pole or not_duplicate_cell) { + if (not is_node_SB[ii]) { + ++nnodes; + is_node_SB[ii] = true; + } + bool is_pentagon = (pole_elements == "pentagons" and (iy == 1 or iy == ny - 2)); + if (is_pentagon) { + ++npents; + // mark the pentagon right node + idx_t iil = get_local_id(pentagon_right_idx(ix, iy, ns)); + if (not is_node_SB[iil]) { + ++nnodes; + is_node_SB[iil] = true; + } + } + else { + ++nquads; + } + // mark upper corner + idx_t iil = get_local_id(up_idx(ix, iy, ns)); + if (not is_node_SB[iil]) { + ++nnodes; + is_node_SB[iil] = true; + } + // mark lower corner + iil = get_local_id(down_idx(ix, iy, ns)); + if (not is_node_SB[iil]) { + ++nnodes; + is_node_SB[iil] = true; + } + // mark right corner + iil = get_local_id(right_idx(ix, iy, ns)); + if (not is_node_SB[iil]) { + ++nnodes; + is_node_SB[iil] = true; + } } } ++ii; } } - - // periodic points are always needed, even if they don't belong to a cell - ii_ghost = nnodes_SB - (iy_max - iy_min + 1); - for (ii = 0; ii < iy_max - iy_min + 1; ii++) { - if (!is_node_SB[ii_ghost + ii]) { - // is_node_SB[ii_ghost + ii] = true; - // ++nnodes; - } - } + int ncells = nquads + npents; + ATLAS_ASSERT(ncells > 0); #if DEBUG_OUTPUT Log::info() << "[" << mypart << "] : " - << "nnodes = " << nnodes << ", ncells = " << ncells << ", parts_sidx = " << parts_sidx << std::endl; + << "nnodes = " << nnodes << ", nquads = " << nquads << ", npents = " << npents + << ", parts_sidx = " << parts_sidx << std::endl; Log::info() << "[" << mypart << "] : " << "iy_min = " << iy_min << ", iy_max = " << iy_max << std::endl; #endif @@ -736,14 +804,17 @@ void HealpixMeshGenerator::generate_mesh(const StructuredGrid& grid, const grid: auto flags = array::make_view(nodes.flags()); // define cells and associated properties - mesh.cells().add(new mesh::temporary::Quadrilateral(), ncells); + mesh.cells().add(new mesh::temporary::Quadrilateral(), nquads); + mesh.cells().add(new mesh::temporary::Pentagon(), npents); int quad_begin = mesh.cells().elements(0).begin(); + int pent_begin = mesh.cells().elements(1).begin(); auto cells_part = array::make_view(mesh.cells().partition()); auto cells_glb_idx = array::make_view(mesh.cells().global_index()); auto& node_connectivity = mesh.cells().node_connectivity(); - idx_t quad_nodes[4]; - int jcell = quad_begin; + idx_t cell_nodes[5]; + int jquadcell = quad_begin; + int jpentcell = pent_begin; int inode_nonghost, inode_ghost; @@ -768,18 +839,18 @@ void HealpixMeshGenerator::generate_mesh(const StructuredGrid& grid, const grid: // flags Topology::reset(flags(inode)); - glb_idx(inode) = 1 + match_idx(idx_xy_to_x(ix, iy, ns), ns); + glb_idx(inode) = 1 + match_node_idx(idx_xy_to_x(ix, iy, ns), ns); // grid coordinates double _xy[2]; double xy1[2], xy2[2]; if (iy == 0) { - _xy[0] = (nb_pole_nodes == 8 ? 45. * ix : 180.); + _xy[0] = (nb_pole_nodes == 8 ? 45. * ix : (nb_pole_nodes == 4 ? 90. * ix : 180.)); _xy[1] = 90.; Topology::set(flags(inode), Topology::BC); } else if (iy == ny - 1) { - _xy[0] = (nb_pole_nodes == 8 ? 45. * ix : 180.); + _xy[0] = (nb_pole_nodes == 8 ? 45. * ix : (nb_pole_nodes == 4 ? 90. * ix : 180.)); _xy[1] = -90.; Topology::set(flags(inode), Topology::BC); } @@ -855,59 +926,92 @@ void HealpixMeshGenerator::generate_mesh(const StructuredGrid& grid, const grid: } } + auto get_local_idx_SB = [&local_idx_SB, &get_local_id](gidx_t gidx) { return local_idx_SB[get_local_id(gidx)]; }; + ii = 0; // index inside SB (surrounding belt) int jcell_offset = 0; // global index offset due to extra points at the north pole + gidx_t jcell = 0; // global cell counter for (iy = iy_min; iy <= iy_max; iy++) { int nx = nb_lat_nodes(iy) + 1; for (ix = 0; ix < nx; ix++) { - int not_duplicate_cell = (iy == 0 ? ix % 2 : 1) * (iy == ny - 1 ? ix % 2 : 1); - int iil = idx_xy_to_x(ix, iy, ns); + const bool at_pole = (iy == 0 or iy == ny - 1); + int not_duplicate_cell = (at_pole ? ix % 2 : 1); + if (at_pole and (not not_duplicate_cell or nb_pole_nodes < 8)) { + // if nb_pole_nodes = 4 : the pole nodes do not own any pentagons + // if nb_pole_nodes = 1 : the pole nodes do not own any quads + continue; + } + bool pentagon = (iy == 1 or iy == 4 * ns - 1) and pole_elements == "pentagons"; + int iil = idx_xy_to_x(ix, iy, ns); iil -= (iil < nb_nodes_ ? parts_sidx : -nnodes_SB + iy_max + nb_nodes_ + 1); - if (!is_ghost_SB[iil] && not_duplicate_cell) { - // define cell corners (local indices) - quad_nodes[0] = local_idx_SB[iil]; - - quad_nodes[1] = down_idx(ix, iy, ns); // point to the right, global idx - quad_nodes[1] -= (quad_nodes[1] < nb_nodes_ ? parts_sidx : -nnodes_SB + iy_max + nb_nodes_ + 1); - quad_nodes[1] = local_idx_SB[quad_nodes[1]]; // now, local idx - - quad_nodes[2] = right_idx(ix, iy, ns); // point above right, global idx - quad_nodes[2] -= (quad_nodes[2] < nb_nodes_ ? parts_sidx : -nnodes_SB + iy_max + nb_nodes_ + 1); - quad_nodes[2] = local_idx_SB[quad_nodes[2]]; // now, local idx - - quad_nodes[3] = up_idx(ix, iy, ns); // point above, global idx - quad_nodes[3] -= (quad_nodes[3] < nb_nodes_ ? parts_sidx : -nnodes_SB + iy_max + nb_nodes_ + 1); - quad_nodes[3] = local_idx_SB[quad_nodes[3]]; // now, local idx + if (not is_ghost_SB[iil]) { + jcell++; // a cell will be added + // define cell vertices (in local indices) in cell_nodes + int node_id = 0; + cell_nodes[node_id++] = get_local_idx_SB(idx_xy_to_x(ix, iy, ns)); + cell_nodes[node_id++] = get_local_idx_SB(down_idx(ix, iy, ns)); + bool south_hemisphere = (iy > 2 * ns); + if (pentagon && not south_hemisphere) { + cell_nodes[node_id++] = get_local_idx_SB(pentagon_right_idx(ix, iy, ns)); + } + cell_nodes[node_id++] = get_local_idx_SB(right_idx(ix, iy, ns)); + if (pentagon && south_hemisphere) { + // in the south hemisphere the pentagon point comes in a different clock-wise ordering + cell_nodes[node_id++] = get_local_idx_SB(pentagon_right_idx(ix, iy, ns)); + } + cell_nodes[node_id++] = get_local_idx_SB(up_idx(ix, iy, ns)); - node_connectivity.set(jcell, quad_nodes); - if (iy == 0) { - if (nb_pole_nodes == 8) { - cells_glb_idx(jcell) = 12 * ns * ns + 1 + ix / 2; + // match global cell indexing for the three healpix versions + if (nb_pole_nodes == 1) { + cells_glb_idx(jquadcell) = jquadcell + 1; + } + else if (nb_pole_nodes == 8) { + if (iy == 0) { + cells_glb_idx(jquadcell) = 12 * ns * ns + 1 + ix / 2; jcell_offset++; } - } - else if (iy == ny - 1) { - if (nb_pole_nodes == 8) { - cells_glb_idx(jcell) = 12 * ns * ns + 5 + ix / 2; + else if (iy == ny - 1) { + cells_glb_idx(jquadcell) = 12 * ns * ns + 5 + ix / 2; jcell_offset++; } + else { + cells_glb_idx(jquadcell) = parts_sidx + iil - 3 - (mypart != 0 ? 4 : jcell_offset); + } } - else { - if (nb_pole_nodes == 8) { - cells_glb_idx(jcell) = parts_sidx + iil - 3 - (mypart != 0 ? 4 : jcell_offset); + else if (nb_pole_nodes == 4) { + if (pentagon) { + cells_glb_idx(jpentcell) = jcell; } else { - cells_glb_idx(jcell) = parts_sidx + iil; + cells_glb_idx(jquadcell) = jcell; } } - cells_part(jcell) = mypart; #if DEBUG_OUTPUT_DETAIL - Log::info() << "[" << mypart << "] : " - << "New quad: loc-idx " << jcell << ", glb-idx " << cells_glb_idx(jcell) << ": " - << glb_idx(quad_nodes[0]) << "," << glb_idx(quad_nodes[1]) << "," << glb_idx(quad_nodes[2]) - << "," << glb_idx(quad_nodes[3]) << std::endl; + std::cout << "[" << mypart << "] : "; + if (pentagon) { + std::cout << "New pent: loc-idx " << jpentcell << ", glb-idx " << cells_glb_idx(jpentcell) << ": "; + } + else { + std::cout << "New quad: loc-idx " << jquadcell << ", glb-idx " << cells_glb_idx(jquadcell) << ": "; + } + std::cout << glb_idx(cell_nodes[0]) << "," << glb_idx(cell_nodes[1]) << "," << glb_idx(cell_nodes[2]) + << "," << glb_idx(cell_nodes[3]); + if (pentagon) { + std::cout << "," << glb_idx(cell_nodes[4]); + } + std::cout << std::endl; #endif - ++jcell; + // add cell to the node connectivity table + if (pentagon) { + cells_part(jpentcell) = mypart; + node_connectivity.set(jpentcell, cell_nodes); + ++jpentcell; + } + else { + cells_part(jquadcell) = mypart; + node_connectivity.set(jquadcell, cell_nodes); + ++jquadcell; + } } ii += (ix != nx - 1 ? 1 : 0); } @@ -915,19 +1019,25 @@ void HealpixMeshGenerator::generate_mesh(const StructuredGrid& grid, const grid: #if DEBUG_OUTPUT_DETAIL // list nodes + Log::info() << "Listing nodes ..."; for (inode = 0; inode < nnodes; inode++) { - Log::info() << "[" << mypart << "] : " - << " node " << inode << ": ghost = " << ghost(inode) << ", glb_idx = " << glb_idx(inode) - << ", part = " << part(inode) << ", lon = " << lonlat(inode, 0) << ", lat = " << lonlat(inode, 1) - << ", remote_idx = " << remote_idx(inode) << std::endl; - } - - int* cell_nodes; - for (jcell = 0; jcell < ncells; jcell++) { - Log::info() << "[" << mypart << "] : " - << " cell " << jcell << ", glb-idx " << cells_glb_idx(jcell) << ": " - << glb_idx(node_connectivity(jcell, 0)) << "," << glb_idx(node_connectivity(jcell, 1)) << "," - << glb_idx(node_connectivity(jcell, 2)) << "," << glb_idx(node_connectivity(jcell, 3)) << std::endl; + std::cout << "[" << mypart << "] : " + << " node " << inode << ": ghost = " << ghost(inode) << ", glb_idx = " << glb_idx(inode) + << ", part = " << part(inode) << ", lon = " << lonlat(inode, 0) << ", lat = " << lonlat(inode, 1) + << ", remote_idx = " << remote_idx(inode) << std::endl; + } + + for (gidx_t jcell = quad_begin; jcell < nquads; jcell++) { + std::cout << "[" << mypart << "] : " + << " cell " << jcell << ", glb-idx " << cells_glb_idx(jcell) << ": " + << glb_idx(node_connectivity(jcell, 0)) << "," << glb_idx(node_connectivity(jcell, 1)) << "," + << glb_idx(node_connectivity(jcell, 2)) << "," << glb_idx(node_connectivity(jcell, 3)) << std::endl; + } + for (gidx_t jcell = pent_begin; jcell < nquads; jcell++) { + std::cout << "[" << mypart << "] : " + << " cell " << jcell << ", glb-idx " << cells_glb_idx(jcell) << ": " + << glb_idx(node_connectivity(jcell, 0)) << "," << glb_idx(node_connectivity(jcell, 1)) << "," + << glb_idx(node_connectivity(jcell, 2)) << "," << glb_idx(node_connectivity(jcell, 3)) << std::endl; } #endif diff --git a/src/atlas/meshgenerator/detail/HealpixMeshGenerator.h b/src/atlas/meshgenerator/detail/HealpixMeshGenerator.h index a8feeab19..749ad885e 100644 --- a/src/atlas/meshgenerator/detail/HealpixMeshGenerator.h +++ b/src/atlas/meshgenerator/detail/HealpixMeshGenerator.h @@ -58,11 +58,12 @@ class HealpixMeshGenerator : public MeshGenerator::Implementation { void generate_mesh(const StructuredGrid&, const grid::Distribution&, Mesh& m) const; - gidx_t match_idx(gidx_t gidx, const int ns) const; + gidx_t match_node_idx(const gidx_t& gidx, const int ns) const; gidx_t idx_xy_to_x(const int xidx, const int yidx, const int ns) const; gidx_t up_idx(const int xidx, const int yidx, const int ns) const; gidx_t down_idx(const int xidx, const int yidx, const int ns) const; gidx_t right_idx(const int xidx, const int yidx, const int ns) const; + gidx_t pentagon_right_idx(const int xidx, const int yidx, const int ns) const; private: util::Metadata options; diff --git a/src/atlas/meshgenerator/detail/StructuredMeshGenerator.cc b/src/atlas/meshgenerator/detail/StructuredMeshGenerator.cc index 41ba410ce..c0fe3d4bf 100644 --- a/src/atlas/meshgenerator/detail/StructuredMeshGenerator.cc +++ b/src/atlas/meshgenerator/detail/StructuredMeshGenerator.cc @@ -15,6 +15,7 @@ #include #include +#include "eckit/log/ProgressTimer.h" #include "eckit/types/FloatCompare.h" #include "eckit/utils/Hash.h" @@ -271,39 +272,88 @@ void StructuredMeshGenerator::generate_region(const StructuredGrid& rg, const gr bool unique_pole = options.getBool("unique_pole") && three_dimensional && has_north_pole && has_south_pole; bool periodic_east_west = rg.periodic(); - int n; - /* + + idx_t lat_north = -1; + idx_t lat_south = -1; + + if (distribution.nb_partitions() == 1) { + lat_north = 0; + lat_south = rg.ny() - 1; + } + else { + ATLAS_TRACE_SCOPE("compute bounds") { + // This part does not scale + /* Find min and max latitudes used by this part. */ - n = 0; - idx_t lat_north = -1; - for (idx_t jlat = 0; jlat < rg.ny(); ++jlat) { - for (idx_t jlon = 0; jlon < rg.nx(jlat); ++jlon) { - if (distribution.partition(n) == mypart) { - lat_north = jlat; - goto end_north; + + size_t num_threads = atlas_omp_get_max_threads(); + if (num_threads == 1) { + int n = 0; + for (idx_t jlat = 0; jlat < rg.ny(); ++jlat) { + for (idx_t jlon = 0; jlon < rg.nx(jlat); ++jlon) { + if (distribution.partition(n) == mypart) { + lat_north = jlat; + goto end_north; + } + ++n; + } + } + end_north:; + + n = rg.size() - 1; + for (idx_t jlat = rg.ny() - 1; jlat >= 0; --jlat) { + for (idx_t jlon = rg.nx(jlat) - 1; jlon >= 0; --jlon) { + if (distribution.partition(n) == mypart) { + lat_south = jlat; + goto end_south; + } + --n; + } + } + end_south:; } - ++n; - } - } -end_north: + else { + std::vector thread_reduce_lat_north(num_threads, std::numeric_limits::max()); + std::vector thread_reduce_lat_south(num_threads, std::numeric_limits::min()); + + atlas_omp_parallel { + const idx_t thread_num = atlas_omp_get_thread_num(); + const idx_t jbegin = static_cast(size_t(thread_num) * size_t(rg.ny()) / size_t(num_threads)); + const idx_t jend = static_cast(size_t(thread_num + 1) * size_t(rg.ny()) / size_t(num_threads)); + + auto& thread_lat_north = thread_reduce_lat_north[thread_num]; + auto& thread_lat_south = thread_reduce_lat_south[thread_num]; + + for (idx_t jlat = jbegin; jlat < jend && thread_lat_north > jend; ++jlat) { + int n = rg.index(0, jlat); + for (idx_t jlon = 0; jlon < rg.nx(jlat); ++jlon, ++n) { + if (distribution.partition(n) == mypart) { + thread_lat_north = jlat; + break; + } + } + } - n = rg.size() - 1; - idx_t lat_south = -1; - for (idx_t jlat = rg.ny() - 1; jlat >= 0; --jlat) { - for (idx_t jlon = rg.nx(jlat) - 1; jlon >= 0; --jlon) { - if (distribution.partition(n) == mypart) { - lat_south = jlat; - goto end_south; + for (idx_t jlat = jend - 1; jlat >= jbegin && thread_lat_south < jbegin; --jlat) { + int n = rg.index(0, jlat); + for (idx_t jlon = 0; jlon < rg.nx(jlat); ++jlon, ++n) { + if (distribution.partition(n) == mypart) { + thread_lat_south = jlat; + break; + } + } + } + } + lat_north = *std::min_element(thread_reduce_lat_north.begin(), thread_reduce_lat_north.end()); + lat_south = *std::max_element(thread_reduce_lat_south.begin(), thread_reduce_lat_south.end()); } - --n; } } -end_south: std::vector offset(rg.ny(), 0); - n = 0; + int n = 0; for (idx_t jlat = 0; jlat < rg.ny(); ++jlat) { offset.at(jlat) = n; n += rg.nx(jlat); @@ -335,411 +385,419 @@ We need to connect to next region array::ArrayView elemview = array::make_view(*region.elems); elemview.assign(-1); - for (idx_t jlat = lat_north; jlat < lat_south; ++jlat) { - // std::stringstream filename; filename << "/tmp/debug/"< 1280 ? Log::trace() : blackhole); - ilat = jlat - region.north; + for (idx_t jlat = lat_north; jlat < lat_south; ++jlat, ++progress) { + // std::stringstream filename; filename << "/tmp/debug/"< dS1N2) { - if (ipN1 != ipN2) { - try_make_triangle_down = true; + else if (dN1S2 > dS1N2) { + if (ipN1 != ipN2) { + try_make_triangle_down = true; + } + else { + try_make_triangle_up = true; + } } else { - try_make_triangle_up = true; + throw_Exception("Should not be here", Here()); } } - else { - throw_Exception("Should not be here", Here()); - } - } - } - else { - if (ipN1 == ipN2) { - try_make_triangle_up = true; - } - else if (ipS1 == ipS2) { - try_make_triangle_down = true; } else { - try_make_quad = true; - } + if (ipN1 == ipN2) { + try_make_triangle_up = true; + } + else if (ipS1 == ipS2) { + try_make_triangle_down = true; + } + else { + try_make_quad = true; + } - // try_make_quad = true; - } - } - else { - dN1S2 = std::abs(xN1 - xS2); - dS1N2 = std::abs(xS1 - xN2); - // dN2S2 = std::abs(xN2-xS2); - // Log::info() << " dN1S2 " << dN1S2 << " dS1N2 " << dS1N2 << " - // dN2S2 " << dN2S2 << std::endl; - if ((dN1S2 <= dS1N2) && (ipS1 != ipS2)) { - try_make_triangle_up = true; - } - else if ((dN1S2 >= dS1N2) && (ipN1 != ipN2)) { - try_make_triangle_down = true; + // try_make_quad = true; + } } else { - if (ipN1 == ipN2) { + dN1S2 = std::abs(xN1 - xS2); + dS1N2 = std::abs(xS1 - xN2); + // dN2S2 = std::abs(xN2-xS2); + // Log::info() << " dN1S2 " << dN1S2 << " dS1N2 " << dS1N2 << " + // dN2S2 " << dN2S2 << std::endl; + if ((dN1S2 <= dS1N2) && (ipS1 != ipS2)) { try_make_triangle_up = true; } - else if (ipS1 == ipS2) { + else if ((dN1S2 >= dS1N2) && (ipN1 != ipN2)) { try_make_triangle_down = true; } else { - ATLAS_DEBUG_VAR(dN1S2); - ATLAS_DEBUG_VAR(dS1N2); - ATLAS_DEBUG_VAR(jlat); - Log::info() << ipN1 << "(" << xN1 << ") " << ipN2 << "(" << xN2 << ") " << std::endl; - Log::info() << ipS1 << "(" << xS1 << ") " << ipS2 << "(" << xS2 << ") " << std::endl; - throw_Exception("Should not try to make a quadrilateral!", Here()); + if (ipN1 == ipN2) { + try_make_triangle_up = true; + } + else if (ipS1 == ipS2) { + try_make_triangle_down = true; + } + else { + ATLAS_DEBUG_VAR(dN1S2); + ATLAS_DEBUG_VAR(dS1N2); + ATLAS_DEBUG_VAR(jlat); + Log::info() << ipN1 << "(" << xN1 << ") " << ipN2 << "(" << xN2 << ") " << std::endl; + Log::info() << ipS1 << "(" << xS1 << ") " << ipS2 << "(" << xS2 << ") " << std::endl; + throw_Exception("Should not try to make a quadrilateral!", Here()); + } } } - } - // ------------------------------------------------ - // END RULES - // ------------------------------------------------ + // ------------------------------------------------ + // END RULES + // ------------------------------------------------ #if DEBUG_OUTPUT - ATLAS_DEBUG_VAR(jelem); + ATLAS_DEBUG_VAR(jelem); #endif - auto elem = lat_elems_view.slice(jelem, Range::all()); + auto elem = lat_elems_view.slice(jelem, Range::all()); - if (try_make_quad) { + if (try_make_quad) { // add quadrilateral #if DEBUG_OUTPUT - Log::info() << " " << ipN1 << " " << ipN2 << '\n'; - Log::info() << " " << ipS1 << " " << ipS2 << '\n'; + Log::info() << " " << ipN1 << " " << ipN2 << '\n'; + Log::info() << " " << ipS1 << " " << ipS2 << '\n'; #endif - elem(0) = ipN1; - elem(1) = ipS1; - elem(2) = ipS2; - elem(3) = ipN2; - add_quad = false; - std::array np{pN1, pN2, pS1, pS2}; - std::array pcnts; - for (int j = 0; j < 4; ++j) { - pcnts[j] = static_cast(std::count(np.begin(), np.end(), np[j])); - } - if (pcnts[0] > 2) { // 3 or more of pN1 - pE = pN1; - if (latS == rg.ny() - 1) { - pE = pS1; - } - } - else if (pcnts[2] > 2) { // 3 or more of pS1 - pE = pS1; - if (latN == 0) { - pE = pN1; - } - } - else { - std::array::iterator p_max = std::max_element(pcnts.begin(), pcnts.end()); - if (*p_max > 2) { // 3 or 4 points belong to same part - pE = np[std::distance(np.begin(), p_max)]; + elem(0) = ipN1; + elem(1) = ipS1; + elem(2) = ipS2; + elem(3) = ipN2; + add_quad = false; + std::array np{pN1, pN2, pS1, pS2}; + std::array pcnts; + for (int j = 0; j < 4; ++j) { + pcnts[j] = static_cast(std::count(np.begin(), np.end(), np[j])); } - else { // 3 or 4 points don't belong to mypart + if (pcnts[0] > 2) { // 3 or more of pN1 pE = pN1; if (latS == rg.ny() - 1) { pE = pS1; } } - } - add_quad = (pE == mypart); - if (add_quad) { - ++region.nquads; - ++jelem; - ++nelems; - - if (region.lat_begin.at(latN) == -1) { - region.lat_begin.at(latN) = ipN1; + else if (pcnts[2] > 2) { // 3 or more of pS1 + pE = pS1; + if (latN == 0) { + pE = pN1; + } + } + else { + std::array::iterator p_max = std::max_element(pcnts.begin(), pcnts.end()); + if (*p_max > 2) { // 3 or 4 points belong to same part + pE = np[std::distance(np.begin(), p_max)]; + } + else { // 3 or 4 points don't belong to mypart + pE = pN1; + if (latS == rg.ny() - 1) { + pE = pS1; + } + } } - if (region.lat_begin.at(latS) == -1) { - region.lat_begin.at(latS) = ipS1; + add_quad = (pE == mypart); + if (add_quad) { + ++region.nquads; + ++jelem; + ++nelems; + + if (region.lat_begin.at(latN) == -1) { + region.lat_begin.at(latN) = ipN1; + } + if (region.lat_begin.at(latS) == -1) { + region.lat_begin.at(latS) = ipS1; + } + region.lat_begin.at(latN) = std::min(region.lat_begin.at(latN), ipN1); + region.lat_begin.at(latS) = std::min(region.lat_begin.at(latS), ipS1); + region.lat_end.at(latN) = std::max(region.lat_end.at(latN), ipN2); + region.lat_end.at(latS) = std::max(region.lat_end.at(latS), ipS2); } - region.lat_begin.at(latN) = std::min(region.lat_begin.at(latN), ipN1); - region.lat_begin.at(latS) = std::min(region.lat_begin.at(latS), ipS1); - region.lat_end.at(latN) = std::max(region.lat_end.at(latN), ipN2); - region.lat_end.at(latS) = std::max(region.lat_end.at(latS), ipS2); - } - else { + else { #if DEBUG_OUTPUT - Log::info() << "Quad belongs to other partition" << std::endl; + Log::info() << "Quad belongs to other partition" << std::endl; #endif + } + ipN1 = ipN2; + ipS1 = ipS2; } - ipN1 = ipN2; - ipS1 = ipS2; - } - else if (try_make_triangle_down) // make triangle down - { + else if (try_make_triangle_down) // make triangle down + { // triangle without ip3 #if DEBUG_OUTPUT - Log::info() << " " << ipN1 << " " << ipN2 << '\n'; - Log::info() << " " << ipS1 << '\n'; + Log::info() << " " << ipN1 << " " << ipN2 << '\n'; + Log::info() << " " << ipS1 << '\n'; #endif - elem(0) = ipN1; - elem(1) = ipS1; - elem(2) = -1; - elem(3) = ipN2; - - pE = pN1; - if (latS == rg.ny() - 1) { - pE = pS1; - } - add_triag = (mypart == pE); + elem(0) = ipN1; + elem(1) = ipS1; + elem(2) = -1; + elem(3) = ipN2; - if (add_triag) { - ++region.ntriags; - ++jelem; - ++nelems; - - if (region.lat_begin.at(latN) == -1) { - region.lat_begin.at(latN) = ipN1; + pE = pN1; + if (latS == rg.ny() - 1) { + pE = pS1; } - if (region.lat_begin.at(latS) == -1) { - region.lat_begin.at(latS) = ipS1; + add_triag = (mypart == pE); + + if (add_triag) { + ++region.ntriags; + ++jelem; + ++nelems; + + if (region.lat_begin.at(latN) == -1) { + region.lat_begin.at(latN) = ipN1; + } + if (region.lat_begin.at(latS) == -1) { + region.lat_begin.at(latS) = ipS1; + } + region.lat_begin.at(latN) = std::min(region.lat_begin.at(latN), ipN1); + region.lat_begin.at(latS) = std::min(region.lat_begin.at(latS), ipS1); + region.lat_end.at(latN) = std::max(region.lat_end.at(latN), ipN2); + region.lat_end.at(latS) = std::max(region.lat_end.at(latS), ipS1); } - region.lat_begin.at(latN) = std::min(region.lat_begin.at(latN), ipN1); - region.lat_begin.at(latS) = std::min(region.lat_begin.at(latS), ipS1); - region.lat_end.at(latN) = std::max(region.lat_end.at(latN), ipN2); - region.lat_end.at(latS) = std::max(region.lat_end.at(latS), ipS1); - } - else { + else { #if DEBUG_OUTPUT - Log::info() << "Downward Triag belongs to other partition" << std::endl; + Log::info() << "Downward Triag belongs to other partition" << std::endl; #endif + } + ipN1 = ipN2; + // and ipS1=ipS1; } - ipN1 = ipN2; - // and ipS1=ipS1; - } - else if (try_make_triangle_up) // make triangle up - { + else if (try_make_triangle_up) // make triangle up + { // triangle without ip4 #if DEBUG_OUTPUT - Log::info() << " " << ipN1 << " (" << pN1 << ")" << '\n'; - Log::info() << " " << ipS1 << " (" << pS1 << ")" - << " " << ipS2 << " (" << pS2 << ")" << '\n'; + Log::info() << " " << ipN1 << " (" << pN1 << ")" << '\n'; + Log::info() << " " << ipS1 << " (" << pS1 << ")" + << " " << ipS2 << " (" << pS2 << ")" << '\n'; #endif - elem(0) = ipN1; - elem(1) = ipS1; - elem(2) = ipS2; - elem(3) = -1; - - if (pS1 == pE && pN1 != pE) { - if (xN1 < 0.5 * (xS1 + xS2)) { + elem(0) = ipN1; + elem(1) = ipS1; + elem(2) = ipS2; + elem(3) = -1; + + if (pS1 == pE && pN1 != pE) { + if (xN1 < 0.5 * (xS1 + xS2)) { + pE = pN1; + } // else pE of previous element + } + else { pE = pN1; - } // else pE of previous element - } - else { - pE = pN1; - } - if (ipN1 == rg.nx(latN)) { - pE = pS1; - } - if (latS == rg.ny() - 1) { - pE = pS1; - } + } + if (ipN1 == rg.nx(latN)) { + pE = pS1; + } + if (latS == rg.ny() - 1) { + pE = pS1; + } - add_triag = (mypart == pE); + add_triag = (mypart == pE); - if (add_triag) { - ++region.ntriags; - ++jelem; - ++nelems; + if (add_triag) { + ++region.ntriags; + ++jelem; + ++nelems; - if (region.lat_begin.at(latN) == -1) { - region.lat_begin.at(latN) = ipN1; + if (region.lat_begin.at(latN) == -1) { + region.lat_begin.at(latN) = ipN1; + } + if (region.lat_begin.at(latS) == -1) { + region.lat_begin.at(latS) = ipS1; + } + region.lat_begin.at(latN) = std::min(region.lat_begin.at(latN), ipN1); + region.lat_begin.at(latS) = std::min(region.lat_begin.at(latS), ipS1); + region.lat_end.at(latN) = std::max(region.lat_end.at(latN), ipN1); + region.lat_end.at(latS) = std::max(region.lat_end.at(latS), ipS2); } - if (region.lat_begin.at(latS) == -1) { - region.lat_begin.at(latS) = ipS1; + else { +#if DEBUG_OUTPUT + Log::info() << "Upward Triag belongs to other partition" << std::endl; +#endif } - region.lat_begin.at(latN) = std::min(region.lat_begin.at(latN), ipN1); - region.lat_begin.at(latS) = std::min(region.lat_begin.at(latS), ipS1); - region.lat_end.at(latN) = std::max(region.lat_end.at(latN), ipN1); - region.lat_end.at(latS) = std::max(region.lat_end.at(latS), ipS2); + ipS1 = ipS2; + // and ipN1=ipN1; } else { -#if DEBUG_OUTPUT - Log::info() << "Upward Triag belongs to other partition" << std::endl; -#endif + throw_Exception("Could not detect which element to create", Here()); } - ipS1 = ipS2; - // and ipN1=ipN1; - } - else { - throw_Exception("Could not detect which element to create", Here()); + ipN2 = std::min(endN, ipN1 + 1); + ipS2 = std::min(endS, ipS1 + 1); } - ipN2 = std::min(endN, ipN1 + 1); - ipS2 = std::min(endS, ipS1 + 1); - } - region.nb_lat_elems.at(jlat) = jelem; + region.nb_lat_elems.at(jlat) = jelem; #if DEBUG_OUTPUT - ATLAS_DEBUG_VAR(region.nb_lat_elems.at(jlat)); + ATLAS_DEBUG_VAR(region.nb_lat_elems.at(jlat)); #endif - if (region.nb_lat_elems.at(jlat) == 0 && latN == region.north) { - ++region.north; - } - if (region.nb_lat_elems.at(jlat) == 0 && latS == region.south) { - --region.south; - } - // region.lat_end.at(latN) = std::min(region.lat_end.at(latN), - // int(rg.nx(latN)-1)); - // region.lat_end.at(latS) = std::min(region.lat_end.at(latS), - // int(rg.nx(latS)-1)); - if (yN == 90 && unique_pole) { - region.lat_end.at(latN) = rg.nx(latN) - 1; - } - if (yS == -90 && unique_pole) { - region.lat_end.at(latS) = rg.nx(latS) - 1; - } + if (region.nb_lat_elems.at(jlat) == 0 && latN == region.north) { + ++region.north; + } + if (region.nb_lat_elems.at(jlat) == 0 && latS == region.south) { + --region.south; + } + // region.lat_end.at(latN) = std::min(region.lat_end.at(latN), + // int(rg.nx(latN)-1)); + // region.lat_end.at(latS) = std::min(region.lat_end.at(latS), + // int(rg.nx(latS)-1)); + if (yN == 90 && unique_pole) { + region.lat_end.at(latN) = rg.nx(latN) - 1; + } + if (yS == -90 && unique_pole) { + region.lat_end.at(latS) = rg.nx(latS) - 1; + } - if (region.nb_lat_elems.at(jlat) > 0) { - region.lat_end.at(latN) = std::max(region.lat_end.at(latN), region.lat_begin.at(latN)); - region.lat_end.at(latS) = std::max(region.lat_end.at(latS), region.lat_begin.at(latS)); - } - } // for jlat + if (region.nb_lat_elems.at(jlat) > 0) { + region.lat_end.at(latN) = std::max(region.lat_end.at(latN), region.lat_begin.at(latN)); + region.lat_end.at(latS) = std::max(region.lat_end.at(latS), region.lat_begin.at(latS)); + } + } // for jlat + } // Log::info() << "nb_triags = " << region.ntriags << std::endl; // Log::info() << "nb_quads = " << region.nquads << std::endl; @@ -1205,6 +1263,10 @@ void StructuredMeshGenerator::generate_mesh(const StructuredGrid& rg, const grid nodes[1] = tmp; }; + bool regular_cells_glb_idx = atlas::RegularLonLatGrid(rg); + if( options.getBool("triangulate") ) { + regular_cells_glb_idx = false; + } for (idx_t jlat = region.north; jlat < region.south; ++jlat) { idx_t ilat = jlat - region.north; @@ -1239,6 +1301,13 @@ void StructuredMeshGenerator::generate_mesh(const StructuredGrid& rg, const grid node_connectivity.set(jcell, quad_nodes); cells_glb_idx(jcell) = jcell + 1; cells_part(jcell) = mypart; + if( regular_cells_glb_idx ) { + gidx_t nx = rg.nx(jlatN) - 1; + if (periodic_east_west) { + ++nx; + } + cells_glb_idx(jcell) = jlatN * nx + region.lat_begin[jlat] + jelem; + } } else // This is a triag { @@ -1499,7 +1568,9 @@ void StructuredMeshGenerator::generate_mesh(const StructuredGrid& rg, const grid } } } - generateGlobalElementNumbering(mesh); + if (not regular_cells_glb_idx) { + generateGlobalElementNumbering(mesh); + } } namespace { diff --git a/src/atlas/parallel/GatherScatter.cc b/src/atlas/parallel/GatherScatter.cc index 10b300b0f..4439fe8be 100644 --- a/src/atlas/parallel/GatherScatter.cc +++ b/src/atlas/parallel/GatherScatter.cc @@ -14,6 +14,8 @@ #include #include +#include "eckit/log/Bytes.h" + #include "atlas/array.h" #include "atlas/array/ArrayView.h" #include "atlas/parallel/GatherScatter.h" @@ -21,6 +23,8 @@ #include "atlas/runtime/Log.h" #include "atlas/runtime/Trace.h" +#include "atlas/parallel/omp/sort.h" + namespace atlas { namespace parallel { @@ -80,70 +84,135 @@ GatherScatter::GatherScatter(const std::string& name): name_(name), is_setup_(fa void GatherScatter::setup(const int part[], const idx_t remote_idx[], const int base, const gidx_t glb_idx[], const int mask[], const idx_t parsize) { ATLAS_TRACE("GatherScatter::setup"); - parsize_ = parsize; glbcounts_.resize(nproc); glbcounts_.assign(nproc, 0); glbdispls_.resize(nproc); glbdispls_.assign(nproc, 0); - const idx_t nvar = 3; - std::vector sendnodes(parsize_ * nvar); + std::vector sendnodes_gidx(parsize_); + std::vector sendnodes_part(parsize_); + std::vector sendnodes_ridx(parsize_); loccnt_ = 0; for (idx_t n = 0; n < parsize_; ++n) { if (!mask[n]) { - sendnodes[loccnt_++] = glb_idx[n]; - sendnodes[loccnt_++] = part[n]; - sendnodes[loccnt_++] = remote_idx[n] - base; + sendnodes_gidx[loccnt_] = glb_idx[n]; + sendnodes_part[loccnt_] = part[n]; + sendnodes_ridx[loccnt_] = remote_idx[n] - base; + ++loccnt_; } } - ATLAS_TRACE_MPI(ALLGATHER) { mpi::comm().allGather(loccnt_, glbcounts_.begin(), glbcounts_.end()); } + ATLAS_TRACE_MPI(ALLGATHER) { + mpi::comm().allGather(loccnt_, glbcounts_.begin(), glbcounts_.end()); + } - glbcnt_ = std::accumulate(glbcounts_.begin(), glbcounts_.end(), 0); + size_t glbcnt_size_t = std::accumulate(glbcounts_.begin(), glbcounts_.end(), size_t(0)); + if (glbcnt_size_t > std::numeric_limits::max()) { + ATLAS_THROW_EXCEPTION("Due to limitation of MPI we cannot use larger counts"); + } + + glbcnt_ = std::accumulate(glbcounts_.begin(), glbcounts_.end(), size_t(0)); glbdispls_[0] = 0; for (idx_t jproc = 1; jproc < nproc; ++jproc) // start at 1 { glbdispls_[jproc] = glbcounts_[jproc - 1] + glbdispls_[jproc - 1]; } - std::vector recvnodes(glbcnt_); - ATLAS_TRACE_MPI(ALLGATHER) { - mpi::comm().allGatherv(sendnodes.begin(), sendnodes.begin() + loccnt_, recvnodes.data(), glbcounts_.data(), - glbdispls_.data()); + idx_t nb_recv_nodes = glbcnt_; + std::vector node_sort; + + try { + node_sort.resize(nb_recv_nodes); + } + catch (const std::bad_alloc& e) { + ATLAS_THROW_EXCEPTION("Could not allocate node_sort with size " << eckit::Bytes(nb_recv_nodes * sizeof(Node))); + } + + { + std::vector recvnodes_gidx; + try { + recvnodes_gidx.resize(glbcnt_); + } + catch (const std::bad_alloc& e) { + ATLAS_THROW_EXCEPTION("Could not allocate recvnodes_gidx with size " + << eckit::Bytes(glbcnt_ * sizeof(gidx_t))); + } + ATLAS_TRACE_MPI(ALLGATHER) { + mpi::comm().allGatherv(sendnodes_gidx.begin(), sendnodes_gidx.begin() + loccnt_, recvnodes_gidx.data(), + glbcounts_.data(), glbdispls_.data()); + } + atlas_omp_parallel_for(idx_t n = 0; n < nb_recv_nodes; ++n) { + node_sort[n].g = recvnodes_gidx[n]; + } + sendnodes_gidx.clear(); } - // Load recvnodes in sorting structure - idx_t nb_recv_nodes = glbcnt_ / nvar; - std::vector node_sort(nb_recv_nodes); - for (idx_t n = 0; n < nb_recv_nodes; ++n) { - node_sort[n].g = recvnodes[n * nvar + 0]; - node_sort[n].p = recvnodes[n * nvar + 1]; - node_sort[n].i = recvnodes[n * nvar + 2]; + + { + std::vector recvnodes_part; + try { + recvnodes_part.resize(glbcnt_); + } + catch (const std::bad_alloc& e) { + ATLAS_THROW_EXCEPTION("Could not allocate recvnodes_part with size " + << eckit::Bytes(glbcnt_ * sizeof(int))); + } + ATLAS_TRACE_MPI(ALLGATHER) { + mpi::comm().allGatherv(sendnodes_part.begin(), sendnodes_part.begin() + loccnt_, recvnodes_part.data(), + glbcounts_.data(), glbdispls_.data()); + } + atlas_omp_parallel_for(idx_t n = 0; n < nb_recv_nodes; ++n) { + node_sort[n].p = recvnodes_part[n]; + } + sendnodes_part.clear(); + } + + { + std::vector recvnodes_ridx; + try { + recvnodes_ridx.resize(glbcnt_); + } + catch (const std::bad_alloc& e) { + ATLAS_THROW_EXCEPTION("Could not allocate recvnodes_ridx with size " + << eckit::Bytes(glbcnt_ * sizeof(idx_t))); + } + ATLAS_TRACE_MPI(ALLGATHER) { + mpi::comm().allGatherv(sendnodes_ridx.begin(), sendnodes_ridx.begin() + loccnt_, recvnodes_ridx.data(), + glbcounts_.data(), glbdispls_.data()); + } + atlas_omp_parallel_for(idx_t n = 0; n < nb_recv_nodes; ++n) { + node_sort[n].i = recvnodes_ridx[n]; + } + sendnodes_ridx.clear(); } - recvnodes.clear(); // Sort on "g" member, and remove duplicates ATLAS_TRACE_SCOPE("sorting") { + // omp::sort(node_sort.begin(), node_sort.end()); std::sort(node_sort.begin(), node_sort.end()); node_sort.erase(std::unique(node_sort.begin(), node_sort.end()), node_sort.end()); } + glbcounts_.assign(nproc, 0); glbdispls_.assign(nproc, 0); + for (size_t n = 0; n < node_sort.size(); ++n) { ++glbcounts_[node_sort[n].p]; } + glbdispls_[0] = 0; for (idx_t jproc = 1; jproc < nproc; ++jproc) // start at 1 { glbdispls_[jproc] = glbcounts_[jproc - 1] + glbdispls_[jproc - 1]; } - glbcnt_ = std::accumulate(glbcounts_.begin(), glbcounts_.end(), 0); + + glbcnt_ = std::accumulate(glbcounts_.begin(), glbcounts_.end(), size_t(0)); loccnt_ = glbcounts_[myproc]; glbmap_.clear(); @@ -152,7 +221,7 @@ void GatherScatter::setup(const int part[], const idx_t remote_idx[], const int locmap_.resize(loccnt_); std::vector idx(nproc, 0); - int n{0}; + size_t n{0}; for (const auto& node : node_sort) { idx_t jproc = node.p; glbmap_[glbdispls_[jproc] + idx[jproc]] = n++; diff --git a/src/atlas/parallel/GatherScatter.h b/src/atlas/parallel/GatherScatter.h index 78a147528..d882acd89 100644 --- a/src/atlas/parallel/GatherScatter.h +++ b/src/atlas/parallel/GatherScatter.h @@ -193,7 +193,7 @@ class GatherScatter : public util::Object { idx_t parsize_; friend class Checksum; - int glb_cnt(idx_t root) const { return myproc == root ? glbcnt_ : 0; } + size_t glb_cnt(idx_t root) const { return myproc == root ? glbcnt_ : 0; } }; //-------------------------------------------------------------------------------------------------- @@ -212,8 +212,8 @@ void GatherScatter::gather(parallel::Field lfields[], parallel: const idx_t gvar_size = std::accumulate(gfields[jfield].var_shape.data(), gfields[jfield].var_shape.data() + gfields[jfield].var_rank, 1, std::multiplies()); - const int loc_size = loccnt_ * lvar_size; - const int glb_size = glb_cnt(root) * gvar_size; + const size_t loc_size = loccnt_ * lvar_size; + const size_t glb_size = glb_cnt(root) * gvar_size; std::vector loc_buffer(loc_size); std::vector glb_buffer(glb_size); std::vector glb_displs(nproc); @@ -230,7 +230,9 @@ void GatherScatter::gather(parallel::Field lfields[], parallel: /// Gather - ATLAS_TRACE_MPI(GATHER) { mpi::comm().gatherv(loc_buffer, glb_buffer, glb_counts, glb_displs, root); } + ATLAS_TRACE_MPI(GATHER) { + mpi::comm().gatherv(loc_buffer, glb_buffer, glb_counts, glb_displs, root); + } /// Unpack if (myproc == root) @@ -261,8 +263,8 @@ void GatherScatter::scatter(parallel::Field gfields[], parallel const int gvar_size = std::accumulate(gfields[jfield].var_shape.data(), gfields[jfield].var_shape.data() + gfields[jfield].var_rank, 1, std::multiplies()); - const int loc_size = loccnt_ * lvar_size; - const int glb_size = glb_cnt(root) * gvar_size; + const size_t loc_size = loccnt_ * lvar_size; + const size_t glb_size = glb_cnt(root) * gvar_size; std::vector loc_buffer(loc_size); std::vector glb_buffer(glb_size); std::vector glb_displs(nproc); @@ -351,7 +353,7 @@ void GatherScatter::unpack_recv_buffer(const std::vector& recvmap, const DA const parallel::Field& field) const { const idx_t recvcnt = static_cast(recvmap.size()); - int ibuf = 0; + size_t ibuf = 0; const idx_t recv_stride = field.var_strides[0] * field.var_shape[0]; switch (field.var_rank) { diff --git a/src/atlas/runtime/trace/CallStack.cc b/src/atlas/runtime/trace/CallStack.cc index 8813f429d..7e32208d7 100644 --- a/src/atlas/runtime/trace/CallStack.cc +++ b/src/atlas/runtime/trace/CallStack.cc @@ -18,11 +18,33 @@ namespace atlas { namespace runtime { namespace trace { +namespace { +static size_t hash_cstr(const char *s) { + // http://www.cse.yorku.ca/~oz/hash.html + size_t h = 5381; + int c; + while ((c = *s++)) { + h = ((h << 5) + h) + c; + } + return h; +} +static size_t hash_combine(size_t h1, size_t h2) { + return h1 ^ (h2 << 1); +}; +static size_t hash_codelocation( const CodeLocation& loc ) { + return hash_combine( hash_cstr(loc.file()), std::hash{}(loc.line()) ); +} +} + void CallStack::push(const CodeLocation& loc, const std::string& id) { if (stack_.size() == size_) { stack_.resize(2 * size_); } - stack_[size_++] = std::hash{}(loc.asString() + id); + auto hash = hash_codelocation(loc); + if( ! id.empty() ) { + hash = hash_combine( hash, std::hash{}(id) ); + } + stack_[size_++] = hash; } void CallStack::pop() { diff --git a/src/atlas/trans/Trans.cc b/src/atlas/trans/Trans.cc index 7780eb78b..d3d4636c3 100644 --- a/src/atlas/trans/Trans.cc +++ b/src/atlas/trans/Trans.cc @@ -103,6 +103,19 @@ void Trans::dirtrans_wind2vordiv(const Field& gpwind, Field& spvor, Field& spdiv get()->dirtrans_wind2vordiv(gpwind, spvor, spdiv, options(config)); } +void Trans::dirtrans_adj(const Field& spfield, Field& gpfield, const eckit::Configuration& config) const { + get()->dirtrans_adj(spfield, gpfield, options(config)); +} + +void Trans::dirtrans_adj(const FieldSet& spfields, FieldSet& gpfields, const eckit::Configuration& config) const { + get()->dirtrans_adj(spfields, gpfields, options(config)); +} + +void Trans::dirtrans_wind2vordiv_adj(const Field& spvor, const Field& spdiv, Field& gpwind, + const eckit::Configuration& config) const { + get()->dirtrans_wind2vordiv_adj(spvor, spdiv, gpwind, options(config)); +} + void Trans::invtrans(const Field& spfield, Field& gpfield, const eckit::Configuration& config) const { get()->invtrans(spfield, gpfield, options(config)); } diff --git a/src/atlas/trans/Trans.h b/src/atlas/trans/Trans.h index fe35a4cb6..75e70fa94 100644 --- a/src/atlas/trans/Trans.h +++ b/src/atlas/trans/Trans.h @@ -74,6 +74,13 @@ class Trans : DOXYGEN_HIDE(public util::ObjectHandle) { void dirtrans_wind2vordiv(const Field& gpwind, Field& spvor, Field& spdiv, const eckit::Configuration& = util::NoConfig()) const; + void dirtrans_adj(const Field& spfield, Field& gpfield, const eckit::Configuration& = util::NoConfig()) const; + + void dirtrans_adj(const FieldSet& spfields, FieldSet& gpfields, const eckit::Configuration& = util::NoConfig()) const; + + void dirtrans_wind2vordiv_adj(const Field& spvor, const Field& spdiv, Field& gpwind, + const eckit::Configuration& = util::NoConfig()) const; + void invtrans(const Field& spfield, Field& gpfield, const eckit::Configuration& = util::NoConfig()) const; void invtrans(const FieldSet& spfields, FieldSet& gpfields, const eckit::Configuration& = util::NoConfig()) const; diff --git a/src/atlas/trans/detail/TransImpl.h b/src/atlas/trans/detail/TransImpl.h index 274be36f4..cfa024882 100644 --- a/src/atlas/trans/detail/TransImpl.h +++ b/src/atlas/trans/detail/TransImpl.h @@ -60,6 +60,15 @@ class TransImpl : public util::Object { virtual void dirtrans_wind2vordiv(const Field& gpwind, Field& spvor, Field& spdiv, const eckit::Configuration& = util::NoConfig()) const = 0; + virtual void dirtrans_adj(const Field& spfield, Field& gpfield, + const eckit::Configuration& = util::NoConfig()) const = 0; + + virtual void dirtrans_adj(const FieldSet& spfields, FieldSet& gpfields, + const eckit::Configuration& = util::NoConfig()) const = 0; + + virtual void dirtrans_wind2vordiv_adj(const Field& spvor, const Field& spdiv, Field& gpwind, + const eckit::Configuration& = util::NoConfig()) const = 0; + virtual void invtrans(const Field& spfield, Field& gpfield, const eckit::Configuration& = util::NoConfig()) const = 0; diff --git a/src/atlas/trans/ifs/TransIFS.cc b/src/atlas/trans/ifs/TransIFS.cc index 3f15c2efe..1e4076da1 100644 --- a/src/atlas/trans/ifs/TransIFS.cc +++ b/src/atlas/trans/ifs/TransIFS.cc @@ -174,6 +174,78 @@ void TransIFS::dirtrans(const FieldSet& gpfields, FieldSet& spfields, const ecki } } +void TransIFS::dirtrans_wind2vordiv(const Field& gpwind, Field& spvor, Field& spdiv, + const eckit::Configuration& config) const { + ATLAS_ASSERT(Spectral(spvor.functionspace())); + ATLAS_ASSERT(Spectral(spdiv.functionspace())); + if (StructuredColumns(gpwind.functionspace())) { + __dirtrans_wind2vordiv(StructuredColumns(gpwind.functionspace()), gpwind, Spectral(spvor.functionspace()), + spvor, spdiv, config); + } + else if (NodeColumns(gpwind.functionspace())) { + __dirtrans_wind2vordiv(NodeColumns(gpwind.functionspace()), gpwind, Spectral(spvor.functionspace()), spvor, + spdiv, config); + } + else { + ATLAS_NOTIMPLEMENTED; + } +} + + +void TransIFS::dirtrans_adj(const Field& spfield, Field& gpfield, const eckit::Configuration& config) const { + ATLAS_ASSERT(functionspace::Spectral(spfield.functionspace())); + if (functionspace::StructuredColumns(gpfield.functionspace())) { + __dirtrans_adj(functionspace::Spectral(spfield.functionspace()), spfield, + functionspace::StructuredColumns(gpfield.functionspace()), gpfield, + config); + } + else if (functionspace::NodeColumns(gpfield.functionspace())) { + __dirtrans_adj(functionspace::Spectral(spfield.functionspace()), spfield, + functionspace::NodeColumns(gpfield.functionspace()), gpfield, + config); + } + else { + ATLAS_NOTIMPLEMENTED; + } +} + +void TransIFS::dirtrans_adj(const FieldSet& spfields, FieldSet& gpfields, const eckit::Configuration& config) const { + ATLAS_TRACE(); + assert_spectral_functionspace(spfields); + std::string functionspace(fieldset_functionspace(gpfields)); + if (functionspace == StructuredColumns::type()) { + __dirtrans_adj(Spectral(spfields[0].functionspace()), spfields, + StructuredColumns(gpfields[0].functionspace()), gpfields, + config); + } + else if (functionspace == NodeColumns::type()) { + // __dirtrans_adj(NodeColumns(gpfields[0].functionspace()), gpfields, Spectral(spfields[0].functionspace()), spfields, + // config); + } + else { + ATLAS_NOTIMPLEMENTED; + } +} + +void TransIFS::dirtrans_wind2vordiv_adj(const Field& spvor, const Field& spdiv, Field& gpwind, + const eckit::Configuration& config) const { + ATLAS_ASSERT(Spectral(spvor.functionspace())); + ATLAS_ASSERT(Spectral(spdiv.functionspace())); + if (StructuredColumns(gpwind.functionspace())) { + __dirtrans_wind2vordiv_adj(Spectral(spvor.functionspace()), spvor, spdiv, + StructuredColumns(gpwind.functionspace()), gpwind, + config); + } + else if (NodeColumns(gpwind.functionspace())) { +// __dirtrans_wind2vordiv_adj(Spectral(spvor.functionspace()), spvor, spdiv, +// NodeColumns(gpwind.functionspace()), gpwind, +// config); + } + else { + ATLAS_NOTIMPLEMENTED; + } +} + void TransIFS::invtrans(const Field& spfield, Field& gpfield, const eckit::Configuration& config) const { ATLAS_TRACE(); ATLAS_ASSERT(Spectral(spfield.functionspace())); @@ -311,24 +383,6 @@ void TransIFS::invtrans_grad_adj(const FieldSet& gradfields, FieldSet& spfields, } } -// -------------------------------------------------------------------------------------------- - -void TransIFS::dirtrans_wind2vordiv(const Field& gpwind, Field& spvor, Field& spdiv, - const eckit::Configuration& config) const { - ATLAS_ASSERT(Spectral(spvor.functionspace())); - ATLAS_ASSERT(Spectral(spdiv.functionspace())); - if (StructuredColumns(gpwind.functionspace())) { - __dirtrans_wind2vordiv(StructuredColumns(gpwind.functionspace()), gpwind, Spectral(spvor.functionspace()), - spvor, spdiv, config); - } - else if (NodeColumns(gpwind.functionspace())) { - __dirtrans_wind2vordiv(NodeColumns(gpwind.functionspace()), gpwind, Spectral(spvor.functionspace()), spvor, - spdiv, config); - } - else { - ATLAS_NOTIMPLEMENTED; - } -} void TransIFS::invtrans_vordiv2wind(const Field& spvor, const Field& spdiv, Field& gpwind, @@ -353,7 +407,6 @@ void TransIFS::invtrans_vordiv2wind_adj(const Field& gpwind, Field& spvor, Field const eckit::Configuration& config) const { ATLAS_ASSERT(Spectral(spvor.functionspace())); ATLAS_ASSERT(Spectral(spdiv.functionspace())); - ATLAS_ASSERT(NodeColumns(gpwind.functionspace())); if (StructuredColumns(gpwind.functionspace())) { __invtrans_vordiv2wind_adj(Spectral(spvor.functionspace()), spvor, spdiv, StructuredColumns(gpwind.functionspace()), gpwind, config); @@ -1473,6 +1526,369 @@ void TransIFS::__dirtrans_wind2vordiv(const functionspace::NodeColumns& gp, cons unpack_div(spdiv); } +// -------------------------------------------------------------------------------------------- + +void TransIFS::__dirtrans_adj(const Spectral& sp, const FieldSet& spfields, + const StructuredColumns& gp, FieldSet& gpfields, + const eckit::Configuration&) const { +#if ATLAS_HAVE_ECTRANS + assertCompatibleDistributions(gp, sp); + + // Count total number of fields and do sanity checks + const idx_t nfld = compute_nfld(gpfields); + for (idx_t jfld = 0; jfld < gpfields.size(); ++jfld) { + const Field& f = gpfields[jfld]; + ATLAS_ASSERT(f.functionspace() == 0 || functionspace::StructuredColumns(f.functionspace())); + } + + const int trans_sp_nfld = compute_nfld(spfields); + + if (nfld != trans_sp_nfld) { + throw_Exception("dirtrans_adj: different number of gridpoint fields than spectral fields", Here()); + } + // Arrays Trans expects + std::vector rgp(nfld * ngptot()); + std::vector rsp(nspec2() * nfld); + auto rgpview = LocalView(rgp.data(), make_shape(nfld, ngptot())); + auto rspview = LocalView(rsp.data(), make_shape(nspec2(), nfld)); + + // Pack spectral fields + { + PackSpectral pack(rspview); + for (idx_t jfld = 0; jfld < spfields.size(); ++jfld) { + pack(spfields[jfld]); + } + } + + // Do transform + { + struct ::DirTransAdj_t transform = ::new_dirtrans_adj(trans_.get()); + transform.nscalar = int(nfld); + transform.rgp = rgp.data(); + transform.rspscalar = rsp.data(); + + TRANS_CHECK(::trans_dirtrans_adj(&transform)); + } + + // Unpack the gridpoint fields + { + UnpackStructuredColumns unpack(rgpview); + for (idx_t jfld = 0; jfld < gpfields.size(); ++jfld) { + unpack(gp, gpfields[jfld]); + } + } +#else + ATLAS_NOTIMPLEMENTED; +#endif +} + +// -------------------------------------------------------------------------------------------- + +void TransIFS::__dirtrans_adj( const Spectral& sp, const Field& spfield, + const NodeColumns& gp, Field& gpfield, + const eckit::Configuration& ) const { +#if ATLAS_HAVE_ECTRANS + ATLAS_ASSERT(gpfield.functionspace() == 0 || functionspace::NodeColumns(gpfield.functionspace())); + ATLAS_ASSERT(spfield.functionspace() == 0 || functionspace::Spectral(spfield.functionspace())); + + assertCompatibleDistributions(gp, sp); + + if (compute_nfld(gpfield) != compute_nfld(spfield)) { + throw_Exception("dirtrans: different number of gridpoint fields than spectral fields", Here()); + } + if ((int)gpfield.shape(0) < ngptot()) { + throw_Exception("dirtrans: slowest moving index must be >= ngptot", Here()); + } + const int nfld = compute_nfld(gpfield); + + // Arrays Trans expects + std::vector rgp(nfld * ngptot()); + std::vector rsp(nspec2() * nfld); + auto rgpview = LocalView(rgp.data(), make_shape(nfld, ngptot())); + auto rspview = LocalView(rsp.data(), make_shape(nspec2(), nfld)); + + // Pack spectral fields + { + PackSpectral pack(rspview); + pack(spfield); + } + + // Do transform + { + struct ::DirTransAdj_t transform = ::new_dirtrans_adj(trans_.get()); + transform.nscalar = nfld; + transform.rgp = rgp.data(); + transform.rspscalar = rsp.data(); + TRANS_CHECK(::trans_dirtrans_adj(&transform)); + } + + // Unpack gridpoint fields + { + UnpackNodeColumns unpack(rgpview, gp); + unpack(gpfield); + } +#else + ATLAS_NOTIMPLEMENTED; +#endif +} + +// -------------------------------------------------------------------------------------------- + +void TransIFS::__dirtrans_adj( const Spectral& sp, const Field& spfield, + const StructuredColumns& gp, Field& gpfield, + const eckit::Configuration& ) const { +#if ATLAS_HAVE_ECTRANS + ATLAS_ASSERT(gpfield.functionspace() == 0 || functionspace::StructuredColumns(gpfield.functionspace())); + ATLAS_ASSERT(spfield.functionspace() == 0 || functionspace::Spectral(spfield.functionspace())); + + assertCompatibleDistributions(gp, sp); + + if (compute_nfld(gpfield) != compute_nfld(spfield)) { + throw_Exception("dirtrans: different number of gridpoint fields than spectral fields", Here()); + } + if ((int)gpfield.shape(0) < ngptot()) { + throw_Exception("dirtrans: slowest moving index must be >= ngptot", Here()); + } + const int nfld = compute_nfld(gpfield); + + // Arrays Trans expects + std::vector rgp(nfld * ngptot()); + std::vector rsp(nspec2() * nfld); + auto rgpview = LocalView(rgp.data(), make_shape(nfld, ngptot())); + auto rspview = LocalView(rsp.data(), make_shape(nspec2(), nfld)); + + // Pack spectral fields + { + PackSpectral pack(rspview); + pack(spfield); + } + + // Do transform + { + struct ::DirTransAdj_t transform = ::new_dirtrans_adj(trans_.get()); + transform.nscalar = nfld; + transform.rgp = rgp.data(); + transform.rspscalar = rsp.data(); + TRANS_CHECK(::trans_dirtrans_adj(&transform)); + } + + // Unpack gridpoint fields + { + UnpackStructuredColumns unpack(rgpview); + unpack(gp, gpfield); + } +#else + ATLAS_NOTIMPLEMENTED; +#endif +} + + +// -------------------------------------------------------------------------------------------- + +void TransIFS::__dirtrans_adj(const Spectral& sp, const FieldSet& spfields, + const NodeColumns& gp, FieldSet& gpfields, + const eckit::Configuration&) const { +#if ATLAS_HAVE_ECTRANS + assertCompatibleDistributions(gp, sp); + + // Count total number of fields and do sanity checks + const idx_t nfld = compute_nfld(gpfields); + for (idx_t jfld = 0; jfld < gpfields.size(); ++jfld) { + const Field& f = gpfields[jfld]; + ATLAS_ASSERT(f.functionspace() == 0 || functionspace::NodeColumns(f.functionspace())); + } + + const int trans_sp_nfld = compute_nfld(spfields); + + if (nfld != trans_sp_nfld) { + throw_Exception("dirtrans_adj: different number of gridpoint fields than spectral fields", Here()); + } + // Arrays Trans expects + std::vector rgp(nfld * ngptot()); + std::vector rsp(nspec2() * nfld); + auto rgpview = LocalView(rgp.data(), make_shape(nfld, ngptot())); + auto rspview = LocalView(rsp.data(), make_shape(nspec2(), nfld)); + + // Pack spectral fields + { + PackSpectral pack(rspview); + for (idx_t jfld = 0; jfld < spfields.size(); ++jfld) { + pack(spfields[jfld]); + } + } + + // Do transform + { + struct ::DirTransAdj_t transform = ::new_dirtrans_adj(trans_.get()); + transform.nscalar = int(nfld); + transform.rgp = rgp.data(); + transform.rspscalar = rsp.data(); + + TRANS_CHECK(::trans_dirtrans_adj(&transform)); + } + + // Unpack the gridpoint fields + { + UnpackNodeColumns unpack(rgpview, gp); + for (idx_t jfld = 0; jfld < gpfields.size(); ++jfld) { + unpack(gpfields[jfld]); + } + } +#else + ATLAS_NOTIMPLEMENTED; +#endif +} + +// -------------------------------------------------------------------------------------------- + +void TransIFS::__dirtrans_wind2vordiv_adj(const Spectral& sp, const Field& spvor, const Field& spdiv, + const functionspace::StructuredColumns& gp, Field& gpwind, + const eckit::Configuration&) const { +#if ATLAS_HAVE_ECTRANS + assertCompatibleDistributions(gp, sp); + + // Count total number of fields and do sanity checks + const size_t nfld = compute_nfld(spvor); + if (spdiv.shape(0) != spvor.shape(0)) { + throw_Exception("dirtrans_adj:vorticity not compatible with divergence for 1st dimension", Here()); + } + if (spdiv.shape(1) != spvor.shape(1)) { + throw_Exception("dirtrans_adj:vorticity not compatible with divergence for 2st dimension", Here()); + } + const size_t nwindfld = compute_nfld(gpwind); + if (nwindfld != 2 * nfld && nwindfld != 3 * nfld) { + throw_Exception("dirtrans_adj:wind field is not compatible with vorticity, divergence.", Here()); + } + + if (spdiv.shape(0) != nspec2()) { + std::stringstream msg; + msg << "dirtrans_adj: Spectral vorticity and divergence have wrong dimension: " + "nspec2 " + << spdiv.shape(0) << " should be " << nspec2(); + throw_Exception(msg.str(), Here()); + } + + if (spvor.size() == 0) { + throw_Exception("dirtrans_adj: spectral vorticity field is empty."); + } + if (spdiv.size() == 0) { + throw_Exception("dirtrans_adj: spectral divergence field is empty."); + } + + // Arrays Trans expects + std::vector rgp(2 * nfld * ngptot()); + std::vector rspvor(nspec2() * nfld); + std::vector rspdiv(nspec2() * nfld); + auto rgpview = LocalView(rgp.data(), make_shape(2 * nfld, ngptot())); + auto rspvorview = LocalView(rspvor.data(), make_shape(nspec2(), nfld)); + auto rspdivview = LocalView(rspdiv.data(), make_shape(nspec2(), nfld)); + + // Pack spectral fields + PackSpectral pack_vor(rspvorview); + pack_vor(spvor); + PackSpectral pack_div(rspdivview); + pack_div(spdiv); + + // Do transform + { + struct ::DirTransAdj_t transform = ::new_dirtrans_adj(trans_.get()); + transform.nvordiv = int(nfld); + transform.rgp = rgp.data(); + transform.rspvor = rspvor.data(); + transform.rspdiv = rspdiv.data(); + + ATLAS_ASSERT(transform.rspvor); + ATLAS_ASSERT(transform.rspdiv); + TRANS_CHECK(::trans_dirtrans_adj(&transform)); + } + + // Unpack the gridpoint fields + { + UnpackStructuredColumns unpack(rgpview); + int wind_components = 2; + unpack(gp, gpwind, wind_components); + } +#else + ATLAS_NOTIMPLEMENTED; +#endif +} + +// -------------------------------------------------------------------------------------------- + +void TransIFS::__dirtrans_wind2vordiv_adj(const Spectral& sp, const Field& spvor, const Field& spdiv, + const functionspace::NodeColumns& gp, Field& gpwind, + const eckit::Configuration&) const { +#if ATLAS_HAVE_ECTRANS + assertCompatibleDistributions(gp, sp); + + // Count total number of fields and do sanity checks + const size_t nfld = compute_nfld(spvor); + if (spdiv.shape(0) != spvor.shape(0)) { + throw_Exception("dirtrans_adj:vorticity not compatible with divergence for 1st dimension", Here()); + } + if (spdiv.shape(1) != spvor.shape(1)) { + throw_Exception("dirtrans_adj:vorticity not compatible with divergence for 2st dimension", Here()); + } + const size_t nwindfld = compute_nfld(gpwind); + if (nwindfld != 2 * nfld && nwindfld != 3 * nfld) { + throw_Exception("dirtrans_adj:wind field is not compatible with vorticity, divergence.", Here()); + } + + if (spdiv.shape(0) != nspec2()) { + std::stringstream msg; + msg << "dirtrans_adj: Spectral vorticity and divergence have wrong dimension: " + "nspec2 " + << spdiv.shape(0) << " should be " << nspec2(); + throw_Exception(msg.str(), Here()); + } + + if (spvor.size() == 0) { + throw_Exception("dirtrans_adj: spectral vorticity field is empty."); + } + if (spdiv.size() == 0) { + throw_Exception("dirtrans_adj: spectral divergence field is empty."); + } + + // Arrays Trans expects + std::vector rgp(2 * nfld * ngptot()); + std::vector rspvor(nspec2() * nfld); + std::vector rspdiv(nspec2() * nfld); + auto rgpview = LocalView(rgp.data(), make_shape(2 * nfld, ngptot())); + auto rspvorview = LocalView(rspvor.data(), make_shape(nspec2(), nfld)); + auto rspdivview = LocalView(rspdiv.data(), make_shape(nspec2(), nfld)); + + // Pack spectral fields + PackSpectral pack_vor(rspvorview); + pack_vor(spvor); + PackSpectral pack_div(rspdivview); + pack_div(spdiv); + + // Do transform + { + struct ::DirTransAdj_t transform = ::new_dirtrans_adj(trans_.get()); + transform.nvordiv = int(nfld); + transform.rgp = rgp.data(); + transform.rspvor = rspvor.data(); + transform.rspdiv = rspdiv.data(); + + ATLAS_ASSERT(transform.rspvor); + ATLAS_ASSERT(transform.rspdiv); + TRANS_CHECK(::trans_dirtrans_adj(&transform)); + } + + // Unpack the gridpoint fields + { + UnpackNodeColumns unpack(rgpview, gp); + int wind_components = 2; + unpack(gpwind, wind_components); + } +#else + ATLAS_NOTIMPLEMENTED; +#endif +} + + + // -------------------------------------------------------------------------------------------- void TransIFS::__invtrans(const functionspace::Spectral& sp, const Field& spfield, diff --git a/src/atlas/trans/ifs/TransIFS.h b/src/atlas/trans/ifs/TransIFS.h index 7fff8eed7..0f8b4f5d2 100644 --- a/src/atlas/trans/ifs/TransIFS.h +++ b/src/atlas/trans/ifs/TransIFS.h @@ -108,6 +108,15 @@ class TransIFS : public trans::TransImpl { virtual void dirtrans_wind2vordiv(const Field& gpwind, Field& spvor, Field& spdiv, const eckit::Configuration& = util::NoConfig()) const override; + virtual void dirtrans_adj(const Field& spfield, Field& gpfield, + const eckit::Configuration& = util::NoConfig()) const override; + + virtual void dirtrans_adj(const FieldSet& spfields, FieldSet& gpfields, + const eckit::Configuration& = util::NoConfig()) const override; + + virtual void dirtrans_wind2vordiv_adj(const Field& spvor, const Field& spdiv, Field& gpwind, + const eckit::Configuration& = util::NoConfig()) const override; + virtual void invtrans(const Field& spfield, Field& gpfield, const eckit::Configuration& = util::NoConfig()) const override; @@ -249,6 +258,27 @@ class TransIFS : public trans::TransImpl { void __dirtrans_wind2vordiv(const functionspace::NodeColumns&, const Field& gpwind, const functionspace::Spectral&, Field& spvor, Field& spdiv, const eckit::Configuration& = util::NoConfig()) const; + void __dirtrans_adj(const functionspace::Spectral&, const Field& spfield, + const functionspace::StructuredColumns&, Field& gpfield, + const eckit::Configuration& = util::NoConfig()) const; + void __dirtrans_adj(const functionspace::Spectral&, const Field& spfield, + const functionspace::NodeColumns&, Field& gpfield, + const eckit::Configuration& = util::NoConfig()) const; + + void __dirtrans_adj(const functionspace::Spectral&, const FieldSet& spfields, + const functionspace::StructuredColumns&, FieldSet& gpfields, + const eckit::Configuration& = util::NoConfig()) const; + void __dirtrans_adj(const functionspace::Spectral&, const FieldSet& spfields, + const functionspace::NodeColumns&, FieldSet& gpfields, + const eckit::Configuration& = util::NoConfig()) const; + + void __dirtrans_wind2vordiv_adj(const functionspace::Spectral&, const Field& spvor, const Field& spdiv, + const functionspace::StructuredColumns&, Field& gpwind, + const eckit::Configuration& = util::NoConfig()) const; + void __dirtrans_wind2vordiv_adj(const functionspace::Spectral&, const Field& spvor, const Field& spdiv, + const functionspace::NodeColumns&, Field& gpwind, + const eckit::Configuration& = util::NoConfig()) const; + void __invtrans(const functionspace::Spectral&, const Field& spfield, const functionspace::StructuredColumns&, Field& gpfield, const eckit::Configuration& = util::NoConfig()) const; void __invtrans(const functionspace::Spectral&, const Field& spfield, const functionspace::NodeColumns&, diff --git a/src/atlas/trans/local/TransLocal.cc b/src/atlas/trans/local/TransLocal.cc index 4eb5c790b..b448b486b 100644 --- a/src/atlas/trans/local/TransLocal.cc +++ b/src/atlas/trans/local/TransLocal.cc @@ -1674,6 +1674,29 @@ void TransLocal::dirtrans_wind2vordiv(const Field& gpwind, Field& spvor, Field& // Use the TransIFS implementation instead. } + +void TransLocal::dirtrans_adj(const Field& spfield, Field& gpfield, + const eckit::Configuration& config) const { + ATLAS_NOTIMPLEMENTED; + // Not implemented and not planned. + // Use the TransIFS implementation instead. +} + +void TransLocal::dirtrans_adj(const FieldSet& spfields, FieldSet& gpfields, + const eckit::Configuration& config) const { + ATLAS_NOTIMPLEMENTED; + // Not implemented and not planned. + // Use the TransIFS implementation instead. +} + +void TransLocal::dirtrans_wind2vordiv_adj(const Field& spvor, const Field& spdiv, Field& gpwind, + const eckit::Configuration& config) const { + ATLAS_NOTIMPLEMENTED; + // Not implemented and not planned. + // Use the TransIFS implementation instead. +} + + // -------------------------------------------------------------------------------------------------------------------- void TransLocal::dirtrans(const int nb_fields, const double scalar_fields[], double scalar_spectra[], diff --git a/src/atlas/trans/local/TransLocal.h b/src/atlas/trans/local/TransLocal.h index bc3ac0cf9..7d1883406 100644 --- a/src/atlas/trans/local/TransLocal.h +++ b/src/atlas/trans/local/TransLocal.h @@ -156,6 +156,15 @@ class TransLocal : public trans::TransImpl { virtual void dirtrans_wind2vordiv(const Field& gpwind, Field& spvor, Field& spdiv, const eckit::Configuration& = util::NoConfig()) const override; + virtual void dirtrans_adj(const Field& spfield, Field& gpfield, + const eckit::Configuration& = util::NoConfig()) const override; + + virtual void dirtrans_adj(const FieldSet& spfields, FieldSet& gpfields, + const eckit::Configuration& = util::NoConfig()) const override; + + virtual void dirtrans_wind2vordiv_adj(const Field& spvor, const Field& spdiv, Field& gpwind, + const eckit::Configuration& = util::NoConfig()) const override; + virtual void dirtrans(const int nb_fields, const double scalar_fields[], double scalar_spectra[], const eckit::Configuration& = util::NoConfig()) const override; diff --git a/src/atlas_f/functionspace/atlas_FunctionSpace_module.F90 b/src/atlas_f/functionspace/atlas_FunctionSpace_module.F90 index 93b0af2aa..c99896c6b 100644 --- a/src/atlas_f/functionspace/atlas_FunctionSpace_module.F90 +++ b/src/atlas_f/functionspace/atlas_FunctionSpace_module.F90 @@ -55,6 +55,9 @@ module atlas_functionspace_module procedure, private :: halo_exchange_field procedure, private :: halo_exchange_fieldset + procedure, private :: adjoint_halo_exchange_field + procedure, private :: adjoint_halo_exchange_fieldset + generic, public :: create_field => & & create_field_args, & & create_field_template, & @@ -62,6 +65,7 @@ module atlas_functionspace_module & deprecated_create_field_2 generic, public :: halo_exchange => halo_exchange_field, halo_exchange_fieldset + generic, public :: adjoint_halo_exchange => adjoint_halo_exchange_field, adjoint_halo_exchange_fieldset #if FCKIT_FINAL_NOT_INHERITING final :: atlas_FunctionSpace__final_auto @@ -179,6 +183,25 @@ subroutine halo_exchange_field(this,field) call atlas__FunctionSpace__halo_exchange_field(this%CPTR_PGIBUG_A,field%CPTR_PGIBUG_A) end subroutine + +!------------------------------------------------------------------------------ + +subroutine adjoint_halo_exchange_fieldset(this,fieldset) + use atlas_functionspace_c_binding + class(atlas_Functionspace), intent(in) :: this + type(atlas_FieldSet), intent(inout) :: fieldset + call atlas__FunctionSpace__adjoint_halo_exchange_fieldset(this%CPTR_PGIBUG_A,fieldset%CPTR_PGIBUG_A) +end subroutine + +!------------------------------------------------------------------------------ + +subroutine adjoint_halo_exchange_field(this,field) + use atlas_functionspace_c_binding + class(atlas_Functionspace), intent(in) :: this + type(atlas_Field), intent(inout) :: field + call atlas__FunctionSpace__adjoint_halo_exchange_field(this%CPTR_PGIBUG_A,field%CPTR_PGIBUG_A) +end subroutine + !------------------------------------------------------------------------------ !------------------------------------------------------------------------------ diff --git a/src/atlas_f/functionspace/atlas_functionspace_PointCloud_module.F90 b/src/atlas_f/functionspace/atlas_functionspace_PointCloud_module.F90 index f30896ac1..5e9db35ff 100644 --- a/src/atlas_f/functionspace/atlas_functionspace_PointCloud_module.F90 +++ b/src/atlas_f/functionspace/atlas_functionspace_PointCloud_module.F90 @@ -14,6 +14,7 @@ module atlas_functionspace_PointCloud_module use atlas_functionspace_module, only : atlas_FunctionSpace use atlas_Grid_module, only: atlas_Grid use atlas_Field_module, only: atlas_Field +use atlas_FieldSet_module, only: atlas_FieldSet use atlas_kinds_module, only: ATLAS_KIND_GIDX implicit none @@ -21,6 +22,7 @@ module atlas_functionspace_PointCloud_module private :: c_str, c_ptr_to_string, c_ptr_free private :: atlas_FunctionSpace private :: atlas_Field +private :: atlas_FieldSet private :: atlas_Grid private :: ATLAS_KIND_GIDX @@ -54,6 +56,7 @@ module atlas_functionspace_PointCloud_module module procedure ctor_cptr module procedure ctor_lonlat module procedure ctor_lonlat_ghost + module procedure ctor_fieldset module procedure ctor_grid end interface @@ -96,6 +99,16 @@ function ctor_lonlat_ghost(lonlat,ghost) result(this) !------------------------------------------------------------------------------ +function ctor_fieldset(fset) result(this) + use atlas_functionspace_PointCloud_c_binding + type(atlas_functionspace_PointCloud) :: this + class(atlas_FieldSet), intent(in) :: fset + call this%reset_c_ptr( atlas__functionspace__PointCloud__new__fieldset( fset%CPTR_PGIBUG_A ) ) + call this%return() +end function + +!------------------------------------------------------------------------------ + function ctor_grid(grid) result(this) use atlas_functionspace_PointCloud_c_binding type(atlas_functionspace_PointCloud) :: this diff --git a/src/atlas_f/runtime/atlas_Trace_module.F90 b/src/atlas_f/runtime/atlas_Trace_module.F90 index 26595d948..1faa0e8d1 100644 --- a/src/atlas_f/runtime/atlas_Trace_module.F90 +++ b/src/atlas_f/runtime/atlas_Trace_module.F90 @@ -77,6 +77,7 @@ function atlas_Trace__labels_1(file,line,title,label) result(this) character(len=*) , intent(in) :: label call this%reset_c_ptr( new_atlas_Trace_labels_1( c_str(file), line, c_str(title), c_str(label) ), & & fckit_c_deleter(delete_atlas_Trace) ) + call this%return() end function function atlas_Trace__labels_2(file,line,title,label1,label2) result(this) use, intrinsic :: iso_c_binding, only : c_ptr @@ -90,6 +91,7 @@ function atlas_Trace__labels_2(file,line,title,label1,label2) result(this) character(len=*) , intent(in) :: label2 call this%reset_c_ptr( new_atlas_Trace_labels_2( c_str(file), line, c_str(title), c_str(label1), c_str(label2) ), & & fckit_c_deleter(delete_atlas_Trace) ) + call this%return() end function function atlas_Trace__labels_3(file,line,title,label1,label2,label3) result(this) use, intrinsic :: iso_c_binding, only : c_ptr @@ -104,6 +106,7 @@ function atlas_Trace__labels_3(file,line,title,label1,label2,label3) result(this character(len=*) , intent(in) :: label3 call this%reset_c_ptr( new_atlas_Trace_labels_3( c_str(file), line, c_str(title), c_str(label1), c_str(label2), & & c_str(label3) ), fckit_c_deleter(delete_atlas_Trace) ) + call this%return() end function function atlas_Trace__labels_4(file,line,title,label1,label2,label3,label4) result(this) use, intrinsic :: iso_c_binding, only : c_ptr @@ -119,6 +122,7 @@ function atlas_Trace__labels_4(file,line,title,label1,label2,label3,label4) resu character(len=*) , intent(in) :: label4 call this%reset_c_ptr( new_atlas_Trace_labels_4( c_str(file), line, c_str(title), c_str(label1), c_str(label2), & & c_str(label3), c_str(label4) ), fckit_c_deleter(delete_atlas_Trace) ) + call this%return() end function function atlas_Trace__labels_5(file,line,title,label1,label2,label3,label4,label5) result(this) use, intrinsic :: iso_c_binding, only : c_ptr @@ -135,6 +139,7 @@ function atlas_Trace__labels_5(file,line,title,label1,label2,label3,label4,label character(len=*) , intent(in) :: label5 call this%reset_c_ptr( new_atlas_Trace_labels_5( c_str(file), line, c_str(title), c_str(label1), c_str(label2), & & c_str(label3), c_str(label4), c_str(label5) ), fckit_c_deleter(delete_atlas_Trace) ) + call this%return() end function !------------------------------------------------------------------------------- diff --git a/src/sandbox/interpolation/atlas-conservative-interpolation.cc b/src/sandbox/interpolation/atlas-conservative-interpolation.cc index 5242e7321..36f8b67b4 100644 --- a/src/sandbox/interpolation/atlas-conservative-interpolation.cc +++ b/src/sandbox/interpolation/atlas-conservative-interpolation.cc @@ -187,9 +187,10 @@ int AtlasParallelInterpolation::execute(const AtlasTool::Args& args) { timers.target_setup.stop(); timers.source_setup.start(); - auto src_meshgenerator = MeshGenerator{src_grid.meshgenerator() | option::halo(2)}; - auto src_partitioner = grid::MatchingPartitioner{tgt_mesh}; - auto src_mesh = src_meshgenerator.generate(src_grid, src_partitioner); + auto src_meshgenerator = + MeshGenerator{src_grid.meshgenerator() | option::halo(2) | util::Config("pole_elements", "pentagons")}; + auto src_partitioner = grid::MatchingPartitioner{tgt_mesh}; + auto src_mesh = src_meshgenerator.generate(src_grid, src_partitioner); auto src_functionspace = create_functionspace(src_mesh, args.getLong("source.halo", 2), args.getString("source.functionspace", "")); auto src_field = src_functionspace.createField(); diff --git a/src/tests/functionspace/CMakeLists.txt b/src/tests/functionspace/CMakeLists.txt index 53ec29d3a..7d81a7b9d 100644 --- a/src/tests/functionspace/CMakeLists.txt +++ b/src/tests/functionspace/CMakeLists.txt @@ -6,7 +6,6 @@ # granted to it by virtue of its status as an intergovernmental organisation nor # does it submit to any jurisdiction. -ecbuild_warn_var( HAVE_FCTEST ) if( HAVE_FCTEST ) if( NOT HAVE_TRANS ) @@ -14,11 +13,6 @@ if( HAVE_FCTEST ) set( ectrans_HAVE_MPI 1 ) endif() - ecbuild_warn_var( HAVE_TRANS ) - ecbuild_warn_var( eckit_HAVE_MPI ) - ecbuild_warn_var( transi_HAVE_MPI ) - - add_fctest( TARGET atlas_fctest_functionspace MPI 4 CONDITION eckit_HAVE_MPI AND (transi_HAVE_MPI OR ectrans_HAVE_MPI) @@ -76,7 +70,9 @@ ecbuild_add_test( TARGET atlas_test_structuredcolumns ecbuild_add_test( TARGET atlas_test_pointcloud SOURCES test_pointcloud.cc LIBS atlas + MPI 2 ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} + CONDITION eckit_HAVE_MPI AND NOT HAVE_GRIDTOOLS_STORAGE ) ecbuild_add_test( TARGET atlas_test_reduced_halo diff --git a/src/tests/functionspace/fctest_functionspace.F90 b/src/tests/functionspace/fctest_functionspace.F90 index 6eb49b05c..a7203cc2b 100644 --- a/src/tests/functionspace/fctest_functionspace.F90 +++ b/src/tests/functionspace/fctest_functionspace.F90 @@ -587,15 +587,19 @@ module fcta_FunctionSpace_fxt TEST( test_pointcloud ) #if 1 +implicit none type(atlas_StructuredGrid) :: grid type(atlas_functionspace_PointCloud) :: fs type(atlas_functionspace) :: fs_base +type(atlas_trace) :: trace character(len=10) str type(atlas_Field) :: field, field2 type(atlas_Field) :: field_lonlat real(8), pointer :: lonlat(:,:), x(:) +trace = atlas_Trace("fctest_functionspace.F90",__LINE__,"test_pointcloud") + grid = atlas_Grid("O8") fs = atlas_functionspace_PointCloud(grid) @@ -624,8 +628,8 @@ module fcta_FunctionSpace_fxt field2 = fs%create_field(name="field2",kind=atlas_real(8),levels=3,variables=2) -FCTEST_CHECK_EQUAL( field%shape(), ([5,grid%size()]) ) -FCTEST_CHECK_EQUAL( field2%shape(), ([2,3,grid%size()]) ) +FCTEST_CHECK_EQUAL( field%shape(), ([5,int(grid%size())]) ) +FCTEST_CHECK_EQUAL( field2%shape(), ([2,3,int(grid%size())]) ) #ifndef _CRAYFTN FCTEST_CHECK_EQUAL( field%owners(), 1 ) @@ -643,6 +647,123 @@ module fcta_FunctionSpace_fxt #endif END_TEST +TEST( test_pointcloud_partition_remote ) +#if 1 +use fckit_module +implicit none + +type(atlas_functionspace_PointCloud) :: fs +type(atlas_Field) :: fld_points +type(atlas_Field) :: fld_ghost +type(atlas_Field) :: fld_partition +type(atlas_Field) :: fld_remote_index +type(atlas_Field) :: fld_values +type(atlas_Field) :: fld_values_save +type(atlas_FieldSet) :: fset + +type(atlas_functionspace) :: fs_base +type(atlas_trace) :: trace + +integer i, k +integer rank +real space +real(c_double), allocatable :: point_values(:,:) +integer(c_int), dimension(4) :: ghost_values +integer(c_int), dimension(4) :: partition_index +integer(ATLAS_KIND_IDX), dimension(4) :: remote_index +real(c_double), pointer :: field_values(:,:) +real(c_double), pointer :: field_values_save(:,:) +real(c_double), pointer :: field_values_new(:,:) +real(c_double), parameter :: tol = 1.e-12_dp + +trace = atlas_Trace("fctest_functionspace.F90",__LINE__,"test_pointcloud_partition_remote") +rank = fckit_mpi%rank() + +! the idea here is that we have a latitudinal circle around the equator +! that is partitioned equally +! +! test here is for 4 PEs (but could be more!) +! PE 0 PE 1 PE 2 PE 3 +! +! Owned values +! 0 45 90 135 180 225 270 315 +! +! Each PE has owned data followed by ghost data +! Ghost data is 1 before followed by 1 point after owned data ie. +! +! Correct halo_exchange for PE 0 +! 0 45 | 315 90 + + +space = 360.0 /(2.0 * fckit_mpi%size()) +fset = atlas_FieldSet() + +allocate(point_values(2, 4)) + +point_values = reshape((/space * (2 * rank), 0.0, & + space * (2 * rank + 1), 0.0, & + space * (2 * rank - 1), 0.0, & + space * (2 * rank + 2), 0.0/), shape(point_values)) + +if (point_values(1, 3) < 0.0) point_values(1, 3) = 360.0 + point_values(1, 3) +if (point_values(1, 4) > 360.0 - tol) point_values(1, 4) = 0.0d0 +fld_points = atlas_Field("lonlat", point_values(:, :)) +call fset%add(fld_points) + +ghost_values = (/ 0, 0, 1, 1 /) +fld_ghost = atlas_Field("ghost", ghost_values(:)) +call fset%add(fld_ghost) + +partition_index = (/ rank, rank, rank - 1, rank + 1 /) +if (partition_index(3) < 0) partition_index(3) = fckit_mpi%size() - 1 +if (partition_index(4) == fckit_mpi%size()) partition_index(4) = 0 +fld_partition = atlas_Field("partition", partition_index(:)) +call fset%add(fld_partition) + +remote_index = (/ 1, 2, 2, 1 /) +fld_remote_index = atlas_Field("remote_index", remote_index(:)) +call fset%add(fld_remote_index) + +fs = atlas_functionspace_PointCloud(fset) +fld_values = fs%create_field(name="values", kind=atlas_real(c_double), levels=1) +call fld_values%data(field_values) +fld_values_save = fs%create_field(name="values_save", kind=atlas_real(c_double), levels=1) +call fld_values_save%data(field_values_save) + +do i = 1, 4 + do k = 1, 1 + field_values(k, i) = point_values(1, i) + field_values_save(k, i) = point_values(1, i) + end do +end do + +call fs%halo_exchange(fld_values) +call fld_values%data(field_values_new) + +do i = 1, 4 + do k = 1, 1 + FCTEST_CHECK_CLOSE(field_values_new(k, i), field_values_save(k, i), tol) + end do +end do + +call fs%adjoint_halo_exchange(fld_values) + +do i = 1, 4 + do k = 1, 1 + if (i < 3) then + FCTEST_CHECK_CLOSE(field_values_new(k, i), 2.0 * field_values_save(k, i), tol) + else + FCTEST_CHECK_CLOSE(field_values_new(k, i), 0.0_dp, tol) + end if +end do +end do + +#else +#warning test_pointcloud_partition_remote disabled +#endif +END_TEST + + ! ----------------------------------------------------------------------------- diff --git a/src/tests/functionspace/test_pointcloud.cc b/src/tests/functionspace/test_pointcloud.cc index 336be54ba..4115f593c 100644 --- a/src/tests/functionspace/test_pointcloud.cc +++ b/src/tests/functionspace/test_pointcloud.cc @@ -9,8 +9,10 @@ */ #include "atlas/array.h" +#include "atlas/field.h" #include "atlas/functionspace/PointCloud.h" #include "atlas/option.h" +#include "atlas/parallel/mpi/mpi.h" #include "tests/AtlasTestEnvironment.h" @@ -155,6 +157,177 @@ CASE("test_createField") { EXPECT_EQ(f3.shape(2), 5); } +//----------------------------------------------------------------------------- + +double innerproductwithhalo(const atlas::Field& f1, const atlas::Field& f2) { + long sum(0); + + auto view1 = atlas::array::make_view(f1); + auto view2 = atlas::array::make_view(f2); + + for (atlas::idx_t jn = 0; jn < f1.shape(0); ++jn) { + for (atlas::idx_t jl = 0; jl < f1.levels(); ++jl) { + sum += view1(jn, jl) * view2(jn, jl); + } + } + + atlas::mpi::comm().allReduceInPlace(sum, eckit::mpi::sum()); + return sum; +} + + +CASE("test_createFieldSet") { + // Here is some ascii art to describe the test. + // Remote index is the same for both PEs in this case + // + // _1____0__ + // 1| 0 1 | 0 + // 3| 2 3 | 2 + // --------- + // 3 2 + // + // Point order (locally is also the same) + // + // _4____5__ + // 11| 0 1 | 6 + // 10| 2 3 | 7 + // --------- + // 9 8 + // + // Partition index + // PE 0: PE 1: + // + // _1____1__ _0____0__ + // 1| 0 0 | 1 0| 1 1 | 0 + // 1| 0 0 | 1 0| 1 1 | 0 + // --------- --------- + // 1 1 0 0 + // + // Initial Values + // PE 0: PE 1: + // + // _0____0__ _0____0__ + // 0|10 11 | 0 0|20 21 | 0 + // 0|12 13 | 0 0|22 23 | 0 + // --------- --------- + // 0 0 0 0 + // + // Values after halo exchange + // PE 0: PE 1: + // + // 21___20__ 11___10__ + // 21|10 11 | 20 11|20 21 | 10 + // 23|12 13 | 22 13|22 23 | 12 + // --------- --------- + // 23 22 13 12 + // + // Values after halo exchange and adjoint halo exchange + // PE 0: PE 1: + // + // _0____0__ _0___0__ + // 0|30 33 | 0 0|60 63 | 0 + // 0|36 39 | 0 0|66 69 | 0 + // --------- --------- + // 0 0 0 0 + + double tolerance(1e-16); + Field lonlat("lonlat", array::make_datatype(), array::make_shape(12, 2)); + Field ghost("ghost", array::make_datatype(), array::make_shape(12)); + Field remote_index("remote_index", array::make_datatype(), array::make_shape(12)); + Field partition("partition", array::make_datatype(), array::make_shape(12)); + + auto lonlatv = array::make_view(lonlat); + auto ghostv = array::make_view(ghost); + auto remote_indexv = array::make_indexview(remote_index); + auto partitionv = array::make_view(partition); + + ghostv.assign({ 0, 0, 0, 0, + 1, 1, 1, 1, 1, 1, 1, 1}); + + remote_indexv.assign({0, 1, 2, 3, 1, 0, 0, 2, 2, 3, 3, 1}); + + + if (atlas::mpi::rank() == 0) { + // center followed by clockwise halo starting from top left + lonlatv.assign({-45.0, 45.0, + 45.0, 45.0, + -45.0, -45.0, + -45.0, 45.0, // center + 225.0, 45.0, 135.0, 45.0, // up + 135.0, 45.0, 135.0, -45.0, //left + 135.0, -45.0, 135.0, 45.0, // down + 135.0, 45.0, 225.0, 45.0}); // right + + partitionv.assign({0, 0, 0, 0, + 1, 1, 1, 1, 1, 1, 1, 1}); + + + + } else if (atlas::mpi::rank() == 1) { + // center followed by clockwise halo starting from top left + lonlatv.assign({135.0, 45.0, 225.0, 45.0, 135.0, -45.0, 135.0, 45.0, // center + 45.0, 45.0, -45.0, 45.0, // up + -45.0, 45.0, -45.0, -45.0, // left + -45.0, -45.0, -45.0, 45.0, // down + -45.0, 45.0, 45.0, 45.0}); // right + + partitionv.assign({1, 1, 1, 1, + 0, 0, 0, 0, 0, 0, 0, 0}); + + } + + atlas::FieldSet fset; + fset.add(lonlat); + fset.add(ghost); + fset.add(remote_index); + fset.add(partition); + + auto fs2 = functionspace::PointCloud(fset); + Field f1 = fs2.createField(option::name("f1") | option::levels(2)); + Field f2 = fs2.createField(option::name("f2") | option::levels(2)); + auto f1v = array::make_view(f1); + auto f2v = array::make_view(f2); + + f1v.assign(0.0); + f2v.assign(0.0); + for (idx_t i = 0; i < f2v.shape(0); ++i) { + for (idx_t l = 0; l < f2v.shape(1); ++l) { + auto ghostv2 = array::make_view(fs2.ghost()); + if (ghostv2(i) == 0) { + f1v(i, l) = (atlas::mpi::rank() +1) * 10.0 + i; + f2v(i, l) = f1v(i, l); + } + } + } + + f2.haloExchange(); + + // adjoint test + double sum1 = innerproductwithhalo(f2, f2); + + f2.adjointHaloExchange(); + + double sum2 = innerproductwithhalo(f1, f2); + + atlas::mpi::comm().allReduceInPlace(sum1, eckit::mpi::sum()); + atlas::mpi::comm().allReduceInPlace(sum2, eckit::mpi::sum()); + EXPECT(std::abs(sum1 - sum2)/ std::abs(sum1) < tolerance); + atlas::Log::info() << "adjoint test passed :: " + << "sum1 " << sum1 << " sum2 " << sum2 << " normalised difference " + << std::abs(sum1 - sum2)/ std::abs(sum1) << std::endl; + + // In this case the effect of the halo exchange followed by the adjoint halo exchange + // multiples the values by a factor of 3 (see pictures above) + for (idx_t i = 0; i < f2v.shape(0); ++i) { + for (idx_t l = 0; l < f2v.shape(1); ++l) { + EXPECT( std::abs(f2v(i, l) - 3.0 * f1v(i, l)) < tolerance); + } + } + atlas::Log::info() << "values from halo followed by halo adjoint are as expected " + << std::endl; +} + + //----------------------------------------------------------------------------- diff --git a/src/tests/interpolation/test_interpolation_unstructured_bilinear_lonlat.cc b/src/tests/interpolation/test_interpolation_unstructured_bilinear_lonlat.cc index dcd954e3b..83a1c2ee5 100644 --- a/src/tests/interpolation/test_interpolation_unstructured_bilinear_lonlat.cc +++ b/src/tests/interpolation/test_interpolation_unstructured_bilinear_lonlat.cc @@ -32,6 +32,27 @@ namespace test { //----------------------------------------------------------------------------- // +CASE("test_interpolation_O64_to_empty_PointCloud") { + Grid grid("O64"); + Mesh mesh(grid); + NodeColumns fs(mesh); + + atlas::Field points("lonlat", atlas::array::make_datatype(), atlas::array::make_shape(0, 2)); + // Some points at the equator + PointCloud pointcloud(points); + + Interpolation interpolation(option::type("unstructured-bilinear-lonlat"), fs, pointcloud); + + Field field_source = fs.createField(option::name("source")); + Field field_target("target", array::make_datatype(), array::make_shape(pointcloud.size())); + + auto source = array::make_view(field_source); + for (idx_t j = 0; j < fs.nodes().size(); ++j) { + source(j) = 0; + } + + interpolation.execute(field_source, field_target); +} CASE("test_interpolation_O64_to_points_bilinear_remapping") { Grid grid("O64"); diff --git a/src/tests/trans/test_transgeneral.cc b/src/tests/trans/test_transgeneral.cc index 4213c52e8..741575ea7 100644 --- a/src/tests/trans/test_transgeneral.cc +++ b/src/tests/trans/test_transgeneral.cc @@ -1582,7 +1582,7 @@ CASE("test_trans_levels") { #endif -#if ATLAS_HAVE_TRANS && (ATLAS_HAVE_ECTRANS || defined(TRANS_HAVE_INVTRANS_ADJ)) +#if ATLAS_HAVE_ECTRANS || defined(TRANS_HAVE_INVTRANS_ADJ) CASE("test_2level_adjoint_test_with_powerspectrum_convolution") { std::string grid_uid("F64"); // Regular Gaussian F ( 8 N^2) atlas::StructuredGrid g2(grid_uid); @@ -1597,15 +1597,11 @@ CASE("test_2level_adjoint_test_with_powerspectrum_convolution") { atlas::functionspace::Spectral specFS(2 * N - 1, atlas::option::levels(levels)); atlas::functionspace::StructuredColumns gridFS( - g2, atlas::grid::Partitioner(new atlas::grid::detail::partitioner::TransPartitioner()), - atlas::option::levels(levels)); - - auto lonlatview = atlas::array::make_view(gridFS.lonlat()); + g2, atlas::grid::Partitioner(new atlas::grid::detail::partitioner::TransPartitioner()), + atlas::option::levels(levels)); atlas::trans::Trans transIFS(gridFS, specFS); - // Log::info() << "transIFS llatlon = " << dynamic_cast(transIFS.get())->trans()->llatlon << std::endl; - Log::info() << "transIFS backend" << transIFS.backend() << std::endl; std::vector powerSpectrum(2 * N, 0.0); @@ -1617,13 +1613,15 @@ CASE("test_2level_adjoint_test_with_powerspectrum_convolution") { for (std::size_t w = 0; w < powerSpectrum.size(); ++w) { powerSpectrum[w] = powerSpectrum[w] / tot; } - Log::info() << "create a fictitous power spectrum" << atlas::mpi::rank() << " " << powerSpectrum[0] << " " + Log::info() << "create a fictitous power spectrum" << atlas::mpi::rank() + << " " << powerSpectrum[0] << " " << powerSpectrum[1] << " " << tot << std::endl; atlas::Field gpf = gridFS.createField(atlas::option::name("gpf")); atlas::Field gpf2 = gridFS.createField(atlas::option::name("gpf2")); atlas::Field spf = specFS.createField(atlas::option::name("spf")); - atlas::Field spfg = specFS.createField(atlas::option::name("spfg") | atlas::option::global()); + atlas::Field spfdiv = specFS.createField(atlas::option::name("spfdv")); + auto gpfView = atlas::array::make_view(gpf); for (atlas::idx_t j = gridFS.j_begin(); j < gridFS.j_end(); ++j) { @@ -1631,22 +1629,150 @@ CASE("test_2level_adjoint_test_with_powerspectrum_convolution") { atlas::idx_t jn = gridFS.index(i, j); for (atlas::idx_t jl : {0, 1}) { gpfView(jn, jl) = 0.0; - if ((j == gridFS.j_end() - 1) && (i == gridFS.i_end(gridFS.j_end() - 1) - 1) && (jl == 0) && - (atlas::mpi::rank() == 0)) { + if ((j == gridFS.j_end() - 1) && + (i == gridFS.i_end(gridFS.j_end() - 1) - 1) && (jl == 0) && + (atlas::mpi::rank() == 0)) { gpfView(jn, jl) = 1.0; } - if ((j == gridFS.j_end() - 1) && (i == gridFS.i_end(gridFS.j_end() - 1) - 1) && (jl == 1) && - (atlas::mpi::rank() == 1)) { + if ((j == gridFS.j_end() - 1) && + (i == gridFS.i_end(gridFS.j_end() - 1) - 1) && (jl == 1) && + (atlas::mpi::rank() == 1)) { gpfView(jn, jl) = 1.0; } } } } - // transform fields to spectral and view - transIFS.invtrans_adj(gpf, spf); + std::vector test_name{"inverse", "direct"}; + + for (std::size_t test_type = 0; test_type < test_name.size(); ++test_type) { + auto spfView = atlas::array::make_view(spf); + auto gpf2View = atlas::array::make_view(gpf2); + + spfView.assign(0.0); + gpf2View.assign(0.0); + + // transform fields to spectral and view + if (test_name[test_type].compare("inverse") == 0) { + transIFS.invtrans_adj(gpf, spf); + } else if (test_name[test_type].compare("direct") == 0) { + transIFS.dirtrans(gpf, spf); + } + + const auto zonal_wavenumbers = specFS.zonal_wavenumbers(); + const int nb_zonal_wavenumbers = zonal_wavenumbers.size(); + double adj_value(0.0); + int i{0}; + for (int jm = 0; jm < nb_zonal_wavenumbers; ++jm) { + const std::size_t m1 = zonal_wavenumbers(jm); + for (std::size_t n1 = m1; n1 <= static_cast(2 * N - 1); ++n1) { + for (int imag1 : {0, 1}) { + for (int jl : {0, 1}) { + // scale by the square root of the power spectrum + // note here that we need the modal power spectrum + // i.e. divide powerSpectra by (2n +1) + spfView(i, jl) = spfView(i, jl) * + std::sqrt(std::sqrt(static_cast(powerSpectrum[n1]))) / + std::sqrt(static_cast(2 * n1 + 1)); + + // adjoint at the heart + double temp = spfView(i, jl) * spfView(i, jl); + adj_value += (m1 > 0 ? 2 * temp : temp); // to take account of -ve m. + + // scale by the square root of the power spectrum (again!) + spfView(i, jl) = spfView(i, jl) * + std::sqrt(std::sqrt(static_cast(powerSpectrum[n1]))) / + std::sqrt(static_cast(2 * n1 + 1)); + } + ++i; + } + } + } + atlas::mpi::comm().allReduceInPlace(adj_value, eckit::mpi::sum()); + + // transform fields to spectral and view + if (test_name[test_type].compare("inverse") == 0) { + transIFS.invtrans(spf, gpf2); + } else if (test_name[test_type].compare("direct") == 0) { + transIFS.dirtrans_adj(spf, gpf2); + } + Log::info() << "adjoint test transforms " << test_name[test_type] << std::endl; + + + double adj_value2(0.0); + for (atlas::idx_t j = gridFS.j_begin(); j < gridFS.j_end(); ++j) { + for (atlas::idx_t i = gridFS.i_begin(j); i < gridFS.i_end(j); ++i) { + atlas::idx_t jn = gridFS.index(i, j); + for (atlas::idx_t jl : {0, 1}) { + adj_value2 += gpfView(jn, jl) * gpf2View(jn, jl); + } + } + } + atlas::mpi::comm().allReduceInPlace(adj_value2, eckit::mpi::sum()); + + Log::info() << "adjoint test " << test_name[test_type] + << " " << adj_value << " " << adj_value2 << std::endl; + EXPECT(std::abs(adj_value - adj_value2) / adj_value < 1e-12); + } +} + +CASE("test_2level_adjoint_test_with_vortdiv") { + std::string grid_uid("F64"); // Regular Gaussian F ( 8 N^2) + atlas::StructuredGrid g2(grid_uid); + + auto N = atlas::GaussianGrid(g2).N(); // -> cast to Gaussian grid and get the Gaussian number + auto levels = 2; + std::vector imag{false, true}; // imaginary component or not + + atlas::functionspace::Spectral specFS(2 * N - 1, atlas::option::levels(levels)); + + atlas::functionspace::StructuredColumns gridFS( + g2, atlas::grid::Partitioner(new atlas::grid::detail::partitioner::TransPartitioner()), + atlas::option::levels(levels)); + + atlas::trans::Trans transIFS(gridFS, specFS); + + atlas::Field spfvor = specFS.createField(atlas::option::name("spfvor")); + atlas::Field spfdiv = specFS.createField(atlas::option::name("spfdiv")); + + atlas::Field gpfuv = gridFS.createField(atlas::option::name("gpfuv") | + atlas::option::variables(2)); + atlas::Field gpfuv2 = gridFS.createField(atlas::option::name("gpfuv") | + atlas::option::variables(2)); + + + auto gpfuvView = atlas::array::make_view(gpfuv); + for (atlas::idx_t j = gridFS.j_begin(); j < gridFS.j_end(); ++j) { + for (atlas::idx_t i = gridFS.i_begin(j); i < gridFS.i_end(j); ++i) { + atlas::idx_t jn = gridFS.index(i, j); + for (atlas::idx_t jl : {0, 1}) { + gpfuvView(jn, jl, 0) = 0.0; + gpfuvView(jn, jl, 1) = 0.0; + if ((j == gridFS.j_end() - 1) && + (i == gridFS.i_end(gridFS.j_end() - 1) - 1) && (jl == 0) && + (atlas::mpi::rank() == 0)) { + gpfuvView(jn, jl, 0) = 10.0; + gpfuvView(jn, jl, 1) = 1.0; + } + if ((j == gridFS.j_end() - 1) && + (i == gridFS.i_end(gridFS.j_end() - 1) - 1) && (jl == 1) && + (atlas::mpi::rank() == 1)) { + gpfuvView(jn, jl, 0) = 10.0; + gpfuvView(jn, jl, 1) = 1.0; + } + } + } + } + + auto spfdivView = atlas::array::make_view(spfdiv); + auto spfvorView = atlas::array::make_view(spfvor); + auto gpfuv2View = atlas::array::make_view(gpfuv2); + spfdivView.assign(0.0); + spfvorView.assign(0.0); + gpfuv2View.assign(0.0); + + transIFS.dirtrans_wind2vordiv(gpfuv, spfvor, spfdiv); - auto spfView = atlas::array::make_view(spf); const auto zonal_wavenumbers = specFS.zonal_wavenumbers(); const int nb_zonal_wavenumbers = zonal_wavenumbers.size(); double adj_value(0.0); @@ -1656,19 +1782,10 @@ CASE("test_2level_adjoint_test_with_powerspectrum_convolution") { for (std::size_t n1 = m1; n1 <= static_cast(2 * N - 1); ++n1) { for (int imag1 : {0, 1}) { for (int jl : {0, 1}) { - // scale by the square root of the power spectrum - // note here that we need the modal power spectrum - // i.e. divide powerSpectra by (2n +1) - spfView(i, jl) = spfView(i, jl) * std::sqrt(std::sqrt(static_cast(powerSpectrum[n1]))) / - std::sqrt(static_cast(2 * n1 + 1)); - // adjoint at the heart - double temp = spfView(i, jl) * spfView(i, jl); + double temp = spfdivView(i, jl) * spfdivView(i, jl) + + spfvorView(i, jl) * spfvorView(i, jl); adj_value += (m1 > 0 ? 2 * temp : temp); // to take account of -ve m. - - // scale by the square root of the power spectrum (again!) - spfView(i, jl) = spfView(i, jl) * std::sqrt(std::sqrt(static_cast(powerSpectrum[n1]))) / - std::sqrt(static_cast(2 * n1 + 1)); } ++i; } @@ -1676,22 +1793,21 @@ CASE("test_2level_adjoint_test_with_powerspectrum_convolution") { } atlas::mpi::comm().allReduceInPlace(adj_value, eckit::mpi::sum()); - transIFS.invtrans(spf, gpf2); + transIFS.dirtrans_wind2vordiv_adj(spfvor, spfdiv, gpfuv2); - Log::info() << "adjoint test transforms " << std::endl; - auto gpf2View = atlas::array::make_view(gpf2); double adj_value2(0.0); for (atlas::idx_t j = gridFS.j_begin(); j < gridFS.j_end(); ++j) { for (atlas::idx_t i = gridFS.i_begin(j); i < gridFS.i_end(j); ++i) { atlas::idx_t jn = gridFS.index(i, j); for (atlas::idx_t jl : {0, 1}) { - adj_value2 += gpfView(jn, jl) * gpf2View(jn, jl); + adj_value2 += gpfuv2View(jn, jl, 0) * gpfuvView(jn, jl, 0) + + gpfuv2View(jn, jl, 1) * gpfuvView(jn, jl, 1); } } } atlas::mpi::comm().allReduceInPlace(adj_value2, eckit::mpi::sum()); - Log::info() << "adjoint test " + Log::info() << "adjoint test for wind2vordiv" << " " << adj_value << " " << adj_value2 << std::endl; EXPECT(std::abs(adj_value - adj_value2) / adj_value < 1e-12); } diff --git a/tools/apply-clang-format.sh b/tools/apply-clang-format.sh index ed44034ba..521ccb2a3 100755 --- a/tools/apply-clang-format.sh +++ b/tools/apply-clang-format.sh @@ -105,9 +105,9 @@ if [[ $all =~ "yes" ]]; then cd $SCRIPTDIR/../atlas_io if [[ $dryrun =~ "yes" ]]; then - echo "+ find . -iname *.h -o -iname *.cc | xargs clang-format -i -style=file" + echo "+ find . -iname \"*.h\" -o -iname \"*.cc\" | xargs clang-format -i -style=file" else - find . -iname *.h -o -iname *.cc | xargs clang-format -i -style=file + find . -iname "*.h" -o -iname "*.cc" | xargs clang-format -i -style=file fi else @@ -120,10 +120,10 @@ else fi cd $1 if [[ $dryrun =~ "yes" ]]; then - echo "+ find . -iname *.h -o -iname *.cc | xargs clang-format -i -style=file" - find . -iname *.h -o -iname *.cc | xargs echo + echo "+ find . -iname \"*.h\" -o -iname \"*.cc\" | xargs clang-format -i -style=file" + find . -iname "*.h" -o -iname "*.cc" | xargs echo else - find . -iname *.h -o -iname *.cc | xargs clang-format -i -style=file + find . -iname "*.h" -o -iname "*.cc" | xargs clang-format -i -style=file fi shift else