diff --git a/.github/docs-config.yml b/.github/docs-config.yml new file mode 100644 index 000000000..e08272a73 --- /dev/null +++ b/.github/docs-config.yml @@ -0,0 +1,14 @@ +# note: each step is executed in own process +build-steps: + - git clone --depth 1 https://github.com/ecmwf/atlas-docs.git $RUNNER_TEMP/atlas-docs + - sudo apt install -y -q doxygen texlive-full + - | + cd $RUNNER_TEMP/atlas-docs + make PUBLIC=1 WITH_ECKIT=1 WITH_DOXYGEN=1 ATLAS_SOURCE_DIR=$GITHUB_WORKSPACE clean html + echo "DOC_BUILD_PATH=$RUNNER_TEMP/atlas-docs/build/html" >> "$GITHUB_ENV" + +hosts: + ecmwf-sites: + space: docs + name: atlas + diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml new file mode 100644 index 000000000..1c1df2a90 --- /dev/null +++ b/.github/workflows/docs.yml @@ -0,0 +1,15 @@ +name: docs + +on: + push: + tags: + - '**' + + workflow_dispatch: ~ + +jobs: + publish: + uses: ecmwf-actions/reusable-workflows/.github/workflows/cd-docs.yml@v2 + secrets: inherit + with: + config: .github/docs-config.yml diff --git a/CHANGELOG.md b/CHANGELOG.md index 55f29a917..f7a0e8ece 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,19 @@ This project adheres to [Semantic Versioning](http://semver.org/spec/v2.0.0.html ## [Unreleased] +## [0.37.0] - 2024-04-09 +### Added +- Add SphericalVector interpolation method using parallel transport (#163) +- Added arrayForEachDim method + +### Fixed +- Bugfix for regional grids with ny > nx +- Fix build for configuration setting ATLAS_BITS_LOCAL=64 + +### Changed +- Use new LocalConfiguration baseclass functions in util::Config and util::Metadata instead of eckit::Value backdoor +- atlas_io is an adaptor library when eckit_codec is available (#181) + ## [0.36.0] - 2023-12-11 ### Added - Add TriangularMeshBuilder with Fortran API, so far for serial meshes only @@ -514,6 +527,7 @@ Fix StructuredInterpolation2D with retry for failed stencils ## 0.13.0 - 2018-02-16 [Unreleased]: https://github.com/ecmwf/atlas/compare/master...develop +[0.37.0]: https://github.com/ecmwf/atlas/compare/0.36.0...0.37.0 [0.36.0]: https://github.com/ecmwf/atlas/compare/0.35.1...0.36.0 [0.35.1]: https://github.com/ecmwf/atlas/compare/0.35.0...0.35.1 [0.35.0]: https://github.com/ecmwf/atlas/compare/0.34.0...0.35.0 diff --git a/VERSION b/VERSION index 93d4c1ef0..0f1a7dfc7 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -0.36.0 +0.37.0 diff --git a/atlas_io/CMakeLists.txt b/atlas_io/CMakeLists.txt index 26f549808..2a821b761 100644 --- a/atlas_io/CMakeLists.txt +++ b/atlas_io/CMakeLists.txt @@ -28,6 +28,20 @@ ecbuild_debug( " eckit_FEATURES : [${eckit_FEATURES}]" ) ################################################################################ # Features that can be enabled / disabled with -DENABLE_ +ecbuild_add_option( FEATURE ECKIT_DEVELOP + DESCRIPTION "Used to enable new features or API depending on eckit develop branch, not yet in a tagged release" + DEFAULT OFF ) + +set( eckit_HAVE_ECKIT_CODEC 0 ) +if( TARGET eckit_codec ) + set( eckit_HAVE_ECKIT_CODEC 1 ) +endif() + +ecbuild_add_option( FEATURE ECKIT_CODEC + DEFAULT ON + DESCRIPTION "Use eckit::codec with adaptor" + CONDITION eckit_HAVE_ECKIT_CODEC ) + ################################################################################ # sources @@ -47,8 +61,12 @@ else() set( ATLAS_IO_LITTLE_ENDIAN 1 ) endif() - -add_subdirectory( src ) +if( HAVE_ECKIT_CODEC ) + ecbuild_info("atlas_io is configured to be an adaptor library which delegates calls to eckit_codec") + add_subdirectory(eckit_codec_adaptor) +else() + add_subdirectory( src ) +endif() add_subdirectory( tests ) ################################################################################ diff --git a/atlas_io/eckit_codec_adaptor/CMakeLists.txt b/atlas_io/eckit_codec_adaptor/CMakeLists.txt new file mode 100644 index 000000000..febd4f0ab --- /dev/null +++ b/atlas_io/eckit_codec_adaptor/CMakeLists.txt @@ -0,0 +1 @@ +add_subdirectory(src) diff --git a/atlas_io/eckit_codec_adaptor/src/CMakeLists.txt b/atlas_io/eckit_codec_adaptor/src/CMakeLists.txt new file mode 100644 index 000000000..e91bc52fd --- /dev/null +++ b/atlas_io/eckit_codec_adaptor/src/CMakeLists.txt @@ -0,0 +1,2 @@ +add_subdirectory(atlas_io) + diff --git a/atlas_io/eckit_codec_adaptor/src/atlas_io/CMakeLists.txt b/atlas_io/eckit_codec_adaptor/src/atlas_io/CMakeLists.txt new file mode 100644 index 000000000..1794974b1 --- /dev/null +++ b/atlas_io/eckit_codec_adaptor/src/atlas_io/CMakeLists.txt @@ -0,0 +1,19 @@ + +ecbuild_add_library( TARGET atlas_io + + INSTALL_HEADERS ALL + HEADER_DESTINATION include/atlas_io + PUBLIC_LIBS eckit_codec + PUBLIC_INCLUDES + $ + + SOURCES + atlas-io.h + Trace.cc + Trace.h + detail/BlackMagic.h + LINKER_LANGUAGE CXX +) + +target_compile_features( atlas_io PUBLIC cxx_std_17 ) + diff --git a/atlas_io/eckit_codec_adaptor/src/atlas_io/Trace.cc b/atlas_io/eckit_codec_adaptor/src/atlas_io/Trace.cc new file mode 100644 index 000000000..5e52b7dc9 --- /dev/null +++ b/atlas_io/eckit_codec_adaptor/src/atlas_io/Trace.cc @@ -0,0 +1,43 @@ +/* + * (C) Copyright 2020 ECMWF. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + * In applying this licence, ECMWF does not waive the privileges and immunities + * granted to it by virtue of its status as an intergovernmental organisation + * nor does it submit to any jurisdiction. + */ + +#include "Trace.h" + +#include "eckit/log/CodeLocation.h" + +namespace atlas { +namespace io { + +atlas::io::Trace::Trace(const eckit::CodeLocation& loc) { + 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 (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 (size_t id = 0; id < TraceHookRegistry::size(); ++id) { + if (TraceHookRegistry::enabled(id)) { + hooks_.emplace_back(TraceHookRegistry::hook(id)(loc, title)); + } + } +} + +} // namespace io +} // namespace atlas diff --git a/atlas_io/eckit_codec_adaptor/src/atlas_io/Trace.h b/atlas_io/eckit_codec_adaptor/src/atlas_io/Trace.h new file mode 100644 index 000000000..43da3f526 --- /dev/null +++ b/atlas_io/eckit_codec_adaptor/src/atlas_io/Trace.h @@ -0,0 +1,77 @@ +/* + * (C) Copyright 2020 ECMWF. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + * In applying this licence, ECMWF does not waive the privileges and immunities + * granted to it by virtue of its status as an intergovernmental organisation + * nor does it submit to any jurisdiction. + */ + +#pragma once + + +#include +#include +#include +#include +#include + +namespace eckit { +class CodeLocation; +} + +namespace atlas { +namespace io { + +struct TraceHook { + TraceHook() = default; + virtual ~TraceHook() = default; +}; + +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 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]; } + static size_t invalidId() { return std::numeric_limits::max(); } + +private: + TraceHookRegistry() = default; +}; + +struct Trace { + using Labels = std::vector; + Trace(const eckit::CodeLocation& loc); + Trace(const eckit::CodeLocation& loc, const std::string& title); + Trace(const eckit::CodeLocation& loc, const std::string& title, const Labels& labels); + +private: + std::vector> hooks_; +}; + +} // namespace io +} // namespace atlas + +#include "atlas_io/detail/BlackMagic.h" + +#define ATLAS_IO_TRACE(...) __ATLAS_IO_TYPE(::atlas::io::Trace, Here() __ATLAS_IO_COMMA_ARGS(__VA_ARGS__)) +#define ATLAS_IO_TRACE_SCOPE(...) __ATLAS_IO_TYPE_SCOPE(::atlas::io::Trace, Here() __ATLAS_IO_COMMA_ARGS(__VA_ARGS__)) diff --git a/atlas_io/eckit_codec_adaptor/src/atlas_io/atlas-io.h b/atlas_io/eckit_codec_adaptor/src/atlas_io/atlas-io.h new file mode 100644 index 000000000..7c5f49d20 --- /dev/null +++ b/atlas_io/eckit_codec_adaptor/src/atlas_io/atlas-io.h @@ -0,0 +1,99 @@ +/* + * (C) Copyright 2023- ECMWF. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + * In applying this licence, ECMWF does not waive the privileges and immunities + * granted to it by virtue of its status as an intergovernmental organisation + * nor does it submit to any jurisdiction. + */ + +#pragma once + +#include "eckit/codec/codec.h" + +namespace atlas::io { + + // Encoding/Decoding + using Metadata = ::eckit::codec::Metadata; + using DataType = ::eckit::codec::DataType; + using Data = ::eckit::codec::Data; + using Encoder = ::eckit::codec::Encoder; + using Decoder = ::eckit::codec::Decoder; + + // Record + using Record = ::eckit::codec::Record; + using RecordWriter = ::eckit::codec::RecordWriter; + using RecordPrinter = ::eckit::codec::RecordPrinter; + using RecordItemReader = ::eckit::codec::RecordItemReader; + using RecordReader = ::eckit::codec::RecordReader; + + // I/O + using Session = ::eckit::codec::Session; + using Mode = ::eckit::codec::Mode; + using Stream = ::eckit::codec::Stream; + using FileStream = ::eckit::codec::FileStream; + using InputFileStream = ::eckit::codec::InputFileStream; + using OutputFileStream = ::eckit::codec::OutputFileStream; + + // Array + using ArrayReference = ::eckit::codec::ArrayReference; + using ArrayMetadata = ::eckit::codec::ArrayMetadata; + using ArrayShape = ::eckit::codec::ArrayShape; + + // Exceptions + using Exception = ::eckit::codec::Exception; + using NotEncodable = ::eckit::codec::NotEncodable; + using NotDecodable = ::eckit::codec::NotDecodable; + using InvalidRecord = ::eckit::codec::InvalidRecord; + using DataCorruption = ::eckit::codec::DataCorruption; + using WriteError = ::eckit::codec::WriteError; + + // Type traits + template + static constexpr bool is_interpretable() { + return ::eckit::codec::is_interpretable(); + } + template + static constexpr bool is_encodable() { + return ::eckit::codec::is_encodable(); + } + template + static constexpr bool is_decodable() { + return ::eckit::codec::is_decodable(); + } + template + static constexpr bool can_encode_metadata() { + return ::eckit::codec::can_encode_metadata(); + } + template + static constexpr bool can_encode_data() { + return ::eckit::codec::can_encode_metadata(); + } + + namespace tag { + using disable_static_assert = ::eckit::codec::tag::disable_static_assert; + using enable_static_assert = ::eckit::codec::tag::enable_static_assert; + } + + // Functions + using ::eckit::codec::ref; + using ::eckit::codec::copy; + using ::eckit::codec::encode; + using ::eckit::codec::decode; + using ::eckit::codec::interprete; + using ::eckit::codec::link; + using ::eckit::codec::make_datatype; + // template + // using make_datatype = eckit::codec::make_datatype; + + namespace defaults { + using ::eckit::codec::defaults::compression_algorithm; + using ::eckit::codec::defaults::checksum_algorithm; + using ::eckit::codec::defaults::checksum_read; + using ::eckit::codec::defaults::checksum_write; + } +} + +#define ATLAS_IO_ASSERT(X) ASSERT(X) +#define ATLAS_IO_ASSERT_MSG(X, M) ASSERT_MSG(X, M) diff --git a/atlas_io/eckit_codec_adaptor/src/atlas_io/detail/BlackMagic.h b/atlas_io/eckit_codec_adaptor/src/atlas_io/detail/BlackMagic.h new file mode 100644 index 000000000..a889c5e19 --- /dev/null +++ b/atlas_io/eckit_codec_adaptor/src/atlas_io/detail/BlackMagic.h @@ -0,0 +1,101 @@ +/* + * (C) Copyright 2013 ECMWF. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + * In applying this licence, ECMWF does not waive the privileges and immunities + * granted to it by virtue of its status as an intergovernmental organisation + * nor does it submit to any jurisdiction. + */ + +#pragma once + +// This file contains preprocessor black magic. It contains macros that +// can return the number of arguments passed + +//----------------------------------------------------------------------------------------------------------- +// Public + +/// Returns the number of passed arguments +#define __ATLAS_IO_NARG(...) + +/// Splice a and b together +#define __ATLAS_IO_SPLICE(a, b) + +#define __ATLAS_IO_STRINGIFY(a) a + +#define __ATLAS_IO_TYPE(Type, ...) +#define __ATLAS_IO_TYPE_SCOPE(Type, ...) + +//----------------------------------------------------------------------------------------------------------- +// Details + +// Undefine these, to be redefined further down. +#undef __ATLAS_IO_NARG +#undef __ATLAS_IO_SPLICE +#undef __ATLAS_IO_TYPE +#undef __ATLAS_IO_TYPE_SCOPE + +#define __ATLAS_IO_REVERSE 5, 4, 3, 2, 1, 0 +#define __ATLAS_IO_ARGN(_1, _2, _3, _4, _5, N, ...) N +#define __ATLAS_IO_NARG__(dummy, ...) __ATLAS_IO_ARGN(__VA_ARGS__) +#define __ATLAS_IO_NARG_(...) __ATLAS_IO_NARG__(dummy, ##__VA_ARGS__, __ATLAS_IO_REVERSE) +#define __ATLAS_IO_SPLICE(a, b) __ATLAS_IO_SPLICE_1(a, b) +#define __ATLAS_IO_SPLICE_1(a, b) __ATLAS_IO_SPLICE_2(a, b) +#define __ATLAS_IO_SPLICE_2(a, b) a##b + +#define __ATLAS_IO_ARG16(_0, _1, _2, _3, _4, _5, _6, _7, _8, _9, _10, _11, _12, _13, _14, _15, ...) _15 +#define __ATLAS_IO_HAS_COMMA(...) __ATLAS_IO_ARG16(__VA_ARGS__, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0) +#define __ATLAS_IO_TRIGGER_PARENTHESIS(...) , +#define __ATLAS_IO_ISEMPTY(...) \ + __ATLAS_IO_ISEMPTY_(/* test if there is just one argument, eventually an empty \ + one */ \ + __ATLAS_IO_HAS_COMMA(__VA_ARGS__), /* test if \ + _TRIGGER_PARENTHESIS_ \ + together with the \ + argument adds a comma */ \ + __ATLAS_IO_HAS_COMMA( \ + __ATLAS_IO_TRIGGER_PARENTHESIS __VA_ARGS__), /* test if the argument together with \ + a parenthesis adds a comma */ \ + __ATLAS_IO_HAS_COMMA(__VA_ARGS__(/*empty*/)), /* test if placing it between \ + __ATLAS_IO_TRIGGER_PARENTHESIS and the \ + parenthesis adds a comma */ \ + __ATLAS_IO_HAS_COMMA(__ATLAS_IO_TRIGGER_PARENTHESIS __VA_ARGS__(/*empty*/))) + +#define __ATLAS_IO_PASTE5(_0, _1, _2, _3, _4) _0##_1##_2##_3##_4 +#define __ATLAS_IO_ISEMPTY_(_0, _1, _2, _3) \ + __ATLAS_IO_HAS_COMMA(__ATLAS_IO_PASTE5(__ATLAS_IO_IS_EMPTY_CASE_, _0, _1, _2, _3)) +#define __ATLAS_IO_IS_EMPTY_CASE_0001 , + +#define __ATLAS_IO_NARG(...) __ATLAS_IO_SPLICE(__ATLAS_IO_CALL_NARG_, __ATLAS_IO_ISEMPTY(__VA_ARGS__))(__VA_ARGS__) +#define __ATLAS_IO_CALL_NARG_1(...) 0 +#define __ATLAS_IO_CALL_NARG_0 __ATLAS_IO_NARG_ + +#define __ATLAS_IO_COMMA_ARGS(...) \ + __ATLAS_IO_SPLICE(__ATLAS_IO_COMMA_ARGS_, __ATLAS_IO_ISEMPTY(__VA_ARGS__))(__VA_ARGS__) +#define __ATLAS_IO_COMMA_ARGS_1(...) +#define __ATLAS_IO_COMMA_ARGS_0(...) , __VA_ARGS__ + +#define __ATLAS_IO_ARGS_OR_DUMMY(...) \ + __ATLAS_IO_SPLICE(__ATLAS_IO_ARGS_OR_DUMMY_, __ATLAS_IO_ISEMPTY(__VA_ARGS__)) \ + (__VA_ARGS__) +#define __ATLAS_IO_ARGS_OR_DUMMY_0(...) __VA_ARGS__ +#define __ATLAS_IO_ARGS_OR_DUMMY_1(...) 0 + +#define __ATLAS_IO_TYPE(Type, ...) \ + __ATLAS_IO_SPLICE(__ATLAS_IO_TYPE_, __ATLAS_IO_ISEMPTY(__VA_ARGS__)) \ + (Type, __ATLAS_IO_ARGS_OR_DUMMY(__VA_ARGS__)) +#define __ATLAS_IO_TYPE_1(Type, dummy) Type __ATLAS_IO_SPLICE(__variable_, __LINE__) +#define __ATLAS_IO_TYPE_0(Type, ...) Type __ATLAS_IO_SPLICE(__variable_, __LINE__)(__VA_ARGS__) + +#define __ATLAS_IO_TYPE_SCOPE(Type, ...) \ + __ATLAS_IO_SPLICE(__ATLAS_IO_TYPE_SCOPE_, __ATLAS_IO_ISEMPTY(__VA_ARGS__)) \ + (Type, __ATLAS_IO_ARGS_OR_DUMMY(__VA_ARGS__)) +#define __ATLAS_IO_TYPE_SCOPE_1(Type, ...) \ + for (bool __ATLAS_IO_SPLICE(__done_, __LINE__) = false; __ATLAS_IO_SPLICE(__done_, __LINE__) != true;) \ + for (Type __ATLAS_IO_SPLICE(__variable_, __LINE__); __ATLAS_IO_SPLICE(__done_, __LINE__) != true; \ + __ATLAS_IO_SPLICE(__done_, __LINE__) = true) +#define __ATLAS_IO_TYPE_SCOPE_0(Type, ...) \ + for (bool __ATLAS_IO_SPLICE(__done_, __LINE__) = false; __ATLAS_IO_SPLICE(__done_, __LINE__) != true;) \ + for (Type __ATLAS_IO_SPLICE(__variable_, __LINE__)(__VA_ARGS__); __ATLAS_IO_SPLICE(__done_, __LINE__) != true; \ + __ATLAS_IO_SPLICE(__done_, __LINE__) = true) diff --git a/atlas_io/src/atlas_io/CMakeLists.txt b/atlas_io/src/atlas_io/CMakeLists.txt index 877d61ead..6113287b2 100644 --- a/atlas_io/src/atlas_io/CMakeLists.txt +++ b/atlas_io/src/atlas_io/CMakeLists.txt @@ -1,3 +1,14 @@ + + +if( HAVE_ECKIT_DEVELOP ) + set( ATLAS_IO_ECKIT_DEVELOP 1 ) +else() + set( ATLAS_IO_ECKIT_DEVELOP 0 ) +endif() + +ecbuild_parse_version( ${eckit_VERSION} PREFIX ATLAS_IO_ECKIT ) +math( EXPR ATLAS_IO_ECKIT_VERSION_INT "( 10000 * ${ATLAS_IO_ECKIT_VERSION_MAJOR} ) + ( 100 * ${ATLAS_IO_ECKIT_VERSION_MINOR} ) + ${ATLAS_IO_ECKIT_VERSION_PATCH}" ) + configure_file( detail/defines.h.in detail/defines.h ) install( FILES ${CMAKE_CURRENT_BINARY_DIR}/detail/defines.h diff --git a/atlas_io/src/atlas_io/Metadata.h b/atlas_io/src/atlas_io/Metadata.h index b6f2ab123..d6689ba21 100644 --- a/atlas_io/src/atlas_io/Metadata.h +++ b/atlas_io/src/atlas_io/Metadata.h @@ -57,7 +57,11 @@ class Metadata : public eckit::LocalConfiguration { // extended LocalConfiguration: using eckit::LocalConfiguration::set; - Metadata& set(const eckit::LocalConfiguration& other) { + + Metadata& set(const eckit::Configuration& other) { +#if ATLAS_IO_ECKIT_VERSION_AT_LEAST(1, 26, 0) || ATLAS_IO_ECKIT_DEVELOP + LocalConfiguration::set(other); +#else eckit::Value& root = const_cast(get()); auto& other_root = other.get(); std::vector other_keys; @@ -65,6 +69,7 @@ class Metadata : public eckit::LocalConfiguration { for (auto& key : other_keys) { root[key] = other_root[key]; } +#endif return *this; } @@ -79,17 +84,24 @@ class Metadata : public eckit::LocalConfiguration { Metadata& remove(const std::string& name) { +#if ATLAS_IO_ECKIT_VERSION_AT_LEAST(1, 26, 0) || ATLAS_IO_ECKIT_DEVELOP + LocalConfiguration::remove(name); +#else eckit::Value& root = const_cast(get()); root.remove(name); +#endif return *this; } std::vector keys() const { - // Preserves order of keys +#if ATLAS_IO_ECKIT_VERSION_AT_LEAST(1, 26, 0) || ATLAS_IO_ECKIT_DEVELOP + return LocalConfiguration::keys(); +#else std::vector result; eckit::fromValue(result, get().keys()); return result; +#endif } }; diff --git a/atlas_io/src/atlas_io/detail/defines.h.in b/atlas_io/src/atlas_io/detail/defines.h.in index 40d8132c8..f3fb79f42 100644 --- a/atlas_io/src/atlas_io/detail/defines.h.in +++ b/atlas_io/src/atlas_io/detail/defines.h.in @@ -18,4 +18,9 @@ #cmakedefine01 ATLAS_IO_LITTLE_ENDIAN #cmakedefine01 ATLAS_IO_BIG_ENDIAN +#define ATLAS_IO_ECKIT_VERSION_INT @ATLAS_IO_ECKIT_VERSION_INT@ +#define ATLAS_IO_ECKIT_DEVELOP @ATLAS_IO_ECKIT_DEVELOP@ + +#define ATLAS_IO_ECKIT_VERSION_AT_LEAST(x, y, z) (ATLAS_IO_ECKIT_VERSION_INT >= x * 10000 + y * 100 + z) + #endif diff --git a/atlas_io/tests/TestEnvironment.h b/atlas_io/tests/TestEnvironment.h index aa2b99e90..7b66cee8f 100644 --- a/atlas_io/tests/TestEnvironment.h +++ b/atlas_io/tests/TestEnvironment.h @@ -27,7 +27,7 @@ #include "eckit/types/Types.h" #include "atlas_io/atlas-io.h" -#include "atlas_io/detail/BlackMagic.h" +#include "atlas_io/Trace.h" namespace atlas { namespace test { diff --git a/atlas_io/tests/test_io_stream.cc b/atlas_io/tests/test_io_stream.cc index 46ed49874..598e26b04 100644 --- a/atlas_io/tests/test_io_stream.cc +++ b/atlas_io/tests/test_io_stream.cc @@ -8,12 +8,10 @@ * nor does it submit to any jurisdiction. */ -#include "atlas_io/FileStream.h" -#include "atlas_io/Session.h" +#include "atlas_io/atlas-io.h" #include "eckit/io/FileHandle.h" #include "eckit/io/PooledHandle.h" - #include "TestEnvironment.h" diff --git a/src/atlas/CMakeLists.txt b/src/atlas/CMakeLists.txt index 02d2ec38d..634a22dce 100644 --- a/src/atlas/CMakeLists.txt +++ b/src/atlas/CMakeLists.txt @@ -524,7 +524,7 @@ functionspace/detail/PointCloudInterface.cc functionspace/detail/CubedSphereStructure.h functionspace/detail/CubedSphereStructure.cc -# for cubedsphere matching mesh partitioner +# for cubedsphere matching mesh partitioner interpolation/method/cubedsphere/CellFinder.cc interpolation/method/cubedsphere/CellFinder.h interpolation/Vector2D.cc @@ -541,8 +541,8 @@ interpolation/element/Triag3D.cc interpolation/element/Triag3D.h interpolation/method/Intersect.cc interpolation/method/Intersect.h -interpolation/method/Ray.cc # For testing Quad -interpolation/method/Ray.h # For testing Quad +interpolation/method/Ray.cc # For testing Quad +interpolation/method/Ray.h # For testing Quad # for BuildConvexHull3D @@ -634,6 +634,11 @@ interpolation/method/knn/KNearestNeighboursBase.cc interpolation/method/knn/KNearestNeighboursBase.h interpolation/method/knn/NearestNeighbour.cc interpolation/method/knn/NearestNeighbour.h +interpolation/method/sphericalvector/ComplexMatrixMultiply.h +interpolation/method/sphericalvector/SparseMatrix.h +interpolation/method/sphericalvector/SphericalVector.cc +interpolation/method/sphericalvector/SphericalVector.h +interpolation/method/sphericalvector/Types.h interpolation/method/structured/Cubic2D.cc interpolation/method/structured/Cubic2D.h interpolation/method/structured/Cubic3D.cc @@ -866,7 +871,7 @@ if( NOT atlas_HAVE_ATLAS_FUNCTIONSPACE ) unset( atlas_parallel_srcs ) unset( atlas_output_srcs ) unset( atlas_redistribution_srcs ) - unset( atlas_linalg_srcs ) # only depends on array + unset( atlas_linalg_srcs ) # only depends on array endif() if( NOT atlas_HAVE_ATLAS_INTERPOLATION ) diff --git a/src/atlas/array/helpers/ArrayForEach.h b/src/atlas/array/helpers/ArrayForEach.h index 38d4d6441..f20b4c80b 100644 --- a/src/atlas/array/helpers/ArrayForEach.h +++ b/src/atlas/array/helpers/ArrayForEach.h @@ -1,5 +1,5 @@ /* - * (C) Crown Copyright 2023 Met Office + * (C) Crown Copyright 2024 Met Office * * This software is licensed under the terms of the Apache Licence Version 2.0 * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. @@ -7,9 +7,10 @@ #pragma once +#include #include #include -#include +#include #include "atlas/array/ArrayView.h" #include "atlas/array/Range.h" @@ -22,18 +23,19 @@ namespace atlas { namespace execution { -// As in C++17 std::execution namespace. Note: unsequenced_policy is a C++20 addition. +// As in C++17 std::execution namespace. Note: unsequenced_policy is a C++20 +// addition. class sequenced_policy {}; class unsequenced_policy {}; class parallel_unsequenced_policy {}; class parallel_policy {}; -// execution policy objects as in C++ std::execution namespace. Note: unseq is a C++20 addition. -inline constexpr sequenced_policy seq{ /*unspecified*/ }; -inline constexpr parallel_policy par{ /*unspecified*/ }; -inline constexpr parallel_unsequenced_policy par_unseq{ /*unspecified*/ }; -inline constexpr unsequenced_policy unseq{ /*unspecified*/ }; - +// execution policy objects as in C++ std::execution namespace. Note: unseq is a +// C++20 addition. +inline constexpr sequenced_policy seq{/*unspecified*/}; +inline constexpr parallel_policy par{/*unspecified*/}; +inline constexpr parallel_unsequenced_policy par_unseq{/*unspecified*/}; +inline constexpr unsequenced_policy unseq{/*unspecified*/}; // Type names for execution policy (Not in C++ standard) template @@ -64,34 +66,31 @@ constexpr std::string_view policy_name(execution_policy) { // Type check for execution policy (Not in C++ standard) template -constexpr auto is_execution_policy() { return - std::is_same_v || - std::is_same_v || - std::is_same_v || - std::is_same_v; +constexpr auto is_execution_policy() { + return std::is_same_v || + std::is_same_v || + std::is_same_v || + std::is_same_v; } template constexpr auto demote_policy() { - if constexpr(std::is_same_v) { + if constexpr (std::is_same_v) { return unseq; - } - else if constexpr (std::is_same_v) { + } else if constexpr (std::is_same_v) { return seq; - } - else { + } else { return ExecutionPolicy{}; } } template -constexpr auto is_omp_policy() { return - std::is_same_v || - std::is_same_v; +constexpr auto is_omp_policy() { + return std::is_same_v || + std::is_same_v; } - -template +template using demote_policy_t = decltype(demote_policy()); } // namespace execution @@ -102,7 +101,8 @@ namespace option { template util::Config execution_policy() { - return util::Config("execution_policy", execution::policy_name>()); + return util::Config("execution_policy", + execution::policy_name>()); } template @@ -110,7 +110,7 @@ util::Config execution_policy(T) { return execution_policy>(); } -} // namespace option +} // namespace option namespace array { namespace helpers { @@ -119,7 +119,9 @@ namespace detail { struct NoMask { template - constexpr bool operator()(Args...) const { return 0; } + constexpr bool operator()(Args...) const { + return 0; + } }; inline constexpr NoMask no_mask; @@ -131,13 +133,11 @@ constexpr auto tuplePushBack(const std::tuple& tuple, T value) { template void forEach(idx_t idxMax, const Functor& functor) { - - if constexpr(execution::is_omp_policy()) { + if constexpr (execution::is_omp_policy()) { atlas_omp_parallel_for(auto idx = idx_t{}; idx < idxMax; ++idx) { functor(idx); } - } - else { + } else { // Simple for-loop for sequenced or unsequenced execution policies. for (auto idx = idx_t{}; idx < idxMax; ++idx) { functor(idx); @@ -147,11 +147,10 @@ void forEach(idx_t idxMax, const Functor& functor) { template constexpr auto argPadding() { - if constexpr(NPad > 0) { + if constexpr (NPad > 0) { return std::tuple_cat(std::make_tuple(Range::all()), argPadding()); - } - else { + } else { return std::make_tuple(); } } @@ -159,29 +158,31 @@ constexpr auto argPadding() { template auto makeSlices(const std::tuple& slicerArgs, ArrayViewTuple&& arrayViews) { - constexpr auto nb_views = std::tuple_size_v; auto&& arrayView = std::get(arrayViews); using ArrayView = std::decay_t; - constexpr auto Dim = sizeof...(SlicerArgs); - constexpr auto Rank = ArrayView::rank(); + constexpr auto Dim = sizeof...(SlicerArgs); + constexpr auto Rank = ArrayView::rank(); - static_assert (Dim <= Rank, "Error: number of slicer arguments exceeds the rank of ArrayView."); + static_assert( + Dim <= Rank, + "Error: number of slicer arguments exceeds the rank of ArrayView."); const auto paddedArgs = std::tuple_cat(slicerArgs, argPadding()); const auto slicer = [&arrayView](const auto&... args) { return std::make_tuple(arrayView.slice(args...)); }; - if constexpr (ViewIdx == nb_views-1) { + if constexpr (ViewIdx == nb_views - 1) { return std::apply(slicer, paddedArgs); - } - else { + } else { // recurse - return std::tuple_cat(std::apply(slicer, paddedArgs), - makeSlices(slicerArgs, std::forward(arrayViews))); + return std::tuple_cat( + std::apply(slicer, paddedArgs), + makeSlices(slicerArgs, + std::forward(arrayViews))); } } @@ -192,80 +193,79 @@ template struct ArrayForEachImpl { template - static void apply(ArrayViewTuple&& arrayViews, - const Mask& mask, + static void apply(ArrayViewTuple&& arrayViews, const Mask& mask, const Function& function, const std::tuple& slicerArgs, const std::tuple& maskArgs) { // Iterate over this dimension. - if constexpr(Dim == ItrDim) { - + if constexpr (Dim == ItrDim) { // Get size of iteration dimenion from first view argument. const auto idxMax = std::get<0>(arrayViews).shape(ItrDim); forEach(idxMax, [&](idx_t idx) { - // Demote parallel execution policy to a non-parallel one in further recursion - ArrayForEachImpl, Dim + 1, ItrDims...>::apply( - std::forward(arrayViews), mask, function, - tuplePushBack(slicerArgs, idx), - tuplePushBack(maskArgs, idx)); + // Demote parallel execution policy to a non-parallel one in further + // recursion + ArrayForEachImpl< + execution::demote_policy_t, Dim + 1, + ItrDims...>::apply(std::forward(arrayViews), mask, + function, tuplePushBack(slicerArgs, idx), + tuplePushBack(maskArgs, idx)); }); } // Add a RangeAll to arguments. else { ArrayForEachImpl::apply( std::forward(arrayViews), mask, function, - tuplePushBack(slicerArgs, Range::all()), - maskArgs); + tuplePushBack(slicerArgs, Range::all()), maskArgs); } } }; template - struct is_applicable : std::false_type {}; +struct is_applicable : std::false_type {}; template -struct is_applicable> : std::is_invocable {}; +struct is_applicable> + : std::is_invocable {}; template -inline constexpr bool is_applicable_v = is_applicable::value; +inline constexpr bool is_applicable_v = is_applicable::value; template struct ArrayForEachImpl { - template - static void apply(ArrayViewTuple&& arrayViews, - const Mask& mask, + static void apply(ArrayViewTuple&& arrayViews, const Mask& mask, const Function& function, const std::tuple& slicerArgs, const std::tuple& maskArgs) { - constexpr auto maskPresent = !std::is_same_v; if constexpr (maskPresent) { - - constexpr auto invocableMask = std::is_invocable_r_v; - static_assert (invocableMask, - "Cannot invoke mask function with given arguments.\n" - "Make sure you arguments are N integers (or auto...) " - "where N == sizeof...(ItrDims). Function must return an int." - ); - - if (std::apply(mask, maskArgs)) { - return; - } - + constexpr auto invocableMask = + std::is_invocable_r_v; + static_assert( + invocableMask, + "Cannot invoke mask function with given arguments.\n" + "Make sure you arguments are N integers (or auto...) " + "where N == sizeof...(ItrDims). Function must return an int."); + + if (std::apply(mask, maskArgs)) { + return; + } } - auto slices = makeSlices(slicerArgs, std::forward(arrayViews)); + auto slices = + makeSlices(slicerArgs, std::forward(arrayViews)); - constexpr auto applicable = is_applicable_v; - static_assert(applicable, "Cannot invoke function with given arguments. " - "Make sure you the arguments are rvalue references (Slice&&) or const references (const Slice&) or regular value (Slice)" ); + constexpr auto applicable = is_applicable_v; + static_assert( + applicable, + "Cannot invoke function with given arguments. " + "Make sure you the arguments are rvalue references (Slice&&) or const " + "references (const Slice&) or regular value (Slice)"); std::apply(function, std::move(slices)); } - }; } // namespace detail @@ -286,41 +286,36 @@ struct ArrayForEach { /// and is executed with signature g(idx_i, idx_j,...), where the idxs /// are indices of ItrDims. /// When a config is supplied containing "execution_policy" = - /// "sequenced_policy" (default). All loops are then executed in sequential - /// (row-major) order. - /// With "execution_policy" = "parallel_unsequenced" the first loop is executed - /// using OpenMP. The remaining loops are executed in serial. - /// Note: The lowest ArrayView.rank() must be greater than or equal - /// to the highest dim in ItrDims. TODO: static checking for this. + /// "sequenced_policy" (default). All loops are then executed in + /// sequential (row-major) order. With "execution_policy" = + /// "parallel_unsequenced" the first loop is executed using OpenMP. + /// The remaining loops are executed in serial. Note: The lowest + /// ArrayView.rank() must be greater than or equal to the highest dim + /// in ItrDims. TODO: static checking for this. template static void apply(const eckit::Parametrisation& conf, - std::tuple&& arrayViews, - const Mask& mask, const Function& function) { - + std::tuple&& arrayViews, const Mask& mask, + const Function& function) { auto execute = [&](auto execution_policy) { apply(execution_policy, std::move(arrayViews), mask, function); }; using namespace execution; std::string execution_policy; - if (conf.get("execution_policy",execution_policy)) { + if (conf.get("execution_policy", execution_policy)) { if (execution_policy == policy_name(par_unseq)) { execute(par_unseq); - } - else if (execution_policy == policy_name(par)) { + } else if (execution_policy == policy_name(par)) { execute(par); - } - else if (execution_policy == policy_name(unseq)) { + } else if (execution_policy == policy_name(unseq)) { execute(unseq); - } - else if (execution_policy == policy_name(seq)) { + } else if (execution_policy == policy_name(seq)) { execute(seq); + } else { + throw_Exception("Unrecognized execution policy " + execution_policy, + Here()); } - else { - throw_Exception("Unrecognized execution policy "+execution_policy, Here()); - } - } - else { + } else { execute(seq); } } @@ -328,35 +323,47 @@ struct ArrayForEach { /// brief Apply "For-Each" method. /// /// details As above, but Execution policy is determined at compile-time. - template ()>> - static void apply(ExecutionPolicy, std::tuple&& arrayViews, const Mask& mask, const Function& function) { - detail::ArrayForEachImpl::apply( - std::move(arrayViews), mask, function, std::make_tuple(), std::make_tuple()); + template ()>> + static void apply(ExecutionPolicy, std::tuple&& arrayViews, + const Mask& mask, const Function& function) { + detail::ArrayForEachImpl::apply( + std::move(arrayViews), mask, function, std::make_tuple(), + std::make_tuple()); } /// brief Apply "For-Each" method /// /// details Apply ForEach with default execution policy. template - static void apply(std::tuple&& arrayViews, const Mask& mask, const Function& function) { - apply(std::move(arrayViews), mask, function); + static void apply(std::tuple&& arrayViews, const Mask& mask, + const Function& function) { + apply(std::move(arrayViews), mask, function); } /// brief Apply "For-Each" method /// - /// details Apply ForEach with run-time determined execution policy and no mask. + /// details Apply ForEach with run-time determined execution policy and no + /// mask. template - static void apply(const eckit::Parametrisation& conf, std::tuple&& arrayViews, const Function& function) { + static void apply(const eckit::Parametrisation& conf, + std::tuple&& arrayViews, + const Function& function) { apply(conf, std::move(arrayViews), detail::no_mask, function); } /// brief Apply "For-Each" method /// - /// details Apply ForEach with compile-time determined execution policy and no mask. + /// details Apply ForEach with compile-time determined execution policy and no + /// mask. template ()>> - static void apply(ExecutionPolicy executionPolicy, std::tuple&& arrayViews, const Function& function) { + typename = std::enable_if_t< + execution::is_execution_policy()>> + static void apply(ExecutionPolicy executionPolicy, + std::tuple&& arrayViews, + const Function& function) { apply(executionPolicy, std::move(arrayViews), detail::no_mask, function); } @@ -364,12 +371,22 @@ struct ArrayForEach { /// /// details Apply ForEach with default execution policy and no mask. template - static void apply(std::tuple&& arrayViews, const Function& function) { + static void apply(std::tuple&& arrayViews, + const Function& function) { apply(execution::seq, std::move(arrayViews), function); } - }; +/// brief Construct ArrayForEach and call apply +/// +/// details Construct an ArrayForEach using std::integer_sequence +/// . Remaining arguments are forwarded to apply +/// method. +template +void arrayForEachDim(std::integer_sequence, Args&&... args) { + ArrayForEach::apply(std::forward(args)...); +} + } // namespace helpers } // namespace array } // namespace atlas diff --git a/src/atlas/functionspace/NodeColumns.cc b/src/atlas/functionspace/NodeColumns.cc index 86259313b..76fe397e3 100644 --- a/src/atlas/functionspace/NodeColumns.cc +++ b/src/atlas/functionspace/NodeColumns.cc @@ -283,6 +283,10 @@ void NodeColumns::set_field_metadata(const eckit::Configuration& config, Field& idx_t variables(0); config.get("variables", variables); field.set_variables(variables); + + if (config.has("type")) { + field.metadata().set("type", config.getString("type")); + } } array::DataType NodeColumns::config_datatype(const eckit::Configuration& config) const { diff --git a/src/atlas/functionspace/PointCloud.cc b/src/atlas/functionspace/PointCloud.cc index 61388a2f5..702e4bef5 100644 --- a/src/atlas/functionspace/PointCloud.cc +++ b/src/atlas/functionspace/PointCloud.cc @@ -1001,7 +1001,7 @@ void PointCloud::setupGatherScatter() { array::make_view(remote_index_).data(), REMOTE_IDX_BASE, array::make_view(global_index_).data(), - array::make_view(ghost_).data(), + array::make_view(ghost_).data(), ghost_.size()); size_global_ = gather_scatter_->glb_dof(); } diff --git a/src/atlas/functionspace/detail/StructuredColumns_setup.cc b/src/atlas/functionspace/detail/StructuredColumns_setup.cc index 01a9a1be4..4bdc0a4e9 100644 --- a/src/atlas/functionspace/detail/StructuredColumns_setup.cc +++ b/src/atlas/functionspace/detail/StructuredColumns_setup.cc @@ -371,7 +371,7 @@ void StructuredColumns::setup(const grid::Distribution& distribution, const ecki idx_t jj_max = j + halo; if (regional) { jj_min = std::max(jj_min, idx_t{0}); - jj_max = std::min(jj_max, grid_->nx(j) - idx_t{1}); + jj_max = std::min(jj_max, grid_->ny() - idx_t{1}); } for (idx_t jj = jj_min; jj <= jj_max; ++jj) { idx_t jjj = compute_j(jj); diff --git a/src/atlas/interpolation/method/Method.h b/src/atlas/interpolation/method/Method.h index 43bd1f4d8..5e5abbfde 100644 --- a/src/atlas/interpolation/method/Method.h +++ b/src/atlas/interpolation/method/Method.h @@ -127,6 +127,8 @@ class Method : public util::Object { virtual void do_setup(const FunctionSpace& source, const Field& target); virtual void do_setup(const FunctionSpace& source, const FieldSet& target); + void check_compatibility(const Field& src, const Field& tgt, const Matrix& W) const; + private: template void interpolate_field(const Field& src, Field& tgt, const Matrix&) const; @@ -152,8 +154,6 @@ class Method : public util::Object { template void adjoint_interpolate_field_rank3(Field& src, const Field& tgt, const Matrix&) const; - void check_compatibility(const Field& src, const Field& tgt, const Matrix& W) const; - private: const Matrix* matrix_ = nullptr; std::shared_ptr matrix_shared_; diff --git a/src/atlas/interpolation/method/MethodFactory.cc b/src/atlas/interpolation/method/MethodFactory.cc index 7c9f6952c..4baa44c37 100644 --- a/src/atlas/interpolation/method/MethodFactory.cc +++ b/src/atlas/interpolation/method/MethodFactory.cc @@ -16,6 +16,7 @@ #include "knn/GridBoxMaximum.h" #include "knn/KNearestNeighbours.h" #include "knn/NearestNeighbour.h" +#include "sphericalvector/SphericalVector.h" #include "structured/Cubic2D.h" #include "structured/Cubic3D.h" #include "structured/Linear2D.h" @@ -47,6 +48,7 @@ void force_link() { MethodBuilder(); MethodBuilder(); MethodBuilder(); + MethodBuilder(); } } link; } diff --git a/src/atlas/interpolation/method/sphericalvector/ComplexMatrixMultiply.h b/src/atlas/interpolation/method/sphericalvector/ComplexMatrixMultiply.h new file mode 100644 index 000000000..34cfdacb2 --- /dev/null +++ b/src/atlas/interpolation/method/sphericalvector/ComplexMatrixMultiply.h @@ -0,0 +1,217 @@ +/* + * (C) Crown Copyright 2024 Met Office + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#pragma once + +#include +#include +#include +#include +#include +#include + +#include "atlas/array/ArrayView.h" +#include "atlas/array/LocalView.h" +#include "atlas/array/Range.h" +#include "atlas/array/helpers/ArrayForEach.h" +#include "atlas/interpolation/method/sphericalvector/Types.h" +#include "atlas/parallel/omp/omp.h" + +namespace atlas { +namespace interpolation { +namespace method { +namespace detail { + +struct TwoVectorTag {}; +struct ThreeVectorTag {}; +constexpr auto twoVector = TwoVectorTag{}; +constexpr auto threeVector = ThreeVectorTag{}; + +// Permits the equivalent of static_assert(false, msg). Will be addressed in +// C++26. +template +constexpr bool always_false_v = false; + +/// @brief Helper class to perform complex matrix multiplications +/// +/// @details Performs matrix multiplication between fields of 2-vectors or +/// 3-vectors. Fields must have Rank >= 2. Here, the assumption is +/// that Dim = 0 is the horizontal dimension, and Dim = (Rank - 1) is +/// the vector element dimension. +template +class ComplexMatrixMultiply { + public: + ComplexMatrixMultiply() = default; + + /// @brief Construct object from sparse matrices. + /// + /// @details complexWeights is a SparseMatrix of weights. realWeights is a + /// SparseMatrix containing the magnitudes of the elements of + /// complexWeights. + ComplexMatrixMultiply(const ComplexMatrix& complexWeights, + const RealMatrix& realWeights) + : complexWeightsPtr_{&complexWeights}, realWeightsPtr_{&realWeights} { + if constexpr (ATLAS_BUILD_TYPE_DEBUG) { + ATLAS_ASSERT(complexWeightsPtr_->rows() == realWeightsPtr_->rows()); + ATLAS_ASSERT(complexWeightsPtr_->cols() == realWeightsPtr_->cols()); + ATLAS_ASSERT(complexWeightsPtr_->nonZeros() == + realWeightsPtr_->nonZeros()); + for (auto rowIndex = Index{0}; rowIndex < complexWeightsPtr_->rows(); + ++rowIndex) { + for (auto [complexRowIter, realRowIter] = rowIters(rowIndex); + complexRowIter; ++complexRowIter, ++realRowIter) { + ATLAS_ASSERT(realRowIter); + ATLAS_ASSERT(realRowIter.row() == complexRowIter.row()); + ATLAS_ASSERT(realRowIter.col() == complexRowIter.col()); + // tinyNum ~= 2.3e-13 for double. + constexpr auto tinyNum = 1024 * std::numeric_limits::epsilon(); + const auto complexMagnitude = std::abs(complexRowIter.value()); + const auto realValue = realRowIter.value(); + const auto error = std::abs(complexMagnitude - realValue); + + const auto printError = [&]() { + auto strStream = std::ostringstream{}; + strStream << "Error complex components: "; + strStream << std::setprecision(19) << error; + return strStream.str(); + }; + ATLAS_ASSERT(error < tinyNum, printError()); + } + } + } + } + + /// @brief Apply complex matrix vector multiplication. + /// + /// @details Multiply weights by the elements in sourceView to give + /// elements in targetView. If VectorType == TwoVectorTag, + /// complexWeights are applied to the horizontal elements of + /// sourceView. If VectorType == ThreeVectorTag, then realWeights + /// are additionally applied to the vertical elements of sourceView. + template + void apply(const array::ArrayView& sourceView, + array::ArrayView& targetView, VectorType) const { + if constexpr (std::is_same_v) { + applyTwoVector(sourceView, targetView); + } else if constexpr (std::is_same_v) { + applyThreeVector(sourceView, targetView); + } else { + static_assert(always_false_v, "Unknown vector type"); + } + } + + private: + /// @brief Apply 2-vector MatMul. + template + void applyTwoVector(const array::ArrayView& sourceView, + array::ArrayView& targetView) const { + // We could probably optimise contiguous arrays using + // reinterpret_cast*>(view.data()). This is fine + // according to the C++ standard! + atlas_omp_parallel_for(auto rowIndex = Index{0}; + rowIndex < complexWeightsPtr_->rows(); ++rowIndex) { + auto targetSlice = sliceColumn(targetView, rowIndex); + if constexpr (InitialiseTarget) { + targetSlice.assign(0.); + } + + for (auto complexRowIter = complexWeightsPtr_->rowIter(rowIndex); + complexRowIter; ++complexRowIter) { + const auto colIndex = complexRowIter.col(); + const auto complexWeight = complexRowIter.value(); + const auto sourceSlice = sliceColumn(sourceView, colIndex); + multiplyAdd(sourceSlice, targetSlice, complexWeight); + } + } + } + + /// @brief Apply 3-vector MatMul. + template + void applyThreeVector(const array::ArrayView& sourceView, + array::ArrayView& targetView) const { + atlas_omp_parallel_for(auto rowIndex = Index{0}; + rowIndex < complexWeightsPtr_->rows(); ++rowIndex) { + auto targetSlice = sliceColumn(targetView, rowIndex); + if constexpr (InitialiseTarget) { + targetSlice.assign(0.); + } + + for (auto [complexRowIter, realRowIter] = rowIters(rowIndex); + complexRowIter; ++complexRowIter, ++realRowIter) { + const auto colIndex = complexRowIter.col(); + const auto complexWeight = complexRowIter.value(); + const auto realWeight = realRowIter.value(); + const auto sourceSlice = sliceColumn(sourceView, colIndex); + multiplyAdd(sourceSlice, targetSlice, complexWeight, realWeight); + } + } + } + + /// @brief Multiply source column by weight(s) and add to target column. + template + void multiplyAdd(const array::LocalView& sourceColumn, + array::LocalView& targetColumn, + Weights... weights) const { + for (auto levelIdx = 0; levelIdx < sourceColumn.shape(0); ++levelIdx) { + const auto sourceElem = sliceColumn(sourceColumn, levelIdx); + auto targetElem = sliceColumn(targetColumn, levelIdx); + multiplyAdd(sourceElem, targetElem, weights...); + } + } + + /// @brief Multiply source element by complex weight and add to target + /// element. + template + void multiplyAdd(const array::LocalView& sourceElem, + array::LocalView& targetElem, + Complex complexWeight) const { + const auto targetVector = + complexWeight * Complex(sourceElem(0), sourceElem(1)); + targetElem(0) += targetVector.real(); + targetElem(1) += targetVector.imag(); + } + + /// @brief Multiply source element by complex and real weights and add to + /// target element. + template + void multiplyAdd(const array::LocalView& sourceElem, + array::LocalView& targetElem, + Complex complexWeight, Real realWeight) const { + multiplyAdd(sourceElem, targetElem, complexWeight); + targetElem(2) += realWeight * sourceElem(2); + } + + /// @brief Return a pair of complex and real row iterators + std::pair rowIters( + Index rowIndex) const { + return std::make_pair(complexWeightsPtr_->rowIter(rowIndex), + realWeightsPtr_->rowIter(rowIndex)); + } + + /// @brief Makes the slice arrayView.slice(index, Range::all()...). + template + static auto sliceColumn(View& arrayView, Index index) { + constexpr auto Rank = std::decay_t::rank(); + using RangeAll = decltype(array::Range::all()); + + const auto slicerArgs = std::tuple_cat(std::make_tuple(index), + std::array{}); + const auto slicer = [&](auto&&... args) { + return arrayView.slice(args...); + }; + + return std::apply(slicer, slicerArgs); + } + + const ComplexMatrix* complexWeightsPtr_{}; + const RealMatrix* realWeightsPtr_{}; +}; + +} // namespace detail +} // namespace method +} // namespace interpolation +} // namespace atlas diff --git a/src/atlas/interpolation/method/sphericalvector/SparseMatrix.h b/src/atlas/interpolation/method/sphericalvector/SparseMatrix.h new file mode 100644 index 000000000..a399d80ca --- /dev/null +++ b/src/atlas/interpolation/method/sphericalvector/SparseMatrix.h @@ -0,0 +1,107 @@ +/* + * (C) Crown Copyright 2024 Met Office + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#pragma once + +#include +#include +#include + +#include "atlas/library/defines.h" +#if ATLAS_HAVE_EIGEN + +#include +#endif + +#include "atlas/runtime/Exception.h" +#include "eckit/log/CodeLocation.h" + +namespace atlas { +namespace interpolation { +namespace method { +namespace detail { + +#if ATLAS_HAVE_EIGEN +/// @brief Wrapper class for Eigen sparse matrix +/// +/// @details Eigen sparse matrix wrapper. Allows preprocessor disabling of class +/// if Eigen library is not present. +template +class SparseMatrix { + public: + using Index = int; + using EigenMatrix = Eigen::SparseMatrix; + using Triplet = Eigen::Triplet; + using Triplets = std::vector; + using RowIter = typename EigenMatrix::InnerIterator; + + SparseMatrix() = default; + + SparseMatrix(Index nRows, Index nCols, const Triplets& triplets) + : eigenMatrix_(nRows, nCols) { + eigenMatrix_.setFromTriplets(triplets.begin(), triplets.end()); + } + + Index nonZeros() const { return eigenMatrix_.nonZeros(); } + Index rows() const { return eigenMatrix_.rows(); } + Index cols() const { return eigenMatrix_.cols(); } + RowIter rowIter(Index rowIndex) const { + return RowIter(eigenMatrix_, rowIndex); + } + SparseMatrix adjoint() const { + auto adjointMatrix = SparseMatrix{}; + adjointMatrix.eigenMatrix_ = eigenMatrix_.adjoint(); + return adjointMatrix; + } + + private: + EigenMatrix eigenMatrix_{}; +}; +#else + +template +class SparseMatrix { + public: + using Index = int; + + class Triplet { + public: + template + constexpr Triplet(const Args&... args) {} + }; + using Triplets = std::vector; + + class RowIter { + public: + template + constexpr RowIter(const Args&... args) {} + constexpr Index row() const { return Index{}; } + constexpr Index col() const { return Index{}; } + constexpr Value value() const { return Value{}; } + constexpr operator bool() const { return false; } + constexpr RowIter& operator++() { return *this; } + }; + + constexpr SparseMatrix() = default; + template + SparseMatrix(const Args&... args) { + throw_Exception("Atlas has been compiled without Eigen", Here()); + } + constexpr Index nonZeros() const { return Index{}; } + constexpr Index rows() const { return Index{}; } + constexpr Index cols() const { return Index{}; } + constexpr RowIter rowIter(Index rowIndex) const { return RowIter{}; } + constexpr SparseMatrix adjoint() const { + return SparseMatrix{}; + } +}; +#endif + +} // namespace detail +} // namespace method +} // namespace interpolation +} // namespace atlas diff --git a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc new file mode 100644 index 000000000..7e29f41c1 --- /dev/null +++ b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc @@ -0,0 +1,245 @@ +/* + * (C) Crown Copyright 2024 Met Office + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#include "atlas/interpolation/method/sphericalvector/SphericalVector.h" + +#include "atlas/array/ArrayView.h" +#include "atlas/field/Field.h" +#include "atlas/field/FieldSet.h" +#include "atlas/interpolation/Cache.h" +#include "atlas/interpolation/Interpolation.h" +#include "atlas/interpolation/method/MethodFactory.h" +#include "atlas/interpolation/method/sphericalvector/ComplexMatrixMultiply.h" +#include "atlas/interpolation/method/sphericalvector/Types.h" +#include "atlas/option/Options.h" +#include "atlas/parallel/omp/omp.h" +#include "atlas/runtime/Exception.h" +#include "atlas/runtime/Trace.h" +#include "atlas/util/Constants.h" +#include "atlas/util/Geometry.h" +#include "eckit/config/LocalConfiguration.h" + +namespace atlas { +namespace interpolation { +namespace method { + +namespace { +MethodBuilder __builder("spherical-vector"); +} + +using namespace detail; + +using WeightsMatMul = detail::ComplexMatrixMultiply; +using WeightsMatMulAdjoint = detail::ComplexMatrixMultiply; + +SphericalVector::SphericalVector(const Config& config) : Method(config) { + const auto* conf = dynamic_cast(&config); + ATLAS_ASSERT(conf, "config must be derived from eckit::LocalConfiguration"); + interpolationScheme_ = conf->getSubConfiguration("scheme"); + adjoint_ = conf->getBool("adjoint", false); +} + +void SphericalVector::do_setup(const Grid& source, const Grid& target, + const Cache&) { + ATLAS_NOTIMPLEMENTED; +} + +void SphericalVector::do_setup(const FunctionSpace& source, + const FunctionSpace& target) { + ATLAS_TRACE("interpolation::method::SphericalVector::do_setup"); + source_ = source; + target_ = target; + + if (target_.size() == 0) { + return; + } + + setMatrix(Interpolation(interpolationScheme_, source_, target_)); + + // Get matrix data. + const auto nRows = static_cast(matrix().rows()); + const auto nCols = static_cast(matrix().cols()); + const auto nNonZeros = static_cast(matrix().nonZeros()); + const auto* outerIndices = matrix().outer(); + const auto* innerIndices = matrix().inner(); + const auto* baseWeights = matrix().data(); + + // Note: need to store copy of weights as Eigen3 sorts compressed rows by j + // whereas eckit does not. + auto complexTriplets = ComplexTriplets(nNonZeros); + auto realTriplets = RealTriplets(nNonZeros); + + const auto sourceLonLatsView = array::make_view(source_.lonlat()); + const auto targetLonLatsView = array::make_view(target_.lonlat()); + const auto unitSphere = geometry::UnitSphere{}; + + atlas_omp_parallel_for(auto rowIndex = Index{0}; rowIndex < nRows; + ++rowIndex) { + for (auto dataIndex = outerIndices[rowIndex]; + dataIndex < outerIndices[rowIndex + 1]; ++dataIndex) { + const auto colIndex = static_cast(innerIndices[dataIndex]); + const auto baseWeight = baseWeights[dataIndex]; + + const auto sourceLonLat = PointLonLat(sourceLonLatsView(colIndex, 0), + sourceLonLatsView(colIndex, 1)); + const auto targetLonLat = PointLonLat(targetLonLatsView(rowIndex, 0), + targetLonLatsView(rowIndex, 1)); + + const auto alpha = + unitSphere.greatCircleCourse(sourceLonLat, targetLonLat); + + const auto deltaAlpha = + (alpha.first - alpha.second) * util::Constants::degreesToRadians(); + + complexTriplets[dataIndex] = + ComplexTriplet{rowIndex, colIndex, + Complex{baseWeight * std::cos(deltaAlpha), + baseWeight * std::sin(deltaAlpha)}}; + realTriplets[dataIndex] = RealTriplet{rowIndex, colIndex, baseWeight}; + } + } + + complexWeights_ = ComplexMatrix(nRows, nCols, complexTriplets); + realWeights_ = RealMatrix(nRows, nCols, realTriplets); + + if (adjoint_) { + complexWeightsAdjoint_ = complexWeights_.adjoint(); + realWeightsAdjoint_ = realWeights_.adjoint(); + } +} + +void SphericalVector::print(std::ostream&) const { ATLAS_NOTIMPLEMENTED; } + +void SphericalVector::do_execute(const FieldSet& sourceFieldSet, + FieldSet& targetFieldSet, + Metadata& metadata) const { + ATLAS_TRACE("atlas::interpolation::method::SphericalVector::do_execute()"); + ATLAS_ASSERT(sourceFieldSet.size() == targetFieldSet.size()); + + for (auto i = 0; i < sourceFieldSet.size(); ++i) { + do_execute(sourceFieldSet[i], targetFieldSet[i], metadata); + } +} + +void SphericalVector::do_execute(const Field& sourceField, Field& targetField, + Metadata&) const { + ATLAS_TRACE("atlas::interpolation::method::SphericalVector::do_execute()"); + const auto fieldType = sourceField.metadata().getString("type", ""); + if (fieldType != "vector") { + auto metadata = Metadata(); + Method::do_execute(sourceField, targetField, metadata); + + return; + } + + Method::check_compatibility(sourceField, targetField, matrix()); + + haloExchange(sourceField); + interpolate_vector_field(sourceField, targetField, + WeightsMatMul(complexWeights_, realWeights_)); + targetField.set_dirty(); +} + +void SphericalVector::do_execute_adjoint(FieldSet& sourceFieldSet, + const FieldSet& targetFieldSet, + Metadata& metadata) const { + ATLAS_TRACE( + "atlas::interpolation::method::SphericalVector::do_execute_adjoint()"); + ATLAS_ASSERT(sourceFieldSet.size() == targetFieldSet.size()); + + for (auto i = 0; i < sourceFieldSet.size(); ++i) { + do_execute_adjoint(sourceFieldSet[i], targetFieldSet[i], metadata); + } +} + +void SphericalVector::do_execute_adjoint(Field& sourceField, + const Field& targetField, + Metadata& metadata) const { + ATLAS_TRACE( + "atlas::interpolation::method::SphericalVector::do_execute_adjoint()"); + + const auto fieldType = sourceField.metadata().getString("type", ""); + if (fieldType != "vector") { + auto metadata = Metadata(); + Method::do_execute_adjoint(sourceField, targetField, metadata); + + return; + } + + Method::check_compatibility(sourceField, targetField, matrix()); + + ATLAS_ASSERT(adjoint_, "\"adjoint\" needs to be set to \"true\" in Config."); + interpolate_vector_field( + targetField, sourceField, + WeightsMatMulAdjoint(complexWeightsAdjoint_, realWeightsAdjoint_)); + adjointHaloExchange(sourceField); +} + +template +void SphericalVector::interpolate_vector_field(const Field& sourceField, + Field& targetField, + const MatMul& matMul) { + if (targetField.size() == 0) { + return; + } + + ATLAS_ASSERT_MSG(sourceField.variables() == 2 || sourceField.variables() == 3, + "Vector field can only have 2 or 3 components."); + + if (sourceField.datatype().kind() == array::DataType::KIND_REAL64) { + interpolate_vector_field(sourceField, targetField, matMul); + return; + } + + if (sourceField.datatype().kind() == array::DataType::KIND_REAL32) { + interpolate_vector_field(sourceField, targetField, matMul); + return; + } + + ATLAS_NOTIMPLEMENTED; +}; + +template +void SphericalVector::interpolate_vector_field(const Field& sourceField, + Field& targetField, + const MatMul& matMul) { + if (sourceField.rank() == 2) { + interpolate_vector_field(sourceField, targetField, matMul); + return; + } + + if (sourceField.rank() == 3) { + interpolate_vector_field(sourceField, targetField, matMul); + return; + } + + ATLAS_NOTIMPLEMENTED; +} + +template +void SphericalVector::interpolate_vector_field(const Field& sourceField, + Field& targetField, + const MatMul& matMul) { + const auto sourceView = array::make_view(sourceField); + auto targetView = array::make_view(targetField); + + if (sourceField.variables() == 2) { + matMul.apply(sourceView, targetView, twoVector); + return; + } + + if (sourceField.variables() == 3) { + matMul.apply(sourceView, targetView, threeVector); + return; + } + + ATLAS_NOTIMPLEMENTED; +} + +} // namespace method +} // namespace interpolation +} // namespace atlas diff --git a/src/atlas/interpolation/method/sphericalvector/SphericalVector.h b/src/atlas/interpolation/method/sphericalvector/SphericalVector.h new file mode 100644 index 000000000..c3514b0d1 --- /dev/null +++ b/src/atlas/interpolation/method/sphericalvector/SphericalVector.h @@ -0,0 +1,88 @@ +/* + * (C) Crown Copyright 2024 Met Office + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#pragma once + +#include "atlas/functionspace/FunctionSpace.h" +#include "atlas/interpolation/method/Method.h" +#include "atlas/interpolation/method/sphericalvector/Types.h" + +namespace atlas { +namespace interpolation { +namespace method { + +class SphericalVector : public Method { + public: + /// @brief Interpolation post-processor for vector field interpolation + /// + /// @details Takes a base interpolation config keyed to "scheme" and creates + /// A set of complex intepolation weights which rotate source vector + /// elements into the target elements' individual frame of reference. + /// Method works by creating a great-circle arc between the source + /// and target locations; the amount of rotation is determined by the + /// difference the in the great-cricle course (relative to north) at + /// the source and target location. + /// Both source and target fields require the "type" metadata to be + /// set to "vector" for this method to be invoked. Otherwise, the + /// base scalar interpolation method is invoked. + /// Note: This method only works with matrix-based interpolation + /// schemes. + /// + SphericalVector(const Config& config); + ~SphericalVector() override {} + + void print(std::ostream&) const override; + const FunctionSpace& source() const override { return source_; } + const FunctionSpace& target() const override { return target_; } + + void do_execute(const FieldSet& sourceFieldSet, FieldSet& targetFieldSet, + Metadata& metadata) const override; + void do_execute(const Field& sourceField, Field& targetField, + Metadata& metadata) const override; + + void do_execute_adjoint(FieldSet& sourceFieldSet, + const FieldSet& targetFieldSet, + Metadata& metadata) const override; + void do_execute_adjoint(Field& sourceField, const Field& targetField, + Metadata& metadata) const override; + + private: + template + static void interpolate_vector_field(const Field& sourceField, + Field& targetField, + const MatMul& matMul); + + template + static void interpolate_vector_field(const Field& sourceField, + Field& targetField, + const MatMul& matMul); + + template + static void interpolate_vector_field(const Field& sourceField, + Field& targetField, + const MatMul& matMul); + + using Method::do_setup; + void do_setup(const FunctionSpace& source, + const FunctionSpace& target) override; + void do_setup(const Grid& source, const Grid& target, const Cache&) override; + + eckit::LocalConfiguration interpolationScheme_{}; + bool adjoint_{}; + + FunctionSpace source_{}; + FunctionSpace target_{}; + + detail::ComplexMatrix complexWeights_{}; + detail::RealMatrix realWeights_{}; + detail::ComplexMatrix complexWeightsAdjoint_{}; + detail::RealMatrix realWeightsAdjoint_{}; +}; + +} // namespace method +} // namespace interpolation +} // namespace atlas diff --git a/src/atlas/interpolation/method/sphericalvector/Types.h b/src/atlas/interpolation/method/sphericalvector/Types.h new file mode 100644 index 000000000..f83e9640b --- /dev/null +++ b/src/atlas/interpolation/method/sphericalvector/Types.h @@ -0,0 +1,32 @@ +/* + * (C) Crown Copyright 2024 Met Office + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#pragma once + +#include + +#include "atlas/interpolation/method/sphericalvector/SparseMatrix.h" + +namespace atlas { +namespace interpolation { +namespace method { +namespace detail { + +using Real = double; +using Complex = std::complex; +using ComplexMatrix = SparseMatrix; +using RealMatrix = SparseMatrix; +using Index = ComplexMatrix::Index; +using ComplexTriplet = ComplexMatrix::Triplet; +using ComplexTriplets = ComplexMatrix::Triplets; +using RealTriplet = RealMatrix::Triplet; +using RealTriplets = RealMatrix::Triplets; + +} // detail +} // method +} // interpolation +} // atlas diff --git a/src/atlas/mesh/Mesh.cc b/src/atlas/mesh/Mesh.cc index 411260617..1526a74c4 100644 --- a/src/atlas/mesh/Mesh.cc +++ b/src/atlas/mesh/Mesh.cc @@ -9,11 +9,13 @@ */ #include "atlas/mesh/Mesh.h" -#include "atlas/mesh/Nodes.h" + #include "atlas/grid/Grid.h" #include "atlas/grid/Partitioner.h" +#include "atlas/mesh/Nodes.h" #include "atlas/meshgenerator/MeshGenerator.h" #include "atlas/parallel/mpi/mpi.h" +#include "atlas/runtime/Exception.h" namespace atlas { @@ -23,13 +25,13 @@ Mesh::Mesh(): Handle(new Implementation()) {} Mesh::Mesh(const Grid& grid, const eckit::Configuration& config): Handle([&]() { - if(config.has("mpi_comm")) { + if (config.has("mpi_comm")) { mpi::push(config.getString("mpi_comm")); } - util::Config cfg = grid.meshgenerator()|util::Config(config); - auto meshgenerator = MeshGenerator{grid.meshgenerator()|config}; - auto mesh = meshgenerator.generate(grid, grid::Partitioner(grid.partitioner()|config)); - if(config.has("mpi_comm")) { + auto cfg = grid.meshgenerator() | util::Config(config); + auto meshgenerator = MeshGenerator{grid.meshgenerator() | config}; + auto mesh = meshgenerator.generate(grid, grid::Partitioner(grid.partitioner() | config)); + if (config.has("mpi_comm")) { mpi::pop(); } mesh.get()->attach(); @@ -40,13 +42,13 @@ Mesh::Mesh(const Grid& grid, const eckit::Configuration& config): Mesh::Mesh(const Grid& grid, const grid::Partitioner& partitioner, const eckit::Configuration& config): Handle([&]() { - std::string mpi_comm = partitioner.mpi_comm(); - if(config.has("mpi_comm")) { + auto mpi_comm = partitioner.mpi_comm(); + if (config.has("mpi_comm")) { mpi_comm = config.getString("mpi_comm"); ATLAS_ASSERT(mpi_comm == partitioner.mpi_comm()); } mpi::Scope mpi_scope(mpi_comm); - auto meshgenerator = MeshGenerator{grid.meshgenerator()|config}; + auto meshgenerator = MeshGenerator{grid.meshgenerator() | config}; auto mesh = meshgenerator.generate(grid, partitioner); mesh.get()->attach(); return mesh.get(); @@ -57,7 +59,9 @@ Mesh::Mesh(const Grid& grid, const grid::Partitioner& partitioner, const eckit:: Mesh::Mesh(eckit::Stream& stream): Handle(new Implementation(stream)) {} -Mesh::operator bool() const { return get()->nodes().size() > 0; } +Mesh::operator bool() const { + return get()->nodes().size() > 0; +} //---------------------------------------------------------------------------------------------------------------------- diff --git a/src/atlas/mesh/actions/BuildEdges.cc b/src/atlas/mesh/actions/BuildEdges.cc index eb6461e91..08a568a17 100644 --- a/src/atlas/mesh/actions/BuildEdges.cc +++ b/src/atlas/mesh/actions/BuildEdges.cc @@ -380,7 +380,9 @@ void build_edges(Mesh& mesh, const eckit::Configuration& config) { for (int halo = 0; halo <= mesh_halo; ++halo) { edge_start = edge_end; - edge_end += (edge_halo_offsets[halo + 1] - edge_halo_offsets[halo]); + if (halo+1 < edge_halo_offsets.size()) { + edge_end += (edge_halo_offsets[halo + 1] - edge_halo_offsets[halo]); + } if (/*sort edges based on lowest node local index = */ sort_edges) { if (sorted_edge_nodes_data.empty()) { @@ -552,6 +554,11 @@ void build_edges(Mesh& mesh, const eckit::Configuration& config) { ss << "nb_edges_including_halo[" << i << "]"; mesh.metadata().set(ss.str(), nb_edges_including_halo[i]); } + for (int i = max_halo + 1; i <= mesh_halo; ++i) { + std::stringstream ss; + ss << "nb_edges_including_halo[" << i << "]"; + mesh.metadata().set(ss.str(), nb_edges_including_halo.back()); + } mesh.metadata().set("built_edges_for_halo", mesh_halo); diff --git a/src/atlas/meshgenerator/detail/CubedSphereDualMeshGenerator.cc b/src/atlas/meshgenerator/detail/CubedSphereDualMeshGenerator.cc index 953a84ed8..9026d92e8 100644 --- a/src/atlas/meshgenerator/detail/CubedSphereDualMeshGenerator.cc +++ b/src/atlas/meshgenerator/detail/CubedSphereDualMeshGenerator.cc @@ -63,9 +63,7 @@ CubedSphereDualMeshGenerator::CubedSphereDualMeshGenerator(const eckit::Parametr // Get partitioner. std::string partitioner; - try { p.get("partitioner.type", partitioner); } catch( std::exception& ) {} - p.get("partitioner", partitioner); - if( partitioner.size() ) { + if (p.get("partitioner", partitioner) && partitioner.size()) { options.set("partitioner", partitioner); } } diff --git a/src/atlas/meshgenerator/detail/CubedSphereMeshGenerator.cc b/src/atlas/meshgenerator/detail/CubedSphereMeshGenerator.cc index c2bd898ab..250d42ed0 100644 --- a/src/atlas/meshgenerator/detail/CubedSphereMeshGenerator.cc +++ b/src/atlas/meshgenerator/detail/CubedSphereMeshGenerator.cc @@ -73,9 +73,7 @@ CubedSphereMeshGenerator::CubedSphereMeshGenerator(const eckit::Parametrisation& // Get partitioner. std::string partitioner; - try { p.get("partitioner.type", partitioner); } catch( std::exception& ) {} - p.get("partitioner", partitioner); - if( partitioner.size() ) { + if (p.get("partitioner", partitioner) && partitioner.size()) { options.set("partitioner", partitioner); } } diff --git a/src/atlas/meshgenerator/detail/DelaunayMeshGenerator.cc b/src/atlas/meshgenerator/detail/DelaunayMeshGenerator.cc index 4ab849ba0..7004d6725 100644 --- a/src/atlas/meshgenerator/detail/DelaunayMeshGenerator.cc +++ b/src/atlas/meshgenerator/detail/DelaunayMeshGenerator.cc @@ -197,6 +197,9 @@ void DelaunayMeshGenerator::generate(const Grid& grid, const grid::Distribution& element_nodes_uncertainty.insert(g_node_connectivity(jelem,2)); element_uncertainty.insert(jelem); } + else if (elem_ownership != GHOST && elem_ownership != CERTAIN) { + ATLAS_NOTIMPLEMENTED; + } } // Log::info() << "element_uncertainty" << std::endl; // for( auto& jelem: element_uncertainty ) { @@ -291,6 +294,9 @@ void DelaunayMeshGenerator::generate(const Grid& grid, const grid::Distribution& else if (e2 != UNCERTAIN) { elem_part = e2; } + if (elem_part == UNCERTAIN) { + elem_part = elem_node_partition(jelem,0); + } if (elem_part == part_) { collect_element(jelem); } diff --git a/src/atlas/output/detail/GmshIO.cc b/src/atlas/output/detail/GmshIO.cc index 305051de4..810c30ec1 100644 --- a/src/atlas/output/detail/GmshIO.cc +++ b/src/atlas/output/detail/GmshIO.cc @@ -143,6 +143,7 @@ void write_level(std::ostream& out, GlobalIndex gidx, const array::LocalView data_vec; + data_vec.fill(static_cast(0)); for (idx_t n = 0; n < ndata; ++n) { if (include(n)) { for (idx_t v = 0; v < nvars; ++v) { @@ -158,6 +159,7 @@ void write_level(std::ostream& out, GlobalIndex gidx, const array::LocalView data_vec; + data_vec.fill(static_cast(0)); if (nvars == 4) { for (int n = 0; n < ndata; ++n) { if (include(n)) { diff --git a/src/atlas/projection/detail/MercatorProjection.cc b/src/atlas/projection/detail/MercatorProjection.cc index d3804523a..783ef74c8 100644 --- a/src/atlas/projection/detail/MercatorProjection.cc +++ b/src/atlas/projection/detail/MercatorProjection.cc @@ -10,6 +10,7 @@ #include #include +#include #include #include "eckit/config/Parametrisation.h" @@ -96,10 +97,18 @@ void MercatorProjectionT::lonlat2xy(double crd[]) const { rotation_.unrotate(crd); // then project - crd[XX] = k_radius_ * (D2R(normalise_mercator_(crd[LON]) - lon0_)); - crd[YY] = k_radius_ * 0.5 * std::log(t(crd[LAT])); - crd[XX] += false_easting_; - crd[YY] += false_northing_; + if (crd[LAT] >= 90. - 1e-3) { + crd[XX] = false_easting_; + crd[YY] = std::numeric_limits::infinity(); + } + else if (crd[LAT] <= -90. + 1e-3) { + crd[XX] = false_easting_; + crd[YY] = -std::numeric_limits::infinity(); + } + else { + crd[XX] = false_easting_ + k_radius_ * (D2R(normalise_mercator_(crd[LON]) - lon0_)); + crd[YY] = false_northing_ + k_radius_ * 0.5 * std::log(t(crd[LAT])); + } } template diff --git a/src/atlas/projection/detail/ProjectionImpl.cc b/src/atlas/projection/detail/ProjectionImpl.cc index 9b8264ac7..bf16eb28b 100644 --- a/src/atlas/projection/detail/ProjectionImpl.cc +++ b/src/atlas/projection/detail/ProjectionImpl.cc @@ -8,25 +8,23 @@ * nor does it submit to any jurisdiction. */ +#include "atlas/projection/detail/ProjectionImpl.h" + #include #include -#include +#include #include -#include +#include #include "eckit/types/FloatCompare.h" #include "eckit/utils/Hash.h" #include "eckit/utils/MD5.h" -#include "atlas/projection/detail/ProjectionImpl.h" - #include "atlas/projection/detail/ProjectionFactory.h" #include "atlas/runtime/Exception.h" #include "atlas/util/Config.h" -namespace atlas { -namespace projection { -namespace detail { +namespace atlas::projection::detail { // -------------------------------------------------------------------------------------------------------------------- @@ -53,31 +51,19 @@ struct DerivateBuilder : public ProjectionImpl::DerivateFactory { struct DerivateForwards final : ProjectionImpl::Derivate { using Derivate::Derivate; - PointLonLat d(PointXY P) const override { - PointLonLat A(xy2lonlat(P)); - PointLonLat B(xy2lonlat(PointXY::add(P, H_))); - return PointLonLat::div(PointLonLat::sub(B, A), normH_); - } + PointLonLat d(PointXY P) const override { return (xy2lonlat(P + H_) - xy2lonlat(P)) * invnH_; } }; struct DerivateBackwards final : ProjectionImpl::Derivate { using Derivate::Derivate; - PointLonLat d(PointXY P) const override { - PointLonLat A(xy2lonlat(PointXY::sub(P, H_))); - PointLonLat B(xy2lonlat(P)); - return PointLonLat::div(PointLonLat::sub(B, A), normH_); - } + PointLonLat d(PointXY P) const override { return (xy2lonlat(P) - xy2lonlat(P - H_)) * invnH_; } }; struct DerivateCentral final : ProjectionImpl::Derivate { DerivateCentral(const ProjectionImpl& p, PointXY A, PointXY B, double h, double refLongitude): - Derivate(p, A, B, h, refLongitude), H2_{PointXY::mul(H_, 0.5)} {} + Derivate(p, A, B, h, refLongitude), H2_{H_ * 0.5} {} const PointXY H2_; - PointLonLat d(PointXY P) const override { - PointLonLat A(xy2lonlat(PointXY::sub(P, H2_))); - PointLonLat B(xy2lonlat(PointXY::add(P, H2_))); - return PointLonLat::div(PointLonLat::sub(B, A), normH_); - } + PointLonLat d(PointXY P) const override { return (xy2lonlat(P + H2_) - xy2lonlat(P - H2_)) * invnH_; } }; } // namespace @@ -85,7 +71,7 @@ struct DerivateCentral final : ProjectionImpl::Derivate { ProjectionImpl::Derivate::Derivate(const ProjectionImpl& p, PointXY A, PointXY B, double h, double refLongitude): projection_(p), H_{PointXY::mul(PointXY::normalize(PointXY::sub(B, A)), h)}, - normH_(PointXY::norm(H_)), + invnH_(1. / PointXY::norm(H_)), refLongitude_(refLongitude) {} ProjectionImpl::Derivate::~Derivate() = default; @@ -107,7 +93,7 @@ ProjectionImpl::Derivate* ProjectionImpl::DerivateFactory::build(const std::stri return new DerivateDegenerate(p, A, B, h, refLongitude); } - auto factory = get(type); + auto* factory = get(type); return factory->make(p, A, B, h); } @@ -184,7 +170,7 @@ ProjectionImpl::Normalise::Normalise(const eckit::Parametrisation& p) { provided = true; } if (provided) { - normalise_.reset(new util::NormaliseLongitude(values_[0], values_[1])); + normalise_ = std::make_unique(values_[0], values_[1]); } } @@ -192,7 +178,7 @@ ProjectionImpl::Normalise::Normalise(double west) { values_.resize(2); values_[0] = west; values_[1] = values_[0] + 360.; - normalise_.reset(new util::NormaliseLongitude(values_[0], values_[1])); + normalise_ = std::make_unique(values_[0], values_[1]); } @@ -246,8 +232,8 @@ RectangularLonLatDomain ProjectionImpl::lonlatBoundingBox(const Domain& domain) ATLAS_ASSERT(rect); // use central longitude as absolute reference (keep points within +-180 longitude range) - const auto centre = lonlat({(rect.xmin() + rect.xmax()) / 2., (rect.ymin() + rect.ymax()) / 2.}); - + const PointXY centre_xy{(rect.xmin() + rect.xmax()) / 2., (rect.ymin() + rect.ymax()) / 2.}; + const auto centre_lon = lonlat(centre_xy).lon(); const std::string derivative = "central"; constexpr double h_deg = 0.5e-6; // precision to microdegrees @@ -256,15 +242,18 @@ RectangularLonLatDomain ProjectionImpl::lonlatBoundingBox(const Domain& domain) const double h = units() == "degrees" ? h_deg : h_meters; + // 1. determine box from projected corners - const std::vector corners{ - {rect.xmin(), rect.ymax()}, {rect.xmax(), rect.ymax()}, {rect.xmax(), rect.ymin()}, {rect.xmin(), rect.ymin()}}; + const std::pair segments[] = {{{rect.xmin(), rect.ymax()}, {rect.xmax(), rect.ymax()}}, + {{rect.xmax(), rect.ymax()}, {rect.xmax(), rect.ymin()}}, + {{rect.xmax(), rect.ymin()}, {rect.xmin(), rect.ymin()}}, + {{rect.xmin(), rect.ymin()}, {rect.xmin(), rect.ymax()}}}; BoundLonLat bounds; - for (auto& p : corners) { + for (const auto [p, dummy] : segments) { auto q = lonlat(p); - longitude_in_range(centre.lon(), q.lon()); + longitude_in_range(centre_lon, q.lon()); bounds.extend(q, PointLonLat{h_deg, h_deg}); } @@ -272,12 +261,12 @@ RectangularLonLatDomain ProjectionImpl::lonlatBoundingBox(const Domain& domain) // 2. locate latitude extrema by checking if poles are included (in the un-projected frame) and if not, find extrema // not at the corners by refining iteratively - for (size_t i = 0; i < corners.size(); ++i) { - if (!bounds.includesNorthPole() || !bounds.includesSouthPole()) { - PointXY A = corners[i]; - PointXY B = corners[(i + 1) % corners.size()]; + bounds.includesNorthPole(bounds.includesNorthPole() || rect.contains(xy(PointLonLat{0, 90 - h_deg}))); + bounds.includesSouthPole(bounds.includesSouthPole() || rect.contains(xy(PointLonLat{0, -90 + h_deg}))); - std::unique_ptr derivate(DerivateFactory::build(derivative, *this, A, B, h, centre.lon())); + for (auto [A, B] : segments) { + if (!bounds.includesNorthPole() || !bounds.includesSouthPole()) { + std::unique_ptr derivate(DerivateFactory::build(derivative, *this, A, B, h, centre_lon)); double dAdy = derivate->d(A).lat(); double dBdy = derivate->d(B).lat(); @@ -307,12 +296,9 @@ RectangularLonLatDomain ProjectionImpl::lonlatBoundingBox(const Domain& domain) // 3. locate longitude extrema not at the corners by refining iteratively - for (size_t i = 0; i < corners.size(); ++i) { + for (auto [A, B] : segments) { if (!bounds.crossesDateLine()) { - PointXY A = corners[i]; - PointXY B = corners[(i + 1) % corners.size()]; - - std::unique_ptr derivate(DerivateFactory::build(derivative, *this, A, B, h, centre.lon())); + std::unique_ptr derivate(DerivateFactory::build(derivative, *this, A, B, h, centre_lon)); double dAdx = derivate->d(A).lon(); double dBdx = derivate->d(B).lon(); @@ -410,6 +396,4 @@ void atlas__Projection__lonlat2xy(const ProjectionImpl* This, const double lon, } // extern "C" -} // namespace detail -} // namespace projection -} // namespace atlas +} // namespace atlas::projection::detail diff --git a/src/atlas/projection/detail/ProjectionImpl.h b/src/atlas/projection/detail/ProjectionImpl.h index 6f777a5ee..559707643 100644 --- a/src/atlas/projection/detail/ProjectionImpl.h +++ b/src/atlas/projection/detail/ProjectionImpl.h @@ -10,13 +10,11 @@ #pragma once -#include #include #include #include #include "atlas/projection/Jacobian.h" -#include "atlas/runtime/Exception.h" #include "atlas/util/Factory.h" #include "atlas/util/NormaliseLongitude.h" #include "atlas/util/Object.h" @@ -36,9 +34,7 @@ class Config; } } // namespace atlas -namespace atlas { -namespace projection { -namespace detail { +namespace atlas::projection::detail { class ProjectionImpl : public util::Object { public: @@ -49,8 +45,7 @@ class ProjectionImpl : public util::Object { static const ProjectionImpl* create(const eckit::Parametrisation& p); static const ProjectionImpl* create(const std::string& type, const eckit::Parametrisation& p); - ProjectionImpl() = default; - virtual ~ProjectionImpl() = default; // destructor should be virtual + ProjectionImpl() = default; virtual std::string type() const = 0; @@ -123,7 +118,7 @@ class ProjectionImpl : public util::Object { protected: const ProjectionImpl& projection_; const PointXY H_; - const double normH_; + const double invnH_; const double refLongitude_; PointLonLat xy2lonlat(const PointXY& p) const; }; @@ -190,10 +185,8 @@ class NotRotated { static std::string classNamePrefix() { return ""; } // deliberately empty static std::string typePrefix() { return ""; } // deliberately empty - void rotate(double*) const { /* do nothing */ - } - void unrotate(double*) const { /* do nothing */ - } + void rotate(double*) const { /* do nothing */ } + void unrotate(double*) const { /* do nothing */ } bool rotated() const { return false; } @@ -211,6 +204,4 @@ void atlas__Projection__xy2lonlat(const ProjectionImpl* This, const double x, co void atlas__Projection__lonlat2xy(const ProjectionImpl* This, const double lon, const double lat, double& x, double& y); } -} // namespace detail -} // namespace projection -} // namespace atlas +} // namespace atlas::projection::detail diff --git a/src/atlas/runtime/trace/TraceT.h b/src/atlas/runtime/trace/TraceT.h index efe98f45c..21def5df8 100644 --- a/src/atlas/runtime/trace/TraceT.h +++ b/src/atlas/runtime/trace/TraceT.h @@ -93,8 +93,8 @@ class TraceT { template inline std::string TraceT::formatTitle(const std::string& _title) { - std::string title = _title; - +(Barriers::state() ? " [b]" : "") + + std::string title = + _title + (Barriers::state() ? " [b]" : "") + (atlas_omp_get_num_threads() > 1 ? " @thread[" + std::to_string(atlas_omp_get_thread_num()) + "]" : ""); return title; } diff --git a/src/atlas/util/Config.cc b/src/atlas/util/Config.cc index 8ccf234c1..b54ba2318 100644 --- a/src/atlas/util/Config.cc +++ b/src/atlas/util/Config.cc @@ -59,13 +59,16 @@ Config::Config(std::istream& stream, const std::string&): eckit::LocalConfigurat Config::Config(const eckit::PathName& path): eckit::LocalConfiguration(yaml_from_path(path)) {} -Config Config::operator|(const Config& other) const { +Config Config::operator|(const eckit::Configuration& other) const { Config config(*this); config.set(other); return config; } -Config& Config::set(const eckit::LocalConfiguration& other) { +Config& Config::set(const eckit::Configuration& other) { +#if ATLAS_ECKIT_VERSION_AT_LEAST(1, 26, 0) || ATLAS_ECKIT_DEVELOP + eckit::LocalConfiguration::set(other); +#else eckit::Value& root = const_cast(get()); auto& other_root = other.get(); std::vector other_keys; @@ -73,12 +76,17 @@ Config& Config::set(const eckit::LocalConfiguration& other) { for (auto& key : other_keys) { root[key] = other_root[key]; } +#endif return *this; } Config& Config::remove(const std::string& name) { +#if ATLAS_ECKIT_VERSION_AT_LEAST(1, 26, 0) || ATLAS_ECKIT_DEVELOP + LocalConfiguration::remove(name); +#else eckit::Value& root = const_cast(get()); root.remove(name); +#endif return *this; } @@ -94,20 +102,22 @@ Config& Config::set(const std::string& name, const std::vector& values) bool Config::get(const std::string& name, std::vector& value) const { bool found = has(name); if (found) { - std::vector properties = getSubConfigurations(name); + auto properties = getSubConfigurations(name); value.resize(properties.size()); for (size_t i = 0; i < value.size(); ++i) { - value[i] = Config(properties[i]); + value[i].set(properties[i]); } } return found; } +#if ! (ATLAS_ECKIT_VERSION_AT_LEAST(1, 26, 0) || ATLAS_ECKIT_DEVELOP) std::vector Config::keys() const { std::vector result; eckit::fromValue(result, get().keys()); return result; } +#endif std::string Config::json(eckit::JSON::Formatting formatting) const { std::stringstream json; diff --git a/src/atlas/util/Config.h b/src/atlas/util/Config.h index 6c9e9a42b..5032a17f2 100644 --- a/src/atlas/util/Config.h +++ b/src/atlas/util/Config.h @@ -16,6 +16,8 @@ #include "eckit/config/LocalConfiguration.h" #include "eckit/log/JSON.h" +#include "atlas/library/config.h" + namespace eckit { class PathName; class Hash; @@ -61,14 +63,14 @@ class Config : public eckit::LocalConfiguration { Config operator()(const std::string& name, std::initializer_list&& value); // Overload operators to merge two Config objects. - Config operator|(const Config& other) const; + Config operator|(const eckit::Configuration& other) const; /// @brief Set a key-value parameter using eckit::LocalConfiguration::set; Config& set(const std::string& name, const std::vector&); - Config& set(const eckit::LocalConfiguration&); + Config& set(const eckit::Configuration&); template Config& set(const std::string& name, std::initializer_list&& value); @@ -80,7 +82,9 @@ class Config : public eckit::LocalConfiguration { using eckit::LocalConfiguration::get; bool get(const std::string& name, std::vector& value) const; +#if ! (ATLAS_ECKIT_VERSION_AT_LEAST(1, 26, 0) || ATLAS_ECKIT_DEVELOP) std::vector keys() const; +#endif std::string json(eckit::JSON::Formatting = eckit::JSON::Formatting::indent()) const; }; diff --git a/src/atlas/util/Geometry.cc b/src/atlas/util/Geometry.cc index 7d8e07797..a274d2ead 100644 --- a/src/atlas/util/Geometry.cc +++ b/src/atlas/util/Geometry.cc @@ -3,17 +3,18 @@ * * This software is licensed under the terms of the Apache Licence Version 2.0 * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. - * In applying this licence, ECMWF does not waive the privileges and immGeometryies + * In applying this licence, ECMWF does not waive the privileges and immunities * granted to it by virtue of its status as an intergovernmental organisation * nor does it submit to any jurisdiction. */ -#include "eckit/geometry/Point2.h" -#include "eckit/geometry/Point3.h" +#include "atlas/util/Geometry.h" + +#include #include "atlas/library/config.h" #include "atlas/runtime/Exception.h" -#include "atlas/util/Geometry.h" +#include "atlas/util/Constants.h" namespace atlas { @@ -21,79 +22,124 @@ namespace geometry { namespace detail { void GeometrySphere::lonlat2xyz(const Point2& lonlat, Point3& xyz) const { #if ATLAS_ECKIT_VERSION_AT_LEAST(1, 24, 0) - Sphere::convertSphericalToCartesian(radius_, lonlat, xyz, 0., true); + Sphere::convertSphericalToCartesian(radius_, lonlat, xyz, 0., true); #else - Sphere::convertSphericalToCartesian(radius_, lonlat, xyz); + Sphere::convertSphericalToCartesian(radius_, lonlat, xyz); #endif } void GeometrySphere::xyz2lonlat(const Point3& xyz, Point2& lonlat) const { - Sphere::convertCartesianToSpherical(radius_, xyz, lonlat); + Sphere::convertCartesianToSpherical(radius_, xyz, lonlat); } -} // namespace detail -} // namespace geometry +/// @brief Calculate great-cricle course between points +/// +/// @details Calculates the direction (clockwise from north) of a great-circle +/// arc between lonLat1 and lonLat2. Returns the direction of the arc +/// at lonLat1 (first) and lonLat2 (second). Angle is normalised to the +/// range of atan2 (usually (-180, 180]). All input and output values +/// are in units of degrees. +/// @ref https://en.wikipedia.org/wiki/Great-circle_navigation +/// +std::pair greatCircleCourse(const Point2& lonLat1, + const Point2& lonLat2) { + if (lonLat1 == lonLat2) { + return std::make_pair(0., 0.); + } + + const auto lambda1 = lonLat1[0] * util::Constants::degreesToRadians(); + const auto lambda2 = lonLat2[0] * util::Constants::degreesToRadians(); + const auto phi1 = lonLat1[1] * util::Constants::degreesToRadians(); + const auto phi2 = lonLat2[1] * util::Constants::degreesToRadians(); + + const auto sinLambda12 = std::sin(lambda2 - lambda1); + const auto cosLambda12 = std::cos(lambda2 - lambda1); + const auto sinPhi1 = std::sin(phi1); + const auto sinPhi2 = std::sin(phi2); + const auto cosPhi1 = std::cos(phi1); + const auto cosPhi2 = std::cos(phi2); + + const auto alpha1 = + std::atan2(cosPhi2 * sinLambda12, + cosPhi1 * sinPhi2 - sinPhi1 * cosPhi2 * cosLambda12); + + const auto alpha2 = + std::atan2(cosPhi1 * sinLambda12, + -cosPhi2 * sinPhi1 + sinPhi2 * cosPhi1 * cosLambda12); + + return std::make_pair(alpha1 * util::Constants::radiansToDegrees(), + alpha2 * util::Constants::radiansToDegrees()); +}; + +} // namespace detail +} // namespace geometry extern "C" { // ------------------------------------------------------------------ // C wrapper interfaces to C++ routines Geometry::Implementation* atlas__Geometry__new_name(const char* name) { - Geometry::Implementation* geometry; - { - Geometry handle{std::string{name}}; - geometry = handle.get(); - geometry->attach(); - } - geometry->detach(); - return geometry; + Geometry::Implementation* geometry; + { + Geometry handle{std::string{name}}; + geometry = handle.get(); + geometry->attach(); + } + geometry->detach(); + return geometry; } -geometry::detail::GeometryBase* atlas__Geometry__new_radius(const double radius) { - Geometry::Implementation* geometry; - { - Geometry handle{radius}; - geometry = handle.get(); - geometry->attach(); - } - geometry->detach(); - return geometry; +geometry::detail::GeometryBase* atlas__Geometry__new_radius( + const double radius) { + Geometry::Implementation* geometry; + { + Geometry handle{radius}; + geometry = handle.get(); + geometry->attach(); + } + geometry->detach(); + return geometry; } void atlas__Geometry__delete(Geometry::Implementation* This) { - ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry"); - delete This; + ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry"); + delete This; } -void atlas__Geometry__xyz2lonlat(Geometry::Implementation* This, const double x, const double y, const double z, - double& lon, double& lat) { - ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry"); - PointLonLat lonlat; - This->xyz2lonlat(PointXYZ{x, y, z}, lonlat); - lon = lonlat.lon(); - lat = lonlat.lat(); +void atlas__Geometry__xyz2lonlat(Geometry::Implementation* This, const double x, + const double y, const double z, double& lon, + double& lat) { + ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry"); + PointLonLat lonlat; + This->xyz2lonlat(PointXYZ{x, y, z}, lonlat); + lon = lonlat.lon(); + lat = lonlat.lat(); } -void atlas__Geometry__lonlat2xyz(Geometry::Implementation* This, const double lon, const double lat, double& x, +void atlas__Geometry__lonlat2xyz(Geometry::Implementation* This, + const double lon, const double lat, double& x, double& y, double& z) { - ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry"); - PointXYZ xyz; - This->lonlat2xyz(PointLonLat{lon, lat}, xyz); - x = xyz.x(); - y = xyz.y(); - z = xyz.z(); + ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry"); + PointXYZ xyz; + This->lonlat2xyz(PointLonLat{lon, lat}, xyz); + x = xyz.x(); + y = xyz.y(); + z = xyz.z(); } -double atlas__Geometry__distance_lonlat(Geometry::Implementation* This, const double lon1, const double lat1, +double atlas__Geometry__distance_lonlat(Geometry::Implementation* This, + const double lon1, const double lat1, const double lon2, const double lat2) { - ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry"); - return This->distance(PointLonLat{lon1, lat1}, PointLonLat{lon2, lat2}); + ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry"); + return This->distance(PointLonLat{lon1, lat1}, PointLonLat{lon2, lat2}); } -double atlas__Geometry__distance_xyz(Geometry::Implementation* This, const double x1, const double y1, const double z1, - const double x2, const double y2, const double z2) { - ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry"); - return This->distance(PointXYZ{x1, y1, z1}, PointXYZ{x2, y2, z2}); +double atlas__Geometry__distance_xyz(Geometry::Implementation* This, + const double x1, const double y1, + const double z1, const double x2, + const double y2, const double z2) { + ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry"); + return This->distance(PointXYZ{x1, y1, z1}, PointXYZ{x2, y2, z2}); } double atlas__Geometry__radius(Geometry::Implementation* This) { - ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry"); - return This->radius(); + ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry"); + return This->radius(); } double atlas__Geometry__area(Geometry::Implementation* This) { - ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry"); - return This->area(); + ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_Geometry"); + return This->area(); } } // ------------------------------------------------------------------ diff --git a/src/atlas/util/Geometry.h b/src/atlas/util/Geometry.h index 3880645e8..90b868ce8 100644 --- a/src/atlas/util/Geometry.h +++ b/src/atlas/util/Geometry.h @@ -13,6 +13,7 @@ #include #include #include +#include #include "atlas/runtime/Exception.h" #include "atlas/util/Earth.h" @@ -27,6 +28,20 @@ namespace geometry { namespace detail { +// TODO: move greatCircleCourse to eckit::geometry::Sphere + +/// @brief Calculate great-cricle course between points +/// +/// @details Calculates the direction (clockwise from north) of a great-circle +/// arc between lonLat1 and lonLat2. Returns the direction of the arc +/// at lonLat1 (first) and lonLat2 (second). Angle is normalised to the +/// range of atan2 (usually (-180, 180]). All input and output values +/// are in units of degrees. +/// @ref https://en.wikipedia.org/wiki/Great-circle_navigation +std::pair greatCircleCourse(const Point2& lonLat1, const Point2& lonLat2); + +//------------------------------------------------------------------------------------------------------ + class GeometryBase : public util::Object { public: virtual ~GeometryBase() = default; @@ -36,6 +51,7 @@ class GeometryBase : public util::Object { virtual double distance(const Point3& p1, const Point3& p2) const = 0; virtual double radius() const = 0; virtual double area() const = 0; + virtual std::pair greatCircleCourse(const Point2& p1, const Point2& p2) const = 0; Point3 xyz(const Point2& lonlat) const { Point3 xyz; @@ -64,6 +80,9 @@ class GeometrySphereT : public GeometryBase { double distance(const Point3& p1, const Point3& p2) const override { return SphereT::distance(p1, p2); } double radius() const override { return SphereT::radius(); } double area() const override { return SphereT::area(); } + std::pair greatCircleCourse(const Point2& p1, const Point2& p2) const override { + return atlas::geometry::detail::greatCircleCourse(p1, p2); + } }; class GeometrySphere : public GeometryBase { @@ -77,6 +96,9 @@ class GeometrySphere : public GeometryBase { double distance(const Point3& p1, const Point3& p2) const override { return Sphere::distance(radius_, p1, p2); } double radius() const override { return radius_; } double area() const override { return Sphere::area(radius_); } + std::pair greatCircleCourse(const Point2& p1, const Point2& p2) const override { + return atlas::geometry::detail::greatCircleCourse(p1, p2); + } private: double radius_; @@ -107,6 +129,7 @@ class Geometry : public util::ObjectHandle { double distance(const Point3& p1, const Point3& p2) const { return get()->distance(p1, p2); } double radius() const { return get()->radius(); } double area() const { return get()->area(); } + std::pair greatCircleCourse(const Point2& p1, const Point2& p2) const { return get()->greatCircleCourse(p1, p2); } protected: template diff --git a/src/atlas/util/Metadata.cc b/src/atlas/util/Metadata.cc index c24e6c897..583332bec 100644 --- a/src/atlas/util/Metadata.cc +++ b/src/atlas/util/Metadata.cc @@ -119,7 +119,10 @@ void Metadata::broadcast(Metadata& dest, idx_t root) const { } -Metadata& Metadata::set(const eckit::LocalConfiguration& other) { +Metadata& Metadata::set(const eckit::Configuration& other) { +#if ATLAS_ECKIT_VERSION_AT_LEAST(1, 26, 0) || ATLAS_ECKIT_DEVELOP + LocalConfiguration::set(other); +#else eckit::Value& root = const_cast(get()); auto& other_root = other.get(); std::vector other_keys; @@ -127,6 +130,7 @@ Metadata& Metadata::set(const eckit::LocalConfiguration& other) { for (auto& key : other_keys) { root[key] = other_root[key]; } +#endif return *this; } diff --git a/src/atlas/util/Metadata.h b/src/atlas/util/Metadata.h index c697e13fa..f4ce34eef 100644 --- a/src/atlas/util/Metadata.h +++ b/src/atlas/util/Metadata.h @@ -35,7 +35,7 @@ class Metadata : public eckit::LocalConfiguration { return *this; } - Metadata& set(const eckit::LocalConfiguration&); + Metadata& set(const eckit::Configuration&); using eckit::LocalConfiguration::get; diff --git a/src/tests/array/test_array_foreach.cc b/src/tests/array/test_array_foreach.cc index f8c5b0b4a..303da5dd6 100644 --- a/src/tests/array/test_array_foreach.cc +++ b/src/tests/array/test_array_foreach.cc @@ -1,5 +1,5 @@ /* - * (C) Crown Copyright 2023 Met Office + * (C) Crown Copyright 2024 Met Office * * This software is licensed under the terms of the Apache Licence Version 2.0 * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. @@ -7,13 +7,13 @@ #include #include +#include #include "atlas/array.h" #include "atlas/array/MakeView.h" #include "atlas/array/helpers/ArrayForEach.h" #include "atlas/array/helpers/ArraySlicer.h" #include "atlas/util/Config.h" - #include "tests/AtlasTestEnvironment.h" using namespace atlas::array; @@ -207,6 +207,82 @@ CASE("test_array_foreach_3_views") { EXPECT_EQ(count, 60); } +CASE("test_array_foreach_integer_sequence") { + + const auto arr1 = ArrayT(2, 3); + const auto view1 = make_view(arr1); + + const auto arr2 = ArrayT(2, 3, 4); + const auto view2 = make_view(arr2); + + const auto arr3 = ArrayT(2, 3, 4, 5); + const auto view3 = make_view(arr3); + + const auto none = std::integer_sequence{}; + const auto zero = std::integer_sequence{}; + const auto one = std::integer_sequence{}; + const auto zeroOneTwoThree = std::make_integer_sequence{}; + + + // Test slice shapes. + + const auto loopFunctorDimNone = [](auto&& slice1, auto&& slice2, + auto&& slice3) { + EXPECT_EQ(slice1.rank(), 2); + EXPECT_EQ(slice1.shape(0), 2); + EXPECT_EQ(slice1.shape(1), 3); + + EXPECT_EQ(slice2.rank(), 3); + EXPECT_EQ(slice2.shape(0), 2); + EXPECT_EQ(slice2.shape(1), 3); + EXPECT_EQ(slice2.shape(2), 4); + + EXPECT_EQ(slice3.rank(), 4); + EXPECT_EQ(slice3.shape(0), 2); + EXPECT_EQ(slice3.shape(1), 3); + EXPECT_EQ(slice3.shape(2), 4); + EXPECT_EQ(slice3.shape(3), 5); + }; + // No iterations. Apply functor directly to unsliced ArrayViews. + arrayForEachDim(none, std::tie(view1, view2, view3), loopFunctorDimNone); + + const auto loopFunctorDim0 = [](auto&& slice1, auto&& slice2, auto&& slice3) { + EXPECT_EQ(slice1.rank(), 1); + EXPECT_EQ(slice1.shape(0), 3); + + EXPECT_EQ(slice2.rank(), 2); + EXPECT_EQ(slice2.shape(0), 3); + EXPECT_EQ(slice2.shape(1), 4); + + EXPECT_EQ(slice3.rank(), 3); + EXPECT_EQ(slice3.shape(0), 3); + EXPECT_EQ(slice3.shape(1), 4); + EXPECT_EQ(slice3.shape(2), 5); + }; + arrayForEachDim(zero, std::tie(view1, view2, view3), loopFunctorDim0); + + const auto loopFunctorDim1 = [](auto&& slice1, auto&& slice2, auto&& slice3) { + EXPECT_EQ(slice1.rank(), 1); + EXPECT_EQ(slice1.shape(0), 2); + + EXPECT_EQ(slice2.rank(), 2); + EXPECT_EQ(slice2.shape(0), 2); + EXPECT_EQ(slice2.shape(1), 4); + + EXPECT_EQ(slice3.rank(), 3); + EXPECT_EQ(slice3.shape(0), 2); + EXPECT_EQ(slice3.shape(1), 4); + EXPECT_EQ(slice3.shape(2), 5); + }; + arrayForEachDim(one, std::tie(view1, view2, view3), loopFunctorDim1); + + // Test that slice resolves to double. + + const auto loopFunctorDimAll = [](auto&& slice3) { + static_assert(std::is_convertible_v); + }; + arrayForEachDim(zeroOneTwoThree, std::tie(view3), loopFunctorDimAll); +} CASE("test_array_foreach_forwarding") { diff --git a/src/tests/interpolation/CMakeLists.txt b/src/tests/interpolation/CMakeLists.txt index 596729209..29d220590 100644 --- a/src/tests/interpolation/CMakeLists.txt +++ b/src/tests/interpolation/CMakeLists.txt @@ -122,4 +122,12 @@ ecbuild_add_test( TARGET atlas_test_interpolation_cubedsphere ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} ) -endif() \ No newline at end of file +ecbuild_add_test( TARGET atlas_test_interpolation_spherical_vector + OMP 4 + SOURCES test_interpolation_spherical_vector.cc + LIBS atlas + ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} + CONDITION atlas_HAVE_EIGEN +) + +endif() diff --git a/src/tests/interpolation/test_interpolation_cubedsphere.cc b/src/tests/interpolation/test_interpolation_cubedsphere.cc index b2d446ebf..a6d930983 100644 --- a/src/tests/interpolation/test_interpolation_cubedsphere.cc +++ b/src/tests/interpolation/test_interpolation_cubedsphere.cc @@ -60,7 +60,7 @@ void gmshOutput(const std::string& fileName, const FieldSet& fieldSet) { // Return (u, v) field with vortex_rollup as the streamfunction. // This has no physical significance, but it makes a nice swirly field. -std::pair vortexField(double lon, double lat) { +std::pair vortexHorizontal(double lon, double lat) { // set hLon and hLat step size. const double hLon = 0.0001; @@ -222,7 +222,7 @@ CASE("cubedsphere_wind_interpolation") { const auto ll = PointLonLat(lonlat(idx, LON), lonlat(idx, LAT)); // Set (u, v) wind - std::tie(u(idx), v(idx)) = vortexField(ll.lon(), ll.lat()); + std::tie(u(idx), v(idx)) = vortexHorizontal(ll.lon(), ll.lat()); // Get wind transform jacobian. auto jac = windTransform(ll, t); @@ -282,7 +282,7 @@ CASE("cubedsphere_wind_interpolation") { std::tie(u(idx), v(idx)) = matMul(jac, vAlpha(idx), vBeta(idx)); // Get error. - const auto uvTarget = vortexField(ll.lon(), ll.lat()); + const auto uvTarget = vortexHorizontal(ll.lon(), ll.lat()); error0(idx) = Point2::distance(Point2(uvTarget.first, uvTarget.second), Point2(uOrig(idx), vOrig(idx))); error1(idx) = Point2::distance(Point2(uvTarget.first, uvTarget.second), Point2(u(idx), v(idx))); diff --git a/src/tests/interpolation/test_interpolation_spherical_vector.cc b/src/tests/interpolation/test_interpolation_spherical_vector.cc new file mode 100644 index 000000000..4717c5a19 --- /dev/null +++ b/src/tests/interpolation/test_interpolation_spherical_vector.cc @@ -0,0 +1,353 @@ +/* + * (C) Crown Copyright 2024 Met Office + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#include +#include +#include + +#include "atlas/array.h" +#include "atlas/array/helpers/ArrayForEach.h" +#include "atlas/field.h" +#include "atlas/functionspace.h" +#include "atlas/grid.h" +#include "atlas/interpolation.h" +#include "atlas/mesh.h" +#include "atlas/meshgenerator.h" +#include "atlas/output/Gmsh.h" +#include "atlas/util/Config.h" +#include "atlas/util/Constants.h" +#include "atlas/util/Point.h" +#include "atlas/util/function/VortexRollup.h" +#include "tests/AtlasTestEnvironment.h" + +namespace atlas { +namespace test { + +using namespace atlas::util; +using namespace atlas::array::helpers; + +constexpr auto Rank2dField = 2; +constexpr auto Rank3dField = 3; + +// Return (u, v) field with vortex_rollup as the streamfunction. +// This has no physical significance, but it makes a nice swirly field. +std::pair vortexHorizontal(double lon, double lat) { + // set hLon and hLat step size. + const double hLon = 0.0001; + const double hLat = 0.0001; + + // Get finite differences. + const double u = (function::vortex_rollup(lon, lat + 0.5 * hLat, 0.1) - + function::vortex_rollup(lon, lat - 0.5 * hLat, 0.1)) / + hLat; + + const double v = -(function::vortex_rollup(lon + 0.5 * hLon, lat, 0.1) - + function::vortex_rollup(lon - 0.5 * hLon, lat, 0.1)) / + (hLon * std::cos(lat * util::Constants::degreesToRadians())); + + return std::make_pair(u, v); +} + +double vortexVertical(double lon, double lat) { + return function::vortex_rollup(lon, lat, 0.1); +} + +void gmshOutput(const std::string& fileName, const FieldSet& fieldSet) { + const auto functionSpace = fieldSet[0].functionspace(); + const auto structuredColums = functionspace::StructuredColumns(functionSpace); + const auto nodeColumns = functionspace::NodeColumns(functionSpace); + const auto mesh = + structuredColums ? Mesh(structuredColums.grid()) : nodeColumns.mesh(); + + const auto gmshConfig = Config("coordinates", "xyz") | Config("ghost", true) | + Config("info", true); + const auto gmsh = output::Gmsh(fileName, gmshConfig); + + gmsh.write(mesh); + gmsh.write(fieldSet, functionSpace); +} + +// Helper function to generate a NodeColumns functionspace +const auto generateNodeColums(const std::string& gridName, + const std::string& meshName) { + const auto grid = Grid(gridName); + const auto mesh = MeshGenerator(meshName).generate(grid); + return functionspace::NodeColumns(mesh); +} + +// Helper struct to key different Functionspaces to strings +struct FunctionSpaceFixtures { + static const FunctionSpace& get(const std::string& fixture) { + static const auto functionSpaces = + std::map{ + {"cubedsphere_mesh", + generateNodeColums("CS-LFR-48", "cubedsphere_dual")}, + {"gaussian_mesh", generateNodeColums("O48", "structured")}, + {"structured_columns", + functionspace::StructuredColumns(Grid("O48"), option::halo(1))}, + {"structured_columns_lowres", + functionspace::StructuredColumns(Grid("O24"), option::halo(1))}, + {"structured_columns_hires", + functionspace::StructuredColumns(Grid("O96"), option::halo(1))}}; + return functionSpaces.at(fixture); + } +}; + +// Helper struct to key different grid configs to strings +struct FieldSpecFixtures { + static const Config& get(const std::string& fixture) { + static const auto fieldSpecs = std::map{ + {"2vector", option::name("test field") | option::variables(2) | + option::type("vector")}, + {"3vector", option::name("test field") | option::variables(3) | + option::type("vector")}}; + return fieldSpecs.at(fixture); + } +}; + +// Helper stcut to key different interpolation schemes to strings +struct InterpSchemeFixtures { + static const Config& get(const std::string& fixture) { + static const auto cubedsphereBilinear = + option::type("cubedsphere-bilinear"); + static const auto finiteElement = option::type("finite-element"); + static const auto structuredLinear = + option::type("structured-linear2D") | option::halo(1); + + static const auto sphericalVector = + option::type("spherical-vector") | Config("adjoint", true); + + static const auto interpSchemes = std::map{ + {"cubedsphere_bilinear", cubedsphereBilinear}, + {"finite_element", finiteElement}, + {"structured_linear", structuredLinear}, + {"cubedsphere_bilinear_spherical", + sphericalVector | Config("scheme", cubedsphereBilinear)}, + {"finite_element_spherical", + sphericalVector | Config("scheme", finiteElement)}, + {"structured_linear_spherical", + sphericalVector | Config("scheme", structuredLinear)}}; + return interpSchemes.at(fixture); + } +}; + +template +double dotProduct(const array::ArrayView& a, + const array::ArrayView& b) { + auto dotProd = 0.; + arrayForEachDim(std::make_integer_sequence{}, std::tie(a, b), + [&](const double& aElem, const double& bElem) { + dotProd += aElem * bElem; + }); + return dotProd; +} + +template +int countNans(const array::ArrayView& view) { + auto nNans = 0; + arrayForEachDim(std::make_integer_sequence{}, std::tie(view), + [&](const double& viewElem) { + if (!std::isfinite(viewElem)) { + ++nNans; + } + }); + return nNans; +} + +template +void testInterpolation(const Config& config) { + const auto& sourceFunctionSpace = + FunctionSpaceFixtures::get(config.getString("source_fixture")); + const auto& targetFunctionSpace = + FunctionSpaceFixtures::get(config.getString("target_fixture")); + + auto sourceFieldSet = FieldSet{}; + auto targetFieldSet = FieldSet{}; + + const auto sourceLonLat = + array::make_view(sourceFunctionSpace.lonlat()); + const auto targetLonLat = + array::make_view(targetFunctionSpace.lonlat()); + + auto fieldSpec = + FieldSpecFixtures::get(config.getString("field_spec_fixture")); + if constexpr (Rank == 3) { + fieldSpec.set("levels", 2); + } + + auto sourceField = + sourceFieldSet.add(sourceFunctionSpace.createField(fieldSpec)); + auto targetField = + targetFieldSet.add(targetFunctionSpace.createField(fieldSpec)); + + auto sourceView = array::make_view(sourceField); + auto targetView = array::make_view(targetField); + + ArrayForEach<0>::apply( + std::tie(sourceLonLat, sourceView), + [](auto&& lonLat, auto&& sourceColumn) { + const auto setElems = [&](auto&& sourceElem) { + std::tie(sourceElem(0), sourceElem(1)) = + vortexHorizontal(lonLat(0), lonLat(1)); + if (sourceElem.size() == 3) { + sourceElem(2) = vortexVertical(lonLat(0), lonLat(1)); + } + }; + if constexpr (Rank == 2) { + setElems(sourceColumn); + } else if constexpr (Rank == 3) { + ArrayForEach<0>::apply(std::tie(sourceColumn), setElems); + } + }); + + const auto interp = Interpolation( + InterpSchemeFixtures::get(config.getString("interp_fixture")), + sourceFunctionSpace, targetFunctionSpace); + + interp.execute(sourceFieldSet, targetFieldSet); + targetFieldSet.haloExchange(); + + auto errorFieldSpec = fieldSpec; + errorFieldSpec.remove("variables").set("name", "error field"); + + auto errorView = array::make_view(targetFieldSet.add( + targetFunctionSpace.createField(errorFieldSpec))); + errorView.assign(0.); + + auto maxError = 0.; + ArrayForEach<0>::apply( + std::tie(targetLonLat, targetView, errorView), + [&](auto&& lonLat, auto&& targetColumn, auto&& errorColumn) { + const auto calcError = [&](auto&& targetElem, auto&& errorElem) { + auto trueValue = std::vector(targetElem.size()); + std::tie(trueValue[0], trueValue[1]) = + vortexHorizontal(lonLat(0), lonLat(1)); + if (targetElem.size() == 3) { + trueValue[2] = vortexVertical(lonLat(0), lonLat(1)); + } + + auto errorSqrd = 0.; + for (auto k = 0; k < targetElem.size(); ++k) { + errorSqrd += + (targetElem(k) - trueValue[k]) * (targetElem(k) - trueValue[k]); + } + + errorElem = std::sqrt(errorSqrd); + maxError = std::max(maxError, static_cast(errorElem)); + }; + + if constexpr (Rank == 2) { + calcError(targetColumn, errorColumn); + } else if constexpr (Rank == 3) { + ArrayForEach<0>::apply(std::tie(targetColumn, errorColumn), + calcError); + } + }); + + EXPECT_APPROX_EQ(maxError, 0., config.getDouble("tol")); + + gmshOutput(config.getString("file_id") + "_source.msh", sourceFieldSet); + gmshOutput(config.getString("file_id") + "_target.msh", targetFieldSet); + + // Adjoint test + auto targetAdjoint = targetFunctionSpace.createField(fieldSpec); + auto targetAdjointView = array::make_view(targetAdjoint); + targetAdjoint.array().copy(targetField); + targetAdjoint.adjointHaloExchange(); + + auto sourceAdjoint = sourceFunctionSpace.createField(fieldSpec); + auto sourceAdjointView = array::make_view(sourceAdjoint); + sourceAdjointView.assign(0.); + + sourceAdjoint.set_dirty(false); + interp.execute_adjoint(sourceAdjoint, targetAdjoint); + + // Check fields for nans or +/-inf + EXPECT_EQ(countNans(targetView), 0); + EXPECT_EQ(countNans(sourceView), 0); + EXPECT_EQ(countNans(targetAdjointView), 0); + EXPECT_EQ(countNans(sourceAdjointView), 0); + + constexpr auto tinyNum = 1e-13; + const auto targetDotTarget = dotProduct(targetView, targetView); + const auto sourceDotSourceAdjoint = dotProduct(sourceView, sourceAdjointView); + const auto dotProdRatio = targetDotTarget / sourceDotSourceAdjoint; + EXPECT_APPROX_EQ(dotProdRatio, 1., tinyNum); +} + +CASE("cubed sphere vector interpolation (3d-field, 2-vector)") { + const auto config = + Config("source_fixture", "cubedsphere_mesh") + .set("target_fixture", "gaussian_mesh") + .set("field_spec_fixture", "2vector") + .set("interp_fixture", "cubedsphere_bilinear_spherical") + .set("file_id", "spherical_vector_cs2") + .set("tol", 0.00018); + + testInterpolation((config)); +} + +CASE("cubed sphere vector interpolation (3d-field, 3-vector)") { + const auto config = + Config("source_fixture", "cubedsphere_mesh") + .set("target_fixture", "gaussian_mesh") + .set("field_spec_fixture", "3vector") + .set("interp_fixture", "cubedsphere_bilinear_spherical") + .set("file_id", "spherical_vector_cs3") + .set("tol", 0.00096); + + testInterpolation((config)); +} + +CASE("finite element vector interpolation (2d-field, 2-vector)") { + const auto config = Config("source_fixture", "gaussian_mesh") + .set("target_fixture", "cubedsphere_mesh") + .set("field_spec_fixture", "2vector") + .set("interp_fixture", "finite_element_spherical") + .set("file_id", "spherical_vector_fe") + .set("tol", 0.00015); + + testInterpolation((config)); +} + +CASE("structured columns vector interpolation (2d-field, 2-vector)") { + const auto config = Config("source_fixture", "structured_columns") + .set("target_fixture", "cubedsphere_mesh") + .set("field_spec_fixture", "2vector") + .set("interp_fixture", "structured_linear_spherical") + .set("file_id", "spherical_vector_sc") + .set("tol", 0.00017); + + testInterpolation((config)); +} + +CASE("structured columns vector interpolation (2d-field, 2-vector, low-res)") { + const auto config = Config("source_fixture", "structured_columns_lowres") + .set("target_fixture", "gaussian_mesh") + .set("field_spec_fixture", "2vector") + .set("interp_fixture", "structured_linear_spherical") + .set("file_id", "spherical_vector_sc_lr") + .set("tol", 0.00056); + + testInterpolation((config)); +} + +CASE("structured columns vector interpolation (2d-field, 2-vector, hi-res)") { + const auto config = Config("source_fixture", "structured_columns_hires") + .set("target_fixture", "gaussian_mesh") + .set("field_spec_fixture", "2vector") + .set("interp_fixture", "structured_linear_spherical") + .set("file_id", "spherical_vector_sc_hr") + .set("tol", 0.000044); + + testInterpolation((config)); +} + +} // namespace test +} // namespace atlas + +int main(int argc, char** argv) { return atlas::test::run(argc, argv); } diff --git a/src/tests/util/CMakeLists.txt b/src/tests/util/CMakeLists.txt index 4e6e0daf0..47d8b6def 100644 --- a/src/tests/util/CMakeLists.txt +++ b/src/tests/util/CMakeLists.txt @@ -88,3 +88,9 @@ ecbuild_add_test( TARGET atlas_test_convexsphericalpolygon ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} ) +ecbuild_add_test( TARGET atlas_test_unitsphere + SOURCES test_unitsphere.cc + LIBS atlas + ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} +) + diff --git a/src/tests/util/test_unitsphere.cc b/src/tests/util/test_unitsphere.cc new file mode 100644 index 000000000..8bc249c6a --- /dev/null +++ b/src/tests/util/test_unitsphere.cc @@ -0,0 +1,57 @@ +/* + * (C) Crown Copyright 2023 Met Office + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ + +#include "atlas/util/Geometry.h" +#include "atlas/util/Point.h" +#include "tests/AtlasTestEnvironment.h" + +namespace atlas { +namespace test { + +using namespace atlas::util; + +CASE("great-circle course") { + geometry::UnitSphere g; + + const auto point1 = PointLonLat(-71.6, -33.0); // Valparaiso + const auto point2 = PointLonLat(121.8, 31.4); // Shanghai + const auto point3 = PointLonLat(0., 89.); + const auto point4 = PointLonLat(180., 89.); + const auto point5 = PointLonLat(0., 90.); + const auto point6 = PointLonLat(180., 90.); + + const auto targetCourse1 = -94.41; + const auto targetCourse2 = -78.42; + const auto targetCourse3 = 0.; + const auto targetCourse4 = 180.; + const auto targetCourse5 = 0.; + const auto targetCourse6 = 180.; + const auto targetCourse7 = 0.; + const auto targetCourse8 = 0.; + + + const auto [course1, course2] = g.greatCircleCourse(point1, point2); + const auto [course3, course4] = g.greatCircleCourse(point3, point4); + + // Colocated points on pole. + const auto [course5, course6] = g.greatCircleCourse(point5, point6); + const auto [course7, course8] = g.greatCircleCourse(point5, point5); + + + EXPECT_APPROX_EQ(course1, targetCourse1, 0.01); + EXPECT_APPROX_EQ(course2, targetCourse2, 0.01); + EXPECT_APPROX_EQ(course3, targetCourse3, 1.e-14); + EXPECT_APPROX_EQ(std::abs(course4), targetCourse4, 1.e-14); + EXPECT_APPROX_EQ(std::abs(course5), targetCourse5, 1.e-14); + EXPECT_APPROX_EQ(std::abs(course6), targetCourse6, 1.e-14); + EXPECT_APPROX_EQ(std::abs(course7), targetCourse7, 1.e-14); + EXPECT_APPROX_EQ(std::abs(course8), targetCourse8, 1.e-14); +} +} // namespace test +} // namespace atlas + +int main(int argc, char** argv) { return atlas::test::run(argc, argv); } diff --git a/tools/install-mpi.sh b/tools/install-mpi.sh index 03a5e3a35..c5c1847ae 100755 --- a/tools/install-mpi.sh +++ b/tools/install-mpi.sh @@ -62,6 +62,8 @@ case "$os" in openmpi) brew ls --versions openmpi || brew install openmpi echo "localhost slots=72" >> $(brew --prefix)/etc/openmpi-default-hostfile + echo "localhost slots=72" >> $(brew --prefix)/etc/prte-default-hostfile + # workaround for open-mpi/omp#7516 echo "setting the mca gds to hash..." echo "gds = hash" >> $(brew --prefix)/etc/pmix-mca-params.conf