diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index 78ff651a1..0bddfc529 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -91,6 +91,7 @@ jobs: compiler_cxx: ~ compiler_fc: gfortran-11 caching: true + cmake_options: -DMPI_SLOTS=4 runs-on: ${{ matrix.os }} steps: @@ -106,6 +107,7 @@ jobs: if [[ "${{ matrix.os }}" =~ macos ]]; then brew install ninja + brew install libomp else sudo apt-get update sudo apt-get install ninja-build diff --git a/CHANGELOG.md b/CHANGELOG.md index e28da8072..3a89a3b52 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,11 +7,32 @@ This project adheres to [Semantic Versioning](http://semver.org/spec/v2.0.0.html ## [Unreleased] -## [0.28.1] - 2021-03-14 +## [0.29.0] - 2022-04-21 +### Added +- MatchingMeshPartitioner "cubedsphere" +- Interpolator "cubedsphere-bilinear" +- Improvements to Interpolation::Cache +- Add support for rank 2 fields when a nonlinear action is added to the interpolator +- Create Array using ArraySpec only + +### Changed +- FieldSet::has(...) replaces FieldSet::has_field(...) +- Metadata return value to Interpolation::execute() +- Rename BilinearRemapping to UnstructuredBilinearLonLat + +### Fixed +- Compatibility with proj version >= 8 +- Compatibility with eckit version <= 1.18.5 +- Compatibility with GridTools backend and using 64bit idx_t +- Wrongly computed Jacobian::transpose() introduced in 0.28.0 +- Fix bug where using ectrans was not enabling adjoint of invtrans +- Avoid segfault when OpenMP tasking is broken, as it is with AppleClang and LLVM libomp + +## [0.28.1] - 2022-03-14 ### Fixed - Fix compilation for GNU 7.3 -## [0.28.0] - 2021-03-02 +## [0.28.0] - 2022-03-02 ### Added - Assignment of ArrayView from ArrayView - Grid "regional_variable_resolution" via a new VariableResolutionProjection @@ -345,6 +366,7 @@ This project adheres to [Semantic Versioning](http://semver.org/spec/v2.0.0.html ## 0.13.0 - 2018-02-16 [Unreleased]: https://github.com/ecmwf/atlas/compare/master...develop +[0.29.0]: https://github.com/ecmwf/atlas/compare/0.28.1...0.29.0 [0.28.1]: https://github.com/ecmwf/atlas/compare/0.28.0...0.28.1 [0.28.0]: https://github.com/ecmwf/atlas/compare/0.27.0...0.28.0 [0.27.0]: https://github.com/ecmwf/atlas/compare/0.26.0...0.27.0 diff --git a/VERSION b/VERSION index 48f7a71df..ae6dd4e20 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -0.28.1 +0.29.0 diff --git a/bamboo/GCC-env.sh b/bamboo/GCC-env.sh index 3397bd06c..a136aa713 100644 --- a/bamboo/GCC-env.sh +++ b/bamboo/GCC-env.sh @@ -16,4 +16,4 @@ module unload metview module unload netcdf4 module load cmake/3.16.5 -module load proj/6.1.1 +module load proj/8.2.1 diff --git a/cmake/FindPROJ.cmake b/cmake/FindPROJ.cmake new file mode 100644 index 000000000..6ac0d60cb --- /dev/null +++ b/cmake/FindPROJ.cmake @@ -0,0 +1,73 @@ +# (C) Copyright 2011- 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. + +# - Try to find the proj library +# Once done this will define +# +# PROJ_FOUND - system has proj +# PROJ_INCLUDE_DIRS - the proj include directory +# PROJ_LIBRARIES - Link these to use proj +# PROJ_VERSION - semantic version of proj +# + +### Set search paths from environment +if ( NOT PROJ_PATH AND PROJ_ROOT ) + set( PROJ_PATH ${PROJ_ROOT} ) +endif() +if ( NOT PROJ_PATH AND NOT "$ENV{PROJ_ROOT}" STREQUAL "" ) + set( PROJ_PATH "$ENV{PROJ_ROOT}" ) +endif() +if ( NOT PROJ_PATH AND NOT "$ENV{PROJ_PATH}" STREQUAL "" ) + set( PROJ_PATH "$ENV{PROJ_PATH}" ) +endif() +if ( NOT PROJ_PATH AND NOT "$ENV{PROJ_DIR}" STREQUAL "" ) + set( PROJ_PATH "$ENV{PROJ_DIR}" ) +endif() + + +### If search paths given, use it, otherwise, use pkg-config +if( PROJ_PATH ) + find_path(PROJ_INCLUDE_DIR NAMES proj.h PATHS ${PROJ_PATH} ${PROJ_PATH}/include PATH_SUFFIXES PROJ NO_DEFAULT_PATH ) + find_library(PROJ_LIBRARY NAMES proj PATHS ${PROJ_PATH} ${PROJ_PATH}/lib PATH_SUFFIXES PROJ NO_DEFAULT_PATH ) +else() + find_package(PkgConfig) + if(PKG_CONFIG_FOUND) + if(PROJ_FIND_VERSION) + pkg_check_modules(PKPROJ ${_pkgconfig_REQUIRED} QUIET PROJ>=${PROJ_FIND_VERSION}) + else() + pkg_check_modules(PKPROJ ${_pkgconfig_REQUIRED} QUIET PROJ) + endif() + + if( PKPROJ_FOUND ) + find_path(PROJ_INCLUDE_DIR proj.h HINTS ${PKPROJ_INCLUDEDIR} ${PKPROJ_INCLUDE_DIRS} PATH_SUFFIXES PROJ NO_DEFAULT_PATH ) + find_library(PROJ_LIBRARY proj HINTS ${PKPROJ_LIBDIR} ${PKPROJ_LIBRARY_DIRS} PATH_SUFFIXES PROJ NO_DEFAULT_PATH ) + endif() + endif() +endif() + +find_path(PROJ_INCLUDE_DIR NAMES proj.h PATHS PATH_SUFFIXES PROJ ) +find_library( PROJ_LIBRARY NAMES proj PATHS PATH_SUFFIXES PROJ ) + +### Detect version +set( PROJ_VERSION 0 ) +file(READ ${PROJ_INCLUDE_DIR}/proj.h proj_version) +string(REGEX REPLACE "^.*PROJ_VERSION_MAJOR +([0-9]+).*$" "\\1" PROJ_VERSION_MAJOR "${proj_version}") +string(REGEX REPLACE "^.*PROJ_VERSION_MINOR +([0-9]+).*$" "\\1" PROJ_VERSION_MINOR "${proj_version}") +string(REGEX REPLACE "^.*PROJ_VERSION_PATCH +([0-9]+).*$" "\\1" PROJ_VERSION_PATCH "${proj_version}") +string(CONCAT PROJ_VERSION ${PROJ_VERSION_MAJOR} "." ${PROJ_VERSION_MINOR} "." ${PROJ_VERSION_PATCH}) + +### Handle the QUIETLY and REQUIRED arguments and set PROJ_FOUND +include(FindPackageHandleStandardArgs) +find_package_handle_standard_args(PROJ + REQUIRED_VARS PROJ_LIBRARY PROJ_INCLUDE_DIR + VERSION_VAR PROJ_VERSION) + +set( PROJ_LIBRARIES ${PROJ_LIBRARY} ) +set( PROJ_INCLUDE_DIRS ${PROJ_INCLUDE_DIR} ) + +mark_as_advanced( PROJ_INCLUDE_DIR PROJ_LIBRARY ) diff --git a/cmake/features/OMP.cmake b/cmake/features/OMP.cmake index 78e30ea67..443239cbe 100644 --- a/cmake/features/OMP.cmake +++ b/cmake/features/OMP.cmake @@ -24,18 +24,18 @@ endif() if( HAVE_OMP_CXX ) if( NOT CMAKE_CXX_COMPILER_ID MATCHES Clang ) - set( ATLAS_OMP_TASK_UNTIED_SUPPORTED 1 ) + set( ATLAS_OMP_TASK_SUPPORTED 1 ) endif() - if( NOT DEFINED ATLAS_OMP_TASK_UNTIED_SUPPORTED ) - try_run( execute_result compile_result + if( NOT DEFINED ATLAS_OMP_TASK_SUPPORTED ) + try_run( execute_result compile_result ${CMAKE_CURRENT_BINARY_DIR} - ${PROJECT_SOURCE_DIR}/cmake/features/OMP/test_omp_untied.cc + ${PROJECT_SOURCE_DIR}/cmake/features/OMP/test_omp_task.cc LINK_LIBRARIES ${OMP_CXX} COMPILE_OUTPUT_VARIABLE compile_output RUN_OUTPUT_VARIABLE execute_output ) - ecbuild_debug("Compiling and running ${PROJECT_SOURCE_DIR}/cmake/features/OMP/test_omp_untied.cc") + ecbuild_debug("Compiling and running ${PROJECT_SOURCE_DIR}/cmake/features/OMP/test_omp_task.cc") ecbuild_debug_var( compile_result ) ecbuild_debug_var( compile_output ) ecbuild_debug_var( execute_result ) @@ -43,15 +43,53 @@ if( HAVE_OMP_CXX ) if( compile_result ) if( execute_result MATCHES 0 ) - set( ATLAS_OMP_TASK_UNTIED_SUPPORTED 1 ) + set( ATLAS_OMP_TASK_SUPPORTED 1 ) else() - ecbuild_info(" Compiler failed to run program with omp pragma with 'untied if' construct." - "Workaround will be enabled.") - set( ATLAS_OMP_TASK_UNTIED_SUPPORTED 0 ) + ecbuild_info(" Compiler failed to correctly run program with 'omp task' pragma." + "Sorting with OMP is disabled.") + set( ATLAS_OMP_TASK_SUPPORTED 0 ) endif() - else() - set( ATLAS_OMP_TASK_UNTIED_SUPPORTED 0 ) + else() + set( ATLAS_OMP_TASK_SUPPORTED 0 ) endif() endif() + + if( ATLAS_OMP_TASK_SUPPORTED ) + if( NOT CMAKE_CXX_COMPILER_ID MATCHES Clang ) + set( ATLAS_OMP_TASK_UNTIED_SUPPORTED 1 ) + endif() + + if( NOT DEFINED ATLAS_OMP_TASK_UNTIED_SUPPORTED ) + try_run( execute_result compile_result + ${CMAKE_CURRENT_BINARY_DIR} + ${PROJECT_SOURCE_DIR}/cmake/features/OMP/test_omp_untied.cc + LINK_LIBRARIES ${OMP_CXX} + COMPILE_OUTPUT_VARIABLE compile_output + RUN_OUTPUT_VARIABLE execute_output ) + + ecbuild_debug("Compiling and running ${PROJECT_SOURCE_DIR}/cmake/features/OMP/test_omp_untied.cc") + ecbuild_debug_var( compile_result ) + ecbuild_debug_var( compile_output ) + ecbuild_debug_var( execute_result ) + ecbuild_debug_var( execute_output ) + + if( compile_result ) + if( execute_result MATCHES 0 ) + set( ATLAS_OMP_TASK_UNTIED_SUPPORTED 1 ) + else() + ecbuild_info(" Compiler failed to run program with omp pragma with 'untied if' construct." + "Workaround will be enabled.") + set( ATLAS_OMP_TASK_UNTIED_SUPPORTED 0 ) + endif() + else() + set( ATLAS_OMP_TASK_UNTIED_SUPPORTED 0 ) + endif() + endif() + else() + set( ATLAS_OMP_TASK_UNTIED_SUPPORTED 0 ) + endif() +else() + set( ATLAS_OMP_TASK_SUPPORTED 0 ) + set( ATLAS_OMP_TASK_UNTIED_SUPPORTED 0 ) endif() diff --git a/cmake/features/OMP/test_omp_task.cc b/cmake/features/OMP/test_omp_task.cc new file mode 100644 index 000000000..744201b1a --- /dev/null +++ b/cmake/features/OMP/test_omp_task.cc @@ -0,0 +1,55 @@ +/* + * (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. + */ + + +// This test mirrors functionality from src/atlas/parallel/omp/sort.h +// It appears that on some installations / compilers (such as apple-clang with llvm omp) +// the tasking does not seem to work properly. +// It all compiles but leads to runtime errors. This executable is +// added to allow compiler introspection for this feature. + +#include +#include +#include + +#include + + +void recursive_task( int begin, int end ) { + int size = end - begin; + + if ( size >= 2 ) { // should be much larger in real case (e.g. 256) + int mid = begin + size / 2; + std::cout << begin << " - " << end << " [" <= 7) + find_package( PROJ CONFIG ) + if( PROJ_FOUND ) + ecbuild_debug("Found PROJ via CONFIG") + else() +# 2) Try to find PROJ4 the CMake way (proj 6) + find_package( PROJ4 CONFIG ) + if( PROJ4_FOUND ) + ecbuild_debug("Found PROJ4 via CONFIG") + set( PROJ_FOUND ${PROJ4_FOUND} ) + set( PROJ_LIBRARIES ${PROJ4_LIBRARIES} ) + set( PROJ_INCLUDE_DIRS ${PROJ4_INCLUDE_DIRS} ) + set( PROJ_VERSION ${PROJ4_VERSION} ) + endif() + endif() + if( NOT PROJ_FOUND ) +# 3) Try to find PROJ via FindPROJ.cmake provided within atlas (proj < 9) + find_package( PROJ MODULE ) + if( PROJ_FOUND ) + ecbuild_debug("Found PROJ via FindPROJ.cmake") + endif() + endif() +else() + ecbuild_debug("Skipping search for PROJ as ENABLE_PROJ=OFF (default)") +endif() + +if( PROJ_FOUND ) + ecbuild_info( "Found PROJ (${PROJ_VERSION}):") + ecbuild_info( " PROJ_LIBRARIES : ${PROJ_LIBRARIES}") + ecbuild_info( " PROJ_INCLUDE_DIRS : ${PROJ_INCLUDE_DIRS}") +endif() + ecbuild_add_option( FEATURE PROJ DESCRIPTION "PROJ-based projections" DEFAULT OFF - CONDITION PROJ4_FOUND ) + CONDITION PROJ_FOUND ) diff --git a/src/atlas/CMakeLists.txt b/src/atlas/CMakeLists.txt index 04fea0bb9..9ff575ddb 100644 --- a/src/atlas/CMakeLists.txt +++ b/src/atlas/CMakeLists.txt @@ -72,6 +72,7 @@ option/TransOptions.cc projection.h projection/Projection.cc projection/Projection.h +projection/Jacobian.h projection/detail/CubedSphereEquiAnglProjection.cc projection/detail/CubedSphereEquiAnglProjection.h projection/detail/CubedSphereEquiDistProjection.cc @@ -203,6 +204,8 @@ grid/detail/partitioner/MatchingMeshPartitioner.h grid/detail/partitioner/MatchingMeshPartitioner.cc grid/detail/partitioner/MatchingMeshPartitionerBruteForce.cc grid/detail/partitioner/MatchingMeshPartitionerBruteForce.h +grid/detail/partitioner/MatchingMeshPartitionerCubedSphere.cc +grid/detail/partitioner/MatchingMeshPartitionerCubedSphere.h grid/detail/partitioner/MatchingMeshPartitionerLonLatPolygon.cc grid/detail/partitioner/MatchingMeshPartitionerLonLatPolygon.h grid/detail/partitioner/MatchingMeshPartitionerSphericalPolygon.cc @@ -555,10 +558,10 @@ interpolation/method/PointSet.cc interpolation/method/PointSet.h interpolation/method/Ray.cc interpolation/method/Ray.h -interpolation/method/fe/FiniteElement.cc -interpolation/method/fe/FiniteElement.h -interpolation/method/bil/BilinearRemapping.cc -interpolation/method/bil/BilinearRemapping.h +interpolation/method/cubedsphere/CellFinder.cc +interpolation/method/cubedsphere/CellFinder.h +interpolation/method/cubedsphere/CubedSphereBilinear.cc +interpolation/method/cubedsphere/CubedSphereBilinear.h interpolation/method/knn/GridBox.cc interpolation/method/knn/GridBox.h interpolation/method/knn/GridBoxAverage.cc @@ -598,6 +601,10 @@ interpolation/method/structured/kernels/LinearVerticalKernel.h interpolation/method/structured/kernels/QuasiCubic3DKernel.cc interpolation/method/structured/kernels/QuasiCubic3DKernel.h interpolation/method/structured/kernels/QuasiCubicHorizontalKernel.h +interpolation/method/unstructured/FiniteElement.cc +interpolation/method/unstructured/FiniteElement.h +interpolation/method/unstructured/UnstructuredBilinearLonLat.cc +interpolation/method/unstructured/UnstructuredBilinearLonLat.h interpolation/nonlinear/Missing.cc interpolation/nonlinear/Missing.h interpolation/nonlinear/NonLinear.cc @@ -916,14 +923,15 @@ if( atlas_HAVE_TRANS ) endif() if( atlas_HAVE_PROJ ) - target_link_libraries( atlas PRIVATE ${PROJ4_LIBRARIES} ) - target_include_directories( atlas PRIVATE ${PROJ4_INCLUDE_DIRS} ) + target_link_libraries( atlas PRIVATE ${PROJ_LIBRARIES} ) + target_include_directories( atlas PRIVATE ${PROJ_INCLUDE_DIRS} ) endif() target_link_libraries( atlas PUBLIC eckit eckit_geometry eckit_linalg + eckit_maths eckit_mpi eckit_option ) diff --git a/src/atlas/array/Array.h b/src/atlas/array/Array.h index c8c263e64..80339b911 100644 --- a/src/atlas/array/Array.h +++ b/src/atlas/array/Array.h @@ -45,6 +45,8 @@ class Array : public util::Object { static Array* create(array::DataType, ArraySpec&&); + static Array* create(ArraySpec&&); + virtual size_t footprint() const = 0; template diff --git a/src/atlas/array/ArraySpec.cc b/src/atlas/array/ArraySpec.cc index c9c389490..afa0d4db3 100644 --- a/src/atlas/array/ArraySpec.cc +++ b/src/atlas/array/ArraySpec.cc @@ -29,11 +29,15 @@ size_t compute_aligned_size(size_t size, size_t alignment) { } } // namespace -ArraySpec::ArraySpec(): size_(), rank_(), allocated_size_(), contiguous_(true), default_layout_(true) {} +ArraySpec::ArraySpec(): + size_(), rank_(), datatype_(DataType::KIND_REAL64), allocated_size_(), contiguous_(true), default_layout_(true) {} ArraySpec::ArraySpec(const ArrayShape& shape): ArraySpec(shape, ArrayAlignment()) {} +ArraySpec::ArraySpec(DataType datatype, const ArrayShape& shape): ArraySpec(shape) { + datatype_ = datatype; +} -ArraySpec::ArraySpec(const ArrayShape& shape, const ArrayAlignment& alignment) { +ArraySpec::ArraySpec(const ArrayShape& shape, const ArrayAlignment& alignment): datatype_(DataType::KIND_REAL64) { ArrayShape aligned_shape = shape; aligned_shape.back() = compute_aligned_size(aligned_shape.back(), alignment); @@ -60,10 +64,16 @@ ArraySpec::ArraySpec(const ArrayShape& shape, const ArrayAlignment& alignment) { #endif }; +ArraySpec::ArraySpec(DataType datatype, const ArrayShape& shape, const ArrayAlignment& alignment): + ArraySpec(shape, alignment) { + datatype_ = datatype; +} + ArraySpec::ArraySpec(const ArrayShape& shape, const ArrayStrides& strides): ArraySpec(shape, strides, ArrayAlignment()) {} -ArraySpec::ArraySpec(const ArrayShape& shape, const ArrayStrides& strides, const ArrayAlignment& alignment) { +ArraySpec::ArraySpec(const ArrayShape& shape, const ArrayStrides& strides, const ArrayAlignment& alignment): + datatype_(DataType::KIND_REAL64) { ATLAS_ASSERT(shape.size() == strides.size(), "dimensions of shape and stride don't match"); rank_ = static_cast(shape.size()); @@ -86,11 +96,24 @@ ArraySpec::ArraySpec(const ArrayShape& shape, const ArrayStrides& strides, const #endif } +ArraySpec::ArraySpec(DataType datatype, const ArrayShape& shape, const ArrayStrides& strides, + const ArrayAlignment& alignment): + ArraySpec(shape, strides, alignment) { + datatype_ = datatype; +} + ArraySpec::ArraySpec(const ArrayShape& shape, const ArrayStrides& strides, const ArrayLayout& layout): ArraySpec(shape, strides, layout, ArrayAlignment()) {} +ArraySpec::ArraySpec(DataType datatype, const ArrayShape& shape, const ArrayStrides& strides, + const ArrayLayout& layout): + ArraySpec(shape, strides, layout) { + datatype_ = datatype; +} + ArraySpec::ArraySpec(const ArrayShape& shape, const ArrayStrides& strides, const ArrayLayout& layout, - const ArrayAlignment& alignment) { + const ArrayAlignment& alignment): + datatype_(DataType::KIND_REAL64) { ATLAS_ASSERT(shape.size() == strides.size(), "dimensions of shape and stride don't match"); rank_ = static_cast(shape.size()); @@ -116,6 +139,11 @@ ArraySpec::ArraySpec(const ArrayShape& shape, const ArrayStrides& strides, const #endif } +ArraySpec::ArraySpec(DataType datatype, const ArrayShape& shape, const ArrayStrides& strides, const ArrayLayout& layout, + const ArrayAlignment& alignment): + ArraySpec(shape, strides, layout, alignment) { + datatype_ = datatype; +} const std::vector& ArraySpec::shapef() const { return shapef_; } diff --git a/src/atlas/array/ArraySpec.h b/src/atlas/array/ArraySpec.h index 49d107d08..55c9330c5 100644 --- a/src/atlas/array/ArraySpec.h +++ b/src/atlas/array/ArraySpec.h @@ -17,6 +17,7 @@ #include "atlas/array/ArrayLayout.h" #include "atlas/array/ArrayShape.h" #include "atlas/array/ArrayStrides.h" +#include "atlas/array/DataType.h" #include "atlas/library/config.h" //------------------------------------------------------------------------------------------------------ @@ -29,6 +30,7 @@ class ArraySpec { size_t size_; idx_t rank_; size_t allocated_size_; + DataType datatype_; ArrayShape shape_; ArrayStrides strides_; ArrayLayout layout_; @@ -46,9 +48,16 @@ class ArraySpec { ArraySpec(const ArrayShape&, const ArrayAlignment&); ArraySpec(const ArrayShape&, const ArrayStrides&, const ArrayAlignment&); ArraySpec(const ArrayShape&, const ArrayStrides&, const ArrayLayout&, const ArrayAlignment&); + ArraySpec(DataType, const ArrayShape&); + ArraySpec(DataType, const ArrayShape&, const ArrayStrides&); + ArraySpec(DataType, const ArrayShape&, const ArrayStrides&, const ArrayLayout&); + ArraySpec(DataType, const ArrayShape&, const ArrayAlignment&); + ArraySpec(DataType, const ArrayShape&, const ArrayStrides&, const ArrayAlignment&); + ArraySpec(DataType, const ArrayShape&, const ArrayStrides&, const ArrayLayout&, const ArrayAlignment&); size_t allocatedSize() const { return allocated_size_; } size_t size() const { return size_; } idx_t rank() const { return rank_; } + DataType datatype() const { return datatype_; } const ArrayShape& shape() const { return shape_; } const ArrayAlignment& alignment() const { return alignment_; } const ArrayStrides& strides() const { return strides_; } diff --git a/src/atlas/array/DataType.h b/src/atlas/array/DataType.h index a18c34035..a8485ef01 100644 --- a/src/atlas/array/DataType.h +++ b/src/atlas/array/DataType.h @@ -86,6 +86,7 @@ class DataType { DataType(const std::string&); DataType(long); DataType(const DataType&); + DataType& operator=(const DataType&); std::string str() const { return kind_to_str(kind_); } kind_t kind() const { return kind_; } size_t size() const { return (kind_ == KIND_UINT64) ? 8 : std::abs(kind_); } @@ -361,6 +362,11 @@ inline bool DataType::kind_valid(kind_t kind) { inline DataType::DataType(const DataType& other): kind_(other.kind_) {} +inline DataType& DataType::operator=(const DataType& other) { + kind_ = other.kind_; + return *this; +} + inline DataType::DataType(const std::string& datatype): kind_(str_to_kind(datatype)) {} inline DataType::DataType(long kind): kind_(kind) {} diff --git a/src/atlas/array/Range.h b/src/atlas/array/Range.h index 5d8b918ec..d2c588a62 100644 --- a/src/atlas/array/Range.h +++ b/src/atlas/array/Range.h @@ -96,8 +96,8 @@ class Range : public helpers::RangeBase { public: static From from(int start) { return From(start); } static To to(int end) { return To(end); } - static All all() { return All(); } - static Dummy dummy() { return Dummy(); } + static constexpr All all() { return All(); } + static constexpr Dummy dummy() { return Dummy(); } public: template diff --git a/src/atlas/array/gridtools/GridToolsArray.cc b/src/atlas/array/gridtools/GridToolsArray.cc index 8fbc4d5f8..27a5ee462 100644 --- a/src/atlas/array/gridtools/GridToolsArray.cc +++ b/src/atlas/array/gridtools/GridToolsArray.cc @@ -410,6 +410,10 @@ Array* Array::create(DataType datatype, ArraySpec&& spec) { } } +Array* Array::create(ArraySpec&& spec) { + return create(spec.datatype(), std::move(spec)); +} + //------------------------------------------------------------------------------ Array::~Array() = default; diff --git a/src/atlas/array/gridtools/GridToolsArrayView.h b/src/atlas/array/gridtools/GridToolsArrayView.h index be0b25e2e..899646d54 100644 --- a/src/atlas/array/gridtools/GridToolsArrayView.h +++ b/src/atlas/array/gridtools/GridToolsArrayView.h @@ -92,16 +92,18 @@ class ArrayView { value_type* data() { return gt_data_view_.data(); } value_type const* data() const { return gt_data_view_.data(); } - template ::type> - ATLAS_HOST_DEVICE value_type& operator()(Coords... c) { - assert(sizeof...(Coords) == Rank); - return gt_data_view_(c...); + template ::type> + ATLAS_HOST_DEVICE value_type& operator()(Ints... c) { + assert(sizeof...(Ints) == Rank); + using common_type = typename std::common_type::type; + return gt_data_view_(static_cast(c)...); } - template ::type> - ATLAS_HOST_DEVICE value_type const& operator()(Coords... c) const { - assert(sizeof...(Coords) == Rank); - return gt_data_view_(c...); + template ::type> + ATLAS_HOST_DEVICE value_type const& operator()(Ints... c) const { + assert(sizeof...(Ints) == Rank); + using common_type = typename std::common_type::type; + return gt_data_view_(static_cast(c)...); } template diff --git a/src/atlas/array/native/NativeArray.cc b/src/atlas/array/native/NativeArray.cc index 746012ec8..e75defeae 100644 --- a/src/atlas/array/native/NativeArray.cc +++ b/src/atlas/array/native/NativeArray.cc @@ -102,6 +102,11 @@ Array* Array::create(DataType datatype, ArraySpec&& spec) { } } +Array* Array::create(ArraySpec&& spec) { + return create(spec.datatype(), std::move(spec)); +} + + template ArrayT::ArrayT(ArrayDataStore* ds, const ArraySpec& spec) { data_store_ = std::unique_ptr(ds); diff --git a/src/atlas/field/FieldSet.cc b/src/atlas/field/FieldSet.cc index f12539d33..9b7a307b5 100644 --- a/src/atlas/field/FieldSet.cc +++ b/src/atlas/field/FieldSet.cc @@ -40,12 +40,13 @@ Field FieldSetImpl::add(const Field& field) { return field; } -bool FieldSetImpl::has_field(const std::string& name) const { +bool FieldSetImpl::has(const std::string& name) const { return index_.count(name); } + Field& FieldSetImpl::field(const std::string& name) const { - if (!has_field(name)) { + if (!has(name)) { const std::string msg("FieldSet" + (name_.length() ? " \"" + name_ + "\"" : "") + ": cannot find field \"" + name + "\""); throw_Exception(msg, Here()); @@ -104,7 +105,7 @@ void atlas__FieldSet__add_field(FieldSetImpl* This, FieldImpl* field) { int atlas__FieldSet__has_field(const FieldSetImpl* This, char* name) { ATLAS_ASSERT(This != nullptr, "Reason: Use of uninitialised atlas_FieldSet"); - return This->has_field(std::string(name)); + return This->has(std::string(name)); } idx_t atlas__FieldSet__size(const FieldSetImpl* This) { diff --git a/src/atlas/field/FieldSet.h b/src/atlas/field/FieldSet.h index 4bde79fed..b93eb43f2 100644 --- a/src/atlas/field/FieldSet.h +++ b/src/atlas/field/FieldSet.h @@ -22,6 +22,8 @@ #include #include +#include "eckit/deprecated.h" + #include "atlas/field/Field.h" #include "atlas/library/config.h" #include "atlas/runtime/Exception.h" @@ -99,7 +101,7 @@ class FieldSetImpl : public util::Object { Field add(const Field&); - bool has_field(const std::string& name) const; + bool has(const std::string& name) const; Field& field(const std::string& name) const; @@ -195,7 +197,7 @@ class FieldSet : DOXYGEN_HIDE(public util::ObjectHandle) { Field add(const Field& field) { return get()->add(field); } - bool has_field(const std::string& name) const { return get()->has_field(name); } + bool has(const std::string& name) const { return get()->has(name); } Field& field(const std::string& name) const { return get()->field(name); } @@ -208,6 +210,9 @@ class FieldSet : DOXYGEN_HIDE(public util::ObjectHandle) { void haloExchange(bool on_device = false) const { get()->haloExchange(on_device); } void set_dirty(bool = true) const; + + // Deprecated API + DEPRECATED("use 'has' instead") bool has_field(const std::string& name) const { return get()->has(name); } }; } // namespace atlas diff --git a/src/atlas/functionspace/CubedSphereColumns.cc b/src/atlas/functionspace/CubedSphereColumns.cc index c2a880b91..7708631fb 100644 --- a/src/atlas/functionspace/CubedSphereColumns.cc +++ b/src/atlas/functionspace/CubedSphereColumns.cc @@ -150,6 +150,11 @@ idx_t CubedSphereColumns::index(idx_t t, idx_t i, idx_t j) co return cubedSphereColumnsHandle_.get()->index(t, i, j); } +template +bool CubedSphereColumns::is_valid_index(idx_t t, idx_t i, idx_t j) const { + return cubedSphereColumnsHandle_.get()->is_valid_index(t, i, j); +} + template Field CubedSphereColumns::tij() const { return cubedSphereColumnsHandle_.get()->tij(); diff --git a/src/atlas/functionspace/CubedSphereColumns.h b/src/atlas/functionspace/CubedSphereColumns.h index 916a6c88d..25bc68104 100644 --- a/src/atlas/functionspace/CubedSphereColumns.h +++ b/src/atlas/functionspace/CubedSphereColumns.h @@ -54,6 +54,9 @@ class CubedSphereColumns : public BaseFunctionSpace { /// Return array_view index for (t, i, j). idx_t index(idx_t t, idx_t i, idx_t j) const; + /// Return true if (t, i, j) is a valid index. + bool is_valid_index(idx_t t, idx_t i, idx_t j) const; + /// Return tij field. Field tij() const; diff --git a/src/atlas/functionspace/NodeColumns.cc b/src/atlas/functionspace/NodeColumns.cc index 8bc2742b5..6ec46973f 100644 --- a/src/atlas/functionspace/NodeColumns.cc +++ b/src/atlas/functionspace/NodeColumns.cc @@ -358,6 +358,27 @@ void dispatch_haloExchange(Field& field, const parallel::HaloExchange& halo_exch } field.set_dirty(false); } + +template +void dispatch_adjointHaloExchange(Field& field, const parallel::HaloExchange& halo_exchange, bool on_device) { + if (field.datatype() == array::DataType::kind()) { + halo_exchange.template execute_adjoint(field.array(), on_device); + } + else if (field.datatype() == array::DataType::kind()) { + halo_exchange.template execute_adjoint(field.array(), on_device); + } + else if (field.datatype() == array::DataType::kind()) { + halo_exchange.template execute_adjoint(field.array(), on_device); + } + else if (field.datatype() == array::DataType::kind()) { + halo_exchange.template execute_adjoint(field.array(), on_device); + } + else { + throw_Exception("datatype not supported", Here()); + } + field.set_dirty(false); +} + } // namespace void NodeColumns::haloExchange(const FieldSet& fieldset, bool on_device) const { @@ -382,11 +403,40 @@ void NodeColumns::haloExchange(const FieldSet& fieldset, bool on_device) const { } } +void NodeColumns::adjointHaloExchange(const FieldSet& fieldset, bool on_device) const { + for (idx_t f = 0; f < fieldset.size(); ++f) { + Field& field = const_cast(fieldset)[f]; + switch (field.rank()) { + case 1: + dispatch_adjointHaloExchange<1>(field, halo_exchange(), on_device); + break; + case 2: + dispatch_adjointHaloExchange<2>(field, halo_exchange(), on_device); + break; + case 3: + dispatch_adjointHaloExchange<3>(field, halo_exchange(), on_device); + break; + case 4: + dispatch_adjointHaloExchange<4>(field, halo_exchange(), on_device); + break; + default: + throw_Exception("Rank not supported", Here()); + } + } +} + void NodeColumns::haloExchange(const Field& field, bool on_device) const { FieldSet fieldset; fieldset.add(field); haloExchange(fieldset, on_device); } + +void NodeColumns::adjointHaloExchange(const Field& field, bool) const { + FieldSet fieldset; + fieldset.add(field); + adjointHaloExchange(fieldset); +} + const parallel::HaloExchange& NodeColumns::halo_exchange() const { if (halo_exchange_) { return *halo_exchange_; diff --git a/src/atlas/functionspace/NodeColumns.h b/src/atlas/functionspace/NodeColumns.h index 01872b615..2f8217fe7 100644 --- a/src/atlas/functionspace/NodeColumns.h +++ b/src/atlas/functionspace/NodeColumns.h @@ -85,6 +85,9 @@ class NodeColumns : public functionspace::FunctionSpaceImpl { void haloExchange(const Field&, bool on_device = false) const override; const parallel::HaloExchange& halo_exchange() const; + void adjointHaloExchange(const FieldSet&, bool on_device = false) const override; + void adjointHaloExchange(const Field&, bool on_device = false) const override; + void gather(const FieldSet&, FieldSet&) const; void gather(const Field&, Field&) const; const parallel::GatherScatter& gather() const; diff --git a/src/atlas/functionspace/detail/CubedSphereStructure.cc b/src/atlas/functionspace/detail/CubedSphereStructure.cc index 612825851..dc369d1af 100644 --- a/src/atlas/functionspace/detail/CubedSphereStructure.cc +++ b/src/atlas/functionspace/detail/CubedSphereStructure.cc @@ -120,6 +120,25 @@ idx_t CubedSphereStructure::index(idx_t t, idx_t i, idx_t j) const { return tijToIdx_[static_cast(t)][vecIndex(t, i, j)]; } +bool CubedSphereStructure::is_valid_index(idx_t t, idx_t i, idx_t j) const { + // Check if t is in range. + if (t < 0 || t > 5) { + return false; + } + + // Check if i and j are in range in index method. + if (i < i_begin(t) || i >= i_end(t) || j < j_begin(t) || j >= j_end(t)) { + return false; + } + + // Check if (t, i, j) is a valid index. + if (index(t, i, j) == invalid_index()) { + return false; + } + + return true; +} + Field CubedSphereStructure::tij() const { return tij_; } diff --git a/src/atlas/functionspace/detail/CubedSphereStructure.h b/src/atlas/functionspace/detail/CubedSphereStructure.h index aaa8354fc..4f21ad926 100644 --- a/src/atlas/functionspace/detail/CubedSphereStructure.h +++ b/src/atlas/functionspace/detail/CubedSphereStructure.h @@ -53,6 +53,10 @@ class CubedSphereStructure : public util::Object { /// Return array_view index for (t, i, j). idx_t index(idx_t t, idx_t i, idx_t j) const; + /// Return true if (t, i, j) is a valid index. + bool is_valid_index(idx_t t, idx_t i, idx_t j) const; + + /// Return ijt field. Field tij() const; diff --git a/src/atlas/grid/CubedSphereGrid.h b/src/atlas/grid/CubedSphereGrid.h index 659abf64b..7cd758bea 100644 --- a/src/atlas/grid/CubedSphereGrid.h +++ b/src/atlas/grid/CubedSphereGrid.h @@ -263,7 +263,9 @@ class IterateTIJ { class CubedSphereGrid : public Grid { public: - using grid_t = grid::detail::grid::CubedSphere; + using grid_t = grid::detail::grid::CubedSphere; + using CubedSphereProjection = projection::detail::CubedSphereProjectionBase; + using CubedSphereTiles = grid::CubedSphereTiles; public: CubedSphereGrid(); @@ -295,8 +297,13 @@ class CubedSphereGrid : public Grid { // Return the size of the cubed sphere grid, where N is the number of grid boxes along the edge of a tile inline int N() const { return grid_->N(); } - // Return the number of tiles - inline atlas::grid::CubedSphereTiles tiles() const { return grid_->tiles(); } + /// @brief return tiles object. + inline CubedSphereTiles tiles() const { return grid_->tiles(); } + + /// @brief return cubed sphere projection object. + inline const CubedSphereProjection& cubedSphereProjection() const { + return dynamic_cast(*projection().get()); + }; temporary::IterateTIJ tij() const { return temporary::IterateTIJ(*grid_); } diff --git a/src/atlas/grid/Tiles.cc b/src/atlas/grid/Tiles.cc index bc8665a14..d0e24af8c 100644 --- a/src/atlas/grid/Tiles.cc +++ b/src/atlas/grid/Tiles.cc @@ -12,6 +12,7 @@ #include "atlas/grid/detail/tiles/FV3Tiles.h" #include "atlas/grid/detail/tiles/LFRicTiles.h" #include "atlas/grid/detail/tiles/Tiles.h" +#include "atlas/projection/Jacobian.h" namespace atlas { @@ -47,10 +48,18 @@ idx_t CubedSphereTiles::indexFromXY(const double xy[]) const { return get()->indexFromXY(xy); } +idx_t CubedSphereTiles::indexFromXY(const PointXY& xy) const { + return get()->indexFromXY(xy); +} + idx_t CubedSphereTiles::indexFromLonLat(const double lonlat[]) const { return get()->indexFromLonLat(lonlat); } +idx_t CubedSphereTiles::indexFromLonLat(const PointLonLat& lonlat) const { + return get()->indexFromLonLat(lonlat.data()); +} + idx_t CubedSphereTiles::size() const { return get()->size(); } @@ -72,5 +81,13 @@ std::ostream& operator<<(std::ostream& os, const CubedSphereTiles& t) { return os; } +const PointXY& CubedSphereTiles::tileCentre(size_t t) const { + return get()->tileCentre(t); +} + +const projection::Jacobian& CubedSphereTiles::tileJacobian(size_t t) const { + return get()->tileJacobian(t); +} + } // namespace grid } // namespace atlas diff --git a/src/atlas/grid/Tiles.h b/src/atlas/grid/Tiles.h index 723e932d1..5b6fe528b 100644 --- a/src/atlas/grid/Tiles.h +++ b/src/atlas/grid/Tiles.h @@ -35,6 +35,10 @@ namespace util { class Config; } // namespace util +namespace projection { +class Jacobian; +} + namespace grid { #ifndef DOXYGEN_SHOULD_SKIP_THIS @@ -70,9 +74,11 @@ class CubedSphereTiles : DOXYGEN_HIDE(public util::ObjectHandle::epsilon() * N; + const auto edgeEpsilon = epsilon; + const size_t listSize = 4; + + // Loop over grid and set partioning[]. + auto lonlatIt = grid.lonlat().begin(); + for (gidx_t i = 0; i < grid.size(); ++i) { + // This is probably more expensive than it needs to be, as it performs + // a dry run of the cubedsphere interpolation method. + const auto& lonlat = *lonlatIt; + partitioning[i] = finder.getCell(lonlat, listSize, edgeEpsilon, epsilon).isect ? mpi::rank() : -1; + ++lonlatIt; + } + + // AllReduce to get full partitioning array. + mpi::comm().allReduceInPlace(partitioning, grid.size(), eckit::mpi::Operation::MAX); + const auto misses = std::count_if(partitioning, partitioning + grid.size(), [](int elem) { return elem < 0; }); + if (misses > 0) { + throw_Exception( + "Could not find partition for target node (source " + "mesh does not contain all target grid points)\n" + + std::to_string(misses) + " misses.", + Here()); + } +} + +} // namespace partitioner +} // namespace detail +} // namespace grid +} // namespace atlas diff --git a/src/atlas/grid/detail/partitioner/MatchingMeshPartitionerCubedSphere.h b/src/atlas/grid/detail/partitioner/MatchingMeshPartitionerCubedSphere.h new file mode 100644 index 000000000..d6cf3be36 --- /dev/null +++ b/src/atlas/grid/detail/partitioner/MatchingMeshPartitionerCubedSphere.h @@ -0,0 +1,38 @@ +/* + * (C) Crown Copyright 2022 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/grid/detail/partitioner/MatchingMeshPartitioner.h" + +namespace atlas { +namespace grid { +namespace detail { +namespace partitioner { + + +class MatchingMeshPartitionerCubedSphere : public MatchingMeshPartitioner { +public: + static std::string static_type() { return "cubedsphere"; } + +public: + MatchingMeshPartitionerCubedSphere(): MatchingMeshPartitioner() {} + MatchingMeshPartitionerCubedSphere(const idx_t nb_partitions): MatchingMeshPartitioner(nb_partitions) {} + MatchingMeshPartitionerCubedSphere(const idx_t nb_partitions, const eckit::Parametrisation&): + MatchingMeshPartitioner(nb_partitions) {} + MatchingMeshPartitionerCubedSphere(const Mesh& mesh): MatchingMeshPartitioner(mesh) {} + + using MatchingMeshPartitioner::partition; + virtual void partition(const Grid& grid, int partitioning[]) const; + + virtual std::string type() const { return static_type(); } +}; + +} // namespace partitioner +} // namespace detail +} // namespace grid +} // namespace atlas diff --git a/src/atlas/grid/detail/partitioner/Partitioner.cc b/src/atlas/grid/detail/partitioner/Partitioner.cc index d8b040c13..9211c63d5 100644 --- a/src/atlas/grid/detail/partitioner/Partitioner.cc +++ b/src/atlas/grid/detail/partitioner/Partitioner.cc @@ -29,6 +29,7 @@ #include "atlas/grid/detail/partitioner/MatchingFunctionSpacePartitionerLonLatPolygon.h" #include "atlas/grid/detail/partitioner/MatchingMeshPartitioner.h" #include "atlas/grid/detail/partitioner/MatchingMeshPartitionerBruteForce.h" +#include "atlas/grid/detail/partitioner/MatchingMeshPartitionerCubedSphere.h" #include "atlas/grid/detail/partitioner/MatchingMeshPartitionerLonLatPolygon.h" #include "atlas/grid/detail/partitioner/MatchingMeshPartitionerSphericalPolygon.h" #include "atlas/grid/detail/partitioner/RegularBandsPartitioner.h" @@ -202,6 +203,9 @@ Partitioner* MatchingPartitionerFactory::build(const std::string& type, const Me else if (type == MatchingMeshPartitionerBruteForce::static_type()) { return new MatchingMeshPartitionerBruteForce(partitioned); } + else if (type == MatchingMeshPartitionerCubedSphere::static_type()) { + return new MatchingMeshPartitionerCubedSphere(partitioned); + } else { ATLAS_NOTIMPLEMENTED; } diff --git a/src/atlas/grid/detail/spacing/Spacing.cc b/src/atlas/grid/detail/spacing/Spacing.cc index 6086ecd85..d3dfccb8a 100644 --- a/src/atlas/grid/detail/spacing/Spacing.cc +++ b/src/atlas/grid/detail/spacing/Spacing.cc @@ -8,9 +8,10 @@ * nor does it submit to any jurisdiction. */ +#include "Spacing.h" + #include "eckit/config/Parametrisation.h" -#include "atlas/grid/detail/spacing/Spacing.h" #include "atlas/grid/detail/spacing/SpacingFactory.h" #include "atlas/runtime/Exception.h" diff --git a/src/atlas/grid/detail/spacing/Spacing.h b/src/atlas/grid/detail/spacing/Spacing.h index bc973e10b..4d5685ca3 100644 --- a/src/atlas/grid/detail/spacing/Spacing.h +++ b/src/atlas/grid/detail/spacing/Spacing.h @@ -11,6 +11,7 @@ #pragma once #include +#include #include #include "atlas/util/Object.h" diff --git a/src/atlas/grid/detail/tiles/FV3Tiles.cc b/src/atlas/grid/detail/tiles/FV3Tiles.cc index ce9ef91f5..da1cf5628 100644 --- a/src/atlas/grid/detail/tiles/FV3Tiles.cc +++ b/src/atlas/grid/detail/tiles/FV3Tiles.cc @@ -453,6 +453,16 @@ void FV3CubedSphereTiles::print(std::ostream& os) const { } +const PointXY& FV3CubedSphereTiles::tileCentre(size_t t) const { + throw_NotImplemented("tileCentre not implemented for FV3Tiles", Here()); +} + +const FV3CubedSphereTiles::Jacobian& FV3CubedSphereTiles::tileJacobian(size_t t) const { + throw_NotImplemented("tileJacobian not implemented for FV3Tiles", Here()); + ; +} + + namespace { static CubedSphereTilesBuilder register_builder(FV3CubedSphereTiles::static_type()); } diff --git a/src/atlas/grid/detail/tiles/FV3Tiles.h b/src/atlas/grid/detail/tiles/FV3Tiles.h index 9a248f54e..8232ac82b 100644 --- a/src/atlas/grid/detail/tiles/FV3Tiles.h +++ b/src/atlas/grid/detail/tiles/FV3Tiles.h @@ -46,6 +46,11 @@ class FV3CubedSphereTiles : public CubedSphereTiles { virtual void print(std::ostream&) const override; + virtual const PointXY& tileCentre(size_t t) const override; + + virtual const Jacobian& tileJacobian(size_t t) const override; + + private: }; diff --git a/src/atlas/grid/detail/tiles/LFRicTiles.cc b/src/atlas/grid/detail/tiles/LFRicTiles.cc index 8705eb244..10ace5374 100644 --- a/src/atlas/grid/detail/tiles/LFRicTiles.cc +++ b/src/atlas/grid/detail/tiles/LFRicTiles.cc @@ -26,68 +26,49 @@ namespace detail { namespace { -static constexpr bool debug = false; // constexpr so compiler can optimize `if ( debug ) { ... }` out +static constexpr bool debug = false; // constexpr so compiler can optimize `if ( debug ) { ... }` out +static constexpr double epsilon = 1.e-12; -using atlas::projection::detail::ProjectionUtilities; +using projection::detail::ProjectionUtilities; static bool is_tiny(const double& x) { - constexpr double epsilon = 1.e-12; return (std::abs(x) < epsilon); } static bool is_same(const double& x, const double& y, const double& tol = 1.0) { - constexpr double epsilon = 1.e-12; return (std::abs(x - y) < epsilon * tol); } static bool is_less(const double& lhs, const double& rhs) { - constexpr double epsilon = 1.e-12; return lhs < rhs - epsilon; } static bool is_geq(const double& lhs, const double& rhs) { - constexpr double epsilon = 1.e-12; return lhs >= rhs - epsilon; } -void sphericalToCartesian(const double lonlat[], double xyz[]) { +static void sphericalToCartesian(const double lonlat[], double xyz[]) { auto crd_sys = ProjectionUtilities::CoordinateSystem::LEFT_HAND; constexpr double radius = 1.; ProjectionUtilities::sphericalToCartesian(lonlat, xyz, crd_sys, radius); } -atlas::PointXY rotatePlus90AboutPt(const atlas::PointXY& xy, const atlas::PointXY& origin) { - return atlas::PointXY{-xy.y() + origin.x() + origin.y(), xy.x() - origin.x() + origin.y()}; +static PointXY rotatePlus90AboutPt(const PointXY& xy, const PointXY& origin) { + return PointXY{-xy.y() + origin.x() + origin.y(), xy.x() - origin.x() + origin.y()}; } -atlas::PointXY rotateMinus90AboutPt(const atlas::PointXY& xy, const atlas::PointXY& origin) { - return atlas::PointXY{xy.y() + origin.x() - origin.y(), -xy.x() + origin.x() + origin.y()}; +static PointXY rotateMinus90AboutPt(const PointXY& xy, const PointXY& origin) { + return PointXY{xy.y() + origin.x() - origin.y(), -xy.x() + origin.x() + origin.y()}; } -atlas::PointXY rotatePlus180AboutPt(const atlas::PointXY& xy, const atlas::PointXY& origin) { - return atlas::PointXY{2.0 * origin.x() - xy.x(), 2.0 * origin.y() - xy.y()}; +static PointXY rotatePlus180AboutPt(const PointXY& xy, const PointXY& origin) { + return PointXY{2.0 * origin.x() - xy.x(), 2.0 * origin.y() - xy.y()}; } } // anonymous namespace // constructor -LFRicCubedSphereTiles::LFRicCubedSphereTiles(const eckit::Parametrisation&) { - botLeftTile_ = - std::array{atlas::PointXY{0., -45.}, atlas::PointXY{90, -45}, atlas::PointXY{180., -45.}, - atlas::PointXY{270, -45}, atlas::PointXY{0., 45.}, atlas::PointXY{0, -135.}}; - - botRightTile_ = - std::array{atlas::PointXY{90., -45.}, atlas::PointXY{180., -45}, atlas::PointXY{270., -45.}, - atlas::PointXY{360., -45}, atlas::PointXY{90., 45.}, atlas::PointXY{90., -135.}}; - - topLeftTile_ = - std::array{atlas::PointXY{0., 45.}, atlas::PointXY{90, 45}, atlas::PointXY{180., 45.}, - atlas::PointXY{270, 45}, atlas::PointXY{0., 135.}, atlas::PointXY{0, -45.}}; - - topRightTile_ = - std::array{atlas::PointXY{90., 45.}, atlas::PointXY{180., 45}, atlas::PointXY{270., 45.}, - atlas::PointXY{360., 45}, atlas::PointXY{90., 135.}, atlas::PointXY{90., -45.}}; -} +LFRicCubedSphereTiles::LFRicCubedSphereTiles(const eckit::Parametrisation&) {} std::array, 2> LFRicCubedSphereTiles::xy2abOffsets() const { return {{{0., 1., 2., 3., 0., 0.}, {1., 1., 1., 1., 2., 0.}}}; @@ -97,78 +78,78 @@ std::array, 2> LFRicCubedSphereTiles::ab2xyOffsets() const return {{{0., 90., 180., 270., 0., 0.}, {-45., -45., -45., -45., 45., -135.}}}; } -static void tile0Rotate(double xyz[]) { +void tile0Rotate(double xyz[]) { // Face 0, no rotation. } -static void tile1Rotate(double xyz[]) { +void tile1Rotate(double xyz[]) { double xyz_in[3]; std::copy(xyz, xyz + 3, xyz_in); xyz[XX] = -xyz_in[YY]; xyz[YY] = xyz_in[XX]; } -static void tile2Rotate(double xyz[]) { +void tile2Rotate(double xyz[]) { double xyz_in[3]; std::copy(xyz, xyz + 3, xyz_in); xyz[XX] = -xyz_in[XX]; xyz[YY] = -xyz_in[YY]; } -static void tile3Rotate(double xyz[]) { +void tile3Rotate(double xyz[]) { double xyz_in[3]; std::copy(xyz, xyz + 3, xyz_in); xyz[XX] = xyz_in[YY]; xyz[YY] = -xyz_in[XX]; } -static void tile4Rotate(double xyz[]) { +void tile4Rotate(double xyz[]) { double xyz_in[3]; std::copy(xyz, xyz + 3, xyz_in); xyz[XX] = xyz_in[ZZ]; xyz[ZZ] = -xyz_in[XX]; } -static void tile5Rotate(double xyz[]) { +void tile5Rotate(double xyz[]) { double xyz_in[3]; std::copy(xyz, xyz + 3, xyz_in); xyz[XX] = -xyz_in[ZZ]; xyz[ZZ] = xyz_in[XX]; } -static void tile0RotateInverse(double xyz[]) { +void tile0RotateInverse(double xyz[]) { // Face 0, no rotation. } -static void tile1RotateInverse(double xyz[]) { +void tile1RotateInverse(double xyz[]) { double xyz_in[3]; std::copy(xyz, xyz + 3, xyz_in); xyz[XX] = xyz_in[YY]; xyz[YY] = -xyz_in[XX]; } -static void tile2RotateInverse(double xyz[]) { +void tile2RotateInverse(double xyz[]) { double xyz_in[3]; std::copy(xyz, xyz + 3, xyz_in); xyz[XX] = -xyz_in[XX]; xyz[YY] = -xyz_in[YY]; } -static void tile3RotateInverse(double xyz[]) { +void tile3RotateInverse(double xyz[]) { double xyz_in[3]; std::copy(xyz, xyz + 3, xyz_in); xyz[XX] = -xyz_in[YY]; xyz[YY] = xyz_in[XX]; } -static void tile4RotateInverse(double xyz[]) { +void tile4RotateInverse(double xyz[]) { double xyz_in[3]; std::copy(xyz, xyz + 3, xyz_in); xyz[XX] = -xyz_in[ZZ]; xyz[ZZ] = xyz_in[XX]; } -static void tile5RotateInverse(double xyz[]) { +void tile5RotateInverse(double xyz[]) { double xyz_in[3]; std::copy(xyz, xyz + 3, xyz_in); xyz[XX] = xyz_in[ZZ]; @@ -380,8 +361,7 @@ void LFRicCubedSphereTiles::enforceXYdomain(double xy[]) const { // input is the xy value as PointXY that is a continous "cross " extension in terms of xy space from the tile in question // the output is an xy value that lives on the standard "|---" shape -atlas::PointXY LFRicCubedSphereTiles::tileCubePeriodicity(const atlas::PointXY& xyExtended, - const atlas::idx_t t) const { +PointXY LFRicCubedSphereTiles::tileCubePeriodicity(const PointXY& xyExtended, const idx_t t) const { // xy space is a function of tiles--Tile 0) // xy space for Tile 1 // // // y ^ // y ^ @@ -484,24 +464,23 @@ atlas::PointXY LFRicCubedSphereTiles::tileCubePeriodicity(const atlas::PointXY& << " is not in the cross extension of tile " << t); } - atlas::PointXY withinRange = xyExtended; + PointXY withinRange = xyExtended; enforceWrapAround(t, withinRange); - atlas::PointXY finalXY = withinRange; - atlas::PointXY tempXY = withinRange; + PointXY finalXY = withinRange; + PointXY tempXY = withinRange; switch (t) { case 0: - finalXY = (withinRange.y() > 135.0) ? rotatePlus180AboutPt(withinRange, atlas::PointXY{135.0, 90.0}) - : withinRange; + finalXY = (withinRange.y() > 135.0) ? rotatePlus180AboutPt(withinRange, PointXY{135.0, 90.0}) : withinRange; break; case 1: if ((withinRange.x() >= 90.0) && (withinRange.x() <= 180.0)) { if (withinRange.y() >= 45.0) { - tempXY = rotatePlus90AboutPt(withinRange, atlas::PointXY{90.0, 45.0}); + tempXY = rotatePlus90AboutPt(withinRange, PointXY{90.0, 45.0}); if (withinRange.y() > 135.0) { - finalXY = rotatePlus90AboutPt(tempXY, atlas::PointXY{0.0, 45.0}); + finalXY = rotatePlus90AboutPt(tempXY, PointXY{0.0, 45.0}); finalXY.x() += 360.0; } else { @@ -509,22 +488,22 @@ atlas::PointXY LFRicCubedSphereTiles::tileCubePeriodicity(const atlas::PointXY& } } else if (withinRange.y() < -45.0) { - finalXY = rotateMinus90AboutPt(withinRange, atlas::PointXY{90.0, -45.0}); + finalXY = rotateMinus90AboutPt(withinRange, PointXY{90.0, -45.0}); } } break; case 2: if ((withinRange.x() >= 180.0) && (withinRange.x() <= 270.0)) { if (withinRange.y() > 135.0) { - finalXY = rotatePlus180AboutPt(tempXY, atlas::PointXY{135.0, 90.0}); + finalXY = rotatePlus180AboutPt(tempXY, PointXY{135.0, 90.0}); } else if (withinRange.y() >= 45.0) { - tempXY = rotatePlus90AboutPt(withinRange, atlas::PointXY{180.0, 45.0}); - finalXY = rotatePlus90AboutPt(tempXY, atlas::PointXY{90.0, 45.0}); + tempXY = rotatePlus90AboutPt(withinRange, PointXY{180.0, 45.0}); + finalXY = rotatePlus90AboutPt(tempXY, PointXY{90.0, 45.0}); } else if (withinRange.y() < -45.0) { - tempXY = rotateMinus90AboutPt(withinRange, atlas::PointXY{180.0, -45.0}); - finalXY = rotateMinus90AboutPt(tempXY, atlas::PointXY{90.0, -45.0}); + tempXY = rotateMinus90AboutPt(withinRange, PointXY{180.0, -45.0}); + finalXY = rotateMinus90AboutPt(tempXY, PointXY{90.0, -45.0}); } } break; @@ -532,47 +511,47 @@ atlas::PointXY LFRicCubedSphereTiles::tileCubePeriodicity(const atlas::PointXY& case 3: if (((withinRange.x() >= 270.0) && (withinRange.x() <= 360.0)) || withinRange.x() == 0.0) { if (withinRange.y() > 135.0) { - finalXY = rotatePlus180AboutPt(tempXY, atlas::PointXY{225.0, 90.0}); + finalXY = rotatePlus180AboutPt(tempXY, PointXY{225.0, 90.0}); } else if (withinRange.y() >= 45.0) { - finalXY = rotateMinus90AboutPt(tempXY, atlas::PointXY{360.0, 45.0}); + finalXY = rotateMinus90AboutPt(tempXY, PointXY{360.0, 45.0}); finalXY.x() -= 360.0; } else if (withinRange.y() < -45.0) { - finalXY = rotatePlus90AboutPt(tempXY, atlas::PointXY{360.0, -45.0}); + finalXY = rotatePlus90AboutPt(tempXY, PointXY{360.0, -45.0}); finalXY.x() -= 360.0; } } break; case 4: if (withinRange.y() > 135.0) { - finalXY = rotatePlus180AboutPt(tempXY, atlas::PointXY{135.0, 90.0}); + finalXY = rotatePlus180AboutPt(tempXY, PointXY{135.0, 90.0}); } else if ((withinRange.y() >= 45.0) && (withinRange.y() <= 135.0)) { if ((withinRange.x() > 90.0) && (withinRange.x() <= 180.0)) { - finalXY = rotateMinus90AboutPt(withinRange, atlas::PointXY{90.0, 45.0}); + finalXY = rotateMinus90AboutPt(withinRange, PointXY{90.0, 45.0}); } if ((withinRange.x() > 180.0) && (withinRange.x() <= 270.0)) { - finalXY = rotatePlus180AboutPt(withinRange, atlas::PointXY{135.0, 0.0}); + finalXY = rotatePlus180AboutPt(withinRange, PointXY{135.0, 0.0}); } if ((withinRange.x() > 270.0) && (withinRange.x() <= 360.0)) { - finalXY = rotatePlus90AboutPt(withinRange, atlas::PointXY{360.0, 45.0}); + finalXY = rotatePlus90AboutPt(withinRange, PointXY{360.0, 45.0}); } } break; case 5: if (withinRange.y() > 135.0) { - finalXY = rotatePlus180AboutPt(tempXY, atlas::PointXY{135.0, 90.0}); + finalXY = rotatePlus180AboutPt(tempXY, PointXY{135.0, 90.0}); } else if ((withinRange.y() <= -45.0) && (withinRange.y() >= -135.0)) { if ((withinRange.x() > 90.0) && (withinRange.x() <= 180.0)) { - finalXY = rotatePlus90AboutPt(withinRange, atlas::PointXY{90.0, -45.0}); + finalXY = rotatePlus90AboutPt(withinRange, PointXY{90.0, -45.0}); } if ((withinRange.x() > 180.0) && (withinRange.x() <= 270.0)) { - finalXY = rotatePlus180AboutPt(withinRange, atlas::PointXY{135.0, 0.0}); + finalXY = rotatePlus180AboutPt(withinRange, PointXY{135.0, 0.0}); } if ((withinRange.x() > 270.0) && (withinRange.x() <= 360.0)) { - finalXY = rotateMinus90AboutPt(withinRange, atlas::PointXY{360.0, -45.0}); + finalXY = rotateMinus90AboutPt(withinRange, PointXY{360.0, -45.0}); } } break; @@ -628,39 +607,39 @@ void LFRicCubedSphereTiles::print(std::ostream& os) const { } -bool LFRicCubedSphereTiles::withinCross(const atlas::idx_t tiles, const atlas::PointXY& withinRange) const { +bool LFRicCubedSphereTiles::withinCross(const idx_t tiles, const PointXY& withinRange) const { std::size_t t = static_cast(tiles); - return !((withinRange.x() < botLeftTile_[t].x() && withinRange.y() < botLeftTile_[t].y()) || - (withinRange.x() > botRightTile_[t].x() && withinRange.y() < botRightTile_[t].y()) || - (withinRange.x() > topRightTile_[t].x() && withinRange.y() > topRightTile_[t].y()) || - (withinRange.x() < topLeftTile_[t].x() && withinRange.y() > topLeftTile_[t].y())); + return !((withinRange.x() < botLeftTile(t).x() && withinRange.y() < botLeftTile(t).y()) || + (withinRange.x() > botRightTile(t).x() && withinRange.y() < botRightTile(t).y()) || + (withinRange.x() > topRightTile(t).x() && withinRange.y() > topRightTile(t).y()) || + (withinRange.x() < topLeftTile(t).x() && withinRange.y() > topLeftTile(t).y())); } -void LFRicCubedSphereTiles::enforceWrapAround(const atlas::idx_t t, atlas::PointXY& withinRange) const { +void LFRicCubedSphereTiles::enforceWrapAround(const idx_t t, PointXY& withinRange) const { if (withinRange.x() < 0.0) { - atlas::PointXY temp = withinRange; + PointXY temp = withinRange; temp.x() += 360; if (withinCross(t, temp)) { withinRange = temp; } } if (withinRange.x() > 360.0) { - atlas::PointXY temp = withinRange; + PointXY temp = withinRange; temp.x() -= 360; if (withinCross(t, temp)) { withinRange = temp; } } if (withinRange.y() <= -135.0) { - atlas::PointXY temp = withinRange; + PointXY temp = withinRange; temp.y() += 360; if (withinCross(t, temp)) { withinRange = temp; } } if (withinRange.y() > 225.0) { - atlas::PointXY temp = withinRange; + PointXY temp = withinRange; temp.y() -= 360; if (withinCross(t, temp)) { withinRange = temp; @@ -670,6 +649,41 @@ void LFRicCubedSphereTiles::enforceWrapAround(const atlas::idx_t t, atlas::Point return; } +const PointXY& LFRicCubedSphereTiles::tileCentre(size_t t) const { + return tileCentres_[t]; +} + +const LFRicCubedSphereTiles::Jacobian& LFRicCubedSphereTiles::tileJacobian(size_t t) const { + return tileJacobians_[t]; +} + +PointXY LFRicCubedSphereTiles::botLeftTile(size_t t) { + return tileCentres_[t] + PointXY{-45., -45.}; +} + +PointXY LFRicCubedSphereTiles::botRightTile(size_t t) { + return tileCentres_[t] + PointXY{45., -45.}; +} + +PointXY LFRicCubedSphereTiles::topLeftTile(size_t t) { + return tileCentres_[t] + PointXY{-45., 45.}; +} + +PointXY LFRicCubedSphereTiles::topRightTile(size_t t) { + return tileCentres_[t] + PointXY{45., 45.}; +} + + +// Centre of each tile in xy space. +const std::array LFRicCubedSphereTiles::tileCentres_{ + PointXY{45., 0.}, PointXY{135., 0.}, PointXY{225., 0.}, PointXY{315., 0.}, PointXY{45., 90.}, PointXY{45., -90.}}; + +// Jacobian of xy space with respect to curvilinear coordinates for each tile. +const std::array LFRicCubedSphereTiles::tileJacobians_{ + Jacobian{{1., 0.}, {0., 1.}}, Jacobian{{1., 0.}, {0., 1.}}, Jacobian{{0., -1.}, {1., 0.}}, + Jacobian{{0., -1.}, {1., 0.}}, Jacobian{{1., 0.}, {0., 1.}}, Jacobian{{0., 1.}, {-1., 0.}}}; + + namespace { static CubedSphereTilesBuilder register_builder(LFRicCubedSphereTiles::static_type()); } diff --git a/src/atlas/grid/detail/tiles/LFRicTiles.h b/src/atlas/grid/detail/tiles/LFRicTiles.h index c1c3aef81..a6c957387 100644 --- a/src/atlas/grid/detail/tiles/LFRicTiles.h +++ b/src/atlas/grid/detail/tiles/LFRicTiles.h @@ -43,20 +43,29 @@ class LFRicCubedSphereTiles : public CubedSphereTiles { virtual void enforceXYdomain(double xy[]) const override; - virtual atlas::PointXY tileCubePeriodicity(const atlas::PointXY& xyExtended, - const atlas::idx_t tile) const override; + virtual PointXY tileCubePeriodicity(const PointXY& xyExtended, const idx_t tile) const override; virtual void print(std::ostream&) const override; + virtual const PointXY& tileCentre(size_t t) const override; + + virtual const Jacobian& tileJacobian(size_t t) const override; + private: - std::array botLeftTile_; - std::array botRightTile_; - std::array topLeftTile_; - std::array topRightTile_; + static PointXY botLeftTile(size_t t); + static PointXY botRightTile(size_t t); + static PointXY topLeftTile(size_t t); + static PointXY topRightTile(size_t t); + + bool withinCross(const idx_t t, const PointXY& withinRange) const; + + void enforceWrapAround(const idx_t t, PointXY& withinRange) const; - bool withinCross(const atlas::idx_t t, const atlas::PointXY& withinRange) const; - void enforceWrapAround(const atlas::idx_t t, atlas::PointXY& withinRange) const; + // Centre of each tile in xy-space. + static const std::array tileCentres_; + // Jacobian of xy with respect to tile curvilinear coordinates. + static const std::array tileJacobians_; }; } // namespace detail diff --git a/src/atlas/grid/detail/tiles/Tiles.cc b/src/atlas/grid/detail/tiles/Tiles.cc index 3d824452a..0a4702bbd 100644 --- a/src/atlas/grid/detail/tiles/Tiles.cc +++ b/src/atlas/grid/detail/tiles/Tiles.cc @@ -22,7 +22,7 @@ namespace grid { namespace detail { const CubedSphereTiles* CubedSphereTiles::create() { - // default: FV3 version (for now) + // default: LFRic version (for now) util::Config params; params.set("type", "cubedsphere_lfric"); return CubedSphereTiles::create(params); @@ -49,6 +49,14 @@ const CubedSphereTiles* CubedSphereTiles::create(const eckit::Parametrisation& p } } +idx_t CubedSphereTiles::indexFromXY(const PointXY& xy) const { + return indexFromXY(xy.data()); +} + +idx_t CubedSphereTiles::indexFromLonLat(const PointLonLat& lonlat) const { + return indexFromXY(lonlat.data()); +} + } // namespace detail } // namespace grid } // namespace atlas diff --git a/src/atlas/grid/detail/tiles/Tiles.h b/src/atlas/grid/detail/tiles/Tiles.h index 74c9f304f..033749574 100644 --- a/src/atlas/grid/detail/tiles/Tiles.h +++ b/src/atlas/grid/detail/tiles/Tiles.h @@ -14,21 +14,27 @@ #include #include "atlas/library/config.h" +#include "atlas/projection/detail/ProjectionImpl.h" #include "atlas/util/Config.h" #include "atlas/util/Object.h" -#include "atlas/util/Point.h" namespace eckit { class Parametrisation; } + +namespace atlas { +class PointXY; +} // namespace atlas + namespace atlas { namespace grid { namespace detail { class CubedSphereTiles : public util::Object { public: - using Spec = util::Config; + using Spec = util::Config; + using Jacobian = projection::Jacobian; public: static const CubedSphereTiles* create(); @@ -51,8 +57,16 @@ class CubedSphereTiles : public util::Object { virtual idx_t indexFromLonLat(const double lonlat[]) const = 0; + idx_t indexFromXY(const PointXY& xy) const; + + idx_t indexFromLonLat(const PointLonLat& lonlat) const; + virtual void enforceXYdomain(double xy[]) const = 0; + virtual const PointXY& tileCentre(size_t t) const = 0; + + virtual const Jacobian& tileJacobian(size_t t) const = 0; + idx_t size() const { return 6; } virtual atlas::PointXY tileCubePeriodicity(const atlas::PointXY& xyExtended, const atlas::idx_t tile) const = 0; diff --git a/src/atlas/interpolation/Cache.cc b/src/atlas/interpolation/Cache.cc index e43c0b1df..bae12b89c 100644 --- a/src/atlas/interpolation/Cache.cc +++ b/src/atlas/interpolation/Cache.cc @@ -24,6 +24,20 @@ Cache::Cache(const Cache& other) { add(other); } +Cache::Cache(const Cache& other, const std::string& filter): Cache(other) { + std::shared_ptr filtered; + for (auto& entry : cache_) { + if (entry.first == filter) { + filtered = entry.second; + } + } + cache_.clear(); + if (filtered) { + cache_[filtered->type()] = filtered; + } +} + + Cache::Cache(const Interpolation& interpolation): Cache(interpolation.createCache()) {} Cache::~Cache() = default; @@ -49,7 +63,8 @@ class MatrixCacheEntryOwned : public MatrixCacheEntry { class MatrixCacheEntryShared : public MatrixCacheEntry { public: - MatrixCacheEntryShared(std::shared_ptr matrix): MatrixCacheEntry(matrix.get()), matrix_{matrix} {} + MatrixCacheEntryShared(std::shared_ptr matrix, const std::string& uid = ""): + MatrixCacheEntry(matrix.get(), uid), matrix_{matrix} {} private: const std::shared_ptr matrix_; @@ -57,11 +72,13 @@ class MatrixCacheEntryShared : public MatrixCacheEntry { MatrixCache::MatrixCache(const Cache& c): - Cache(c), matrix_{dynamic_cast(c.get(MatrixCacheEntry::static_type()))} {} + Cache(c, MatrixCacheEntry::static_type()), + matrix_{dynamic_cast(c.get(MatrixCacheEntry::static_type()))} {} MatrixCache::MatrixCache(Matrix&& m): MatrixCache(std::make_shared(std::move(m))) {} -MatrixCache::MatrixCache(std::shared_ptr m): MatrixCache(std::make_shared(m)) {} +MatrixCache::MatrixCache(std::shared_ptr m, const std::string& uid): + MatrixCache(std::make_shared(m, uid)) {} MatrixCache::MatrixCache(const Matrix* m): MatrixCache(std::make_shared(m)) {} @@ -76,6 +93,11 @@ const MatrixCache::Matrix& MatrixCache::matrix() const { return matrix_->matrix(); } +const std::string& MatrixCache::uid() const { + ATLAS_ASSERT(matrix_); + return matrix_->uid(); +} + size_t MatrixCache::footprint() const { if (matrix_) { return matrix_->footprint(); @@ -95,7 +117,8 @@ const IndexKDTreeCacheEntry::IndexKDTree& IndexKDTreeCacheEntry::tree() const { } IndexKDTreeCache::IndexKDTreeCache(const Cache& c): - Cache(c), tree_{dynamic_cast(c.get(IndexKDTreeCacheEntry::static_type()))} {} + Cache(c, IndexKDTreeCacheEntry::static_type()), + tree_{dynamic_cast(c.get(IndexKDTreeCacheEntry::static_type()))} {} IndexKDTreeCache::IndexKDTreeCache(const IndexKDTree& tree): Cache(std::make_shared(tree)) { tree_ = dynamic_cast(get(IndexKDTreeCacheEntry::static_type())); diff --git a/src/atlas/interpolation/Cache.h b/src/atlas/interpolation/Cache.h index ba0c5d3fa..3fdaebdfd 100644 --- a/src/atlas/interpolation/Cache.h +++ b/src/atlas/interpolation/Cache.h @@ -49,6 +49,7 @@ class Cache { public: Cache() = default; Cache(const Cache& other); + Cache(const Cache& other, const std::string& filter); Cache(const Interpolation&); operator bool() const { return not cache_.empty(); } virtual ~Cache(); @@ -83,8 +84,11 @@ class MatrixCacheEntry : public InterpolationCacheEntry { public: using Matrix = eckit::linalg::SparseMatrix; ~MatrixCacheEntry() override; - MatrixCacheEntry(const Matrix* matrix): matrix_{matrix} { ATLAS_ASSERT(matrix_ != nullptr); } + MatrixCacheEntry(const Matrix* matrix, const std::string& uid = ""): matrix_{matrix}, uid_(uid) { + ATLAS_ASSERT(matrix_ != nullptr); + } const Matrix& matrix() const { return *matrix_; } + const std::string& uid() const { return uid_; } size_t footprint() const override { return matrix_->footprint(); } operator bool() const { return not matrix_->empty(); } static std::string static_type() { return "Matrix"; } @@ -92,6 +96,7 @@ class MatrixCacheEntry : public InterpolationCacheEntry { private: const Matrix* matrix_; + const std::string uid_; }; //----------------------------------------------------------------------------- @@ -104,11 +109,12 @@ class MatrixCache final : public Cache { MatrixCache() = default; MatrixCache(const Cache& c); MatrixCache(Matrix&& m); - MatrixCache(std::shared_ptr m); + MatrixCache(std::shared_ptr m, const std::string& uid = ""); MatrixCache(const Matrix* m); MatrixCache(const Interpolation&); operator bool() const; const Matrix& matrix() const; + const std::string& uid() const; size_t footprint() const; private: diff --git a/src/atlas/interpolation/Interpolation.cc b/src/atlas/interpolation/Interpolation.cc index ae995d245..7a26b95b2 100644 --- a/src/atlas/interpolation/Interpolation.cc +++ b/src/atlas/interpolation/Interpolation.cc @@ -79,20 +79,20 @@ Interpolation::Interpolation(const Interpolation::Config& config, const Function } } -void Interpolation::execute(const FieldSet& source, FieldSet& target) const { - get()->execute(source, target); +Interpolation::Metadata Interpolation::execute(const FieldSet& source, FieldSet& target) const { + return get()->execute(source, target); } -void Interpolation::execute(const Field& source, Field& target) const { - get()->execute(source, target); +Interpolation::Metadata Interpolation::execute(const Field& source, Field& target) const { + return get()->execute(source, target); } -void Interpolation::execute_adjoint(FieldSet& source, const FieldSet& target) const { - get()->execute_adjoint(source, target); +Interpolation::Metadata Interpolation::execute_adjoint(FieldSet& source, const FieldSet& target) const { + return get()->execute_adjoint(source, target); } -void Interpolation::execute_adjoint(Field& source, const Field& target) const { - get()->execute_adjoint(source, target); +Interpolation::Metadata Interpolation::execute_adjoint(Field& source, const Field& target) const { + return get()->execute_adjoint(source, target); } void Interpolation::print(std::ostream& out) const { diff --git a/src/atlas/interpolation/Interpolation.h b/src/atlas/interpolation/Interpolation.h index 1aac5ad6f..50144b4b9 100644 --- a/src/atlas/interpolation/Interpolation.h +++ b/src/atlas/interpolation/Interpolation.h @@ -34,8 +34,9 @@ namespace atlas { class Interpolation : DOXYGEN_HIDE(public util::ObjectHandle) { public: - using Config = eckit::Parametrisation; - using Cache = interpolation::Cache; + using Config = eckit::Parametrisation; + using Cache = interpolation::Cache; + using Metadata = interpolation::Method::Metadata; using Handle::Handle; Interpolation() = default; @@ -52,13 +53,13 @@ class Interpolation : DOXYGEN_HIDE(public util::ObjectHandle bool { - if (std::abs(a) > epsilon) { - // if a is non-zero, we need to solve a quadratic equation for the weight - // ax**2 + bx + x = 0 - double det = b * b - 4. * a * c; - if (det >= 0.) { - double inv_two_a = 1. / (2. * a); - double sqrt_det = std::sqrt(det); - double root_a = (-b + sqrt_det) * inv_two_a; - if ((root_a > -epsilon) && (root_a < (1. + epsilon))) { - wght = root_a; - return true; - } - double root_b = (-b - sqrt_det) * inv_two_a; - if ((root_b > -epsilon) && (root_b < (1. + epsilon))) { - wght = root_b; + auto solve_weight = [&](double a, double b, double c, double& weight) -> bool { + auto checkWeight = [](double weight, double tol) -> bool { return ((weight > -tol) && (weight < (1. + tol))); }; + + // Quadratic equation ax^2 + bx + c = 0. + if (std::abs(a) >= areaEpsilon) { + // Solve numerically stable form of quadratic formula: + // x1 = (-b - sign(b) sqrt(b^2 - 4ac)) / 2a + // x2 = c / ax1 + + // Kahan's algorithm for accurate difference of products. + const auto prodDiff = [](double a, double b, double c, double d) { + // return ab - cd + const double w = d * c; + const double e = std::fma(-d, c, w); + const double f = std::fma(a, b, -w); + return f + e; + }; + + double discriminant = prodDiff(b, b, 4. * a, c); + + if (discriminant > -areaEpsilon * 2. * quadArea) { + // Solution is real. + + double sqrtDiscriminant = std::sqrt(std::max(0., discriminant)); + + // "Classic" solution to quadratic formula with no cancelation + // on numerator. + weight = (-b - std::copysign(sqrtDiscriminant, b)) / (2. * a); + if (checkWeight(weight, edgeEpsilon)) { return true; } + + // Use Vieta's formula x1 * x2 = c / a; + weight = c / (a * weight); + return checkWeight(weight, edgeEpsilon); } - return false; } - if (std::abs(b) > epsilon) { - // solve ax + b = 0 - wght = -c / b; - return true; + else if (std::abs(b) >= areaEpsilon) { + // Linear case bx + c = 0. + weight = -c / b; + return checkWeight(weight, edgeEpsilon); } + + // No real solutions to equation. return false; }; @@ -91,8 +115,8 @@ method::Intersect Quad2D::localRemap(const PointXY& p, double edgeEpsilon, doubl // w3 = u * v // w4 = ( 1 - u ) * v - Vector2D ray(p.x(), p.y()); - Vector2D vA = v00 - ray; + Vector2D point(p.x(), p.y()); + Vector2D vA = v00 - point; Vector2D vB = v10 - v00; Vector2D vC = v01 - v00; Vector2D vD = v00 - v10 - v01 + v11; @@ -164,10 +188,10 @@ double Quad2D::area() const { return std::abs(0.5 * cross2d((v01 - v00), (v11 - v00))) + std::abs(0.5 * cross2d((v10 - v11), (v01 - v11))); } -bool Quad2D::inQuadrilateral(const Vector2D& p) const { +bool Quad2D::inQuadrilateral(const Vector2D& p, double tolerance) const { // point p must be on the inside of all quad edges to be inside the quad. - return cross2d(p - v00, p - v10) >= 0. && cross2d(p - v10, p - v11) >= 0. && cross2d(p - v11, p - v01) >= 0. && - cross2d(p - v01, p - v00) >= 0; + return cross2d(p - v00, p - v10) > -tolerance && cross2d(p - v10, p - v11) > -tolerance && + cross2d(p - v11, p - v01) > -tolerance && cross2d(p - v01, p - v00) > -tolerance; } void Quad2D::print(std::ostream& s) const { diff --git a/src/atlas/interpolation/element/Quad2D.h b/src/atlas/interpolation/element/Quad2D.h index 46c44daca..1e0640ee6 100644 --- a/src/atlas/interpolation/element/Quad2D.h +++ b/src/atlas/interpolation/element/Quad2D.h @@ -12,6 +12,7 @@ #include "atlas/interpolation/Vector2D.h" #include "atlas/interpolation/method/Intersect.h" +#include "atlas/runtime/Exception.h" #include "atlas/util/Point.h" namespace atlas { @@ -28,22 +29,38 @@ class Quad2D { Quad2D(const double* x0, const double* x1, const double* x2, const double* x3): v00(x0), v10(x1), v11(x2), v01(x3) {} - Quad2D(const PointXY& x0, const PointXY& x1, const PointXY& x2, const PointXY& x3): + Quad2D(const Point2& x0, const Point2& x1, const Point2& x2, const Point2& x3): Quad2D(x0.data(), x1.data(), x2.data(), x3.data()) {} Quad2D(const Vector2D& x0, const Vector2D& x1, const Vector2D& x2, const Vector2D& x3): Quad2D(x0.data(), x1.data(), x2.data(), x3.data()) {} - method::Intersect intersects(const PointXY& r, double edgeEpsilon = 5 * std::numeric_limits::epsilon(), + method::Intersect intersects(const Point2& r, double edgeEpsilon = 5 * std::numeric_limits::epsilon(), double epsilon = 5 * std::numeric_limits::epsilon()) const; - method::Intersect localRemap(const PointXY& r, double edgeEpsilon = 5 * std::numeric_limits::epsilon(), + method::Intersect localRemap(const Point2& r, double edgeEpsilon = 5 * std::numeric_limits::epsilon(), double epsilon = 5 * std::numeric_limits::epsilon()) const; bool validate() const; double area() const; + const Vector2D& p(int i) { + if (i == 0) { + return v00; + } + if (i == 1) { + return v10; + } + if (i == 2) { + return v11; + } + if (i == 3) { + return v01; + } + throw_OutOfRange("Quad2D::p(i)", i, 4, Here()); + } + void print(std::ostream&) const; friend std::ostream& operator<<(std::ostream& s, const Quad2D& p) { @@ -57,9 +74,9 @@ class Quad2D { Vector2D v11; // aka v2 Vector2D v01; // aka v3 - static double cross2d(const Vector2D& a, const Vector2D& b) { return a.x() * b.y() - a.y() * b.x(); } + bool inQuadrilateral(const Vector2D& p, double tolerance = 5 * std::numeric_limits::epsilon()) const; - bool inQuadrilateral(const Vector2D& p) const; + static double cross2d(const Vector2D& a, const Vector2D& b) { return a.x() * b.y() - a.y() * b.x(); } }; //---------------------------------------------------------------------------------------------------------------------- diff --git a/src/atlas/interpolation/element/Triag2D.cc b/src/atlas/interpolation/element/Triag2D.cc index 870ba3262..9210762e6 100644 --- a/src/atlas/interpolation/element/Triag2D.cc +++ b/src/atlas/interpolation/element/Triag2D.cc @@ -32,7 +32,8 @@ method::Intersect Triag2D::intersects(const PointXY& r, double edgeEpsilon, doub Vector2D rvec{r.data()}; - if (!inTriangle(rvec)) { + + if (!inTriangle(rvec, epsilon * area())) { return isect.fail(); } @@ -45,11 +46,13 @@ method::Intersect Triag2D::intersects(const PointXY& r, double edgeEpsilon, doub isect.u = (pvec.x() * e2.y() - e2.x() * pvec.y()) * invDet; isect.v = (e1.x() * pvec.y() - pvec.x() * e1.y()) * invDet; - // clamp values between 0 and 1 - isect.u = std::max(0.0, std::min(isect.u, 1.0)); - isect.v = std::max(0.0, std::min(isect.v, 1.0)); - - return isect.success(); + // make sure weights are valid + if (isect.u > -edgeEpsilon || isect.u < 1. + edgeEpsilon || isect.v > -edgeEpsilon || isect.v < 1. + edgeEpsilon) { + return isect.success(); + } + else { + return isect.fail(); + } } bool Triag2D::validate() const { @@ -89,20 +92,10 @@ double Triag2D::area() const { return std::abs(0.5 * cross2d((v10 - v00), (v11 - v00))); } -bool Triag2D::inTriangle(const Vector2D& p) const { +bool Triag2D::inTriangle(const Vector2D& p, double tolerance) const { // point p must be on the inside of all triangle edges to be inside the triangle. - double zst1 = cross2d(p - v00, p - v10); - if (zst1 >= 0.0) { - double zst2 = cross2d(p - v10, p - v11); - if (zst2 >= 0.0) { - double zst3 = cross2d(p - v11, p - v00); - if (zst3 >= 0.0) { - return true; - } - } - } - - return false; + return cross2d(p - v00, p - v10) > -tolerance && cross2d(p - v10, p - v11) > -tolerance && + cross2d(p - v11, p - v00) > -tolerance; } void Triag2D::print(std::ostream& s) const { diff --git a/src/atlas/interpolation/element/Triag2D.h b/src/atlas/interpolation/element/Triag2D.h index 9d76f2532..84b0586c5 100644 --- a/src/atlas/interpolation/element/Triag2D.h +++ b/src/atlas/interpolation/element/Triag2D.h @@ -62,7 +62,7 @@ class Triag2D { static double cross2d(const Vector2D& a, const Vector2D& b) { return a.x() * b.y() - a.y() * b.x(); } - bool inTriangle(const Vector2D& p) const; + bool inTriangle(const Vector2D& p, double tolerance = 5 * std::numeric_limits::epsilon()) const; }; //---------------------------------------------------------------------------------------------------------------------- diff --git a/src/atlas/interpolation/method/Method.cc b/src/atlas/interpolation/method/Method.cc index 3f0907d8d..5c402cb07 100644 --- a/src/atlas/interpolation/method/Method.cc +++ b/src/atlas/interpolation/method/Method.cc @@ -97,24 +97,61 @@ void set_missing_values(Field& tgt, const std::vector& missing) { } // anonymous namespace + template void Method::interpolate_field_rank1(const Field& src, Field& tgt, const Matrix& W) const { - auto src_v = array::make_view(src); - auto tgt_v = array::make_view(tgt); - if (std::is_same::value) { - sparse_matrix_multiply(W, src_v, tgt_v, sparse::backend::openmp()); + auto backend = std::is_same::value ? sparse::backend::openmp() : sparse::Backend{linalg_backend_}; + auto src_v = array::make_view(src); + auto tgt_v = array::make_view(tgt); + + if (nonLinear_(src)) { + Matrix W_nl(W); // copy (a big penalty -- copy-on-write would definitely be better) + nonLinear_->execute(W_nl, src); + sparse_matrix_multiply(W_nl, src_v, tgt_v, backend); } else { - sparse_matrix_multiply(W, src_v, tgt_v, sparse::Backend{linalg_backend_}); + sparse_matrix_multiply(W, src_v, tgt_v, backend); } } + template void Method::interpolate_field_rank2(const Field& src, Field& tgt, const Matrix& W) const { sparse::Backend backend{linalg_backend_}; auto src_v = array::make_view(src); auto tgt_v = array::make_view(tgt); - sparse_matrix_multiply(W, src_v, tgt_v, sparse::backend::openmp()); + + if (nonLinear_(src)) { + // We cannot apply the same matrix to full columns as e.g. missing values could be present in only certain parts. + + // Allocate temporary rank-1 fields corresponding to one horizontal level + auto src_slice = Field("s", array::make_datatype(), {src.shape(0)}); + auto tgt_slice = Field("t", array::make_datatype(), {tgt.shape(0)}); + + // Copy metadata to the source rank-1 field + src_slice.metadata() = src.metadata(); + + auto src_slice_v = array::make_view(src_slice); + auto tgt_slice_v = array::make_view(tgt_slice); + + for (idx_t lev = 0; lev < src_v.shape(1); ++lev) { + // Copy this level to temporary rank-1 field + for (idx_t i = 0; i < src.shape(0); ++i) { + src_slice_v(i) = src_v(i, lev); + } + + // Interpolate between rank-1 fields + interpolate_field_rank1(src_slice, tgt_slice, W); + + // Copy rank-1 field to this level in the rank-2 field + for (idx_t i = 0; i < tgt.shape(0); ++i) { + tgt_v(i, lev) = tgt_slice_v(i); + } + } + } + else { + sparse_matrix_multiply(W, src_v, tgt_v, sparse::backend::openmp()); + } } @@ -123,6 +160,9 @@ void Method::interpolate_field_rank3(const Field& src, Field& tgt, const Matrix& sparse::Backend backend{linalg_backend_}; auto src_v = array::make_view(src); auto tgt_v = array::make_view(tgt); + if (not W.empty() && nonLinear_(src)) { + ATLAS_ASSERT(false, "nonLinear interpolation not supported for rank-3 fields."); + } sparse_matrix_multiply(W, src_v, tgt_v, sparse::backend::openmp()); } @@ -243,9 +283,6 @@ Method::Method(const Method::Config& config) { if (config.get("non_linear", non_linear)) { nonLinear_ = NonLinear(non_linear, config); } - matrix_shared_ = std::make_shared(); - matrix_cache_ = interpolation::MatrixCache(matrix_shared_); - matrix_ = matrix_shared_.get(); config.get("adjoint", adjoint_); } @@ -284,24 +321,32 @@ void Method::setup(const Grid& source, const Grid& target, const Cache& cache) { this->do_setup(source, target, cache); } -void Method::execute(const FieldSet& source, FieldSet& target) const { +Method::Metadata Method::execute(const FieldSet& source, FieldSet& target) const { ATLAS_TRACE("atlas::interpolation::method::Method::execute(FieldSet, FieldSet)"); - this->do_execute(source, target); + Metadata metadata; + this->do_execute(source, target, metadata); + return metadata; } -void Method::execute(const Field& source, Field& target) const { +Method::Metadata Method::execute(const Field& source, Field& target) const { ATLAS_TRACE("atlas::interpolation::method::Method::execute(Field, Field)"); - this->do_execute(source, target); + Metadata metadata; + this->do_execute(source, target, metadata); + return metadata; } -void Method::execute_adjoint(FieldSet& source, const FieldSet& target) const { +Method::Metadata Method::execute_adjoint(FieldSet& source, const FieldSet& target) const { ATLAS_TRACE("atlas::interpolation::method::Method::execute(FieldSet, FieldSet)"); - this->do_execute_adjoint(source, target); + Metadata metadata; + this->do_execute_adjoint(source, target, metadata); + return metadata; } -void Method::execute_adjoint(Field& source, const Field& target) const { +Method::Metadata Method::execute_adjoint(Field& source, const Field& target) const { ATLAS_TRACE("atlas::interpolation::method::Method::execute(Field, Field)"); - this->do_execute_adjoint(source, target); + Metadata metadata; + this->do_execute_adjoint(source, target, metadata); + return metadata; } void Method::do_setup(const FunctionSpace& /*source*/, const Field& /*target*/) { @@ -312,36 +357,27 @@ void Method::do_setup(const FunctionSpace& /*source*/, const FieldSet& /*target* ATLAS_NOTIMPLEMENTED; } -void Method::do_execute(const FieldSet& fieldsSource, FieldSet& fieldsTarget) const { +void Method::do_execute(const FieldSet& fieldsSource, FieldSet& fieldsTarget, Metadata& metadata) const { ATLAS_TRACE("atlas::interpolation::method::Method::do_execute()"); const idx_t N = fieldsSource.size(); ATLAS_ASSERT(N == fieldsTarget.size()); for (idx_t i = 0; i < fieldsSource.size(); ++i) { - Method::do_execute(fieldsSource[i], fieldsTarget[i]); + Method::do_execute(fieldsSource[i], fieldsTarget[i], metadata); } } -void Method::do_execute(const Field& src, Field& tgt) const { +void Method::do_execute(const Field& src, Field& tgt, Metadata&) const { ATLAS_TRACE("atlas::interpolation::method::Method::do_execute()"); haloExchange(src); - // non-linearities: a non-empty M matrix contains the corrections applied to matrix_ - Matrix M; - if (not matrix_->empty() && nonLinear_(src)) { - Matrix W(*matrix_); // copy (a big penalty -- copy-on-write would definitely be better) - if (nonLinear_->execute(W, src)) { - M.swap(W); - } - } - if (src.datatype().kind() == array::DataType::KIND_REAL64) { - interpolate_field(src, tgt, M.empty() ? *matrix_ : M); + interpolate_field(src, tgt, *matrix_); } else if (src.datatype().kind() == array::DataType::KIND_REAL32) { - interpolate_field(src, tgt, M.empty() ? *matrix_ : M); + interpolate_field(src, tgt, *matrix_); } else { ATLAS_NOTIMPLEMENTED; @@ -368,18 +404,18 @@ void Method::do_execute(const Field& src, Field& tgt) const { tgt.set_dirty(); } -void Method::do_execute_adjoint(FieldSet& fieldsSource, const FieldSet& fieldsTarget) const { +void Method::do_execute_adjoint(FieldSet& fieldsSource, const FieldSet& fieldsTarget, Metadata& metadata) const { ATLAS_TRACE("atlas::interpolation::method::Method::do_execute_adjoint()"); const idx_t N = fieldsSource.size(); ATLAS_ASSERT(N == fieldsTarget.size()); for (idx_t i = 0; i < fieldsSource.size(); ++i) { - Method::do_execute_adjoint(fieldsSource[i], fieldsTarget[i]); + Method::do_execute_adjoint(fieldsSource[i], fieldsTarget[i], metadata); } } -void Method::do_execute_adjoint(Field& src, const Field& tgt) const { +void Method::do_execute_adjoint(Field& src, const Field& tgt, Metadata&) const { ATLAS_TRACE("atlas::interpolation::method::Method::do_execute_adjoint()"); if (nonLinear_(src)) { diff --git a/src/atlas/interpolation/method/Method.h b/src/atlas/interpolation/method/Method.h index 4820b1bd0..5a0f6eac1 100644 --- a/src/atlas/interpolation/method/Method.h +++ b/src/atlas/interpolation/method/Method.h @@ -16,6 +16,7 @@ #include "atlas/interpolation/Cache.h" #include "atlas/interpolation/NonLinear.h" +#include "atlas/util/Metadata.h" #include "atlas/util/Object.h" #include "eckit/config/Configuration.h" #include "eckit/linalg/SparseMatrix.h" @@ -38,7 +39,8 @@ namespace interpolation { class Method : public util::Object { public: - typedef eckit::Parametrisation Config; + using Config = eckit::Parametrisation; + using Metadata = util::Metadata; Method(const Config&); virtual ~Method() {} @@ -54,8 +56,8 @@ class Method : public util::Object { void setup(const FunctionSpace& source, const FieldSet& target); void setup(const Grid& source, const Grid& target, const Cache&); - void execute(const FieldSet& source, FieldSet& target) const; - void execute(const Field& source, Field& target) const; + Metadata execute(const FieldSet& source, FieldSet& target) const; + Metadata execute(const Field& source, Field& target) const; /** * @brief execute_adjoint @@ -66,8 +68,8 @@ class Method : public util::Object { * to zero. This is not done for efficiency reasons and * because in most cases it is not necessary. */ - void execute_adjoint(FieldSet& source, const FieldSet& target) const; - void execute_adjoint(Field& source, const Field& target) const; + Metadata execute_adjoint(FieldSet& source, const FieldSet& target) const; + Metadata execute_adjoint(Field& source, const Field& target) const; virtual void print(std::ostream&) const = 0; @@ -77,11 +79,11 @@ class Method : public util::Object { virtual interpolation::Cache createCache() const; protected: - virtual void do_execute(const FieldSet& source, FieldSet& target) const; - virtual void do_execute(const Field& source, Field& target) const; + virtual void do_execute(const FieldSet& source, FieldSet& target, Metadata&) const; + virtual void do_execute(const Field& source, Field& target, Metadata&) const; - virtual void do_execute_adjoint(FieldSet& source, const FieldSet& target) const; - virtual void do_execute_adjoint(Field& source, const Field& target) const; + virtual void do_execute_adjoint(FieldSet& source, const FieldSet& target, Metadata&) const; + virtual void do_execute_adjoint(Field& source, const Field& target, Metadata&) const; using Triplet = eckit::linalg::Triplet; using Triplets = std::vector; @@ -98,18 +100,28 @@ class Method : public util::Object { // NOTE : Matrix-free or non-linear interpolation operators do not have matrices, so do not expose here friend class atlas::test::Access; friend class interpolation::MatrixCache; - const Matrix* matrix_; - std::shared_ptr matrix_shared_; - interpolation::MatrixCache matrix_cache_; - NonLinear nonLinear_; - std::string linalg_backend_; - bool allow_halo_exchange_{true}; - std::vector missing_; - bool adjoint_{false}; - Matrix matrix_transpose_; - protected: + void setMatrix(Matrix& m, const std::string& uid = "") { + if (not matrix_shared_) { + matrix_shared_ = std::make_shared(); + } + matrix_shared_->swap(m); + matrix_cache_ = interpolation::MatrixCache(matrix_shared_, uid); + matrix_ = &matrix_cache_.matrix(); + } + + void setMatrix(interpolation::MatrixCache matrix_cache) { + ATLAS_ASSERT(matrix_cache); + matrix_cache_ = matrix_cache; + matrix_ = &matrix_cache_.matrix(); + matrix_shared_.reset(); + } + + bool matrixAllocated() const { return matrix_shared_.use_count(); } + + const Matrix& matrix() const { return *matrix_; } + virtual void do_setup(const FunctionSpace& source, const FunctionSpace& target) = 0; virtual void do_setup(const Grid& source, const Grid& target, const Cache&) = 0; virtual void do_setup(const FunctionSpace& source, const Field& target); @@ -141,6 +153,19 @@ class Method : public util::Object { 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_; + std::shared_ptr matrix_shared_; + interpolation::MatrixCache matrix_cache_; + NonLinear nonLinear_; + std::string linalg_backend_; + bool adjoint_{false}; + Matrix matrix_transpose_; + +protected: + bool allow_halo_exchange_{true}; + std::vector missing_; }; } // namespace interpolation diff --git a/src/atlas/interpolation/method/MethodFactory.cc b/src/atlas/interpolation/method/MethodFactory.cc index 8e3fd8970..7c9f6952c 100644 --- a/src/atlas/interpolation/method/MethodFactory.cc +++ b/src/atlas/interpolation/method/MethodFactory.cc @@ -11,8 +11,7 @@ #include "MethodFactory.h" // for static linking -#include "bil/BilinearRemapping.h" -#include "fe/FiniteElement.h" +#include "cubedsphere/CubedSphereBilinear.h" #include "knn/GridBoxAverage.h" #include "knn/GridBoxMaximum.h" #include "knn/KNearestNeighbours.h" @@ -23,6 +22,9 @@ #include "structured/Linear3D.h" #include "structured/QuasiCubic2D.h" #include "structured/QuasiCubic3D.h" +#include "unstructured/FiniteElement.h" +#include "unstructured/UnstructuredBilinearLonLat.h" + namespace atlas { namespace interpolation { @@ -32,7 +34,7 @@ namespace { void force_link() { static struct Link { Link() { - MethodBuilder(); + MethodBuilder(); MethodBuilder(); MethodBuilder(); MethodBuilder(); @@ -44,6 +46,7 @@ void force_link() { MethodBuilder(); MethodBuilder(); MethodBuilder(); + MethodBuilder(); } } link; } diff --git a/src/atlas/interpolation/method/cubedsphere/CellFinder.cc b/src/atlas/interpolation/method/cubedsphere/CellFinder.cc new file mode 100644 index 000000000..2c9c4477b --- /dev/null +++ b/src/atlas/interpolation/method/cubedsphere/CellFinder.cc @@ -0,0 +1,134 @@ +/* + * (C) Crown Copyright 2022 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/cubedsphere/CellFinder.h" +#include "atlas/grid/CubedSphereGrid.h" +#include "atlas/interpolation/Vector2D.h" +#include "atlas/interpolation/element/Quad2D.h" +#include "atlas/interpolation/element/Triag2D.h" +#include "atlas/mesh/HybridElements.h" +#include "atlas/mesh/Nodes.h" +#include "atlas/util/Topology.h" + +namespace atlas { +namespace interpolation { +namespace method { +namespace cubedsphere { + +CellFinder::CellFinder(const Mesh& mesh, const util::Config& config): mesh_{mesh} { + // Check mesh and get projection. + const auto csGrid = CubedSphereGrid(mesh_.grid()); + ATLAS_ASSERT_MSG(csGrid, "cubedsphere::CellFinder requires a cubed sphere mesh."); + projection_ = &(csGrid.cubedSphereProjection()); + + // Get views to cell data. + const auto xyView = array::make_view(mesh_.cells().field("xy")); + const auto tijView = array::make_view(mesh_.cells().field("tij")); + const auto haloView = array::make_view(mesh_.cells().halo()); + + // make points and payloads vectors. + auto points = std::array, 6>{}; + auto payloads = std::array, 6>{}; + + // Iterate over cells. + auto halo = config.getInt("halo", 0); + for (idx_t i = 0; i < mesh_.cells().size(); ++i) { + if (haloView(i) <= halo) { + const auto t = tijView(i, 0); + const auto xy = PointXY(xyView(i, XX), xyView(i, YY)); + + points[size_t(t)].push_back(projection_->xy2alphabeta(xy, t)); + payloads[size_t(t)].push_back(i); + } + } + + // build cell kd-trees. + for (size_t t = 0; t < 6; ++t) { + trees_[t].build(points[t], payloads[t]); + } +} + +CellFinder::Cell CellFinder::getCell(const PointXY& xy, size_t listSize, double edgeEpsilon, double epsilon) const { + // Convert xy to alphabeta and t; + const auto& tiles = projection_->getCubedSphereTiles(); + const auto t = tiles.indexFromXY(xy.data()); + const auto alphabeta = projection_->xy2alphabeta(xy, t); + + + // Get mesh nodes and connectivity table. + const auto nodeXyView = array::make_view(mesh_.nodes().xy()); + const auto& nodeConnectivity = mesh_.cells().node_connectivity(); + + // Get four nearest cell-centres to xy. + const auto& tree = trees_[size_t(t)]; + if (tree.size() == 0) { + return Cell{{}, Intersect{}.fail()}; + } + + const auto valueList = tree.closestPoints(alphabeta, std::min(listSize, tree.size())); + + // Get view to cell flags. + const auto cellFlagsView = array::make_view(mesh_.cells().flags()); + + for (const auto& value : valueList) { + const auto i = value.payload(); + + // Ignore invalid cells. + if (!(cellFlagsView(i) & Topology::INVALID)) { + const auto row = nodeConnectivity.row(i); + + if (row.size() == 4) { + auto quadAlphabeta = std::array{}; + for (size_t k = 0; k < 4; ++k) { + const auto xyNode = PointXY(nodeXyView(row(k), XX), nodeXyView(row(k), YY)); + quadAlphabeta[k] = Vector2D(projection_->xy2alphabeta(xyNode, t).data()); + } + + const auto quad = + element::Quad2D(quadAlphabeta[0], quadAlphabeta[1], quadAlphabeta[2], quadAlphabeta[3]); + + const auto isect = quad.localRemap(alphabeta, edgeEpsilon, epsilon); + + if (isect) { + return Cell{{row(0), row(1), row(2), row(3)}, isect}; + } + } + else { + // Cell is triangle. + auto triagAlphabeta = std::array{}; + for (size_t k = 0; k < 3; ++k) { + const auto xyNode = PointXY(nodeXyView(row(k), XX), nodeXyView(row(k), YY)); + triagAlphabeta[k] = Vector2D(projection_->xy2alphabeta(xyNode, t).data()); + } + + const auto triag = element::Triag2D(triagAlphabeta[0], triagAlphabeta[1], triagAlphabeta[2]); + + const auto isect = triag.intersects(alphabeta, edgeEpsilon, epsilon); + + if (isect) { + return Cell{{row(0), row(1), row(2)}, isect}; + } + } + } + } + + // Couldn't find a cell. + return Cell{{}, Intersect{}.fail()}; +} + +CellFinder::Cell CellFinder::getCell(const PointLonLat& lonlat, size_t listSize, double edgeEpsilon, + double epsilon) const { + auto xy = PointXY(lonlat.data()); + projection_->lonlat2xy(xy.data()); + return getCell(xy, listSize, edgeEpsilon, epsilon); +} + + +} // namespace cubedsphere +} // namespace method +} // namespace interpolation +} // namespace atlas diff --git a/src/atlas/interpolation/method/cubedsphere/CellFinder.h b/src/atlas/interpolation/method/cubedsphere/CellFinder.h new file mode 100644 index 000000000..76932ba99 --- /dev/null +++ b/src/atlas/interpolation/method/cubedsphere/CellFinder.h @@ -0,0 +1,66 @@ +/* + * (C) Crown Copyright 2022 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/Intersect.h" +#include "atlas/mesh/Mesh.h" +#include "atlas/util/Config.h" +#include "atlas/util/KDTree.h" +#include "atlas/util/Point.h" + + +namespace atlas { +namespace projection { +namespace detail { +class CubedSphereProjectionBase; +} // namespace detail +} // namespace projection +} // namespace atlas + + +namespace atlas { +namespace interpolation { +namespace method { +namespace cubedsphere { + +using namespace util; + +/// @brief class to find points within cells of cubedsphere mesh. +class CellFinder { +public: + struct Cell { + std::vector nodes; + Intersect isect; + }; + + /// @brief Constructor. + CellFinder(const Mesh& mesh, const util::Config& config = util::Config("halo", 0)); + + /// @brief Find a cell which encompasses an xy point. + Cell getCell(const PointXY& xy, size_t listSize = 4, + double edgeEpsilon = 5. * std::numeric_limits::epsilon(), + double epsilon = 5. * std::numeric_limits::epsilon()) const; + + /// @brief Find a cell which encompasses a lonlat point. + Cell getCell(const PointLonLat& lonlat, size_t listSize = 4, + double edgeEpsilon = 5. * std::numeric_limits::epsilon(), + double epsilon = 5. * std::numeric_limits::epsilon()) const; + + +private: + Mesh mesh_{}; + const projection::detail::CubedSphereProjectionBase* projection_{}; + std::array trees_{}; +}; + +} // namespace cubedsphere +} // namespace method +} // namespace interpolation +} // namespace atlas diff --git a/src/atlas/interpolation/method/cubedsphere/CubedSphereBilinear.cc b/src/atlas/interpolation/method/cubedsphere/CubedSphereBilinear.cc new file mode 100644 index 000000000..da39e1ecf --- /dev/null +++ b/src/atlas/interpolation/method/cubedsphere/CubedSphereBilinear.cc @@ -0,0 +1,99 @@ +/* + * (C) Crown Copyright 2022 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/cubedsphere/CubedSphereBilinear.h" +#include "atlas/functionspace/NodeColumns.h" +#include "atlas/grid/CubedSphereGrid.h" +#include "atlas/interpolation/method/MethodFactory.h" +#include "atlas/interpolation/method/cubedsphere/CellFinder.h" +#include "atlas/util/CoordinateEnums.h" + +namespace atlas { +namespace interpolation { +namespace method { + +namespace { +MethodBuilder __builder("cubedsphere-bilinear"); +} // namespace + +void CubedSphereBilinear::do_setup(const Grid& source, const Grid& target, const Cache&) { + ATLAS_NOTIMPLEMENTED; +} + +void CubedSphereBilinear::do_setup(const FunctionSpace& source, const FunctionSpace& target) { + source_ = source; + target_ = target; + + const auto ncSource = functionspace::NodeColumns(source); + + ATLAS_ASSERT(ncSource); + ATLAS_ASSERT(target_); + + const auto finder = cubedsphere::CellFinder(ncSource.mesh(), util::Config("halo", halo_)); + + // Numeric tolerance should scale with N. + const auto N = CubedSphereGrid(ncSource.mesh().grid()).N(); + const auto tolerance = 2. * std::numeric_limits::epsilon() * N; + + // Enable or disable halo exchange. + this->allow_halo_exchange_ = halo_exchange_; + + // Loop over target at calculate interpolation weights. + auto weights = std::vector{}; + const auto ghostView = array::make_view(target_.ghost()); + const auto lonlatView = array::make_view(target_.lonlat()); + + for (idx_t i = 0; i < target_.size(); ++i) { + if (!ghostView(i)) { + const auto cell = + finder.getCell(PointLonLat(lonlatView(i, LON), lonlatView(i, LAT)), listSize_, tolerance, tolerance); + + if (!cell.isect) { + ATLAS_THROW_EXCEPTION( + "Cannot find a cell surrounding target" + "point " + + std::to_string(i) + "."); + } + + const auto& isect = cell.isect; + const auto& j = cell.nodes; + + switch (cell.nodes.size()) { + case (3): { + // Cell is a triangle. + weights.emplace_back(i, j[0], 1. - isect.u - isect.v); + weights.emplace_back(i, j[1], isect.u); + weights.emplace_back(i, j[2], isect.v); + break; + } + case (4): { + // Cell is quad. + weights.emplace_back(i, j[0], (1. - isect.u) * (1. - isect.v)); + weights.emplace_back(i, j[1], isect.u * (1. - isect.v)); + weights.emplace_back(i, j[2], isect.u * isect.v); + weights.emplace_back(i, j[3], (1. - isect.u) * isect.v); + break; + } + default: { + ATLAS_THROW_EXCEPTION("Unknown cell type with " + std::to_string(cell.nodes.size()) + " nodes."); + } + } + } + + // fill sparse matrix and return. + Matrix A(target_.size(), source_.size(), weights); + setMatrix(A); + } +} + +void CubedSphereBilinear::print(std::ostream&) const { + ATLAS_NOTIMPLEMENTED; +} + +} // namespace method +} // namespace interpolation +} // namespace atlas diff --git a/src/atlas/interpolation/method/cubedsphere/CubedSphereBilinear.h b/src/atlas/interpolation/method/cubedsphere/CubedSphereBilinear.h new file mode 100644 index 000000000..3cf75c87f --- /dev/null +++ b/src/atlas/interpolation/method/cubedsphere/CubedSphereBilinear.h @@ -0,0 +1,55 @@ +/* + * (C) Crown Copyright 2022 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" + +namespace atlas { +namespace interpolation { +namespace method { + +/// @brief Cubed sphere bilinear interpolation method. +/// +/// @details Performs bilinear (quads) and baycentric (triangles) interpolation +/// accross cubed sphere tiles in (alpha, beta) coordinates. +/// Adding int "halo" (default 0) to the config controls the amount of +/// halo to consider when seraching for interpolation polygons. Adding +/// int "list_size" (defualt 4) will change the number of cells +/// returned by the internal kd-tree search. Increasing both values +/// may resolve errors if setup method fails to find cells. +/// The automatic halo exchange in the execute method can be disabled +/// by setting "halo_exchange" to false. +class CubedSphereBilinear : public Method { +public: + CubedSphereBilinear(const Config& config): Method(config) { + config.get("halo", halo_); + config.get("list_size", listSize_); + config.get("halo_exchange", halo_exchange_); + } + virtual ~CubedSphereBilinear() override {} + + virtual void print(std::ostream&) const override; + virtual const FunctionSpace& source() const override { return source_; } + virtual const FunctionSpace& target() const override { return target_; } + +private: + virtual void do_setup(const FunctionSpace& source, const FunctionSpace& target) override; + virtual void do_setup(const Grid& source, const Grid& target, const Cache&) override; + + FunctionSpace source_; + FunctionSpace target_; + + int halo_{0}; + int listSize_{4}; + bool halo_exchange_{true}; +}; + +} // namespace method +} // namespace interpolation +} // namespace atlas diff --git a/src/atlas/interpolation/method/knn/GridBoxAverage.cc b/src/atlas/interpolation/method/knn/GridBoxAverage.cc index 6997c190c..fce550c02 100644 --- a/src/atlas/interpolation/method/knn/GridBoxAverage.cc +++ b/src/atlas/interpolation/method/knn/GridBoxAverage.cc @@ -31,28 +31,28 @@ MethodBuilder __builder("grid-box-average"); } -void GridBoxAverage::do_execute(const FieldSet& source, FieldSet& target) const { +void GridBoxAverage::do_execute(const FieldSet& source, FieldSet& target, Metadata& metadata) const { ATLAS_ASSERT(source.size() == target.size()); // Matrix-based interpolation is handled by base (Method) class // TODO: exploit sparse/dense matrix multiplication for (idx_t i = 0; i < source.size(); ++i) { if (matrixFree_) { - GridBoxAverage::do_execute(source[i], target[i]); + GridBoxAverage::do_execute(source[i], target[i], metadata); } else { - Method::do_execute(source[i], target[i]); + Method::do_execute(source[i], target[i], metadata); } } } -void GridBoxAverage::do_execute(const Field& source, Field& target) const { +void GridBoxAverage::do_execute(const Field& source, Field& target, Metadata& metadata) const { ATLAS_TRACE("atlas::interpolation::method::GridBoxAverage::do_execute()"); // Matrix-based interpolation is handled by base (Method) class if (!matrixFree_) { - Method::do_execute(source, target); + Method::do_execute(source, target, metadata); return; } diff --git a/src/atlas/interpolation/method/knn/GridBoxAverage.h b/src/atlas/interpolation/method/knn/GridBoxAverage.h index 612863b16..1887712ec 100644 --- a/src/atlas/interpolation/method/knn/GridBoxAverage.h +++ b/src/atlas/interpolation/method/knn/GridBoxAverage.h @@ -24,8 +24,8 @@ class GridBoxAverage final : public GridBoxMethod { using GridBoxMethod::GridBoxMethod; private: - void do_execute(const FieldSet& source, FieldSet& target) const override; - void do_execute(const Field& source, Field& target) const override; + void do_execute(const FieldSet& source, FieldSet& target, Metadata&) const override; + void do_execute(const Field& source, Field& target, Metadata&) const override; }; diff --git a/src/atlas/interpolation/method/knn/GridBoxMaximum.cc b/src/atlas/interpolation/method/knn/GridBoxMaximum.cc index 4c3444b66..8ca86672a 100644 --- a/src/atlas/interpolation/method/knn/GridBoxMaximum.cc +++ b/src/atlas/interpolation/method/knn/GridBoxMaximum.cc @@ -34,16 +34,16 @@ MethodBuilder __builder("grid-box-maximum"); } -void GridBoxMaximum::do_execute(const FieldSet& source, FieldSet& target) const { +void GridBoxMaximum::do_execute(const FieldSet& source, FieldSet& target, Metadata& metadata) const { ATLAS_ASSERT(source.size() == target.size()); for (idx_t i = 0; i < source.size(); ++i) { - do_execute(source[i], target[i]); + do_execute(source[i], target[i], metadata); } } -void GridBoxMaximum::do_execute(const Field& source, Field& target) const { +void GridBoxMaximum::do_execute(const Field& source, Field& target, Metadata&) const { ATLAS_TRACE("atlas::interpolation::method::GridBoxMaximum::do_execute()"); // set arrays @@ -60,7 +60,7 @@ void GridBoxMaximum::do_execute(const Field& source, Field& target) const { if (!matrixFree_) { - const Matrix& m = *matrix_; + const Matrix& m = matrix(); Matrix::const_iterator k(m); for (decltype(m.rows()) i = 0, j = 0; i < m.rows(); ++i) { diff --git a/src/atlas/interpolation/method/knn/GridBoxMaximum.h b/src/atlas/interpolation/method/knn/GridBoxMaximum.h index 220fbb2c7..ca907e9c2 100644 --- a/src/atlas/interpolation/method/knn/GridBoxMaximum.h +++ b/src/atlas/interpolation/method/knn/GridBoxMaximum.h @@ -24,8 +24,8 @@ class GridBoxMaximum final : public GridBoxMethod { using GridBoxMethod::GridBoxMethod; private: - void do_execute(const FieldSet& source, FieldSet& target) const override; - void do_execute(const Field& source, Field& target) const override; + void do_execute(const FieldSet& source, FieldSet& target, Metadata&) const override; + void do_execute(const Field& source, Field& target, Metadata&) const override; }; diff --git a/src/atlas/interpolation/method/knn/GridBoxMethod.cc b/src/atlas/interpolation/method/knn/GridBoxMethod.cc index ebd1fe4eb..5ac390f47 100644 --- a/src/atlas/interpolation/method/knn/GridBoxMethod.cc +++ b/src/atlas/interpolation/method/knn/GridBoxMethod.cc @@ -107,10 +107,9 @@ void GridBoxMethod::do_setup(const Grid& source, const Grid& target, const Cache target_ = tgt; if (not matrixFree_ && interpolation::MatrixCache(cache)) { - matrix_cache_ = cache; - matrix_ = &matrix_cache_.matrix(); - ATLAS_ASSERT(matrix_cache_.matrix().rows() == target.size()); - ATLAS_ASSERT(matrix_cache_.matrix().cols() == source.size()); + setMatrix(cache); + ATLAS_ASSERT(matrix().rows() == target.size()); + ATLAS_ASSERT(matrix().cols() == source.size()); return; } @@ -125,8 +124,6 @@ void GridBoxMethod::do_setup(const Grid& source, const Grid& target, const Cache failures_.clear(); if (matrixFree_) { - Matrix A; - matrix_shared_->swap(A); return; } @@ -158,7 +155,7 @@ void GridBoxMethod::do_setup(const Grid& source, const Grid& target, const Cache { ATLAS_TRACE("GridBoxMethod::setup: build interpolant matrix"); Matrix A(targetBoxes_.size(), sourceBoxes_.size(), allTriplets); - matrix_shared_->swap(A); + setMatrix(A); } } @@ -184,7 +181,7 @@ void GridBoxMethod::giveUp(const std::forward_list& failures) { Cache GridBoxMethod::createCache() const { Cache cache; cache.add(interpolation::IndexKDTreeCache(pTree_)); - if (not matrix_->empty()) { + if (not matrix().empty()) { cache.add(Method::createCache()); } return cache; diff --git a/src/atlas/interpolation/method/knn/GridBoxMethod.h b/src/atlas/interpolation/method/knn/GridBoxMethod.h index 5658a86dd..170c80b7e 100644 --- a/src/atlas/interpolation/method/knn/GridBoxMethod.h +++ b/src/atlas/interpolation/method/knn/GridBoxMethod.h @@ -46,8 +46,8 @@ class GridBoxMethod : public KNearestNeighboursBase { bool intersect(size_t i, const GridBox& iBox, const util::IndexKDTree::ValueList&, std::vector&) const; - virtual void do_execute(const FieldSet& source, FieldSet& target) const override = 0; - virtual void do_execute(const Field& source, Field& target) const override = 0; + virtual void do_execute(const FieldSet& source, FieldSet& target, Metadata&) const override = 0; + virtual void do_execute(const Field& source, Field& target, Metadata&) const override = 0; virtual Cache createCache() const override; diff --git a/src/atlas/interpolation/method/knn/KNearestNeighbours.cc b/src/atlas/interpolation/method/knn/KNearestNeighbours.cc index c50ef6aa3..d39ffbf88 100644 --- a/src/atlas/interpolation/method/knn/KNearestNeighbours.cc +++ b/src/atlas/interpolation/method/knn/KNearestNeighbours.cc @@ -130,7 +130,7 @@ void KNearestNeighbours::do_setup(const FunctionSpace& source, const FunctionSpa // fill sparse matrix and return Matrix A(out_npts, inp_npts, weights_triplets); - matrix_shared_->swap(A); + setMatrix(A); } } // namespace method diff --git a/src/atlas/interpolation/method/knn/NearestNeighbour.cc b/src/atlas/interpolation/method/knn/NearestNeighbour.cc index 2ca8b24fb..fb500ba83 100644 --- a/src/atlas/interpolation/method/knn/NearestNeighbour.cc +++ b/src/atlas/interpolation/method/knn/NearestNeighbour.cc @@ -102,7 +102,7 @@ void NearestNeighbour::do_setup(const FunctionSpace& source, const FunctionSpace // fill sparse matrix and return Matrix A(out_npts, inp_npts, weights_triplets); - matrix_shared_->swap(A); + setMatrix(A); } } // namespace method diff --git a/src/atlas/interpolation/method/structured/StructuredInterpolation2D.h b/src/atlas/interpolation/method/structured/StructuredInterpolation2D.h index 1918145c1..9cd9ffb0c 100644 --- a/src/atlas/interpolation/method/structured/StructuredInterpolation2D.h +++ b/src/atlas/interpolation/method/structured/StructuredInterpolation2D.h @@ -56,9 +56,9 @@ class StructuredInterpolation2D : public Method { virtual void do_setup(const FunctionSpace& source, const FieldSet& target) override; - virtual void do_execute(const Field& src, Field& tgt) const override; + virtual void do_execute(const Field& src, Field& tgt, Metadata&) const override; - virtual void do_execute(const FieldSet& src, FieldSet& tgt) const override; + virtual void do_execute(const FieldSet& src, FieldSet& tgt, Metadata&) const override; template void execute_impl(const Kernel& kernel, const FieldSet& src, FieldSet& tgt) const; diff --git a/src/atlas/interpolation/method/structured/StructuredInterpolation2D.tcc b/src/atlas/interpolation/method/structured/StructuredInterpolation2D.tcc index 8b2c52774..2ee6cc2c0 100644 --- a/src/atlas/interpolation/method/structured/StructuredInterpolation2D.tcc +++ b/src/atlas/interpolation/method/structured/StructuredInterpolation2D.tcc @@ -189,23 +189,23 @@ void StructuredInterpolation2D::setup( const FunctionSpace& source ) { } // fill sparse matrix and return Matrix A( out_npts, inp_npts, triplets ); - matrix_shared_->swap( A ); + setMatrix(A); } } } template -void StructuredInterpolation2D::do_execute( const Field& src_field, Field& tgt_field ) const { +void StructuredInterpolation2D::do_execute( const Field& src_field, Field& tgt_field, Metadata& metadata ) const { FieldSet tgt( tgt_field ); - do_execute( FieldSet( src_field ), tgt ); + do_execute( FieldSet( src_field ), tgt, metadata ); } template -void StructuredInterpolation2D::do_execute( const FieldSet& src_fields, FieldSet& tgt_fields ) const { +void StructuredInterpolation2D::do_execute( const FieldSet& src_fields, FieldSet& tgt_fields, Metadata& metadata ) const { if ( not matrix_free_ ) { - Method::do_execute( src_fields, tgt_fields ); + Method::do_execute( src_fields, tgt_fields, metadata ); return; } diff --git a/src/atlas/interpolation/method/structured/StructuredInterpolation3D.h b/src/atlas/interpolation/method/structured/StructuredInterpolation3D.h index 1e63c774c..60347aac1 100644 --- a/src/atlas/interpolation/method/structured/StructuredInterpolation3D.h +++ b/src/atlas/interpolation/method/structured/StructuredInterpolation3D.h @@ -54,9 +54,9 @@ class StructuredInterpolation3D : public Method { virtual void do_setup(const FunctionSpace& source, const FieldSet& target) override; - virtual void do_execute(const Field& src, Field& tgt) const override; + virtual void do_execute(const Field& src, Field& tgt, Metadata&) const override; - virtual void do_execute(const FieldSet& src, FieldSet& tgt) const override; + virtual void do_execute(const FieldSet& src, FieldSet& tgt, Metadata&) const override; template void execute_impl(const Kernel& kernel, const FieldSet& src, FieldSet& tgt) const; diff --git a/src/atlas/interpolation/method/structured/StructuredInterpolation3D.tcc b/src/atlas/interpolation/method/structured/StructuredInterpolation3D.tcc index 911ec3ff4..359749e1a 100644 --- a/src/atlas/interpolation/method/structured/StructuredInterpolation3D.tcc +++ b/src/atlas/interpolation/method/structured/StructuredInterpolation3D.tcc @@ -136,16 +136,16 @@ void StructuredInterpolation3D::setup( const FunctionSpace& source ) { template -void StructuredInterpolation3D::do_execute( const Field& src_field, Field& tgt_field ) const { +void StructuredInterpolation3D::do_execute( const Field& src_field, Field& tgt_field, Metadata& metadata ) const { FieldSet tgt( tgt_field ); - do_execute( FieldSet( src_field ), tgt ); + do_execute( FieldSet( src_field ), tgt, metadata ); } template -void StructuredInterpolation3D::do_execute( const FieldSet& src_fields, FieldSet& tgt_fields ) const { +void StructuredInterpolation3D::do_execute( const FieldSet& src_fields, FieldSet& tgt_fields, Metadata& metadata ) const { if ( not matrix_free_ ) { - Method::do_execute( src_fields, tgt_fields ); + Method::do_execute( src_fields, tgt_fields, metadata ); return; } diff --git a/src/atlas/interpolation/method/fe/FiniteElement.cc b/src/atlas/interpolation/method/unstructured/FiniteElement.cc similarity index 98% rename from src/atlas/interpolation/method/fe/FiniteElement.cc rename to src/atlas/interpolation/method/unstructured/FiniteElement.cc index aaf9fbe68..2ec6ad994 100644 --- a/src/atlas/interpolation/method/fe/FiniteElement.cc +++ b/src/atlas/interpolation/method/unstructured/FiniteElement.cc @@ -60,10 +60,9 @@ void FiniteElement::do_setup(const Grid& source, const Grid& target, const Cache // no halo_exchange because we don't have any halo with delaunay or 3d structured meshgenerator if (interpolation::MatrixCache(cache)) { - matrix_cache_ = cache; - matrix_ = &matrix_cache_.matrix(); - ATLAS_ASSERT(matrix_cache_.matrix().rows() == target.size()); - ATLAS_ASSERT(matrix_cache_.matrix().cols() == source.size()); + setMatrix(cache); + ATLAS_ASSERT(matrix().rows() == target.size()); + ATLAS_ASSERT(matrix().cols() == source.size()); return; } if (mpi::size() > 1) { @@ -142,7 +141,7 @@ void FiniteElement::print(std::ostream& out) const { out << ", NodeColumns to NodeColumns stencil weights: " << std::endl; auto gidx_src = array::make_view(src.nodes().global_index()); - ATLAS_ASSERT(tgt.nodes().size() == idx_t(matrix_->rows())); + ATLAS_ASSERT(tgt.nodes().size() == idx_t(matrix().rows())); auto field_stencil_points_loc = tgt.createField(option::variables(Stencil::max_stencil_size)); @@ -154,7 +153,7 @@ void FiniteElement::print(std::ostream& out) const { auto stencil_size_loc = array::make_view(field_stencil_size_loc); stencil_size_loc.assign(0); - for (Matrix::const_iterator it = matrix_->begin(); it != matrix_->end(); ++it) { + for (auto it = matrix().begin(); it != matrix().end(); ++it) { idx_t p = idx_t(it.row()); idx_t& i = stencil_size_loc(p); stencil_points_loc(p, i) = gidx_src(it.col()); @@ -162,7 +161,6 @@ void FiniteElement::print(std::ostream& out) const { ++i; } - gidx_t global_size = tgt.gather().glb_dof(); auto field_stencil_points_glb = @@ -309,7 +307,7 @@ void FiniteElement::setup(const FunctionSpace& source) { // fill sparse matrix and return Matrix A(out_npts, inp_npts, weights_triplets); - matrix_shared_->swap(A); + setMatrix(A); } struct ElementEdge { diff --git a/src/atlas/interpolation/method/fe/FiniteElement.h b/src/atlas/interpolation/method/unstructured/FiniteElement.h similarity index 100% rename from src/atlas/interpolation/method/fe/FiniteElement.h rename to src/atlas/interpolation/method/unstructured/FiniteElement.h diff --git a/src/atlas/interpolation/method/bil/BilinearRemapping.cc b/src/atlas/interpolation/method/unstructured/UnstructuredBilinearLonLat.cc similarity index 64% rename from src/atlas/interpolation/method/bil/BilinearRemapping.cc rename to src/atlas/interpolation/method/unstructured/UnstructuredBilinearLonLat.cc index 571f1de8d..d84bad6c2 100644 --- a/src/atlas/interpolation/method/bil/BilinearRemapping.cc +++ b/src/atlas/interpolation/method/unstructured/UnstructuredBilinearLonLat.cc @@ -9,7 +9,7 @@ #include #include -#include "BilinearRemapping.h" +#include "UnstructuredBilinearLonLat.h" #include "eckit/log/Plural.h" #include "eckit/log/ProgressTimer.h" @@ -24,7 +24,8 @@ #include "atlas/interpolation/method/Ray.h" #include "atlas/mesh/ElementType.h" #include "atlas/mesh/Nodes.h" -#include "atlas/mesh/actions/Build2DCellCentres.h" +#include "atlas/mesh/actions/BuildCellCentres.h" +#include "atlas/mesh/actions/BuildXYZField.h" #include "atlas/meshgenerator.h" #include "atlas/parallel/GatherScatter.h" #include "atlas/parallel/mpi/Buffer.h" @@ -43,7 +44,7 @@ namespace method { namespace { -MethodBuilder __builder("bilinear-remapping"); +MethodBuilder __builder_2("unstructured-bilinear-lonlat"); // epsilon used to scale edge tolerance when projecting ray to intesect element static const double parametricEpsilon = 1e-15; @@ -51,15 +52,14 @@ static const double parametricEpsilon = 1e-15; } // namespace -void BilinearRemapping::do_setup(const Grid& source, const Grid& target, const Cache& cache) { +void UnstructuredBilinearLonLat::do_setup(const Grid& source, const Grid& target, const Cache& cache) { allow_halo_exchange_ = false; // no halo_exchange because we don't have any halo with delaunay or 3d structured meshgenerator if (interpolation::MatrixCache(cache)) { - matrix_cache_ = cache; - matrix_ = &matrix_cache_.matrix(); - ATLAS_ASSERT(matrix_cache_.matrix().rows() == target.size()); - ATLAS_ASSERT(matrix_cache_.matrix().cols() == source.size()); + setMatrix(cache); + ATLAS_ASSERT(matrix().rows() == target.size()); + ATLAS_ASSERT(matrix().cols() == source.size()); return; } if (mpi::size() > 1) { @@ -79,7 +79,7 @@ void BilinearRemapping::do_setup(const Grid& source, const Grid& target, const C do_setup(make_nodecolumns(source), functionspace::PointCloud{target}); } -void BilinearRemapping::do_setup(const FunctionSpace& source, const FunctionSpace& target) { +void UnstructuredBilinearLonLat::do_setup(const FunctionSpace& source, const FunctionSpace& target) { ATLAS_TRACE("atlas::interpolation::method::BilinearRemapping::do_setup()"); source_ = source; @@ -89,13 +89,25 @@ void BilinearRemapping::do_setup(const FunctionSpace& source, const FunctionSpac if (functionspace::NodeColumns tgt = target) { Mesh meshTarget = tgt.mesh(); + target_xyz_ = mesh::actions::BuildXYZField("xyz")(meshTarget); target_ghost_ = meshTarget.nodes().ghost(); target_lonlat_ = meshTarget.nodes().lonlat(); } else if (functionspace::PointCloud tgt = target) { const idx_t N = tgt.size(); + target_xyz_ = Field("xyz", array::make_datatype(), array::make_shape(N, 3)); target_ghost_ = tgt.ghost(); target_lonlat_ = tgt.lonlat(); + auto lonlat = array::make_view(tgt.lonlat()); + auto xyz = array::make_view(target_xyz_); + PointXYZ p2; + for (idx_t n = 0; n < N; ++n) { + const PointLonLat p1(lonlat(n, 0), lonlat(n, 1)); + util::Earth::convertSphericalToCartesian(p1, p2); + xyz(n, 0) = p2.x(); + xyz(n, 1) = p2.y(); + xyz(n, 2) = p2.z(); + } } else { ATLAS_NOTIMPLEMENTED; @@ -112,7 +124,7 @@ struct Stencil { }; }; -void BilinearRemapping::print(std::ostream& out) const { +void UnstructuredBilinearLonLat::print(std::ostream& out) const { functionspace::NodeColumns src(source_); functionspace::NodeColumns tgt(target_); out << "atlas::interpolation::method::BilinearRemapping{" << std::endl; @@ -125,7 +137,7 @@ void BilinearRemapping::print(std::ostream& out) const { out << ", NodeColumns to NodeColumns stencil weights: " << std::endl; auto gidx_src = array::make_view(src.nodes().global_index()); - ATLAS_ASSERT(tgt.nodes().size() == idx_t(matrix_->rows())); + ATLAS_ASSERT(tgt.nodes().size() == idx_t(matrix().rows())); auto field_stencil_points_loc = tgt.createField(option::variables(Stencil::max_stencil_size)); @@ -137,7 +149,7 @@ void BilinearRemapping::print(std::ostream& out) const { auto stencil_size_loc = array::make_view(field_stencil_size_loc); stencil_size_loc.assign(0); - for (Matrix::const_iterator it = matrix_->begin(); it != matrix_->end(); ++it) { + for (auto it = matrix().begin(); it != matrix().end(); ++it) { idx_t p = idx_t(it.row()); idx_t& i = stencil_size_loc(p); stencil_points_loc(p, i) = gidx_src(it.col()); @@ -180,7 +192,7 @@ void BilinearRemapping::print(std::ostream& out) const { out << "}" << std::endl; } -void BilinearRemapping::setup(const FunctionSpace& source) { +void UnstructuredBilinearLonLat::setup(const FunctionSpace& source) { const functionspace::NodeColumns src = source; ATLAS_ASSERT(src); @@ -189,21 +201,22 @@ void BilinearRemapping::setup(const FunctionSpace& source) { auto trace_setup_source = atlas::Trace{Here(), "Setup source"}; - // 2D point coordinates - auto source_lonlat = array::make_view(src.lonlat()); + // 3D point coordinates + Field source_xyz = mesh::actions::BuildXYZField("xyz")(meshSource); // generate barycenters of each triangle & insert them on a kd-tree util::Config config; config.set("name", "centre "); config.set("flatten_virtual_elements", false); - Field cell_centres = mesh::actions::Build2DCellCentres(config)(meshSource); + Field cell_centres = mesh::actions::BuildCellCentres(config)(meshSource); - std::unique_ptr eTree(create_element2D_kdtree(meshSource, cell_centres)); + std::unique_ptr eTree(create_element_kdtree(meshSource, cell_centres)); trace_setup_source.stop(); ilonlat_.reset(new array::ArrayView(array::make_view(meshSource.nodes().lonlat()))); olonlat_.reset(new array::ArrayView(array::make_view(target_lonlat_))); + oxyz_.reset(new array::ArrayView(array::make_view(target_xyz_))); igidx_.reset(new array::ArrayView(array::make_view(src.nodes().global_index()))); connectivity_ = &meshSource.cells().node_connectivity(); const mesh::Nodes& i_nodes = meshSource.nodes(); @@ -224,9 +237,8 @@ void BilinearRemapping::setup(const FunctionSpace& source) { // search nearest k cell centres const idx_t maxNbElemsToTry = std::max(8, idx_t(Nelements * max_fraction_elems_to_try_)); - idx_t max_neighbours = 0; - std::vector failures; + std::vector failures; ATLAS_TRACE_SCOPE("Computing interpolation matrix") { eckit::ProgressTimer progress("Computing interpolation weights", out_npts, "point", double(5), Log::debug()); @@ -235,30 +247,18 @@ void BilinearRemapping::setup(const FunctionSpace& source) { continue; } - double lon{(*olonlat_)(ip, 0)}; - while (lon > 360.0) { - lon -= 360.0; - } - while (lon < 0.0) { - lon += 360.0; - } - PointXY p{lon, (*olonlat_)(ip, 1)}; // lookup point + PointXYZ p{(*oxyz_)(ip, XX), (*oxyz_)(ip, YY), (*oxyz_)(ip, ZZ)}; // lookup point - idx_t kpts = 1; + idx_t kpts = 8; bool success = false; std::ostringstream failures_log; - while (!success && kpts <= maxNbElemsToTry) { - max_neighbours = std::max(kpts, max_neighbours); + ElemIndex3::NodeList cs = eTree->kNearestNeighbours(p, kpts); + Triplets triplets = projectPointToElements(ip, cs, failures_log); - ElemIndex2::NodeList cs = eTree->kNearestNeighbours(p, kpts); - Triplets triplets = projectPointToElements(ip, cs, failures_log); - - if (triplets.size()) { - std::copy(triplets.begin(), triplets.end(), std::back_inserter(weights_triplets)); - success = true; - } - kpts *= 2; + if (triplets.size()) { + std::copy(triplets.begin(), triplets.end(), std::back_inserter(weights_triplets)); + success = true; } if (!success) { @@ -266,13 +266,13 @@ void BilinearRemapping::setup(const FunctionSpace& source) { if (not treat_failure_as_missing_value_) { Log::debug() << "------------------------------------------------------" "---------------------\n"; - Log::debug() << "Failed to project point (lon,lat)=" << p << '\n'; + Log::debug() << "Failed to project point (lon,lat)=" << (*olonlat_)(ip, LON) << " " + << (*olonlat_)(ip, LAT) << '\n'; Log::debug() << failures_log.str(); } } } } - Log::debug() << "Maximum neighbours searched was " << eckit::Plural(max_neighbours, "element") << std::endl; if (failures.size()) { if (treat_failure_as_missing_value_) { @@ -283,8 +283,8 @@ void BilinearRemapping::setup(const FunctionSpace& source) { // If this fails, consider lowering atlas::grid::parametricEpsilon std::ostringstream msg; msg << "Rank " << mpi::rank() << " failed to project points:\n"; - for (std::vector::const_iterator i = failures.begin(); i != failures.end(); ++i) { - const PointLonLat pll{(*olonlat_)(*i, (size_t)0), (*olonlat_)(*i, (size_t)1)}; // lookup point + for (auto i : failures) { + const PointLonLat pll{(*olonlat_)(i, LON), (*olonlat_)(i, LAT)}; // lookup point msg << "\t(lon,lat) = " << pll << "\n"; } @@ -295,30 +295,48 @@ void BilinearRemapping::setup(const FunctionSpace& source) { // fill sparse matrix and return Matrix A(out_npts, inp_npts, weights_triplets); - matrix_shared_->swap(A); + setMatrix(A); } -Method::Triplets BilinearRemapping::projectPointToElements(size_t ip, const ElemIndex2::NodeList& elems, - std::ostream& /* failures_log */) const { +Method::Triplets UnstructuredBilinearLonLat::projectPointToElements(size_t ip, const ElemIndex3::NodeList& elems, + std::ostream& /* failures_log */) const { ATLAS_ASSERT(elems.begin() != elems.end()); const size_t inp_points = ilonlat_->shape(0); std::array idx; std::array w; + std::array inv_dist_w; Triplets triplets; triplets.reserve(4); - double lon{(*olonlat_)(ip, 0)}; - while (lon > 360.0) { - lon -= 360.0; + double o_lon{(*olonlat_)(ip, LON)}; + while (o_lon >= 360.0) { + o_lon -= 360.0; } - while (lon < 0.0) { - lon += 360.0; + while (o_lon < 0.0) { + o_lon += 360.0; } - PointXY ob_loc{lon, (*olonlat_)(ip, 1)}; // lookup point + PointLonLat o_loc{o_lon, (*olonlat_)(ip, LAT)}; // lookup point + + auto inv_dist_weight = [](element::Quad2D& q, const PointXY& loc, std::array& w) { + double d[4]; + d[0] = util::Earth::distance({q.p(0).data()}, loc); + d[1] = util::Earth::distance({q.p(1).data()}, loc); + d[2] = util::Earth::distance({q.p(2).data()}, loc); + d[3] = util::Earth::distance({q.p(3).data()}, loc); + w[0] = d[1] * d[2] * d[3]; + w[1] = d[0] * d[2] * d[3]; + w[2] = d[1] * d[0] * d[3]; + w[3] = d[1] * d[0] * d[2]; + + double suminv = 1. / (w[0] + w[1] + w[2] + w[3]); + for (size_t i = 0; i < 4; ++i) { + w[i] *= suminv; + } + }; - for (ElemIndex2::NodeList::const_iterator itc = elems.begin(); itc != elems.end(); ++itc) { + for (ElemIndex3::NodeList::const_iterator itc = elems.begin(); itc != elems.end(); ++itc) { const idx_t elem_id = idx_t((*itc).value().payload()); ATLAS_ASSERT(elem_id < connectivity_->rows()); @@ -332,16 +350,16 @@ Method::Triplets BilinearRemapping::projectPointToElements(size_t ip, const Elem if (nb_cols == 3) { /* triangle */ - element::Triag2D triag(PointXY{(*ilonlat_)(idx[0], 0), (*ilonlat_)(idx[0], 1)}, - PointXY{(*ilonlat_)(idx[1], 0), (*ilonlat_)(idx[1], 1)}, - PointXY{(*ilonlat_)(idx[2], 0), (*ilonlat_)(idx[2], 1)}); + element::Triag2D triag(PointLonLat{(*ilonlat_)(idx[0], LON), (*ilonlat_)(idx[0], LAT)}, + PointLonLat{(*ilonlat_)(idx[1], LON), (*ilonlat_)(idx[1], LAT)}, + PointLonLat{(*ilonlat_)(idx[2], LON), (*ilonlat_)(idx[2], LAT)}); // pick an epsilon based on a characteristic length (sqrt(area)) // (this scales linearly so it better compares with linear weights u,v,w) const double edgeEpsilon = parametricEpsilon * std::sqrt(triag.area()); ATLAS_ASSERT(edgeEpsilon >= 0); - Intersect is = triag.intersects(ob_loc, edgeEpsilon); + Intersect is = triag.intersects(o_loc, edgeEpsilon); if (is) { // weights are the linear Lagrange function evaluated at u,v (aka @@ -359,36 +377,72 @@ Method::Triplets BilinearRemapping::projectPointToElements(size_t ip, const Elem } else { /* quadrilateral */ - element::Quad2D quad(PointXY{(*ilonlat_)(idx[0], 0), (*ilonlat_)(idx[0], 1)}, - PointXY{(*ilonlat_)(idx[1], 0), (*ilonlat_)(idx[1], 1)}, - PointXY{(*ilonlat_)(idx[2], 0), (*ilonlat_)(idx[2], 1)}, - PointXY{(*ilonlat_)(idx[3], 0), (*ilonlat_)(idx[3], 1)}); + double lons[4]; + lons[0] = (*ilonlat_)(idx[0], LON); + lons[1] = (*ilonlat_)(idx[1], LON); + lons[2] = (*ilonlat_)(idx[2], LON); + lons[3] = (*ilonlat_)(idx[3], LON); + // adjust quadrilaterals to lie within [0,360] + for (idx_t i = 0; i < 4; ++i) { + while (lons[i] > 360.0) { + lons[i] -= 360.0; + } + while (lons[i] < 0.0) { + lons[i] += 360.0; + } + } + // shift cells on the east-west periodic boundary from the east to the west + // so that the quad surrounds a point with output longitude in [0,360) + if (lons[0] > lons[1] || lons[3] > lons[2]) { + lons[0] -= 360; + lons[3] -= 360; + } + + element::Quad2D quad( + PointLonLat{lons[0], (*ilonlat_)(idx[0], LAT)}, PointLonLat{lons[1], (*ilonlat_)(idx[1], LAT)}, + PointLonLat{lons[2], (*ilonlat_)(idx[2], LAT)}, PointLonLat{lons[3], (*ilonlat_)(idx[3], LAT)}); + + if (itc == elems.begin()) { + inv_dist_weight(quad, o_loc, inv_dist_w); + } // pick an epsilon based on a characteristic length (sqrt(area)) // (this scales linearly so it better compares with linear weights u,v,w) const double edgeEpsilon = parametricEpsilon * std::sqrt(quad.area()); ATLAS_ASSERT(edgeEpsilon >= 0); - Intersect is = quad.localRemap(ob_loc, edgeEpsilon); + Intersect is = quad.localRemap(o_loc, edgeEpsilon); + + if (!is) { + // repeat the calculation to catch points which are near the east boundary of the grid. + is = quad.localRemap({o_lon - 360, (*olonlat_)(ip, LAT)}, edgeEpsilon); + } if (is) { - //std::cout << "intersection found for p: "; - //std::cout << ob_loc << " and " << quad << std::endl; // weights are the bilinear Lagrange function evaluated at u,v w[0] = (1. - is.u) * (1. - is.v); w[1] = is.u * (1. - is.v); w[2] = is.u * is.v; w[3] = (1. - is.u) * is.v; - for (size_t i = 0; i < 4; ++i) { triplets.emplace_back(ip, idx[i], w[i]); } break; // stop looking for elements } } - } // loop over nearest elements + if (triplets.empty()) { + // Crude inverse distance weighting to catch cells that are difficult + // to identify and interpolate in a lon/lat projection, near the north + // pole + const idx_t elem_id = idx_t((*elems.begin()).value().payload()); + for (size_t i = 0; i < 4; ++i) { + idx[i] = (*connectivity_)(elem_id, i); + triplets.emplace_back(ip, idx[i], inv_dist_w[i]); + } + } + if (!triplets.empty()) { normalise(triplets); } diff --git a/src/atlas/interpolation/method/bil/BilinearRemapping.h b/src/atlas/interpolation/method/unstructured/UnstructuredBilinearLonLat.h similarity index 86% rename from src/atlas/interpolation/method/bil/BilinearRemapping.h rename to src/atlas/interpolation/method/unstructured/UnstructuredBilinearLonLat.h index 5bb0baba3..53f022835 100644 --- a/src/atlas/interpolation/method/bil/BilinearRemapping.h +++ b/src/atlas/interpolation/method/unstructured/UnstructuredBilinearLonLat.h @@ -16,20 +16,20 @@ #include "atlas/array/ArrayView.h" #include "atlas/functionspace/FunctionSpace.h" -#include "atlas/interpolation/method/PointIndex2.h" +#include "atlas/interpolation/method/PointIndex3.h" #include "atlas/mesh/Elements.h" namespace atlas { namespace interpolation { namespace method { -class BilinearRemapping : public Method { +class UnstructuredBilinearLonLat : public Method { public: - BilinearRemapping(const Config& config): Method(config) { + UnstructuredBilinearLonLat(const Config& config): Method(config) { config.get("max_fraction_elems_to_try", max_fraction_elems_to_try_); } - virtual ~BilinearRemapping() override {} + virtual ~UnstructuredBilinearLonLat() override {} virtual void print(std::ostream&) const override; @@ -53,7 +53,7 @@ class BilinearRemapping : public Method { * point to the nearest element(s), returning the (normalized) interpolation * weights */ - Triplets projectPointToElements(size_t ip, const ElemIndex2::NodeList& elems, std::ostream& failures_log) const; + Triplets projectPointToElements(size_t ip, const ElemIndex3::NodeList& elems, std::ostream& failures_log) const; virtual const FunctionSpace& source() const override { return source_; } virtual const FunctionSpace& target() const override { return target_; } @@ -68,6 +68,7 @@ class BilinearRemapping : public Method { mesh::MultiBlockConnectivity* connectivity_; std::unique_ptr> ilonlat_; std::unique_ptr> olonlat_; + std::unique_ptr> oxyz_; std::unique_ptr> igidx_; Field target_lonlat_; diff --git a/src/atlas/interpolation/nonlinear/Missing.h b/src/atlas/interpolation/nonlinear/Missing.h index 78b6f4bdc..aa8eadceb 100644 --- a/src/atlas/interpolation/nonlinear/Missing.h +++ b/src/atlas/interpolation/nonlinear/Missing.h @@ -35,9 +35,10 @@ struct MissingIfAllMissing : Missing { field::MissingValue mv(field); auto& missingValue = mv.ref(); - // NOTE only for scalars (for now) + ATLAS_ASSERT(field.rank() == 1); + auto values = make_view_field_values(field); - ATLAS_ASSERT(idx_t(W.cols()) == values.size()); + ATLAS_ASSERT(idx_t(W.cols()) == values.shape(0)); auto data = const_cast(W.data()); bool modif = false; diff --git a/src/atlas/library/defines.h.in b/src/atlas/library/defines.h.in index 650f56fb3..3316a349e 100644 --- a/src/atlas/library/defines.h.in +++ b/src/atlas/library/defines.h.in @@ -15,6 +15,7 @@ #define atlas_library_defines_h #define ATLAS_HAVE_OMP @atlas_HAVE_OMP_CXX@ +#define ATLAS_OMP_TASK_SUPPORTED @ATLAS_OMP_TASK_SUPPORTED@ #define ATLAS_OMP_TASK_UNTIED_SUPPORTED @ATLAS_OMP_TASK_UNTIED_SUPPORTED@ #define ATLAS_HAVE_ACC @atlas_HAVE_ACC@ #define ATLAS_HAVE_TESSELATION @atlas_HAVE_TESSELATION@ diff --git a/src/atlas/meshgenerator/detail/CubedSphereDualMeshGenerator.cc b/src/atlas/meshgenerator/detail/CubedSphereDualMeshGenerator.cc index b6d3eea46..235b94e4e 100644 --- a/src/atlas/meshgenerator/detail/CubedSphereDualMeshGenerator.cc +++ b/src/atlas/meshgenerator/detail/CubedSphereDualMeshGenerator.cc @@ -123,14 +123,6 @@ void CubedSphereDualMeshGenerator::generate(const Grid& grid, const grid::Distri throw_Exception("CubedSphereDualMeshGenerator will only work with a cell-centroid grid.", Here()); } - // Enforce compatible halo size. - if (options.get("halo") != 0) { - throw_Exception( - "Halo size CubedSphereDualMeshGenerator is currently " - "limited to 0.", - Here()); - } - // Clone some grid properties. setGrid(mesh, csGrid, distribution); @@ -141,6 +133,20 @@ void CubedSphereDualMeshGenerator::generate(const Grid& grid, const grid::Distri namespace { +// (i, j) pair +class IJ { +public: + IJ(idx_t i, idx_t j): i_(i), j_(j) {} + idx_t i() const { return i_; } + idx_t j() const { return j_; } + IJ operator+(const IJ& ij) const { return IJ{i() + ij.i(), j() + ij.j()}; } + IJ operator-(const IJ& ij) const { return IJ{i() - ij.i(), j() - ij.j()}; } + +private: + idx_t i_{}; + idx_t j_{}; +}; + // Helper function to copy fields. template void copyField(const Field& sourceField, Field& targetField) { @@ -148,12 +154,109 @@ void copyField(const Field& sourceField, Field& targetField) { array::make_view(targetField).assign(array::make_view(sourceField)); } +// Get the surrounding node (i, j) pairs from a cell (i, j) pair. +std::vector getIjNodes(const IJ& ijCell, idx_t N) { + // Rotate ij 90 degrees anitclockwise about ijPivot. + auto rotateAnticlockwise = [&](const IJ& ij, const IJ& ijPivot) { + const auto ijTemp = ij - ijPivot; + return IJ{-ijTemp.j(), ijTemp.i()} + ijPivot; + }; + + // Rotate ij 90 degrees clockwise about ijPivot. + auto rotateClockwise = [&](const IJ& ij, const IJ& ijPivot) { + const auto ijTemp = ij - ijPivot; + return IJ{ijTemp.j(), -ijTemp.i()} + ijPivot; + }; + + // Set standard surrounding nodes. + auto ijNodes = std::vector{{ijCell.i() - 1, ijCell.j() - 1}, + {ijCell.i(), ijCell.j() - 1}, + {ijCell.i(), ijCell.j()}, + {ijCell.i() - 1, ijCell.j()}}; + + // Modify nodes that lie in invalid corners of ij space. + // Either remove a node to make cell triangular, or rotate two of the + // nodes out of the forbidden space. + + // Bottom-left corner. + if (ijCell.i() <= 0 && ijCell.j() <= 0) { + // Triangle. + if (ijCell.i() == 0 && ijCell.j() == 0) { + ijNodes.erase(ijNodes.begin()); + } + // Quad (i) + else if (ijCell.i() == 0) { + ijNodes[0] = rotateClockwise(ijNodes[1], IJ{0, 0}); + ijNodes[3] = rotateClockwise(ijNodes[2], IJ{0, 0}); + } + else { + // Quad (ii) + ijNodes[0] = rotateAnticlockwise(ijNodes[3], IJ{0, 0}); + ijNodes[1] = rotateAnticlockwise(ijNodes[2], IJ{0, 0}); + } + } + // Bottom-right corner. + else if (ijCell.i() >= N && ijCell.j() <= 0) { + // Triangle. + if (ijCell.i() == N && ijCell.j() == 0) { + ijNodes.erase(ijNodes.begin() + 1); + } + // Quad (i) + else if (ijCell.j() == 0) { + ijNodes[0] = rotateClockwise(ijNodes[3], IJ{N - 1, 0}); + ijNodes[1] = rotateClockwise(ijNodes[2], IJ{N - 1, 0}); + } + // Quad (ii) + else { + ijNodes[1] = rotateAnticlockwise(ijNodes[0], IJ{N - 1, 0}); + ijNodes[2] = rotateAnticlockwise(ijNodes[3], IJ{N - 1, 0}); + } + } + // Top-right corner. + else if (ijCell.i() >= N && ijCell.j() >= N) { + // Triangle. + if (ijCell.i() == N && ijCell.j() == N) { + ijNodes.erase(ijNodes.begin() + 2); + } + // Quad (i) + else if (ijCell.i() == N) { + ijNodes[1] = rotateClockwise(ijNodes[0], IJ{N - 1, N - 1}); + ijNodes[2] = rotateClockwise(ijNodes[3], IJ{N - 1, N - 1}); + } + // Quad (ii) + else { + ijNodes[2] = rotateAnticlockwise(ijNodes[1], IJ{N - 1, N - 1}); + ijNodes[3] = rotateAnticlockwise(ijNodes[0], IJ{N - 1, N - 1}); + } + } + // Top-left corner. + else if (ijCell.i() <= 0 && ijCell.j() >= N) { + // Triangle. + if (ijCell.i() == 0 && ijCell.j() == N) { + ijNodes.erase(ijNodes.begin() + 3); + } + // Quad (i) + else if (ijCell.j() == N) { + ijNodes[2] = rotateClockwise(ijNodes[1], IJ{0, N - 1}); + ijNodes[3] = rotateClockwise(ijNodes[0], IJ{0, N - 1}); + } + // Quad (ii) + else { + ijNodes[0] = rotateAnticlockwise(ijNodes[1], IJ{0, N - 1}); + ijNodes[3] = rotateAnticlockwise(ijNodes[2], IJ{0, N - 1}); + } + } + + return ijNodes; +} + } // namespace void CubedSphereDualMeshGenerator::generate_mesh(const CubedSphereGrid& csGrid, const grid::Distribution& distribution, Mesh& mesh) const { ATLAS_TRACE("CubedSphereDualMeshGenerator::generate"); + using Topology = atlas::mesh::Nodes::Topology; using namespace detail::cubedsphere; const idx_t N = csGrid.N(); @@ -163,36 +266,36 @@ void CubedSphereDualMeshGenerator::generate_mesh(const CubedSphereGrid& csGrid, // Create a cubed-sphere mesh. //-------------------------------------------------------------------------- - // Generate cubed sphere mesh. - auto csOptions = options; - csOptions.set("halo", nHalo + 1); - const auto csMesh = MeshGenerator("cubedsphere", csOptions).generate(csGrid, distribution); + // Generate cubed sphere primal mesh. + auto primalOptions = options; + primalOptions.set("halo", nHalo + 1); + const auto primalMesh = MeshGenerator("cubedsphere", primalOptions).generate(csGrid, distribution); - // Generate fucntionspaces cubed sphere mesh. - const auto csCellsFunctionSpace = functionspace::CubedSphereCellColumns(csMesh); - const auto csNodesFunctionSpace = functionspace::CubedSphereNodeColumns(csMesh); - const auto& csCells = csCellsFunctionSpace.cells(); - const auto& csNodes = csNodesFunctionSpace.nodes(); + // Generate fucntionspaces for cubed sphere primal mesh. + const auto primalCellsFunctionSpace = functionspace::CubedSphereCellColumns(primalMesh); + const auto primalNodesFunctionSpace = functionspace::CubedSphereNodeColumns(primalMesh); + const auto& primalCells = primalCellsFunctionSpace.cells(); + const auto& primalNodes = primalNodesFunctionSpace.nodes(); //-------------------------------------------------------------------------- // Set dual mesh nodes (easy part). //-------------------------------------------------------------------------- auto& nodes = mesh.nodes(); - nodes.resize(csCellsFunctionSpace.size()); + nodes.resize(primalCellsFunctionSpace.size()); nodes.add(Field("tij", array::make_datatype(), array::make_shape(nodes.size(), 3))); // Copy mesh fields to dual mesh. - copyField(csCells.global_index(), nodes.global_index()); - copyField(csCells.remote_index(), nodes.remote_index()); - copyField(csCells.partition(), nodes.partition()); - copyField(csCells.halo(), nodes.halo()); - copyField(csCells.flags(), nodes.flags()); - copyField(csCells.field("ghost"), nodes.ghost()); - copyField(csCells.field("xy"), nodes.xy()); - copyField(csCells.field("lonlat"), nodes.lonlat()); - copyField(csCells.field("tij"), nodes.field("tij")); + copyField(primalCells.global_index(), nodes.global_index()); + copyField(primalCells.remote_index(), nodes.remote_index()); + copyField(primalCells.partition(), nodes.partition()); + copyField(primalCells.halo(), nodes.halo()); + copyField(primalCells.flags(), nodes.flags()); + copyField(primalCells.field("ghost"), nodes.ghost()); + copyField(primalCells.field("xy"), nodes.xy()); + copyField(primalCells.field("lonlat"), nodes.lonlat()); + copyField(primalCells.field("tij"), nodes.field("tij")); // Need to decrement halo by one. auto nodesHalo = array::make_view(nodes.halo()); @@ -206,37 +309,65 @@ void CubedSphereDualMeshGenerator::generate_mesh(const CubedSphereGrid& csGrid, //-------------------------------------------------------------------------- // Make views to cubed sphere nodes. - const auto csNodesHalo = array::make_view(csNodes.halo()); - const auto csNodesTij = array::make_view(csNodes.field("tij")); + const auto primalNodesHalo = array::make_view(primalNodes.halo()); + const auto primalNodesTij = array::make_view(primalNodes.field("tij")); - // Figure out element types. - enum struct ElemType - { - tri, - quad - }; - const auto getType = [&](idx_t idx) -> ElemType { - // Nodes of csMesh are cells of dual mesh. - const idx_t i = csNodesTij(idx, Coordinates::I); - const idx_t j = csNodesTij(idx, Coordinates::J); - auto cornerCell = (i == 0 && j == 0) || (i == N && j == 0) || (i == N && j == N) || (i == 0 && j == N); - return cornerCell ? ElemType::tri : ElemType::quad; + // Loop over all nodes of primal mesh, excluding outermost halo. + // Find dual mesh nodes around primal mesh nodes. + // Note: some halo cells near corners may be incomplete. + + // Set of nodes around a cell. + struct NodeList { + std::vector nodes{}; // Node indices. + bool incomplete{}; // True if nodes are missing. }; - // Set first element type. ( first = type, second = count ) - auto typeCounts = std::vector >{std::make_pair(getType(0), 1)}; + auto nodeLists = std::vector{}; - // Count the number of consecutive triangles and quadtrilaterals in dual mesh. - // This is an attempt to keep dual mesh cells in the same order as mesh nodes. - // Otherwise, the halo exchange bookkeeping is invalidated. - for (idx_t idx = 1; idx < csNodes.size(); ++idx) { + for (idx_t idx = 0; idx < primalNodes.size(); ++idx) { // Exclude outer ring of cubed sphere mesh halo. - if (csNodesHalo(idx) == nHalo + 1) { + if (primalNodesHalo(idx) == nHalo + 1) { break; } + nodeLists.emplace_back(); + auto& nodeList = nodeLists.back(); + + // Get tij of cell. + const auto tCell = primalNodesTij(idx, Coordinates::T); + const auto ijCell = IJ{primalNodesTij(idx, Coordinates::I), primalNodesTij(idx, Coordinates::J)}; + + // Get ij of surrounding nodes. + auto ijNodes = getIjNodes(ijCell, N); + + // Add indices to nodes vector. + for (const auto& ijNode : ijNodes) { + if (primalCellsFunctionSpace.is_valid_index(tCell, ijNode.i(), ijNode.j())) { + nodeList.nodes.push_back(primalCellsFunctionSpace.index(tCell, ijNode.i(), ijNode.j())); + } + else { + nodeList.incomplete = true; + } + } + } + + // Figure out element types. + // Set first element type. ( first = type, second = count ) + enum struct ElemType : size_t + { + LINE = 2, + TRIANGLE = 3, + QUADRILATERAL = 4 + }; + auto typeCounts = + std::vector>{std::make_pair(static_cast(nodeLists[0].nodes.size()), 1)}; + + // Count the number of consecutive lines, triangles or quadtrilaterals in dual mesh. + // This is an attempt to keep dual mesh cells in the same order as mesh nodes. + // Otherwise, the halo exchange bookkeeping is invalidated. + for (size_t idx = 1; idx < nodeLists.size(); ++idx) { // Get the element type. - const auto elemType = getType(idx); + const auto elemType = static_cast(nodeLists[idx].nodes.size()); // Increment counter if this elemType is the same as last one if (elemType == typeCounts.back().first) { @@ -255,16 +386,26 @@ void CubedSphereDualMeshGenerator::generate_mesh(const CubedSphereGrid& csGrid, for (const auto& typeCount : typeCounts) { // Select element type. switch (typeCount.first) { - // Add a batch of triangles. - case ElemType::tri: { + // Add a block of lines. + case ElemType::LINE: { + cells.add(new mesh::temporary::Line(), typeCount.second); + break; + } + // Add a block of triangles. + case ElemType::TRIANGLE: { cells.add(new mesh::temporary::Triangle(), typeCount.second); break; } - // Add a batch quadrilaterals. - case ElemType::quad: { + // Add a block of quadrilaterals. + case ElemType::QUADRILATERAL: { cells.add(new mesh::temporary::Quadrilateral(), typeCount.second); break; } + default: { + ATLAS_THROW_EXCEPTION("Unknown element type with " + + std::to_string(static_cast(typeCount.first)) + " nodes."); + break; + } } // Increment the total number of cells. nCells += typeCount.second; @@ -272,85 +413,92 @@ void CubedSphereDualMeshGenerator::generate_mesh(const CubedSphereGrid& csGrid, // Add extra fields to cells. cells.add(Field("tij", array::make_datatype(), array::make_shape(cells.size(), 3))); - cells.add(Field("xy", array::make_datatype(), array::make_shape(cells.size(), 2))); - cells.add(Field("lonlat", array::make_datatype(), array::make_shape(cells.size(), 2))); - cells.add(Field("ghost", array::make_datatype(), array::make_shape(cells.size()))); // Copy dual cells fields from nodes. - copyField(csNodes.global_index(), cells.global_index()); - copyField(csNodes.remote_index(), cells.remote_index()); - copyField(csNodes.partition(), cells.partition()); - copyField(csNodes.ghost(), cells.field("ghost")); - copyField(csNodes.halo(), cells.halo()); - copyField(csNodes.flags(), cells.flags()); - copyField(csNodes.field("tij"), cells.field("tij")); - copyField(csNodes.xy(), cells.field("xy")); - copyField(csNodes.lonlat(), cells.field("lonlat")); + copyField(primalNodes.global_index(), cells.global_index()); + copyField(primalNodes.remote_index(), cells.remote_index()); + copyField(primalNodes.partition(), cells.partition()); + copyField(primalNodes.ghost(), cells.field("ghost")); + copyField(primalNodes.halo(), cells.halo()); + copyField(primalNodes.flags(), cells.flags()); + copyField(primalNodes.field("tij"), cells.field("tij")); + copyField(primalNodes.xy(), cells.field("xy")); + copyField(primalNodes.lonlat(), cells.field("lonlat")); + + // Get view of flags field. + auto dualCellsFlags = array::make_view(cells.flags()); // Get node connectivity. auto& nodeConnectivity = cells.node_connectivity(); // Loop over dual mesh cells and set connectivity. for (idx_t idx = 0; idx < nCells; ++idx) { - const idx_t t = csNodesTij(idx, Coordinates::T); - const idx_t i = csNodesTij(idx, Coordinates::I); - const idx_t j = csNodesTij(idx, Coordinates::J); - - - // Declare vector of surrounding nodes. - auto cellNodes = std::vector{}; - - // Get number of nodes per element. - const auto nNodes = cells.node_connectivity().row(idx).size(); - // Element is a triangle. - if (nNodes == 3) { - // Bottom-left corner. - if (i == 0 && j == 0) { - cellNodes = {csCellsFunctionSpace.index(t, 0, -1), csCellsFunctionSpace.index(t, 0, 0), - csCellsFunctionSpace.index(t, -1, 0)}; - } - // Bottom-right corner. - else if (i == N && j == 0) { - cellNodes = {csCellsFunctionSpace.index(t, N - 1, -1), csCellsFunctionSpace.index(t, N, 0), - csCellsFunctionSpace.index(t, N - 1, 0)}; - } - // Top-right corner. - else if (i == N && j == N) { - cellNodes = {csCellsFunctionSpace.index(t, N - 1, N - 1), csCellsFunctionSpace.index(t, N, N - 1), - csCellsFunctionSpace.index(t, N - 1, N)}; - } - // Top-left corner. - else { - cellNodes = {csCellsFunctionSpace.index(t, -1, N - 1), csCellsFunctionSpace.index(t, 0, N - 1), - csCellsFunctionSpace.index(t, 0, N)}; - } - } - // Element is a quadtrilateral. - else if (nNodes == 4) { - cellNodes = {csCellsFunctionSpace.index(t, i - 1, j - 1), csCellsFunctionSpace.index(t, i, j - 1), - csCellsFunctionSpace.index(t, i, j), csCellsFunctionSpace.index(t, i - 1, j)}; - } - // Couldn't determine element type. - else { - ATLAS_THROW_EXCEPTION("Could not determine element type for cell " + std::to_string(idx) + ".";); - } + // Set connectivity. + nodeConnectivity.set(idx, nodeLists[idx].nodes.data()); - nodeConnectivity.set(idx, cellNodes.data()); + // Set invalid flag if cell is incomplete. + if (nodeLists[idx].incomplete) { + Topology::set(dualCellsFlags(idx), Topology::INVALID); + } } // Set metadata. + set_metadata(mesh); + + return; +} + +// ----------------------------------------------------------------------------- + +void CubedSphereDualMeshGenerator::set_metadata(Mesh& mesh) const { + const auto nHalo = options.get("halo"); + + // Set basic halo metadata. mesh.metadata().set("halo", nHalo); mesh.metadata().set("halo_locked", true); mesh.nodes().metadata().set("parallel", true); mesh.cells().metadata().set("parallel", true); - return; -} + // Loop over nodes and count number of halo elements. + auto nNodes = std::vector(nHalo + 2, 0); + const auto nodeHalo = array::make_view(mesh.nodes().halo()); + for (idx_t i = 0; i < mesh.nodes().size(); ++i) { + ++nNodes[static_cast(nodeHalo(i))]; + } + std::partial_sum(nNodes.begin(), nNodes.end(), nNodes.begin()); + + // Set node halo metadata. + for (size_t i = 0; i < nNodes.size(); ++i) { + const auto str = "nb_nodes_including_halo[" + std::to_string(i) + "]"; + mesh.metadata().set(str, nNodes[i]); + } + // Loop over cells and count number of halo elements. + auto nCells = std::vector>(mesh.cells().nb_types(), std::vector(nHalo + 1, 0)); + const auto cellHalo = array::make_view(mesh.cells().halo()); + + for (idx_t i = 0; i < mesh.cells().nb_types(); ++i) { + const auto& elems = mesh.cells().elements(i); + for (idx_t j = elems.begin(); j < elems.end(); ++j) { + ++nCells[static_cast(i)][static_cast(cellHalo(j))]; + } + std::partial_sum(nCells[static_cast(i)].begin(), nCells[static_cast(i)].end(), + nCells[static_cast(i)].begin()); + } + + // Set cell halo metadata. + for (size_t i = 0; i < nCells.size(); ++i) { + for (size_t j = 0; j < nCells[i].size(); ++j) { + const auto str = "nb_cells_including_halo[" + std::to_string(i) + "][" + std::to_string(j) + "]"; + mesh.metadata().set(str, nCells[i][j]); + } + } +} + // ----------------------------------------------------------------------------- void CubedSphereDualMeshGenerator::hash(eckit::Hash& h) const { diff --git a/src/atlas/meshgenerator/detail/CubedSphereDualMeshGenerator.h b/src/atlas/meshgenerator/detail/CubedSphereDualMeshGenerator.h index 8a1f6cda9..9a0369786 100644 --- a/src/atlas/meshgenerator/detail/CubedSphereDualMeshGenerator.h +++ b/src/atlas/meshgenerator/detail/CubedSphereDualMeshGenerator.h @@ -56,6 +56,8 @@ class CubedSphereDualMeshGenerator : public MeshGenerator::Implementation { void generate_mesh(const CubedSphereGrid&, const grid::Distribution&, Mesh&) const; + void set_metadata(Mesh&) const; + private: util::Metadata options; }; diff --git a/src/atlas/meshgenerator/detail/CubedSphereMeshGenerator.cc b/src/atlas/meshgenerator/detail/CubedSphereMeshGenerator.cc index 93430896b..bae86c1bd 100644 --- a/src/atlas/meshgenerator/detail/CubedSphereMeshGenerator.cc +++ b/src/atlas/meshgenerator/detail/CubedSphereMeshGenerator.cc @@ -84,7 +84,7 @@ void CubedSphereMeshGenerator::configure_defaults() { options.set("part", mpi::rank()); // This options sets the number of halo elements around each partition. - options.set("halo", 1); + options.set("halo", 0); // This options sets the default partitioner. options.set("partitioner", "cubedsphere"); @@ -274,8 +274,8 @@ void CubedSphereMeshGenerator::generate_mesh(const CubedSphereGrid& csGrid, cons const idx_t nCellsTotal = nCellsArray - 6 * 4 * nHalo * nHalo; // Projection and jacobian. - const auto* const csProjection = castProjection(csGrid.projection().get()); - const auto jacobian = NeighbourJacobian(csGrid); + const auto& csProjection = csGrid.cubedSphereProjection(); + const auto jacobian = NeighbourJacobian(csGrid); // Get partition information. const int nParts = options.get("nb_parts"); @@ -487,7 +487,7 @@ void CubedSphereMeshGenerator::generate_mesh(const CubedSphereGrid& csGrid, cons // This will only determine if tGlobal does not match t. // This is cheaper than determining the correct tGlobal. - idx_t tGlobal = csProjection->getCubedSphereTiles().indexFromXY(xy.data()); + idx_t tGlobal = csProjection.getCubedSphereTiles().indexFromXY(xy.data()); if (tGlobal == t) { // Node is an owner. @@ -773,7 +773,7 @@ void CubedSphereMeshGenerator::generate_mesh(const CubedSphereGrid& csGrid, cons nodesXy(nodeLocalIdx, YY) = xyLocal.y(); // Set lon-lat. - const PointLonLat lonLat = csProjection->lonlat(xyGlobal); + const PointLonLat lonLat = csProjection.lonlat(xyGlobal); nodesLonLat(nodeLocalIdx, LON) = lonLat.lon(); nodesLonLat(nodeLocalIdx, LAT) = lonLat.lat(); @@ -884,7 +884,7 @@ void CubedSphereMeshGenerator::generate_mesh(const CubedSphereGrid& csGrid, cons cellsXy(cellLocalIdx, YY) = xyLocal.y(); // Set lon-lat. - const PointLonLat lonLat = csProjection->lonlat(xyGlobal); + const PointLonLat lonLat = csProjection.lonlat(xyGlobal); cellsLonLat(cellLocalIdx, LON) = lonLat.lon(); cellsLonLat(cellLocalIdx, LAT) = lonLat.lat(); diff --git a/src/atlas/meshgenerator/detail/NodalCubedSphereMeshGenerator.cc b/src/atlas/meshgenerator/detail/NodalCubedSphereMeshGenerator.cc index c2bd5db9c..2ed0290d2 100644 --- a/src/atlas/meshgenerator/detail/NodalCubedSphereMeshGenerator.cc +++ b/src/atlas/meshgenerator/detail/NodalCubedSphereMeshGenerator.cc @@ -98,9 +98,9 @@ void NodalCubedSphereMeshGenerator::generate(const Grid& grid, const grid::Distr } // Get tiles - auto csprojection = castProjection(csgrid.projection().get()); + const auto& csprojection = csgrid.cubedSphereProjection(); // grid must use FV3Tiles class. - if (csprojection->getCubedSphereTiles().type() != "cubedsphere_fv3") { + if (csprojection.getCubedSphereTiles().type() != "cubedsphere_fv3") { throw_Exception("NodalCubedSphereMeshGenerator only works with FV3 tiles", Here()); } diff --git a/src/atlas/meshgenerator/detail/cubedsphere/CubedSphereUtility.cc b/src/atlas/meshgenerator/detail/cubedsphere/CubedSphereUtility.cc index 733e37018..f778d537d 100644 --- a/src/atlas/meshgenerator/detail/cubedsphere/CubedSphereUtility.cc +++ b/src/atlas/meshgenerator/detail/cubedsphere/CubedSphereUtility.cc @@ -8,71 +8,15 @@ #include "atlas/meshgenerator/detail/cubedsphere/CubedSphereUtility.h" #include "atlas/grid/CubedSphereGrid.h" #include "atlas/grid/Iterator.h" +#include "atlas/projection/Jacobian.h" #include "atlas/projection/detail/CubedSphereProjectionBase.h" + namespace atlas { namespace meshgenerator { namespace detail { namespace cubedsphere { - -// ----------------------------------------------------------------------------- -// Projection cast -// ----------------------------------------------------------------------------- - -const CubedSphereProjectionBase* castProjection(const ProjectionImpl* projectionPtr) { - const auto* const csProjectionPtr = dynamic_cast(projectionPtr); - - if (!csProjectionPtr) { - throw_Exception("Cannot cast ProjectionImpl* to CubedSphereProjectionBase*.", Here()); - } - return csProjectionPtr; -} - -// ----------------------------------------------------------------------------- -// Jacobian2 class -// ----------------------------------------------------------------------------- - -Jacobian2::Jacobian2(double df0_by_dx0, double df0_by_dx1, double df1_by_dx0, double df1_by_dx1): - df0_by_dx0_(df0_by_dx0), df0_by_dx1_(df0_by_dx1), df1_by_dx0_(df1_by_dx0), df1_by_dx1_(df1_by_dx1) {} - -Jacobian2::Jacobian2(const Point2& f00, const Point2& f10, const Point2& f01, double dx0, double dx1): - Jacobian2((f10[0] - f00[0]) / dx0, (f01[0] - f00[0]) / dx1, (f10[1] - f00[1]) / dx0, (f01[1] - f00[1]) / dx1) {} - -Jacobian2::Jacobian2(const Point2& f00, const Point2& f10, const Point2& f01): Jacobian2(f00, f10, f01, 1., 1.) {} - -double Jacobian2::det() const { - return df0_by_dx0_ * df1_by_dx1_ - df0_by_dx1_ * df1_by_dx0_; -} - -Jacobian2 Jacobian2::operator*(double a) const { - return Jacobian2(df0_by_dx0_ * a, df0_by_dx1_ * a, df1_by_dx0_ * a, df1_by_dx1_ * a); -} - -Point2 Jacobian2::operator*(const Point2& dx) const { - return Point2(dx[0] * df0_by_dx0_ + df0_by_dx1_ * dx[1], dx[0] * df1_by_dx0_ + df1_by_dx1_ * dx[1]); -} - -Jacobian2 Jacobian2::operator*(const Jacobian2& Jb) const { - return Jacobian2(df0_by_dx0_ * Jb.df0_by_dx0_ + df0_by_dx1_ * Jb.df1_by_dx0_, - df0_by_dx0_ * Jb.df0_by_dx1_ + df0_by_dx1_ * Jb.df1_by_dx1_, - df1_by_dx0_ * Jb.df0_by_dx0_ + df1_by_dx1_ * Jb.df1_by_dx0_, - df1_by_dx0_ * Jb.df0_by_dx1_ + df1_by_dx1_ * Jb.df1_by_dx1_); -} - -Jacobian2 Jacobian2::inverse() const { - return Jacobian2(df1_by_dx1_, -df0_by_dx1_, -df1_by_dx0_, df0_by_dx0_) * (1. / det()); -} - -Jacobian2 Jacobian2::sign() const { - const double smallNumber = det() * std::numeric_limits::epsilon(); - const auto signValue = [&](double number) -> double { - return std::abs(number) < smallNumber ? 0. : number < 0. ? -1. : 1.; - }; - return Jacobian2(signValue(df0_by_dx0_), signValue(df0_by_dx1_), signValue(df1_by_dx0_), signValue(df1_by_dx1_)); -} - - // ----------------------------------------------------------------------------- // NeighbourJacobian class // ----------------------------------------------------------------------------- @@ -89,7 +33,7 @@ NeighbourJacobian::NeighbourJacobian(const CubedSphereGrid& csGrid) { } // Get projection. - csProjection_ = castProjection(csGrid.projection().get()); + csProjection_ = &csGrid.cubedSphereProjection(); // Get tiles. const auto& csTiles = csProjection_->getCubedSphereTiles(); @@ -101,78 +45,42 @@ NeighbourJacobian::NeighbourJacobian(const CubedSphereGrid& csGrid) { // Get cell width. const double cellWidth = 90. / N_; - // Get xy of points (i = 0, j = 0), (i = 1, j = 0) and (i = 0, j = 1) on tiles. - std::array xy00; - std::array xy10; - std::array xy01; - - // Loop over grid points. - auto tijIt = csGrid.tij().begin(); - for (const PointXY& xy : csGrid.xy()) { - const auto t = static_cast((*tijIt).t()); - const idx_t i = (*tijIt).i(); - const idx_t j = (*tijIt).j(); - - - if (i == 0 && j == 0) { - xy00[t] = xy; - } - else if (i == 1 && j == 0) { - xy10[t] = xy; - } - else if (i == 0 && j == 1) { - xy01[t] = xy; - } - ++tijIt; - } for (size_t t = 0; t < 6; ++t) { // Calculate tile Jacobians. - dxy_by_dij_[t] = Jacobian2(xy00[t], xy10[t], xy01[t]); - - // Rescale by cell width (gains an extra couple of decimal places of precision). - dxy_by_dij_[t] = dxy_by_dij_[t].sign() * cellWidth; + dxy_by_dij_[t] = csTiles.tileJacobian(t) * cellWidth; // Get inverse. dij_by_dxy_[t] = dxy_by_dij_[t].inverse(); // Set xy00. Grid point needs moving to (i = 0, j = 0). - xy00_[t] = xy00[t] + dxy_by_dij_[t] * PointIJ(-0.5, -0.5); - - // Get other three corners so we can work out xy min/max. - const PointXY xyN0 = xy00_[t] + dxy_by_dij_[t] * PointIJ(N_, 0); - const PointXY xyNN = xy00_[t] + dxy_by_dij_[t] * PointIJ(N_, N_); - const PointXY xy0N = xy00_[t] + dxy_by_dij_[t] * PointIJ(0, N_); + xy00_[t] = csTiles.tileCentre(t) + dxy_by_dij_[t] * PointIJ(-0.5 * N_, -0.5 * N_); // Get xy min/max. - std::tie(xyMin_[t].x(), xyMax_[t].x()) = std::minmax({xy00_[t].x(), xyN0.x(), xyNN.x(), xy0N.x()}); - std::tie(xyMin_[t].y(), xyMax_[t].y()) = std::minmax({xy00_[t].y(), xyN0.y(), xyNN.y(), xy0N.y()}); - - // Round to nearest degree. - xyMin_[t].x() = std::round(xyMin_[t].x()); - xyMax_[t].x() = std::round(xyMax_[t].x()); - xyMin_[t].y() = std::round(xyMin_[t].y()); - xyMax_[t].y() = std::round(xyMax_[t].y()); + xyMin_[t].x() = csTiles.tileCentre(t).x() + 45.; + xyMax_[t].x() = csTiles.tileCentre(t).x() - 45.; + xyMin_[t].y() = csTiles.tileCentre(t).y() + 45.; + xyMax_[t].y() = csTiles.tileCentre(t).y() - 45.; // Neighbour assignment lambda. const auto neighbourAssignment = [&](TileEdge::k k) -> void { // Shift points in to neighbouring tiles. - PointIJ ijDisplacement; + PointIJ Dij; switch (k) { case TileEdge::LEFT: { - ijDisplacement = PointIJ(-2, 0); + Dij = PointIJ(-2, 0); break; } case TileEdge::BOTTOM: { - ijDisplacement = PointIJ(0, -2); + Dij = PointIJ(0, -2); break; } case TileEdge::RIGHT: { - ijDisplacement = PointIJ(N_, 0); + Dij = PointIJ(N_, 0); break; } case TileEdge::TOP: { - ijDisplacement = PointIJ(0, N_); + Dij = PointIJ(0, N_); break; } case TileEdge::UNDEFINED: { @@ -180,13 +88,15 @@ NeighbourJacobian::NeighbourJacobian(const CubedSphereGrid& csGrid) { } } - // Convert displacement from ij to xy. - const PointXY xyDisplacement = dxy_by_dij_[t] * ijDisplacement; + // half-index displacements. + const auto dij00 = PointIJ(0.5, 0.5); + const auto dij10 = PointIJ(1.5, 0.5); + const auto dij01 = PointIJ(0.5, 1.5); // Get neighbour xy points in xy space local to tile. - const PointXY xy00Local = xy00[t] + xyDisplacement; - const PointXY xy10Local = xy10[t] + xyDisplacement; - const PointXY xy01Local = xy01[t] + xyDisplacement; + const PointXY xy00Local = xy00_[t] + dxy_by_dij_[t] * (Dij + dij00); + const PointXY xy10Local = xy00_[t] + dxy_by_dij_[t] * (Dij + dij10); + const PointXY xy01Local = xy00_[t] + dxy_by_dij_[t] * (Dij + dij01); // Convert from local xy to global xy. const PointXY xy00Global = csTiles.tileCubePeriodicity(xy00Local, static_cast(t)); @@ -197,10 +107,16 @@ NeighbourJacobian::NeighbourJacobian(const CubedSphereGrid& csGrid) { neighbours_[t].t_[k] = csTiles.indexFromXY(xy00Global.data()); // Set Jacobian of global xy with respect to local ij. - auto dxyGlobal_by_dij = Jacobian2(xy00Global, xy10Global, xy01Global); + auto dxyGlobal_by_dij = Jacobian{{xy10Global[0] - xy00Global[0], xy01Global[0] - xy00Global[0]}, + {xy10Global[1] - xy00Global[1], xy01Global[1] - xy00Global[1]}}; // Rescale by cell width (gains an extra couple of decimal places of precision). - dxyGlobal_by_dij = dxyGlobal_by_dij.sign() * cellWidth; + const auto sign = [](double num, double tol = 360. * std::numeric_limits::epsilon()) { + return std::abs(num) < tol ? 0 : num > 0 ? 1. : -1.; + }; + dxyGlobal_by_dij = Jacobian{sign(dxyGlobal_by_dij[0][0]), sign(dxyGlobal_by_dij[0][1]), + sign(dxyGlobal_by_dij[1][0]), sign(dxyGlobal_by_dij[1][1])} * + cellWidth; // Chain rule to get Jacobian with respect to local xy. neighbours_[t].dxyGlobal_by_dxyLocal_[k] = dxyGlobal_by_dij * dij_by_dxy_[t]; @@ -222,8 +138,8 @@ NeighbourJacobian::NeighbourJacobian(const CubedSphereGrid& csGrid) { PointXY NeighbourJacobian::xy(const PointIJ& ij, idx_t t) const { // Get jacobian. - const Jacobian2& jac = dxy_by_dij_[static_cast(t)]; - const PointXY& xy00 = xy00_[static_cast(t)]; + const Jacobian& jac = dxy_by_dij_[static_cast(t)]; + const PointXY& xy00 = xy00_[static_cast(t)]; // Return ij return xy00 + jac * ij; @@ -231,8 +147,8 @@ PointXY NeighbourJacobian::xy(const PointIJ& ij, idx_t t) const { PointIJ NeighbourJacobian::ij(const PointXY& xy, idx_t t) const { // Get jacobian. - const Jacobian2& jac = dij_by_dxy_[static_cast(t)]; - const PointXY& xy00 = xy00_[static_cast(t)]; + const Jacobian& jac = dij_by_dxy_[static_cast(t)]; + const PointXY& xy00 = xy00_[static_cast(t)]; // Return ij return jac * (xy - xy00); @@ -290,7 +206,7 @@ PointTXY NeighbourJacobian::xyLocalToGlobal(const PointXY& xyLocal, idx_t tLocal // Get reference points and jacobian. const PointXY& xy00Local_ = neighbours_[static_cast(tLocal)].xy00Local_[k]; const PointXY& xy00Global_ = neighbours_[static_cast(tLocal)].xy00Global_[k]; - const Jacobian2& jac = neighbours_[static_cast(tLocal)].dxyGlobal_by_dxyLocal_[k]; + const Jacobian& jac = neighbours_[static_cast(tLocal)].dxyGlobal_by_dxyLocal_[k]; // Get t. tGlobal = neighbours_[static_cast(tLocal)].t_[k]; diff --git a/src/atlas/meshgenerator/detail/cubedsphere/CubedSphereUtility.h b/src/atlas/meshgenerator/detail/cubedsphere/CubedSphereUtility.h index bcaff0e8a..7fb12f407 100644 --- a/src/atlas/meshgenerator/detail/cubedsphere/CubedSphereUtility.h +++ b/src/atlas/meshgenerator/detail/cubedsphere/CubedSphereUtility.h @@ -8,6 +8,7 @@ #pragma once #include "atlas/library/config.h" +#include "atlas/projection/detail/ProjectionImpl.h" #include "atlas/runtime/Exception.h" #include "atlas/util/Point.h" @@ -20,7 +21,6 @@ namespace atlas {} namespace atlas { namespace projection { namespace detail { -class ProjectionImpl; class CubedSphereProjectionBase; } // namespace detail } // namespace projection @@ -31,8 +31,6 @@ namespace meshgenerator { namespace detail { namespace cubedsphere { -using namespace projection::detail; - /// Enum for (i, j, t) coordinate fields. struct Coordinates { enum k : idx_t @@ -55,9 +53,6 @@ struct TileEdge { }; }; -/// Cast Projection to CubedSphereProjectionBase. -const CubedSphereProjectionBase* castProjection(const ProjectionImpl* projectionPtr); - /// Class to store (i, j) indices as a Point2 coordinate. class PointIJ : public Point2 { public: @@ -115,61 +110,6 @@ class PointTIJ : public std::pair { const PointIJ& ij() const { return second; } }; -/// \brief Jacobian class for 2 dimensional vector fields. -/// -/// \details Wrapper class for an Eigen matrix which stores the partial -/// of f(x). Objects can be constructed directly from the four partial -/// derivatives, or by supplying three Point2 objects with the -/// following relative positions: -/// -/// ^ -/// | *f(X0, X1 + dx1) -/// | -/// x1 | -/// | -/// | *f(X0, X1) *f(X0 + dx0, X1) -/// +----------------------------> -/// x0 -class Jacobian2 { -public: - /// Default constructor. - Jacobian2() = default; - - /// Partial derivative constructor. - Jacobian2(double df0_by_dx0, double df0_by_dx1, double df1_by_dx0, double df1_by_dx1); - - /// Discrete point constructor (explicit dx). - Jacobian2(const Point2& f00, const Point2& f10, const Point2& f01, double dx0, double dx1); - - /// Discrete point contructor (implicit dx). - Jacobian2(const Point2& f00, const Point2& f10, const Point2& f01); - - /// Determinant of Jacobian. - double det() const; - - /// Jacobian-scalar multiplication. - Jacobian2 operator*(double a) const; - - /// Jacobian-vector multiplication. - Point2 operator*(const Point2& dx) const; - - /// Jacobian-Jacobian multiplication. - Jacobian2 operator*(const Jacobian2& Jb) const; - - /// Inverse Jacobian (partial derivatives of x(f)). - Jacobian2 inverse() const; - - /// Get signed elements of matrix (i.e, 0, +1, -1). - Jacobian2 sign() const; - -private: - // Data storage. - double df0_by_dx0_; - double df0_by_dx1_; - double df1_by_dx0_; - double df1_by_dx1_; -}; - /// \brief Class to convert between ij and xy on a tile and its four /// surrounding neighbours. /// @@ -178,6 +118,9 @@ class Jacobian2 { /// is specifcially written to comupute the (x, y) and (t, i, j) /// coordinates of halos that extend across tile boundaries. class NeighbourJacobian { +private: + using Jacobian = projection::Jacobian; + public: /// Default constructor. NeighbourJacobian() = default; @@ -211,16 +154,16 @@ class NeighbourJacobian { private: // Pointer to grid projection. - const CubedSphereProjectionBase* csProjection_{}; + const projection::detail::CubedSphereProjectionBase* csProjection_{}; // Grid size. idx_t N_{}; // Jacobian of xy with respect to ij for each tile. - std::array dxy_by_dij_{}; + std::array dxy_by_dij_{}; // Jacobian of ij with respect to xy for each tile. - std::array dij_by_dxy_{}; + std::array dij_by_dxy_{}; // Lower-left xy position on each tile. std::array xy00_{}; @@ -237,7 +180,7 @@ class NeighbourJacobian { std::array t_{}; // Jacobian of remote xy with respect to local xy. - std::array dxyGlobal_by_dxyLocal_{}; + std::array dxyGlobal_by_dxyLocal_{}; // Lower left most local xy position on neighbour tiles. std::array xy00Local_; diff --git a/src/atlas/parallel/omp/sort.h b/src/atlas/parallel/omp/sort.h index 85eb8a518..2d0bb6b15 100644 --- a/src/atlas/parallel/omp/sort.h +++ b/src/atlas/parallel/omp/sort.h @@ -17,7 +17,7 @@ #include "atlas/library/config.h" #include "atlas/parallel/omp/omp.h" -#if ATLAS_HAVE_OMP +#if ATLAS_HAVE_OMP && ATLAS_OMP_TASK_SUPPORTED #include #define ATLAS_HAVE_OMP_SORTING 1 #else diff --git a/src/atlas/projection/Jacobian.h b/src/atlas/projection/Jacobian.h new file mode 100644 index 000000000..e1003c5b2 --- /dev/null +++ b/src/atlas/projection/Jacobian.h @@ -0,0 +1,118 @@ +/* + * (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 + +#include +#include +#include +#include +#include + +#include "eckit/geometry/Point2.h" + +//--------------------------------------------------------------------------------------------------------------------- + +namespace atlas { + +using Point2 = eckit::geometry::Point2; + +//--------------------------------------------------------------------------------------------------------------------- +namespace projection { + +class Jacobian : public std::array, 2> { +public: + using std::array, 2>::array; + + static Jacobian identity() { return Jacobian{1., 0., 0., 1.}; } + + Jacobian() = default; + + Jacobian(double j00, double j01, double j10, double j11) { + (*this)[0][0] = j00; + (*this)[0][1] = j01; + (*this)[1][0] = j10; + (*this)[1][1] = j11; + } + + Jacobian(std::initializer_list> list): + Jacobian{*(list.begin()->begin()), *(list.begin()->begin() + 1), *((list.begin() + 1)->begin()), + *((list.begin() + 1)->begin() + 1)} {} + + Jacobian operator-(const Jacobian& jac) const { + return Jacobian{(*this)[0][0] - jac[0][0], (*this)[0][1] - jac[0][1], (*this)[1][0] - jac[1][0], + (*this)[1][1] - jac[1][1]}; + } + + Jacobian operator+(const Jacobian& jac) const { + return Jacobian{(*this)[0][0] + jac[0][0], (*this)[0][1] + jac[0][1], (*this)[1][0] + jac[1][0], + (*this)[1][1] + jac[1][1]}; + } + + Jacobian operator*(double a) const { + return Jacobian{(*this)[0][0] * a, (*this)[0][1] * a, (*this)[1][0] * a, (*this)[1][1] * a}; + } + + Jacobian operator*(const Jacobian& jac) const { + return Jacobian{(*this)[0][0] * jac[0][0] + (*this)[0][1] * jac[1][0], + (*this)[0][0] * jac[0][1] + (*this)[0][1] * jac[1][1], + (*this)[1][0] * jac[0][0] + (*this)[1][1] * jac[1][0], + (*this)[1][0] * jac[0][1] + (*this)[1][1] * jac[1][1]}; + } + + Point2 operator*(const Point2& x) const { + return Point2{(*this)[0][0] * x[0] + (*this)[0][1] * x[1], (*this)[1][0] * x[0] + (*this)[1][1] * x[1]}; + } + + double norm() const { + return std::sqrt((*this)[0][0] * (*this)[0][0] + (*this)[0][1] * (*this)[0][1] + (*this)[1][0] * (*this)[1][0] + + (*this)[1][1] * (*this)[1][1]); + } + + double determinant() const { return (*this)[0][0] * (*this)[1][1] - (*this)[0][1] * (*this)[1][0]; } + + Jacobian inverse() const { + return Jacobian{(*this)[1][1], -(*this)[0][1], -(*this)[1][0], (*this)[0][0]} * (1. / determinant()); + } + + Jacobian transpose() const { return Jacobian{(*this)[0][0], (*this)[1][0], (*this)[0][1], (*this)[1][1]}; } + + double dx_dlon() const { return (*this)[JDX][JDLON]; } + double dy_dlon() const { return (*this)[JDY][JDLON]; } + double dx_dlat() const { return (*this)[JDX][JDLAT]; } + double dy_dlat() const { return (*this)[JDY][JDLAT]; } + + double dlon_dx() const { return (*this)[JDLON][JDX]; } + double dlon_dy() const { return (*this)[JDLON][JDY]; } + double dlat_dx() const { return (*this)[JDLAT][JDX]; } + double dlat_dy() const { return (*this)[JDLAT][JDY]; } + + friend std::ostream& operator<<(std::ostream& os, const Jacobian& jac) { + os << jac[0][0] << " " << jac[0][1] << "\n" << jac[1][0] << " " << jac[1][1]; + return os; + } + +private: + enum + { + JDX = 0, + JDY = 1 + }; + enum + { + JDLON = 0, + JDLAT = 1 + }; +}; + +//--------------------------------------------------------------------------------------------------------------------- + +} // namespace projection +} // namespace atlas diff --git a/src/atlas/projection/Projection.h b/src/atlas/projection/Projection.h index 273dcb33f..a42ceb4a4 100644 --- a/src/atlas/projection/Projection.h +++ b/src/atlas/projection/Projection.h @@ -14,6 +14,7 @@ #include "atlas/domain/Domain.h" #include "atlas/library/config.h" +#include "atlas/projection/Jacobian.h" #include "atlas/projection/detail/ProjectionImpl.h" #include "atlas/util/ObjectHandle.h" @@ -49,7 +50,7 @@ class ProjectionImpl; class Projection : DOXYGEN_HIDE(public util::ObjectHandle) { public: using Spec = util::Config; - using Jacobian = projection::detail::ProjectionImpl::Jacobian; + using Jacobian = projection::Jacobian; public: using Handle::Handle; diff --git a/src/atlas/projection/detail/CubedSphereEquiAnglProjection.cc b/src/atlas/projection/detail/CubedSphereEquiAnglProjection.cc index d8748feb3..b7b2a56a8 100644 --- a/src/atlas/projection/detail/CubedSphereEquiAnglProjection.cc +++ b/src/atlas/projection/detail/CubedSphereEquiAnglProjection.cc @@ -27,6 +27,34 @@ namespace { static constexpr bool debug = false; // constexpr so compiler can optimize `if ( debug ) { ... }` out static constexpr double deg2rad = atlas::util::Constants::degreesToRadians(); static constexpr double rad2deg = atlas::util::Constants::radiansToDegrees(); + +// Define small number relative to 360. +constexpr double epsilon = std::numeric_limits::epsilon() * 360.; + +// Define some "fuzzy" comparison operators. + +// a is approximately equal to b. +bool equal(double a, double b) { + return std::abs(a - b) <= epsilon; +} +// a is less than b. +bool lessThan(double a, double b) { + return a < b && !equal(a, b); +} +// a is greater than b. +bool greaterThan(double a, double b) { + return a > b && !equal(a, b); +} +// a is less than or approximately equal to b. +bool lessEqual(double a, double b) { + return a < b || equal(a, b); +} +// a is greater than or approximately equal to b. +bool greaterEqual(double a, double b) { + return a > b || equal(a, b); +} + + } // namespace namespace atlas { @@ -41,6 +69,137 @@ CubedSphereEquiAnglProjection::CubedSphereEquiAnglProjection(const eckit::Parame // ------------------------------------------------------------------------------------------------- +void CubedSphereEquiAnglProjection::xy2alphabeta(double crd[], idx_t t) const { + // Get tile centre. + const auto& xyCentre = getCubedSphereTiles().tileCentre(static_cast(t)); + + // Check that xy coordinate is within valid "+" shaped halo region. + const auto inCross = [&](const double crd[]) -> bool { + return (greaterEqual(crd[XX], xyCentre[XX] - 45.) && lessEqual(crd[XX], xyCentre[XX] + 45.)) || + (greaterEqual(crd[YY], xyCentre[YY] - 45.) && lessEqual(crd[YY], xyCentre[YY] + 45.)); + }; + if (!inCross(crd)) { + auto sStream = std::stringstream(); + sStream << "xy coordinate (" << crd[0] << ", " << crd[1] << ") is not in range for tile " << t << "."; + ATLAS_THROW_EXCEPTION(sStream.str()); + } + + // Get alphaBeta Jacobian. + const auto alphabetaJacobian = getCubedSphereTiles().tileJacobian(static_cast(t)).inverse(); + + // Set (alpha, beta) coord. + const Point2 alphabeta = alphabetaJacobian * (Point2(crd) - xyCentre); + crd[0] = alphabeta[0]; + crd[1] = alphabeta[1]; + + // Define correction. + const auto correction = [](const double crd[]) -> double { + return rad2deg * std::atan(std::tan(crd[0] * deg2rad) * std::tan(crd[1] * deg2rad)); + }; + + // Correct halo (alpha, beta) coord. + if (lessThan(crd[0], -45.)) { + // Left. + crd[1] = -correction(crd); + } + else if (greaterThan(crd[0], 45.)) { + // Right. + crd[1] = correction(crd); + } + else if (lessThan(crd[1], -45.)) { + // Bottom. + crd[0] = -correction(crd); + } + else if (greaterThan(crd[1], 45.)) { + // Top. + crd[0] = correction(crd); + } +} + +// ------------------------------------------------------------------------------------------------- + +void CubedSphereEquiAnglProjection::alphabeta2xy(double crd[], idx_t t) const { + // Define correction. + const auto correction1 = [](const double crd[]) -> double { + return rad2deg * std::atan(std::tan(crd[1] * deg2rad) / std::tan(crd[0] * deg2rad)); + }; + const auto correction2 = [](const double crd[]) -> double { + return rad2deg * std::atan(std::tan(crd[0] * deg2rad) / std::tan(crd[1] * deg2rad)); + }; + + // Correct halo (alpha, beta) coord. + if (lessThan(crd[0], -45.) && greaterThan(crd[1], crd[0]) && lessEqual(crd[1], -crd[0])) { + // Left trapezium. + crd[1] = -correction1(crd); + } + else if (greaterThan(crd[0], 45.) && greaterEqual(crd[1], -crd[0]) && lessThan(crd[1], crd[0])) { + // Right trapezium. + crd[1] = correction1(crd); + } + else if (lessThan(crd[1], -45.) && greaterEqual(crd[0], crd[1]) && lessThan(crd[0], -crd[1])) { + // Bottom trapezium. + crd[0] = -correction2(crd); + } + else if (greaterThan(crd[1], 45.) && greaterThan(crd[0], -crd[1]) && lessEqual(crd[0], crd[1])) { + // Top trapezium. + crd[0] = correction2(crd); + } + + // Get tile centre. + const auto& xyCentre = getCubedSphereTiles().tileCentre(static_cast(t)); + + // Get xy Jacobian. + const auto xyJacobian = getCubedSphereTiles().tileJacobian(static_cast(t)); + + // Set xy coord. + const Point2 xy = xyJacobian * Point2(crd) + xyCentre; + crd[XX] = xy[XX]; + crd[YY] = xy[YY]; +} + +// ------------------------------------------------------------------------------------------------- + +Jacobian CubedSphereEquiAnglProjection::jacobian(const PointLonLat& lonlat, idx_t t) const { + // Note: angular units cancel, so we leave all values in radians. + + // Convert lonlat to xyz on unit sphere. + const double lambda = lonlat.lon() * deg2rad; + const double phi = lonlat.lat() * deg2rad; + auto xyz = PointXYZ{std::cos(lambda) * std::cos(phi), std::sin(lambda) * std::cos(phi), -std::sin(phi)}; + + // Get derivatives of xyz with respect to lambda and phi. + auto dxyz_by_dlambda = PointXYZ{-std::sin(lambda) * std::cos(phi), std::cos(lambda) * std::cos(phi), 0.}; + auto dxyz_by_dphi = PointXYZ{-std::cos(lambda) * std::sin(phi), -std::sin(lambda) * std::sin(phi), -std::cos(phi)}; + + // Rotate vectors. + const auto& tiles = getCubedSphereTiles(); + tiles.unrotate(t, xyz.data()); + tiles.unrotate(t, dxyz_by_dlambda.data()); + tiles.unrotate(t, dxyz_by_dphi.data()); + + // Get derivatives of a and b with respect to xyz. + // Note: a and b are xy displacements from the tile centre, *not* the + // (alpha, beta) coordinates defined in the tileJacobian method. + + const double inv_x2y2 = 1. / (xyz.x() * xyz.x() + xyz.y() * xyz.y()); + const double inv_x2z2 = 1. / (xyz.x() * xyz.x() + xyz.z() * xyz.z()); + + const auto da_by_dxyz = PointXYZ{-xyz.y() * inv_x2y2, xyz.x() * inv_x2y2, 0.}; + const auto db_by_dxyz = PointXYZ{xyz.z() * inv_x2z2, 0., -xyz.x() * inv_x2z2}; + + // Use chain rule to get Jacobian. + const auto& dot = eckit::geometry::Point3::dot; + return Jacobian{{dot(da_by_dxyz, dxyz_by_dlambda), dot(da_by_dxyz, dxyz_by_dphi)}, + {dot(db_by_dxyz, dxyz_by_dlambda), dot(db_by_dxyz, dxyz_by_dphi)}}; +} + +Jacobian CubedSphereEquiAnglProjection::alphabetaJacobian(const PointLonLat& lonlat, idx_t t) const { + const auto& tiles = getCubedSphereTiles(); + return tiles.tileJacobian(t).inverse() * jacobian(lonlat, t); +} + +// ------------------------------------------------------------------------------------------------- + void CubedSphereEquiAnglProjection::lonlat2xy(double crd[]) const { if (debug) { Log::info() << "equiangular lonlat2xy start : lonlat = " << crd[LON] << " " << crd[LAT] << std::endl; @@ -104,8 +263,10 @@ void CubedSphereEquiAnglProjection::xy2lonlat(double crd[]) const { // ------------------------------------------------------------------------------------------------- -ProjectionImpl::Jacobian CubedSphereEquiAnglProjection::jacobian(const PointLonLat&) const { - ATLAS_NOTIMPLEMENTED; +Jacobian CubedSphereEquiAnglProjection::jacobian(const PointLonLat& lonlat) const { + const auto& tiles = getCubedSphereTiles(); + const idx_t t = tiles.indexFromLonLat(lonlat.data()); + return jacobian(lonlat, t); } // ------------------------------------------------------------------------------------------------- diff --git a/src/atlas/projection/detail/CubedSphereEquiAnglProjection.h b/src/atlas/projection/detail/CubedSphereEquiAnglProjection.h index 1d69870b2..c297919df 100644 --- a/src/atlas/projection/detail/CubedSphereEquiAnglProjection.h +++ b/src/atlas/projection/detail/CubedSphereEquiAnglProjection.h @@ -9,6 +9,7 @@ #include "atlas/array.h" #include "atlas/domain.h" +#include "atlas/projection/Jacobian.h" #include "atlas/projection/detail/CubedSphereProjectionBase.h" #include "atlas/projection/detail/ProjectionImpl.h" @@ -27,6 +28,18 @@ class CubedSphereEquiAnglProjection final : public CubedSphereProjectionBase { static std::string static_type() { return "cubedsphere_equiangular"; } std::string type() const override { return static_type(); } + /// @brief Convert (x, y) coordinate to (alpha, beta) on tile t. + void xy2alphabeta(double crd[], idx_t t) const override; + + /// @brief Convert (alpha, beta) coordinate to (x, y) on tile t. + void alphabeta2xy(double crd[], idx_t t) const override; + + /// @brief Jacobian of (x, y) with respect to (lon, lat) on tile t + Jacobian jacobian(const PointLonLat& lonlat, idx_t t) const override; + + /// @brief Jacobian of (alpha, beta) with respect to (lon, lat) on tile t + Jacobian alphabetaJacobian(const PointLonLat& lonlat, idx_t t) const override; + // projection and inverse projection void xy2lonlat(double crd[]) const override; void lonlat2xy(double crd[]) const override; diff --git a/src/atlas/projection/detail/CubedSphereEquiDistProjection.cc b/src/atlas/projection/detail/CubedSphereEquiDistProjection.cc index 925020d5b..0ff91b55c 100644 --- a/src/atlas/projection/detail/CubedSphereEquiDistProjection.cc +++ b/src/atlas/projection/detail/CubedSphereEquiDistProjection.cc @@ -35,6 +35,30 @@ CubedSphereEquiDistProjection::CubedSphereEquiDistProjection(const eckit::Parame // ------------------------------------------------------------------------------------------------- +void CubedSphereEquiDistProjection::xy2alphabeta(double crd[], idx_t t) const { + throw_NotImplemented("xy2alphabeta not implemented for CubedSphereEquiDistProjection", Here()); +} + +// ------------------------------------------------------------------------------------------------- + +void CubedSphereEquiDistProjection::alphabeta2xy(double crd[], idx_t t) const { + throw_NotImplemented("alphabeta2xy not implemented for CubedSphereEquiDistProjection", Here()); +} + +// ------------------------------------------------------------------------------------------------- + +Jacobian CubedSphereEquiDistProjection::jacobian(const PointLonLat& lonlat, idx_t t) const { + throw_NotImplemented("jacobian not implemented for CubedSphereEquiDistProjection", Here()); +} + +// ------------------------------------------------------------------------------------------------- + +Jacobian CubedSphereEquiDistProjection::alphabetaJacobian(const PointLonLat& lonlat, idx_t t) const { + throw_NotImplemented("alphabetaJacobian not implemented for CubedSphereEquiDistProjection", Here()); +} + +// ------------------------------------------------------------------------------------------------- + void CubedSphereEquiDistProjection::lonlat2xy(double crd[]) const { if (debug) { Log::info() << "equidist lonlat2xy start : lonlat = " << crd[LON] << " " << crd[LAT] << std::endl; @@ -93,7 +117,7 @@ void CubedSphereEquiDistProjection::xy2lonlat(double crd[]) const { // ------------------------------------------------------------------------------------------------- -ProjectionImpl::Jacobian CubedSphereEquiDistProjection::jacobian(const PointLonLat&) const { +Jacobian CubedSphereEquiDistProjection::jacobian(const PointLonLat&) const { ATLAS_NOTIMPLEMENTED; } diff --git a/src/atlas/projection/detail/CubedSphereEquiDistProjection.h b/src/atlas/projection/detail/CubedSphereEquiDistProjection.h index 6fc59f1a7..6d7ea0b7f 100644 --- a/src/atlas/projection/detail/CubedSphereEquiDistProjection.h +++ b/src/atlas/projection/detail/CubedSphereEquiDistProjection.h @@ -27,6 +27,18 @@ class CubedSphereEquiDistProjection final : public CubedSphereProjectionBase { static std::string static_type() { return "cubedsphere_equidistant"; } std::string type() const override { return static_type(); } + /// @brief Convert (x, y) coordinate to (alpha, beta) on tile t. + void xy2alphabeta(double crd[], idx_t t) const override; + + /// @brief Convert (alpha, beta) coordinate to (x, y) on tile t. + void alphabeta2xy(double crd[], idx_t t) const override; + + /// @brief Jacobian of (alpha, beta) with respect to (lon, lat) on tile t + Jacobian alphabetaJacobian(const PointLonLat& lonlat, idx_t t) const override; + + /// @brief Jacobian of (x, y) with respect to (lon, lat) on tile t + Jacobian jacobian(const PointLonLat& lonlat, idx_t t) const override; + // projection and inverse projection void xy2lonlat(double crd[]) const override; void lonlat2xy(double crd[]) const override; diff --git a/src/atlas/projection/detail/CubedSphereProjectionBase.cc b/src/atlas/projection/detail/CubedSphereProjectionBase.cc index ba05d72b6..a8295be4b 100644 --- a/src/atlas/projection/detail/CubedSphereProjectionBase.cc +++ b/src/atlas/projection/detail/CubedSphereProjectionBase.cc @@ -232,6 +232,22 @@ void CubedSphereProjectionBase::alphabetat2xy(const idx_t t, const double ab[], // ------------------------------------------------------------------------------------------------- +Point2 CubedSphereProjectionBase::xy2alphabeta(const Point2& xy, idx_t t) const { + auto alphabeta = Point2(xy); + xy2alphabeta(alphabeta.data(), t); + return alphabeta; +} + +// ------------------------------------------------------------------------------------------------- + +Point2 CubedSphereProjectionBase::alphabeta2xy(const Point2& alphabeta, idx_t t) const { + auto xy = Point2(alphabeta); + alphabeta2xy(xy.data(), t); + return xy; +} + +// ------------------------------------------------------------------------------------------------- + } // namespace detail } // namespace projection } // namespace atlas diff --git a/src/atlas/projection/detail/CubedSphereProjectionBase.h b/src/atlas/projection/detail/CubedSphereProjectionBase.h index 17c34a861..cc4ff37e2 100644 --- a/src/atlas/projection/detail/CubedSphereProjectionBase.h +++ b/src/atlas/projection/detail/CubedSphereProjectionBase.h @@ -9,6 +9,7 @@ #include "atlas/grid/Tiles.h" #include "atlas/library/config.h" +#include "atlas/projection/Jacobian.h" #include "atlas/projection/detail/ProjectionImpl.h" #include "eckit/config/Parametrisation.h" #include "eckit/utils/Hash.h" @@ -23,6 +24,8 @@ namespace detail { class CubedSphereProjectionBase : public ProjectionImpl { public: + using ProjectionImpl::jacobian; + // constructor CubedSphereProjectionBase(const eckit::Parametrisation&); @@ -30,6 +33,47 @@ class CubedSphereProjectionBase : public ProjectionImpl { atlas::grid::CubedSphereTiles getCubedSphereTiles() const { return tiles_; }; + /// @brief Convert (x, y) coordinate to (alpha, beta) on tile t. + /// + /// @details Converts the Atlas xy coordinates to the angular coordinates + /// described of tile t, described by by Ronchi et al. (1996, + /// Journal of Computational Physics, 124, 93). + /// Note that the xy coordinate must lie within the domain + /// (x <= xmin && x >= xmax) || (y <= ymin && y >= ymax) where + /// xmin, xmax, ymin and ymax are the boundaries of tile t. + ///@{ + Point2 xy2alphabeta(const Point2& xy, idx_t t) const; + virtual void xy2alphabeta(double crd[], idx_t t) const = 0; + ///@} + + /// @brief Convert (alpha, beta) coordinate to (x, y) on tile t. + /// + /// @details Performs the inverse of xy2alpha beta. Note that the result is + /// degenerate when abs(alpha) > 45 && abs(beta) > 45 && + /// abs(alpha) == abs(beta). In these circumstances, the method + /// will return the anticlockwise-most of the two possible + /// values. + ///@{ + Point2 alphabeta2xy(const Point2& alphabeta, idx_t t) const; + virtual void alphabeta2xy(double crd[], idx_t) const = 0; + ///@} + + /// @brief Jacobian of (x, y) with respect to (lon, lat) on tile t + /// + /// @details Returns the Jacobian + /// ∂x/∂λ, ∂x/∂φ + /// ∂y/∂λ, ∂y/∂φ + /// for tile t. + virtual Jacobian jacobian(const PointLonLat& lonlat, idx_t t) const = 0; + + /// @brief Jacobian of (alpha, beta) with respect to (lon, lat) on tile t + /// + /// @details Returns the Jacobian + /// ∂α/∂λ, ∂α/∂φ + /// ∂β/∂λ, ∂β/∂φ + /// for tile t. + virtual Jacobian alphabetaJacobian(const PointLonLat& lonlat, idx_t t) const = 0; + protected: // projection and inverse projection void xy2lonlat_post(double xyz[], const idx_t t, double crd[]) const; diff --git a/src/atlas/projection/detail/ProjProjection.cc b/src/atlas/projection/detail/ProjProjection.cc index 62e4c8f8a..cf81f9d60 100644 --- a/src/atlas/projection/detail/ProjProjection.cc +++ b/src/atlas/projection/detail/ProjProjection.cc @@ -26,9 +26,19 @@ namespace projection { namespace detail { namespace { +struct pj_t : std::unique_ptr { + using t = std::unique_ptr; + explicit pj_t(PJ* ptr): t(ptr, &proj_destroy) {} + operator PJ*() const { return t::get(); } +}; + +struct ctx_t : std::unique_ptr { + using t = std::unique_ptr; + explicit ctx_t(PJ_CONTEXT* ptr): t(ptr, &proj_context_destroy) {} + operator PJ_CONTEXT*() const { return t::get(); } +}; bool proj_ellipsoid_params(PJ_CONTEXT* ctxt, const std::string& proj_str, double& a, double& b) { - using pj_t = ProjProjection::pj_t; bool success = false; pj_t identity(proj_create_crs_to_crs(ctxt, proj_str.c_str(), proj_str.c_str(), nullptr)); pj_t proj(proj_get_target_crs(ctxt, identity)); @@ -66,32 +76,37 @@ std::string geocentric_str(PJ_CONTEXT* ctxt, const std::string& proj_str) { } // namespace -ProjProjection::ProjProjection(const eckit::Parametrisation& param): - normalise_(param), sourceToTarget_(nullptr), sourceToGeocentric_(nullptr), context_(PJ_DEFAULT_CTX) { - ATLAS_ASSERT(param.get("proj", proj_) && !proj_.empty()); - source_encoded_ = param.get("proj_source", source_ = source_str(context_, proj_)); - geocentric_encoded_ = param.get("proj_geocentric", geocentric_ = geocentric_str(context_, proj_)); +struct Proj { + pj_t sourceToTarget_{nullptr}; + pj_t sourceToGeocentric_{nullptr}; + ctx_t context_{PJ_DEFAULT_CTX}; +}; + +ProjProjection::ProjProjection(const eckit::Parametrisation& param): normalise_(param), proj_(new Proj()) { + ATLAS_ASSERT(param.get("proj", proj_string_) && !proj_string_.empty()); + source_encoded_ = param.get("proj_source", source_ = source_str(proj_->context_, proj_string_)); + geocentric_encoded_ = param.get("proj_geocentric", geocentric_ = geocentric_str(proj_->context_, proj_string_)); // set x/y transformations to/from lon/lat and to/from geocentric coordinates { - pj_t p1(proj_create_crs_to_crs(context_, source_.c_str(), proj_.c_str(), nullptr)); + pj_t p1(proj_create_crs_to_crs(proj_->context_, source_.c_str(), proj_string_.c_str(), nullptr)); ATLAS_ASSERT(p1); - sourceToTarget_.reset(proj_normalize_for_visualization(context_, p1)); - ATLAS_ASSERT(sourceToTarget_); + proj_->sourceToTarget_.reset(proj_normalize_for_visualization(proj_->context_, p1)); + ATLAS_ASSERT(proj_->sourceToTarget_); } { - pj_t p2(proj_create_crs_to_crs(context_, source_.c_str(), geocentric_.c_str(), nullptr)); + pj_t p2(proj_create_crs_to_crs(proj_->context_, source_.c_str(), geocentric_.c_str(), nullptr)); ATLAS_ASSERT(p2); - sourceToGeocentric_.reset(proj_normalize_for_visualization(context_, p2)); - ATLAS_ASSERT(sourceToGeocentric_); + proj_->sourceToGeocentric_.reset(proj_normalize_for_visualization(proj_->context_, p2)); + ATLAS_ASSERT(proj_->sourceToGeocentric_); } // set semi-major/minor axis { double a, b; - if (proj_ellipsoid_params(context_, proj_, a, b)) { + if (proj_ellipsoid_params(proj_->context_, proj_string_, a, b)) { eckit::types::is_approximately_equal(a, b) ? extraSpec_.set("radius", a) : extraSpec_.set("semi_major_axis", a).set("semi_minor_axis", b); } @@ -99,12 +114,12 @@ ProjProjection::ProjProjection(const eckit::Parametrisation& param): // set units { - pj_t target(proj_get_target_crs(context_, sourceToTarget_)); + pj_t target(proj_get_target_crs(proj_->context_, proj_->sourceToTarget_)); ATLAS_ASSERT(target); - pj_t coord(proj_crs_get_coordinate_system(context_, target)); + pj_t coord(proj_crs_get_coordinate_system(proj_->context_, target)); ATLAS_ASSERT(coord); - ATLAS_ASSERT(proj_cs_get_axis_count(context_, coord) > 0); + ATLAS_ASSERT(proj_cs_get_axis_count(proj_->context_, coord) > 0); const char* units_c_str; if (proj_cs_get_axis_info(nullptr, coord, 0, nullptr, nullptr, nullptr, nullptr, &units_c_str, nullptr, @@ -123,7 +138,7 @@ ProjProjection::ProjProjection(const eckit::Parametrisation& param): void ProjProjection::xy2lonlat(double crd[]) const { PJ_COORD P = proj_coord(crd[XX], crd[YY], 0, 0); - P = proj_trans(sourceToTarget_, PJ_INV, P); + P = proj_trans(proj_->sourceToTarget_, PJ_INV, P); // std::memcpy(crd, &P, 2 * sizeof(double)); crd[LON] = P.enu.e; @@ -135,7 +150,7 @@ void ProjProjection::xy2lonlat(double crd[]) const { void ProjProjection::lonlat2xy(double crd[]) const { PJ_COORD P = proj_coord(crd[LON], crd[LAT], 0, 0); - P = proj_trans(sourceToTarget_, PJ_FWD, P); + P = proj_trans(proj_->sourceToTarget_, PJ_FWD, P); // std::memcpy(crd, &P, 2 * sizeof(double)); crd[XX] = P.xy.x; @@ -149,7 +164,7 @@ ProjectionImpl::Jacobian ProjProjection::jacobian(const PointLonLat&) const { PointXYZ ProjProjection::xyz(const PointLonLat& lonlat) const { PJ_COORD P = proj_coord(lonlat.lon(), lonlat.lat(), 0, 0); - P = proj_trans(sourceToGeocentric_, PJ_FWD, P); + P = proj_trans(proj_->sourceToGeocentric_, PJ_FWD, P); return {P.xyz.x, P.xyz.y, P.xyz.z}; } @@ -162,7 +177,7 @@ RectangularLonLatDomain ProjProjection::lonlatBoundingBox(const Domain& domain) ProjProjection::Spec ProjProjection::spec() const { Spec spec; spec.set("type", type()); - spec.set("proj", proj_); + spec.set("proj", proj_string_); if (source_encoded_) { spec.set("proj_source", source_); } @@ -181,7 +196,7 @@ std::string ProjProjection::units() const { void ProjProjection::hash(eckit::Hash& h) const { h.add(type()); - h.add(proj_); + h.add(proj_string_); if (source_encoded_) { h.add(source_); } diff --git a/src/atlas/projection/detail/ProjProjection.h b/src/atlas/projection/detail/ProjProjection.h index 2c75b797d..f0b49328f 100644 --- a/src/atlas/projection/detail/ProjProjection.h +++ b/src/atlas/projection/detail/ProjProjection.h @@ -16,37 +16,17 @@ #include "atlas/projection/detail/ProjectionImpl.h" #include "atlas/util/Config.h" - -extern "C" { -using PJ = struct PJconsts; -PJ* proj_destroy(PJ*); -using PJ_CONTEXT = struct projCtx_t; -PJ_CONTEXT* proj_context_destroy(PJ_CONTEXT*); -} - - namespace atlas { namespace projection { namespace detail { +// Forward declaration of Proj, defined in ProjProjection.cc +struct Proj; class ProjProjection final : public ProjectionImpl { public: - // -- Types - // None - - // -- Exceptions - // None - - // Constructors - ProjProjection(const eckit::Parametrisation&); - // Destructor - // None - - // -- Methods - static std::string static_type() { return "proj"; } // -- Overridden methods @@ -67,98 +47,17 @@ class ProjProjection final : public ProjectionImpl { std::string units() const override; void hash(eckit::Hash&) const override; - // -- Class members - // None - - // -- Class methods - // None - - // -- Members - // None - - // -- Methods - // None - - // -- Overridden methods - // None - - // -- Class members - // None - - // -- Class methods - // None - - // -- Friends - // None - -public: - // -- Types - - struct pj_t : std::unique_ptr { - using t = std::unique_ptr; - explicit pj_t(PJ* ptr): t(ptr, &proj_destroy) {} - operator PJ*() const { return t::get(); } - }; - - struct ctx_t : std::unique_ptr { - using t = std::unique_ptr; - explicit ctx_t(PJ_CONTEXT* ptr): t(ptr, &proj_context_destroy) {} - operator PJ_CONTEXT*() const { return t::get(); } - }; - - // -- Contructors - // None - - // -- Destructor - // None - - // -- Convertors - // None - - // -- Operators - // None - - // -- Methods - // None - - // -- Overridden methods - // None - - // -- Class members - // None - - // -- Class methods - // None - - // -- Members private: Normalise normalise_; - std::string proj_; + std::string proj_string_; std::string source_; std::string geocentric_; bool source_encoded_; bool geocentric_encoded_; - pj_t sourceToTarget_; - pj_t sourceToGeocentric_; - ctx_t context_; + std::unique_ptr proj_; Spec extraSpec_; - - // -- Methods - // None - - // -- Overridden methods - // None - - // -- Class members - // None - - // -- Class methods - // None - - // -- Friends - // None (how sad) }; diff --git a/src/atlas/projection/detail/ProjectionImpl.h b/src/atlas/projection/detail/ProjectionImpl.h index fcb586269..6f777a5ee 100644 --- a/src/atlas/projection/detail/ProjectionImpl.h +++ b/src/atlas/projection/detail/ProjectionImpl.h @@ -15,6 +15,7 @@ #include #include +#include "atlas/projection/Jacobian.h" #include "atlas/runtime/Exception.h" #include "atlas/util/Factory.h" #include "atlas/util/NormaliseLongitude.h" @@ -41,92 +42,8 @@ namespace detail { class ProjectionImpl : public util::Object { public: - using Spec = atlas::util::Config; - class Jacobian : public std::array, 2> { - public: - using std::array, 2>::array; - - static Jacobian identity() { return Jacobian{1., 0., 0., 1.}; } - - Jacobian() = default; - - Jacobian(double j00, double j01, double j10, double j11) { - (*this)[0][0] = j00; - (*this)[0][1] = j01; - (*this)[1][0] = j10; - (*this)[1][1] = j11; - } - - Jacobian(std::initializer_list> list): - Jacobian{*(list.begin()->begin()), *(list.begin()->begin() + 1), *((list.begin() + 1)->begin()), - *((list.begin() + 1)->begin() + 1)} {} - - Jacobian operator-(const Jacobian& jac) const { - return Jacobian{(*this)[0][0] - jac[0][0], (*this)[0][1] - jac[0][1], (*this)[1][0] - jac[1][0], - (*this)[1][1] - jac[1][1]}; - } - - Jacobian operator+(const Jacobian& jac) const { - return Jacobian{(*this)[0][0] + jac[0][0], (*this)[0][1] + jac[0][1], (*this)[1][0] + jac[1][0], - (*this)[1][1] + jac[1][1]}; - } - - Jacobian operator*(double a) const { - return Jacobian{(*this)[0][0] * a, (*this)[0][1] * a, (*this)[1][0] * a, (*this)[1][1] * a}; - } - - Jacobian operator*(const Jacobian& jac) const { - return Jacobian{(*this)[0][0] * jac[0][0] + (*this)[0][1] * jac[1][0], - (*this)[0][0] * jac[0][1] + (*this)[0][1] * jac[1][1], - (*this)[1][0] * jac[0][0] + (*this)[1][1] * jac[1][0], - (*this)[1][0] * jac[0][1] + (*this)[1][1] * jac[1][1]}; - } - - Point2 operator*(const Point2& x) const { - return Point2{(*this)[0][0] * x[0] + (*this)[0][1] * x[1], (*this)[1][0] * x[0] + (*this)[1][1] * x[1]}; - } - - double norm() const { - return std::sqrt((*this)[0][0] * (*this)[0][0] + (*this)[0][1] * (*this)[0][1] + - (*this)[1][0] * (*this)[1][0] + (*this)[1][1] * (*this)[1][1]); - } - - double determinant() const { return (*this)[0][0] * (*this)[1][1] - (*this)[0][1] * (*this)[1][0]; } - - Jacobian inverse() const { - return Jacobian{(*this)[1][1], -(*this)[0][1], -(*this)[1][0], (*this)[0][0]} * (1. / determinant()); - } - - Jacobian transpose() const { return Jacobian{(*this)[1][1], (*this)[0][1], (*this)[1][0], (*this)[0][0]}; } - - double dx_dlon() const { return (*this)[JDX][JDLON]; } - double dy_dlon() const { return (*this)[JDY][JDLON]; } - double dx_dlat() const { return (*this)[JDX][JDLAT]; } - double dy_dlat() const { return (*this)[JDY][JDLAT]; } - - double dlon_dx() const { return (*this)[JDLON][JDX]; } - double dlon_dy() const { return (*this)[JDLON][JDY]; } - double dlat_dx() const { return (*this)[JDLAT][JDX]; } - double dlat_dy() const { return (*this)[JDLAT][JDY]; } - - friend std::ostream& operator<<(std::ostream& os, const Jacobian& jac) { - os << jac[0][0] << " " << jac[0][1] << "\n" << jac[1][0] << " " << jac[1][1]; - return os; - } - - private: - enum - { - JDX = 0, - JDY = 1 - }; - enum - { - JDLON = 0, - JDLAT = 1 - }; - }; - + using Spec = atlas::util::Config; + using Jacobian = projection::Jacobian; public: static const ProjectionImpl* create(const eckit::Parametrisation& p); diff --git a/src/atlas/redistribution/detail/RedistributeGeneric.cc b/src/atlas/redistribution/detail/RedistributeGeneric.cc index 53254341a..e869e9d56 100644 --- a/src/atlas/redistribution/detail/RedistributeGeneric.cc +++ b/src/atlas/redistribution/detail/RedistributeGeneric.cc @@ -215,35 +215,42 @@ std::pair, std::vector> getUidIntersection(const std::ve // Iterate over a field, in the order of an index list, and apply a functor to // each element. -// Rank 1 overload. -template -void iterateField(const std::vector& idxList, array::ArrayView& fieldView, const Functor& f) { - for (const idx_t i : idxList) { - f(fieldView(i)); +// Recursive ForEach to visit all elements of field. +template +struct ForEach { + template + static void apply(const std::vector& idxList, array::ArrayView& fieldView, const Functor& f, + Idxs... idxs) { + // Iterate over dimension Dim of array. + for (idx_t idx = 0; idx < fieldView.shape(Dim); ++idx) { + ForEach::apply(idxList, fieldView, f, idxs..., idx); + } } -} +}; -// Rank 2 overload. -template -void iterateField(const std::vector& idxList, array::ArrayView& fieldView, const Functor f) { - for (const idx_t i : idxList) { - for (idx_t j = 0; j < fieldView.shape(1); ++j) { - f(fieldView(i, j)); +// Beginning of recursion when Dim == 0. +template +struct ForEach { + template + static void apply(const std::vector& idxList, array::ArrayView& fieldView, const Functor& f, + Idxs... idxs) { + // Iterate over dimension 0 of array in order defined by idxList. + for (idx_t idx : idxList) { + ForEach::apply(idxList, fieldView, f, idxs..., idx); } } -} +}; -// Rank 3 overload. -template -void iterateField(const std::vector& idxList, array::ArrayView& fieldView, const Functor f) { - for (const idx_t i : idxList) { - for (idx_t j = 0; j < fieldView.shape(1); ++j) { - for (idx_t k = 0; k < fieldView.shape(2); ++k) { - f(fieldView(i, j, k)); - } - } +// End of recursion when Dim == Rank. +template +struct ForEach { + template + static void apply(const std::vector& idxList, array::ArrayView& fieldView, const Functor& f, + Idxs... idxs) { + // Apply functor. + f(fieldView(idxs...)); } -} +}; } // namespace @@ -297,25 +304,22 @@ void RedistributeGeneric::execute(const FieldSet& sourceFieldSet, FieldSet& targ // Determine datatype. void RedistributeGeneric::do_execute(const Field& sourceField, Field& targetField) const { + // Available datatypes defined in array/LocalView.cc switch (sourceField.datatype().kind()) { case array::DataType::KIND_REAL64: { - do_execute(sourceField, targetField); - break; + return do_execute(sourceField, targetField); } case array::DataType::KIND_REAL32: { - do_execute(sourceField, targetField); - break; + return do_execute(sourceField, targetField); } case array::DataType::KIND_INT64: { - do_execute(sourceField, targetField); - break; + return do_execute(sourceField, targetField); } case array::DataType::KIND_INT32: { - do_execute(sourceField, targetField); - break; + return do_execute(sourceField, targetField); } default: { - throw_NotImplemented("No implementation for data type " + sourceField.datatype().str(), Here()); + ATLAS_THROW_EXCEPTION("No implementation for data type " + sourceField.datatype().str()); } } } @@ -323,21 +327,37 @@ void RedistributeGeneric::do_execute(const Field& sourceField, Field& targetFiel // Determine rank. template void RedistributeGeneric::do_execute(const Field& sourceField, Field& targetField) const { + // Available ranks defined in array/LocalView.cc switch (sourceField.rank()) { case 1: { - do_execute(sourceField, targetField); - break; + return do_execute(sourceField, targetField); } case 2: { - do_execute(sourceField, targetField); - break; + return do_execute(sourceField, targetField); } case 3: { - do_execute(sourceField, targetField); - break; + return do_execute(sourceField, targetField); + } + case 4: { + return do_execute(sourceField, targetField); + } + case 5: { + return do_execute(sourceField, targetField); + } + case 6: { + return do_execute(sourceField, targetField); + } + case 7: { + return do_execute(sourceField, targetField); + } + case 8: { + return do_execute(sourceField, targetField); + } + case 9: { + return do_execute(sourceField, targetField); } default: { - throw_NotImplemented("No implementation for rank " + std::to_string(sourceField.rank()), Here()); + ATLAS_THROW_EXCEPTION("No implementation for rank " + std::to_string(sourceField.rank())); } } } @@ -380,14 +400,14 @@ void RedistributeGeneric::do_execute(const Field& sourceField, Field& targetFiel auto recvBufferIt = recvBuffer.cbegin(); // Copy sourceField to sendBuffer. - iterateField(sourceLocalIdx_, sourceView, [&](const Value& elem) { *sendBufferIt++ = elem; }); + ForEach::apply(sourceLocalIdx_, sourceView, [&](const Value& elem) { *sendBufferIt++ = elem; }); // Perform MPI communication. mpi::comm().allToAllv(sendBuffer.data(), sendCounts.data(), sendDisps.data(), recvBuffer.data(), recvCounts.data(), recvDisps.data()); // Copy recvBuffer to targetField. - iterateField(targetLocalIdx_, targetView, [&](Value& elem) { elem = *recvBufferIt++; }); + ForEach::apply(targetLocalIdx_, targetView, [&](Value& elem) { elem = *recvBufferIt++; }); } namespace { diff --git a/src/atlas/runtime/Exception.cc b/src/atlas/runtime/Exception.cc index 85a4c8a8f..752c2353d 100644 --- a/src/atlas/runtime/Exception.cc +++ b/src/atlas/runtime/Exception.cc @@ -25,10 +25,15 @@ void throw_NotImplemented(const std::string& msg, const eckit::CodeLocation& loc } void throw_AssertionFailed(const std::string& msg) { +#if ATLAS_ECKIT_VERSION_AT_LEAST(1, 18, 5) throw_AssertionFailed(msg, eckit::CodeLocation{}); +#else + throw eckit::AssertionFailed(msg); +#endif } void throw_AssertionFailed(const std::string& msg, const eckit::CodeLocation& loc) { +#if ATLAS_ECKIT_VERSION_AT_LEAST(1, 18, 5) if (loc) { eckit::handle_assert(msg, loc); } @@ -36,6 +41,9 @@ void throw_AssertionFailed(const std::string& msg, const eckit::CodeLocation& lo eckit::handle_assert(msg, eckit::CodeLocation{"unspecified file", 0, "unspecified function"}); } std::abort(); // should never reach here, but makes sure we will never return from this +#else + throw eckit::AssertionFailed(msg, loc); +#endif } void throw_AssertionFailed(const std::string& code, const std::string& msg, const eckit::CodeLocation& loc) { diff --git a/src/atlas/trans/ifs/TransIFS.cc b/src/atlas/trans/ifs/TransIFS.cc index a10be6a09..ba0b874bf 100644 --- a/src/atlas/trans/ifs/TransIFS.cc +++ b/src/atlas/trans/ifs/TransIFS.cc @@ -337,7 +337,7 @@ void TransIFS::invtrans(const int nb_scalar_fields, const double scalar_spectra[ void TransIFS::invtrans_adj(const int nb_scalar_fields, const double gp_fields[], const int nb_vordiv_fields, double vorticity_spectra[], double divergence_spectra[], double scalar_spectra[], const eckit::Configuration& config) const { -#ifdef TRANS_HAVE_INVTRANS_ADJ +#if ATLAS_HAVE_ECTRANS || defined(TRANS_HAVE_INVTRANS_ADJ) ATLAS_TRACE("TransIFS::invtrans_adj"); TransParameters params(*this, config); struct ::InvTransAdj_t args = new_invtrans_adj(trans_.get()); @@ -1344,7 +1344,7 @@ void TransIFS::__invtrans_grad_adj(const Spectral& sp, Field& spfield, const fun void TransIFS::__invtrans_grad_adj(const Spectral& sp, FieldSet& spfields, const functionspace::NodeColumns& gp, const FieldSet& gradfields, const eckit::Configuration& config) const { -#ifdef TRANS_HAVE_INVTRANS_ADJ +#if ATLAS_HAVE_ECTRANS || defined(TRANS_HAVE_INVTRANS_ADJ) assertCompatibleDistributions(gp, sp); // Count total number of fields and do sanity checks @@ -1469,7 +1469,7 @@ void TransIFS::__invtrans(const Spectral& sp, const FieldSet& spfields, const fu void TransIFS::__invtrans_adj(const Spectral& sp, FieldSet& spfields, const functionspace::NodeColumns& gp, const FieldSet& gpfields, const eckit::Configuration& config) const { -#ifdef TRANS_HAVE_INVTRANS_ADJ +#if ATLAS_HAVE_ECTRANS || defined(TRANS_HAVE_INVTRANS_ADJ) assertCompatibleDistributions(gp, sp); @@ -1573,7 +1573,7 @@ void TransIFS::__invtrans(const functionspace::Spectral& sp, const Field& spfiel void TransIFS::__invtrans_adj(const functionspace::Spectral& sp, Field& spfield, const functionspace::StructuredColumns& gp, const Field& gpfield, const eckit::Configuration& config) const { -#ifdef TRANS_HAVE_INVTRANS_ADJ +#if ATLAS_HAVE_ECTRANS || defined(TRANS_HAVE_INVTRANS_ADJ) ATLAS_ASSERT(gpfield.functionspace() == 0 || functionspace::StructuredColumns(gpfield.functionspace())); ATLAS_ASSERT(spfield.functionspace() == 0 || functionspace::Spectral(spfield.functionspace())); @@ -1680,7 +1680,7 @@ void TransIFS::__invtrans(const functionspace::Spectral& sp, const FieldSet& spf void TransIFS::__invtrans_adj(const functionspace::Spectral& sp, FieldSet& spfields, const functionspace::StructuredColumns& gp, const FieldSet& gpfields, const eckit::Configuration& config) const { -#ifdef TRANS_HAVE_INVTRANS_ADJ +#if ATLAS_HAVE_ECTRANS || defined(TRANS_HAVE_INVTRANS_ADJ) assertCompatibleDistributions(gp, sp); @@ -1877,7 +1877,7 @@ void TransIFS::__invtrans_vordiv2wind(const Spectral& sp, const Field& spvor, co void TransIFS::__invtrans_vordiv2wind_adj(const Spectral& sp, Field& spvor, Field& spdiv, const functionspace::NodeColumns& gp, const Field& gpwind, const eckit::Configuration&) const { -#ifdef TRANS_HAVE_INVTRANS_ADJ +#if ATLAS_HAVE_ECTRANS || defined(TRANS_HAVE_INVTRANS_ADJ) assertCompatibleDistributions(gp, sp); diff --git a/src/atlas/util/Config.cc b/src/atlas/util/Config.cc index 1a1c3445b..64955b421 100644 --- a/src/atlas/util/Config.cc +++ b/src/atlas/util/Config.cc @@ -109,6 +109,13 @@ std::vector Config::keys() const { return result; } +std::string Config::json() const { + std::stringstream json; + eckit::JSON js(json); + js << *this; + return json.str(); +} + //================================================================== // ------------------------------------------------------------------ diff --git a/src/atlas/util/Config.h b/src/atlas/util/Config.h index ffb630c65..02170b307 100644 --- a/src/atlas/util/Config.h +++ b/src/atlas/util/Config.h @@ -76,6 +76,8 @@ class Config : public eckit::LocalConfiguration { bool get(const std::string& name, std::vector& value) const; std::vector keys() const; + + std::string json() const; }; // ------------------------------------------------------------------ diff --git a/src/atlas_f/field/atlas_FieldSet_module.F90 b/src/atlas_f/field/atlas_FieldSet_module.F90 index 77f7bddb2..3ecaa4749 100644 --- a/src/atlas_f/field/atlas_FieldSet_module.F90 +++ b/src/atlas_f/field/atlas_FieldSet_module.F90 @@ -41,12 +41,12 @@ module atlas_FieldSet_module !------------------------------------------------------------------------------ contains - procedure, public :: size => FieldSet__size - procedure, public :: has_field + procedure, public :: size => FieldSet__size + procedure, public :: has procedure, private :: field_by_name procedure, private :: field_by_idx_int procedure, private :: field_by_idx_long - procedure, public :: add + procedure, public :: add generic :: field => field_by_name, field_by_idx_int, field_by_idx_long procedure, public :: set_dirty @@ -55,6 +55,7 @@ module atlas_FieldSet_module #if FCKIT_FINAL_NOT_INHERITING final :: atlas_FieldSet__final_auto #endif + procedure, public :: has_field => has ! deprecated ! END TYPE atlas_FieldSet !------------------------------------------------------------------------------ @@ -101,7 +102,7 @@ subroutine add(this,field) call atlas__FieldSet__add_field(this%CPTR_PGIBUG_A, field%CPTR_PGIBUG_A) end subroutine -function has_field(this,name) result(flag) +function has(this,name) result(flag) use, intrinsic :: iso_c_binding, only: c_int use fckit_c_interop_module, only: c_str use atlas_fieldset_c_binding @@ -117,6 +118,7 @@ function has_field(this,name) result(flag) end if end function + function FieldSet__size(this) result(nb_fields) use atlas_fieldset_c_binding class(atlas_FieldSet), intent(in) :: this diff --git a/src/tests/array/test_array_view_util.cc b/src/tests/array/test_array_view_util.cc index ba3034315..b0d2dbcca 100644 --- a/src/tests/array/test_array_view_util.cc +++ b/src/tests/array/test_array_view_util.cc @@ -22,7 +22,7 @@ namespace test { //----------------------------------------------------------------------------- CASE("test_var_size") { - Array* ds = Array::create(4ul, 5ul, 7ul, 9ul); + std::unique_ptr ds{Array::create(4ul, 5ul, 7ul, 9ul)}; auto arrv = make_host_view(*ds); EXPECT(ds->size() == 1260); @@ -31,12 +31,10 @@ CASE("test_var_size") { EXPECT(get_var_size<1>(arrv) == 252); EXPECT(get_var_size<2>(arrv) == 180); EXPECT(get_var_size<3>(arrv) == 140); - - delete ds; } CASE("test_get_parallel_dim") { - Array* ds = Array::create(4ul, 5ul, 7ul, 9ul); + std::unique_ptr ds{Array::create(4ul, 5ul, 7ul, 9ul)}; auto arrv = make_host_view(*ds); EXPECT(ds->size() == 1260); @@ -45,8 +43,14 @@ CASE("test_get_parallel_dim") { EXPECT(get_parallel_dim(arrv) == 3); EXPECT(get_parallel_dim>(arrv) == 1); EXPECT(get_parallel_dim>(arrv) == 2); +} + +CASE("test mixed index types") { + std::unique_ptr ds{Array::create(4ul, 5ul, 7ul, 9ul)}; + + auto arrv = make_host_view(*ds); - delete ds; + arrv(long(0), size_t(0), int(0), unsigned(0)) = 0.; } //----------------------------------------------------------------------------- diff --git a/src/tests/functionspace/test_cubedsphere_functionspace.cc b/src/tests/functionspace/test_cubedsphere_functionspace.cc index 3ce49a710..a1f3e2cf1 100644 --- a/src/tests/functionspace/test_cubedsphere_functionspace.cc +++ b/src/tests/functionspace/test_cubedsphere_functionspace.cc @@ -130,10 +130,8 @@ CASE("cubedsphere_mesh_functionspace") { const auto meshGenCubedSphere = MeshGenerator("cubedsphere", meshConfigCubedSphere); // Set dual mesh generator. - const auto dualMeshGenEqualRegions = - MeshGenerator("cubedsphere_dual", meshConfigEqualRegions | util::Config("halo", 0)); - const auto dualMeshGenCubedSphere = - MeshGenerator("cubedsphere_dual", meshConfigCubedSphere | util::Config("halo", 0)); + const auto dualMeshGenEqualRegions = MeshGenerator("cubedsphere_dual", meshConfigEqualRegions); + const auto dualMeshGenCubedSphere = MeshGenerator("cubedsphere_dual", meshConfigCubedSphere); // Set mesh const auto meshEqualRegions = meshGenEqualRegions.generate(grid); @@ -159,6 +157,8 @@ CASE("cubedsphere_mesh_functionspace") { SECTION("CellColumns: cubedsphere") { testFunctionSpace(cubedSphereCellColumns); } SECTION("NodeColumns: equal_regions") { testFunctionSpace(equalRegionsNodeColumns); } SECTION("NodeColumns: cubedsphere") { testFunctionSpace(cubedSphereNodeColumns); } + SECTION("CellColumns: dual mesh, equal_regions") { testFunctionSpace(equalRegionsDualCellColumns); } + SECTION("CellColumns: dual mesh, cubedsphere") { testFunctionSpace(cubedSphereDualCellColumns); } SECTION("NodeColumns: dual mesh, equal_regions") { testFunctionSpace(equalRegionsDualNodeColumns); } SECTION("NodeColumns: dual mesh, cubedsphere") { testFunctionSpace(cubedSphereDualNodeColumns); } } @@ -211,22 +211,22 @@ CASE("Cubed sphere primal-dual equivalence") { // Create domain decomposed mesh. // We can generate two types of mesh, a primal and a dual. - // The primal mesh can have an arbitarily large halo (default = 1), whereas - // the dual mesh has a fixed halo size. The dual mesh is useful if you - // are interested in node connectivity, otherwise the primal mesh is the - // most general. // The cell-centre (nodes) lonlats of the primal mesh are identical to the // node (cell-centre) lonlats of the dual mesh. - const auto primalMesh = MeshGenerator("cubedsphere").generate(grid); - const auto dualMesh = MeshGenerator("cubedsphere_dual").generate(grid); + const auto primalMesh = + MeshGenerator("cubedsphere", util::Config("halo", 5) | util::Config("partitioner", "equal_regions")) + .generate(grid); + const auto dualMesh = + MeshGenerator("cubedsphere_dual", util::Config("halo", 4) | util::Config("partitioner", "equal_regions")) + .generate(grid); // Create cubed sphere function spaces (these have fancy features, such as // (t, i, j) indexing and parallel_for methods). The halo sizes of the primal // functionspaces are set to match that of the dual functionspaces. - const auto primalNodes = functionspace::CubedSphereNodeColumns(primalMesh, util::Config("halo", 0)); - const auto primalCells = functionspace::CubedSphereCellColumns(primalMesh, util::Config("halo", 1)); - const auto dualNodes = functionspace::CubedSphereNodeColumns(dualMesh); - const auto dualCells = functionspace::CubedSphereCellColumns(dualMesh); + const auto primalNodes = functionspace::CubedSphereNodeColumns(primalMesh, util::Config("halo", 4)); + const auto primalCells = functionspace::CubedSphereCellColumns(primalMesh, util::Config("halo", 5)); + const auto dualNodes = functionspace::CubedSphereNodeColumns(dualMesh, util::Config("halo", 4)); + const auto dualCells = functionspace::CubedSphereCellColumns(dualMesh, util::Config("halo", 4)); // Note, the functionspaces we are usually interested in are primalCells and // dualNodes. The others are there for completeness. @@ -257,7 +257,7 @@ CASE("Cubed sphere primal-dual equivalence") { compareFields(dualCells, primalNodes); } -CASE("Variable halo size functionspaces") { +CASE("Variable halo size functionspaces (primal mesh)") { // Create a mesh with a large halo, and a few functionspaces with different // (smaller) halo sizes. These should create fields with a smaller memory // footprint. @@ -299,6 +299,48 @@ CASE("Variable halo size functionspaces") { checkSize(cellColumns2.sizeOwned()); } +CASE("Variable halo size functionspaces (dual mesh)") { + // Create a mesh with a large halo, and a few functionspaces with different + // (smaller) halo sizes. These should create fields with a smaller memory + // footprint. + + // Set grid. + const auto grid = Grid("CS-LFR-C-12"); + + // Set mesh config. + const auto meshConfig = util::Config("partitioner", "equal_regions") | util::Config("halo", 3); + + // Set mesh. + const auto mesh = MeshGenerator("cubedsphere_dual", meshConfig).generate(grid); + + // Set functionspaces. + const auto nodeColumns0 = functionspace::CubedSphereNodeColumns(mesh, util::Config("halo", 0)); + const auto nodeColumns1 = functionspace::CubedSphereNodeColumns(mesh, util::Config("halo", 1)); + const auto nodeColumns2 = functionspace::CubedSphereNodeColumns(mesh, util::Config("halo", 2)); + + const auto cellColumns0 = functionspace::CubedSphereCellColumns(mesh, util::Config("halo", 0)); + const auto cellColumns1 = functionspace::CubedSphereCellColumns(mesh, util::Config("halo", 1)); + const auto cellColumns2 = functionspace::CubedSphereCellColumns(mesh, util::Config("halo", 2)); + + // Check functionspace sizes. + EXPECT(nodeColumns0.size() < nodeColumns1.size()); + EXPECT(nodeColumns1.size() < nodeColumns2.size()); + EXPECT(nodeColumns2.size() < mesh.nodes().size()); + EXPECT(cellColumns0.size() < cellColumns1.size()); + EXPECT(cellColumns1.size() < cellColumns2.size()); + EXPECT(cellColumns2.size() < mesh.cells().size()); + + // Make sure size of owned cell data matches grid. + auto checkSize = [&](idx_t sizeOwned) { + mpi::comm().allReduceInPlace(sizeOwned, eckit::mpi::Operation::SUM); + EXPECT_EQ(sizeOwned, grid.size()); + }; + + checkSize(nodeColumns0.sizeOwned()); + checkSize(nodeColumns1.sizeOwned()); + checkSize(nodeColumns2.sizeOwned()); +} + } // namespace test } // namespace atlas diff --git a/src/tests/interpolation/CMakeLists.txt b/src/tests/interpolation/CMakeLists.txt index c716e571b..11b9e85be 100644 --- a/src/tests/interpolation/CMakeLists.txt +++ b/src/tests/interpolation/CMakeLists.txt @@ -14,14 +14,12 @@ ecbuild_add_test( TARGET atlas_test_Quad3D ) ecbuild_add_test( TARGET atlas_test_Quad2D - CONDITION eckit_EIGEN_FOUND SOURCES test_Quad2D.cc LIBS atlas ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} ) ecbuild_add_test( TARGET atlas_test_Triag2D - CONDITION eckit_EIGEN_FOUND SOURCES test_Triag2D.cc LIBS atlas ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} @@ -33,8 +31,8 @@ ecbuild_add_test( TARGET atlas_test_interpolation_finite_element ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} ) -ecbuild_add_test( TARGET atlas_test_interpolation_bilinear_remapping - SOURCES test_interpolation_bilinear_remapping.cc +ecbuild_add_test( TARGET atlas_test_interpolation_unstructured_bilinear_lonlat + SOURCES test_interpolation_unstructured_bilinear_lonlat.cc LIBS atlas ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} ) @@ -97,3 +95,11 @@ ecbuild_add_test( TARGET atlas_test_interpolation_structured2D_to_unstructured CONDITION eckit_HAVE_MPI ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} ) + +ecbuild_add_test( TARGET atlas_test_interpolation_cubedsphere + SOURCES test_interpolation_cubedsphere.cc + LIBS atlas + MPI 6 + CONDITION eckit_HAVE_MPI AND MPI_SLOTS GREATER_EQUAL 6 + ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} +) diff --git a/src/tests/interpolation/test_Quad2D.cc b/src/tests/interpolation/test_Quad2D.cc index b9e381803..29e43f404 100644 --- a/src/tests/interpolation/test_Quad2D.cc +++ b/src/tests/interpolation/test_Quad2D.cc @@ -198,6 +198,54 @@ CASE("test_quadrilateral_intersection_corners") { } } +CASE("test_quadrilateral_intersection_arbitrary_corners") { + PointXY v0(338.14, 54.6923); + PointXY v1(340.273, 54.6778); + PointXY v2(340.312, 55.9707); + PointXY v3(338.155, 55.9852); + + Quad2D quad(v0.data(), v1.data(), v2.data(), v3.data()); + + EXPECT(quad.validate()); + + std::array orig{v0, v1, v2, v3}; + std::array uTarget{0., 1., 1., 0.}; + std::array vTarget{0., 0., 1., 1.}; + + for (size_t i = 0; i < 4; ++i) { + Intersect isect = quad.localRemap(orig[i]); + + EXPECT(isect); + EXPECT_APPROX_EQ(isect.u, uTarget[i], relative_error); + EXPECT_APPROX_EQ(isect.v, vTarget[i], relative_error); + } +} + +CASE("test_quadrilateral_intersection_arbitrary_edges") { + PointXY v0(338.14, 54.6923); + PointXY v1(340.273, 54.6778); + PointXY v2(340.312, 55.9707); + PointXY v3(338.155, 55.9852); + + Quad2D quad(v0.data(), v1.data(), v2.data(), v3.data()); + + EXPECT(quad.validate()); + + std::array orig{(v0 + v1) * 0.5, (v1 + v2) * 0.5, (v2 + v3) * 0.5, (v3 + v0) * 0.5}; + + std::array uTarget{0.5, 1., 0.5, 0.}; + std::array vTarget{0., 0.5, 1., 0.5}; + + for (size_t i = 0; i < 4; ++i) { + Intersect isect = quad.localRemap(orig[i]); + + EXPECT(isect); + EXPECT_APPROX_EQ(isect.u, uTarget[i], relative_error); + EXPECT_APPROX_EQ(isect.v, vTarget[i], relative_error); + } +} + + //----------------------------------------------------------------------------- } // namespace test diff --git a/src/tests/interpolation/test_interpolation_cubedsphere.cc b/src/tests/interpolation/test_interpolation_cubedsphere.cc new file mode 100644 index 000000000..04c3b1cb5 --- /dev/null +++ b/src/tests/interpolation/test_interpolation_cubedsphere.cc @@ -0,0 +1,389 @@ +/* + * (C) Crown Copyright 2022 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/field/FieldSet.h" +#include "atlas/functionspace/CellColumns.h" +#include "atlas/functionspace/CubedSphereColumns.h" +#include "atlas/functionspace/NodeColumns.h" +#include "atlas/grid/CubedSphereGrid.h" +#include "atlas/grid/Grid.h" +#include "atlas/grid/Partitioner.h" +#include "atlas/interpolation/Interpolation.h" +#include "atlas/mesh/Mesh.h" +#include "atlas/meshgenerator/MeshGenerator.h" +#include "atlas/output/Gmsh.h" +#include "atlas/parallel/mpi/mpi.h" +#include "atlas/util/Constants.h" +#include "atlas/util/CoordinateEnums.h" +#include "atlas/util/function/VortexRollup.h" + +#include "tests/AtlasTestEnvironment.h" + + +namespace atlas { +namespace test { + + +// 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) { + // set hLon and hLat step size. + const double hLon = 0.0001; + const double hLat = 0.0001; + + // Get finite differences. + + // Set u. + const double u = (util::function::vortex_rollup(lon, lat + 0.5 * hLat, 1.) - + util::function::vortex_rollup(lon, lat - 0.5 * hLat, 1.)) / + hLat; + + const double v = -(util::function::vortex_rollup(lon + 0.5 * hLon, lat, 1.) - + util::function::vortex_rollup(lon - 0.5 * hLon, lat, 1.)) / + (hLon * std::cos(lat * util::Constants::degreesToRadians())); + + return std::make_pair(u, v); +} + + +double dotProd(const Field& a, const Field& b) { + double prod{}; + + const auto aView = array::make_view(a); + const auto bView = array::make_view(b); + + + for (size_t i = 0; i < a.size(); ++i) { + prod += aView(i) * bView(i); + } + mpi::comm().allReduceInPlace(prod, eckit::mpi::Operation::SUM); + return prod; +} + +CASE("cubedsphere_scalar_interpolation") { + // Create a source cubed sphere grid, mesh and functionspace. + const auto sourceGrid = Grid("CS-LFR-24"); + const auto sourceMesh = MeshGenerator("cubedsphere_dual").generate(sourceGrid); + const auto sourceFunctionspace = functionspace::NodeColumns(sourceMesh); + + //-------------------------------------------------------------------------- + // Interpolation test. + //-------------------------------------------------------------------------- + + // Populate analytic source field. + double stDev{}; + auto sourceField = sourceFunctionspace.createField(option::name("test_field")); + { + const auto lonlat = array::make_view(sourceFunctionspace.lonlat()); + const auto ghost = array::make_view(sourceFunctionspace.ghost()); + auto view = array::make_view(sourceField); + for (idx_t i = 0; i < sourceFunctionspace.size(); ++i) { + view(i) = util::function::vortex_rollup(lonlat(i, LON), lonlat(i, LAT), 1.); + if (!ghost(i)) { + stDev += view(i) * view(i); + } + } + } + mpi::comm().allReduceInPlace(stDev, eckit::mpi::Operation::SUM); + stDev = std::sqrt(stDev / sourceGrid.size()); + + + // Create target grid, mesh and functionspace. + const auto partitioner = grid::MatchingPartitioner(sourceMesh, util::Config("type", "cubedsphere")); + const auto targetGrid = Grid("O24"); + const auto targetMesh = MeshGenerator("structured").generate(targetGrid, partitioner); + const auto targetFunctionspace = functionspace::NodeColumns(targetMesh); + + // Set up interpolation object. + const auto scheme = util::Config("type", "cubedsphere-bilinear") | util::Config("adjoint", true); + const auto interp = Interpolation(scheme, sourceFunctionspace, targetFunctionspace); + + // Interpolate from source to target field. + auto targetField = targetFunctionspace.createField(option::name("test_field")); + interp.execute(sourceField, targetField); + targetField.haloExchange(); + + // Make some diagnostic output fields. + auto errorField = targetFunctionspace.createField(option::name("error_field")); + auto partField = targetFunctionspace.createField(option::name("partition")); + { + const auto lonlat = array::make_view(targetFunctionspace.lonlat()); + auto targetView = array::make_view(targetField); + auto errorView = array::make_view(errorField); + auto partView = array::make_view(partField); + for (idx_t i = 0; i < targetFunctionspace.size(); ++i) { + const auto val = util::function::vortex_rollup(lonlat(i, LON), lonlat(i, LAT), 1.); + errorView(i) = std::abs((targetView(i) - val) / stDev); + partView(i) = mpi::rank(); + } + } + partField.haloExchange(); + + // Output source mesh. + const auto gmshConfig = + util::Config("coordinates", "xyz") | util::Config("ghost", true) | util::Config("info", true); + const auto sourceGmsh = output::Gmsh("cubedsphere_source.msh", gmshConfig); + sourceGmsh.write(sourceMesh); + sourceGmsh.write(FieldSet(sourceField), sourceFunctionspace); + + // Output target mesh. + const auto targetGmsh = output::Gmsh("cubedsphere_target.msh", gmshConfig); + targetGmsh.write(targetMesh); + auto targetFields = FieldSet{}; + targetFields.add(targetField); + targetFields.add(errorField); + targetFields.add(partField); + targetGmsh.write(targetFields, targetFunctionspace); + + //-------------------------------------------------------------------------- + // Adjoint test. + //-------------------------------------------------------------------------- + + // Ensure that the adjoint identity relationship holds. + auto targetAdjoint = targetFunctionspace.createField(option::name("target adjoint")); + array::make_view(targetAdjoint).assign(array::make_view(targetField)); + targetAdjoint.adjointHaloExchange(); + + auto sourceAdjoint = sourceFunctionspace.createField(option::name("source adjoint")); + array::make_view(sourceAdjoint).assign(0.); + interp.execute_adjoint(sourceAdjoint, targetAdjoint); + + const auto yDotY = dotProd(targetField, targetField); + const auto xDotXAdj = dotProd(sourceField, sourceAdjoint); + + EXPECT_APPROX_EQ(yDotY / xDotXAdj, 1., 1e-14); +} + +CASE("cubedsphere_wind_interpolation") { + // Create a source cubed sphere grid, mesh and functionspace. + const auto sourceGrid = CubedSphereGrid("CS-LFR-48"); + const auto sourceMesh = MeshGenerator("cubedsphere_dual").generate(sourceGrid); + const auto sourceFunctionspace = functionspace::CubedSphereNodeColumns(sourceMesh); + + // Get projection. + const auto& proj = sourceGrid.cubedSphereProjection(); + + // Set wind transform Jacobian. + const auto windTransform = [&](const PointLonLat& lonlat, idx_t t) { + // Jacobian = [d(alpha, beta) / d(lon, lat)] * [d(lon, lat) / d(u_hat, v_hat)]; + // Note: u_hat and v_hat are increments along the u and v unit vectors. + + const auto secLat = 1. / std::cos(util::Constants::degreesToRadians() * lonlat.lat()); + return proj.alphabetaJacobian(lonlat, t) * projection::Jacobian{{secLat, 0.}, {0., 1.}}; + }; + + // Matrix-vector multiplication helper. y = Ax. + const auto matMul = [](const projection::Jacobian& A, double x0, double x1) { + const Point2 y = A * Point2{x0, x1}; + return std::make_pair(y[0], y[1]); + }; + + //-------------------------------------------------------------------------- + // Interpolation test. + //-------------------------------------------------------------------------- + + // Populate analytic source field. + auto sourceFieldSet = FieldSet{}; + sourceFieldSet.add(sourceFunctionspace.createField(option::name("u_orig"))); + sourceFieldSet.add(sourceFunctionspace.createField(option::name("v_orig"))); + sourceFieldSet.add(sourceFunctionspace.createField(option::name("v_alpha"))); + sourceFieldSet.add(sourceFunctionspace.createField(option::name("v_beta"))); + { + const auto lonlat = array::make_view(sourceFunctionspace.lonlat()); + auto u = array::make_view(sourceFieldSet["u_orig"]); + auto v = array::make_view(sourceFieldSet["v_orig"]); + auto vAlpha = array::make_view(sourceFieldSet["v_alpha"]); + auto vBeta = array::make_view(sourceFieldSet["v_beta"]); + + + // In practice, the (u, v) field needs to be halo exchanged *before* the + // wind transform. Then the transform is applied to the entire field, + // *including* the halo. + + sourceFunctionspace.parallel_for(util::Config("include_halo", true), [&](idx_t idx, idx_t t, idx_t i, idx_t j) { + // Get lonlat + 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()); + + // Get wind transform jacobian. + auto jac = windTransform(ll, t); + + // Transform wind. + std::tie(vAlpha(idx), vBeta(idx)) = matMul(jac, u(idx), v(idx)); + }); + } + + // Create target grid, mesh and functionspace. + const auto partitioner = grid::MatchingPartitioner(sourceMesh, util::Config("type", "cubedsphere")); + const auto targetGrid = Grid("O48"); + const auto targetMesh = MeshGenerator("structured").generate(targetGrid, partitioner); + const auto targetFunctionspace = functionspace::NodeColumns(targetMesh); + + // Set up interpolation object. + // Note: We have to disable the source field halo exhange in the + // interpolation execute and excute_adjoint methods. If left on, the halo + // exchange will corrupt the transformed wind field. + const auto scheme = util::Config("type", "cubedsphere-bilinear") | util::Config("adjoint", true) | + util::Config("halo_exchange", false); + const auto interp = Interpolation(scheme, sourceFunctionspace, targetFunctionspace); + + // Make target fields. + auto targetFieldSet = FieldSet{}; + targetFieldSet.add(targetFunctionspace.createField(option::name("u_orig"))); + targetFieldSet.add(targetFunctionspace.createField(option::name("v_orig"))); + targetFieldSet.add(targetFunctionspace.createField(option::name("v_alpha"))); + targetFieldSet.add(targetFunctionspace.createField(option::name("v_beta"))); + + // Interpolate from source to target fields. + array::make_view(targetFieldSet["v_alpha"]).assign(0.); + array::make_view(targetFieldSet["v_beta"]).assign(0.); + interp.execute(sourceFieldSet, targetFieldSet); + + // Make new (u, v) fields from (v_alpha, v_beta) + targetFieldSet.add(targetFunctionspace.createField(option::name("error_field_0"))); + targetFieldSet.add(targetFunctionspace.createField(option::name("error_field_1"))); + targetFieldSet.add(targetFunctionspace.createField(option::name("u_new"))); + targetFieldSet.add(targetFunctionspace.createField(option::name("v_new"))); + { + const auto lonlat = array::make_view(targetFunctionspace.lonlat()); + const auto ghost = array::make_view(targetFunctionspace.ghost()); + auto u = array::make_view(targetFieldSet["u_new"]); + auto v = array::make_view(targetFieldSet["v_new"]); + const auto uOrig = array::make_view(targetFieldSet["u_orig"]); + const auto vOrig = array::make_view(targetFieldSet["v_orig"]); + const auto vAlpha = array::make_view(targetFieldSet["v_alpha"]); + const auto vBeta = array::make_view(targetFieldSet["v_beta"]); + auto error0 = array::make_view(targetFieldSet["error_field_0"]); + auto error1 = array::make_view(targetFieldSet["error_field_1"]); + for (idx_t idx = 0; idx < targetFunctionspace.size(); ++idx) { + if (!ghost(idx)) { + const auto ll = PointLonLat(lonlat(idx, LON), lonlat(idx, LAT)); + const idx_t t = proj.getCubedSphereTiles().indexFromLonLat(ll); + + // Get inverse wind transform jacobian. + auto jac = windTransform(ll, t).inverse(); + + // Transform wind. + std::tie(u(idx), v(idx)) = matMul(jac, vAlpha(idx), vBeta(idx)); + + // Get error. + const auto uvTarget = vortexField(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))); + } + } + } + targetFieldSet.haloExchange(); + + // Output source mesh. + const auto gmshConfig = + util::Config("coordinates", "xyz") | util::Config("ghost", true) | util::Config("info", true); + const auto sourceGmsh = output::Gmsh("cubedsphere_vec_source.msh", gmshConfig); + sourceGmsh.write(sourceMesh); + sourceGmsh.write(sourceFieldSet, sourceFunctionspace); + + // Output target mesh. + const auto targetGmsh = output::Gmsh("cubedsphere_vec_target.msh", gmshConfig); + targetGmsh.write(targetMesh); + targetGmsh.write(targetFieldSet, targetFunctionspace); + + + //-------------------------------------------------------------------------- + // Adjoint test. + //-------------------------------------------------------------------------- + + // Ensure that the adjoint identity relationship holds. + targetFieldSet.add(targetFunctionspace.createField(option::name("u_adjoint"))); + targetFieldSet.add(targetFunctionspace.createField(option::name("v_adjoint"))); + array::make_view(targetFieldSet["u_adjoint"]) + .assign(array::make_view(targetFieldSet["u_new"])); + array::make_view(targetFieldSet["v_adjoint"]) + .assign(array::make_view(targetFieldSet["v_new"])); + + // Adjoint of target halo exhange. + targetFieldSet["u_adjoint"].adjointHaloExchange(); + targetFieldSet["v_adjoint"].adjointHaloExchange(); + + // Adjoint of inverse wind transform. + targetFieldSet.add(targetFunctionspace.createField(option::name("v_alpha_adjoint"))); + targetFieldSet.add(targetFunctionspace.createField(option::name("v_beta_adjoint"))); + { + const auto lonlat = array::make_view(targetFunctionspace.lonlat()); + const auto ghost = array::make_view(targetFunctionspace.ghost()); + const auto uAdj = array::make_view(targetFieldSet["u_adjoint"]); + const auto vAdj = array::make_view(targetFieldSet["v_adjoint"]); + auto vAlphaAdj = array::make_view(targetFieldSet["v_alpha_adjoint"]); + auto vBetaAdj = array::make_view(targetFieldSet["v_beta_adjoint"]); + vAlphaAdj.assign(0.); + vBetaAdj.assign(0.); + for (idx_t idx = 0; idx < targetFunctionspace.size(); ++idx) { + if (!ghost(idx)) { + const auto ll = PointLonLat(lonlat(idx, LON), lonlat(idx, LAT)); + const idx_t t = proj.getCubedSphereTiles().indexFromLonLat(ll); + + // Get adjoint of inverse wind transform jacobian. + auto jac = windTransform(ll, t).inverse().transpose(); + + // Transform wind. + std::tie(vAlphaAdj(idx), vBetaAdj(idx)) = matMul(jac, uAdj(idx), vAdj(idx)); + } + } + } + + // Adjoint of interpolation. + sourceFieldSet.add(sourceFunctionspace.createField(option::name("v_alpha_adjoint"))); + sourceFieldSet.add(sourceFunctionspace.createField(option::name("v_beta_adjoint"))); + array::make_view(sourceFieldSet["v_alpha_adjoint"]).assign(0.); + array::make_view(sourceFieldSet["v_beta_adjoint"]).assign(0.); + interp.execute_adjoint(sourceFieldSet["v_alpha_adjoint"], targetFieldSet["v_alpha_adjoint"]); + interp.execute_adjoint(sourceFieldSet["v_beta_adjoint"], targetFieldSet["v_beta_adjoint"]); + + // Adjoint of wind transform. + sourceFieldSet.add(sourceFunctionspace.createField(option::name("u_adjoint"))); + sourceFieldSet.add(sourceFunctionspace.createField(option::name("v_adjoint"))); + { + const auto lonlat = array::make_view(sourceFunctionspace.lonlat()); + auto uAdj = array::make_view(sourceFieldSet["u_adjoint"]); + auto vAdj = array::make_view(sourceFieldSet["v_adjoint"]); + uAdj.assign(0.); + vAdj.assign(0.); + const auto vAlphaAdj = array::make_view(sourceFieldSet["v_alpha_adjoint"]); + const auto vBetaAdj = array::make_view(sourceFieldSet["v_beta_adjoint"]); + + sourceFunctionspace.parallel_for(util::Config("include_halo", true), [&](idx_t idx, idx_t t, idx_t i, idx_t j) { + // Get lonlat + const auto ll = PointLonLat(lonlat(idx, LON), lonlat(idx, LAT)); + + // Get adjoint of wind transform jacobian. + auto jac = windTransform(ll, t).transpose(); + + // Transform wind. + std::tie(uAdj(idx), vAdj(idx)) = matMul(jac, vAlphaAdj(idx), vBetaAdj(idx)); + }); + } + + // Check dot products. + const auto yDotY = dotProd(targetFieldSet["u_new"], targetFieldSet["u_new"]) + + dotProd(targetFieldSet["v_new"], targetFieldSet["v_new"]); + + const auto xDotXAdj = dotProd(sourceFieldSet["u_orig"], sourceFieldSet["u_adjoint"]) + + dotProd(sourceFieldSet["v_orig"], sourceFieldSet["v_adjoint"]); + + EXPECT_APPROX_EQ(yDotY / xDotXAdj, 1., 1e-14); +} + + +} // namespace test +} // namespace atlas + +int main(int argc, char** argv) { + return atlas::test::run(argc, argv); +} diff --git a/src/tests/interpolation/test_interpolation_non_linear.cc b/src/tests/interpolation/test_interpolation_non_linear.cc index 8fbe81287..eb560ee67 100644 --- a/src/tests/interpolation/test_interpolation_non_linear.cc +++ b/src/tests/interpolation/test_interpolation_non_linear.cc @@ -136,6 +136,74 @@ CASE("Interpolation with MissingValue") { } } +CASE("Interpolation of rank 2 field with MissingValue") { + RectangularDomain domain({0, 2}, {0, 2}, "degrees"); + Grid gridA("L90", domain); + + Mesh meshA = MeshGenerator("structured").generate(gridA); + + int nlevels = 3; + functionspace::NodeColumns fsA(meshA); + Field fieldA = fsA.createField(option::name("A") | option::levels(nlevels)); + + fieldA.metadata().set("missing_value", missingValue); + fieldA.metadata().set("missing_value_epsilon", missingValueEps); + + auto viewA = array::make_view(fieldA); + for (idx_t j = 0; j < viewA.shape(0); ++j) { + viewA(j, 0) = 10 + j; + viewA(j, 1) = missingValue; + viewA(j, 2) = 30 + j; + } + + const array::ArraySpec spec(array::ArrayShape{fieldA.shape(0)}, array::ArrayStrides{fieldA.shape(1)}); + + // Set output field (2 points) + functionspace::PointCloud fsB({PointLonLat{0.1, 0.1}, PointLonLat{0.9, 0.9}}); + + SECTION("check wrapped data can be indexed with []") { + for (int lev = 0; lev < nlevels; ++lev) { + double* data = const_cast(fieldA.array().data()) + lev; + atlas::Field fieldA_slice = atlas::Field("s", data, spec); + fieldA_slice.metadata().set("missing_value", fieldA.metadata().get("missing_value")); + + auto FlatViewA1 = array::make_view(fieldA_slice); + auto viewA_slice = viewA.slice(array::Range::all(), lev); + EXPECT(FlatViewA1.shape(0) == viewA_slice.shape(0)); + for (idx_t j = 0; j < FlatViewA1.shape(0); ++j) { + EXPECT(FlatViewA1[j] == viewA_slice(j)); + } + } + } + + SECTION("missing-if-all-missing") { + Interpolation interpolation(Config("type", "finite-element").set("non_linear", "missing-if-all-missing"), fsA, + fsB); + + for (std::string type : {"equals", "approximately-equals", "nan"}) { + fieldA.metadata().set("missing_value_type", type); + for (idx_t j = 0; j < viewA.shape(0); ++j) { + viewA(j, 1) = type == "nan" ? nan : missingValue; + } + + EXPECT(MissingValue(fieldA)); + + Field fieldB = fsB.createField(option::name("B") | option::levels(nlevels)); + auto viewB = array::make_view(fieldB); + + interpolation.execute(fieldA, fieldB); + + MissingValue mv(fieldB); + EXPECT(mv(viewB(0, 0)) == false); + EXPECT(mv(viewB(1, 0)) == false); + EXPECT(mv(viewB(0, 1))); + EXPECT(mv(viewB(1, 1))); + EXPECT(mv(viewB(0, 2)) == false); + EXPECT(mv(viewB(1, 2)) == false); + } + } +} + } // namespace test } // namespace atlas diff --git a/src/tests/interpolation/test_interpolation_bilinear_remapping.cc b/src/tests/interpolation/test_interpolation_unstructured_bilinear_lonlat.cc similarity index 93% rename from src/tests/interpolation/test_interpolation_bilinear_remapping.cc rename to src/tests/interpolation/test_interpolation_unstructured_bilinear_lonlat.cc index 8f18c9058..dcd954e3b 100644 --- a/src/tests/interpolation/test_interpolation_bilinear_remapping.cc +++ b/src/tests/interpolation/test_interpolation_unstructured_bilinear_lonlat.cc @@ -44,8 +44,8 @@ CASE("test_interpolation_O64_to_points_bilinear_remapping") { auto func = [](double x) -> double { return std::sin(x * M_PI / 180.); }; - Interpolation interpolation(option::type("bilinear-remapping") | util::Config("max_fraction_elems_to_try", 0.4), fs, - pointcloud); + Interpolation interpolation( + option::type("unstructured-bilinear-lonlat") | util::Config("max_fraction_elems_to_try", 0.4), fs, pointcloud); SECTION("test maximum nearest neighbour settings") { std::stringstream test_stream; @@ -88,8 +88,8 @@ CASE("test_interpolation_N64_to_O32_bilinear_remapping") { Mesh mesh_tgt(grid_tgt); NodeColumns fs_tgt(mesh_tgt); - Interpolation interpolation(option::type("bilinear-remapping") | util::Config("max_fraction_elems_to_try", 0.4), - fs_src, fs_tgt); + Interpolation interpolation( + option::type("unstructured-bilinear-lonlat") | util::Config("max_fraction_elems_to_try", 0.4), fs_src, fs_tgt); Field field_source = fs_src.createField(option::name("source")); Field field_target = fs_tgt.createField(option::name("target")); diff --git a/src/tests/mesh/test_cubedsphere_meshgen.cc b/src/tests/mesh/test_cubedsphere_meshgen.cc index fef7770ad..2f38d73b6 100644 --- a/src/tests/mesh/test_cubedsphere_meshgen.cc +++ b/src/tests/mesh/test_cubedsphere_meshgen.cc @@ -341,8 +341,9 @@ CASE("cubedsphere_dual_mesh_test") { const auto targetPartitioner = grid::MatchingPartitioner(sourceMesh); // Set target grid, mesh and functionspace. - const auto targetGrid = Grid("CS-LFR-C-48"); - const auto targetMesh = MeshGenerator("cubedsphere_dual").generate(targetGrid, targetPartitioner); + const auto targetGrid = Grid("CS-LFR-C-24"); + const auto targetMesh = + MeshGenerator("cubedsphere_dual", util::Config("halo", 3)).generate(targetGrid, targetPartitioner); const auto targetFunctionSpace = functionspace::NodeColumns(targetMesh); auto targetField = targetFunctionSpace.createField(util::Config("name", "targetField") | util::Config("levels", 5)); @@ -393,9 +394,8 @@ CASE("cubedsphere_dual_mesh_test") { const idx_t nLevels = 5; // Set grid, mesh and functionspace. - const auto grid = Grid("CS-LFR-C-48"); - const auto mesh = - MeshGenerator("cubedsphere_dual", util::Config("partitioner", "equal_regions")).generate(grid); + const auto grid = Grid("CS-LFR-C-24"); + const auto mesh = MeshGenerator("cubedsphere_dual", util::Config("halo", 3)).generate(grid); const auto functionSpace = functionspace::CellColumns(mesh); auto field = functionSpace.createField(util::Config("name", "targetField") | util::Config("levels", nLevels)); @@ -432,6 +432,9 @@ CASE("cubedsphere_dual_mesh_test") { } // gmsh output. + auto fieldSet = FieldSet{field}; + fieldSet.add(mesh.cells().halo()); + const auto gmshConfigXy = util::Config("coordinates", "xy") | util::Config("ghost", true) | util::Config("info", true); const auto gmshConfigXyz = @@ -440,8 +443,8 @@ CASE("cubedsphere_dual_mesh_test") { auto gmshXyz = output::Gmsh("dual_cells_xyz.msh", gmshConfigXyz); gmshXy.write(mesh); gmshXyz.write(mesh); - gmshXy.write(FieldSet{field}, functionSpace); - gmshXyz.write(FieldSet{field}, functionSpace); + gmshXy.write(fieldSet, functionSpace); + gmshXyz.write(fieldSet, functionSpace); } } diff --git a/src/tests/projection/CMakeLists.txt b/src/tests/projection/CMakeLists.txt index b44a1ef90..c2c744282 100644 --- a/src/tests/projection/CMakeLists.txt +++ b/src/tests/projection/CMakeLists.txt @@ -9,7 +9,6 @@ foreach(test test_bounding_box test_projection_LAEA - test_projection_cubed_sphere test_projection_variable_resolution test_rotation ) @@ -28,3 +27,11 @@ endif() if( HAVE_FCTEST ) add_fctest( TARGET atlas_fctest_projection SOURCES fctest_projection.F90 LINKER_LANGUAGE Fortran LIBS atlas_f ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} ) endif() + +ecbuild_add_test( TARGET atlas_test_cubedsphere_projection + MPI 6 + CONDITION eckit_HAVE_MPI AND MPI_SLOTS GREATER_EQUAL 6 + SOURCES test_cubedsphere_projection.cc + LIBS atlas + ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} +) diff --git a/src/tests/projection/test_cubedsphere_projection.cc b/src/tests/projection/test_cubedsphere_projection.cc new file mode 100644 index 000000000..24b9b405b --- /dev/null +++ b/src/tests/projection/test_cubedsphere_projection.cc @@ -0,0 +1,282 @@ +/* + * (C) Crown Copyright 2022 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/array/MakeView.h" +#include "atlas/field/FieldSet.h" +#include "atlas/functionspace/NodeColumns.h" +#include "atlas/grid.h" +#include "atlas/grid/CubedSphereGrid.h" +#include "atlas/mesh.h" +#include "atlas/meshgenerator.h" +#include "atlas/output/Gmsh.h" +#include "atlas/projection/detail/CubedSphereProjectionBase.h" +#include "atlas/util/Constants.h" +#include "tests/AtlasTestEnvironment.h" +namespace atlas { +namespace test { + +// Small number relative to 360. +constexpr double epsilon = std::numeric_limits::epsilon() * 360.; + +void testProjection(const std::string& gridType, const std::string& meshType, const std::string& outputID) { + // Create grid. + const auto grid = CubedSphereGrid(gridType); + + // Get projection. + const auto& csProjection = grid.cubedSphereProjection(); + + // Output tile centres and Jacobains. + for (size_t i = 0; i < 6; ++i) { + Log::info() << "Tile " << i << std::endl + << "Centre:" << std::endl + << csProjection.getCubedSphereTiles().tileCentre(i) << std::endl + << "Jacobian:" << std::endl + << csProjection.getCubedSphereTiles().tileJacobian(i) << std::endl + << std::endl; + } + + // Create mesh. + auto mesh = MeshGenerator(meshType, util::Config("halo", 3)).generate(grid); + + // Output mesh. + const auto gmshConfig = + util::Config("coordinates", "xy") | util::Config("ghost", true) | util::Config("info", true); + auto gmsh = output::Gmsh(outputID + "_before.msh", gmshConfig); + gmsh.write(mesh); + + + // "Hack" mesh xy coordinates. + auto xyView = array::make_view(mesh.nodes().xy()); + const auto tijView = array::make_view(mesh.nodes().field("tij")); + for (idx_t i = 0; i < mesh.nodes().size(); ++i) { + const auto t = tijView(i, 0); + + const auto xy = PointXY(xyView(i, XX), xyView(i, YY)); + const auto alphabeta = csProjection.xy2alphabeta(xy, t); + + // Inverse function is degenerate when abs(alpha) == abs(beta) and + // abs(alpha) > 45. + const bool degenerate = + std::abs(alphabeta[0]) > 45. && approx_eq(std::abs(alphabeta[0]), std::abs(alphabeta[1])); + + // Check inverse function. + if (!degenerate) { + const auto newXy = csProjection.alphabeta2xy(alphabeta, t); + EXPECT_APPROX_EQ(xy[XX], newXy[XX], epsilon); + EXPECT_APPROX_EQ(xy[YY], newXy[YY], epsilon); + } + + // overwrite mesh xy field. + xyView(i, XX) = alphabeta[0]; + xyView(i, YY) = alphabeta[1]; + } + + // Output mesh updated mesh. + gmsh = output::Gmsh(outputID + "_after.msh", gmshConfig); + gmsh.write(mesh); +} + +void testJacobian(const std::string& gridType, const std::string& meshType, const std::string& outputID) { + // Create grid. + const auto grid = CubedSphereGrid(gridType); + + // Get projection. + const auto& csProjection = grid.cubedSphereProjection(); + + // Set xy0 values on different tiles. + const auto xy0 = std::array{PointXY{10., -35}, PointXY{100., -35}, PointXY{190., -35}, + PointXY{280., -35}, PointXY{10., 55}, PointXY{10., -125}}; + + // Get equivalent lonlat values. + const auto lonlat0 = std::array{csProjection.lonlat(xy0[0]), csProjection.lonlat(xy0[1]), + csProjection.lonlat(xy0[2]), csProjection.lonlat(xy0[3]), + csProjection.lonlat(xy0[4]), csProjection.lonlat(xy0[5])}; + // Set lonlat step size. + const auto hLonlat0 = PointLonLat{5., 10.}; + + // Perform first order extrapolation test on each tile. + for (size_t t = 0; t < 6; ++t) { + const auto jac = csProjection.jacobian(lonlat0[t]); + + double dxyOld{}; + + for (int i = 0; i < 5; ++i) { + // Sequentially halve lonlat increment. + const auto hLonlat = hLonlat0 * std::pow(0.5, i); + + // Make a first order estimate of xy using Jacobian. + const auto xy = xy0[t] + jac * hLonlat; + + // Get true xy values. + const auto xyTarget = csProjection.xy(lonlat0[t] + hLonlat); + double dxy = Point2::distance(xyTarget, xy); + + // Error should reduce roughly by a factor of four every iteration. + if (i > 0) { + const auto tol = 0.02; + EXPECT(dxy / dxyOld < 0.25 + tol); + } + dxyOld = dxy; + } + } + + // Create mesh. + auto mesh = MeshGenerator(meshType).generate(grid); + + // Create functionspace. + auto functionspace = functionspace::NodeColumns(mesh); + + // Create fields for each component of Jacobian. + auto fieldSet = FieldSet{}; + fieldSet.add(functionspace.createField(option::name("dx_by_dlambda"))); + fieldSet.add(functionspace.createField(option::name("dx_by_dphi"))); + fieldSet.add(functionspace.createField(option::name("dy_by_dlambda"))); + fieldSet.add(functionspace.createField(option::name("dy_by_dphi"))); + + auto dx_dlambda = array::make_view(fieldSet["dx_by_dlambda"]); + auto dx_dphi = array::make_view(fieldSet["dx_by_dphi"]); + auto dy_dlambda = array::make_view(fieldSet["dy_by_dlambda"]); + auto dy_dphi = array::make_view(fieldSet["dy_by_dphi"]); + const auto lonlat = array::make_view(functionspace.lonlat()); + + for (idx_t i = 0; i < functionspace.size(); ++i) { + const auto ll = PointLonLat(lonlat(i, LON), lonlat(i, LAT)); + const auto jac = csProjection.jacobian(ll); + dx_dlambda(i) = jac.dx_dlon(); + dx_dphi(i) = jac.dx_dlat(); + dy_dlambda(i) = jac.dy_dlon(); + dy_dphi(i) = jac.dy_dlat(); + } + + // Output mesh. + const auto gmshConfig = + util::Config("coordinates", "xy") | util::Config("ghost", true) | util::Config("info", true); + auto gmsh = output::Gmsh(outputID + "_jacobian.msh", gmshConfig); + gmsh.write(mesh); + gmsh.write(fieldSet, functionspace); +} + + +CASE("cubedsphere_xy_to_alphabeta_test") { + testProjection("CS-LFR-12", "cubedsphere", "cs_primal"); + testProjection("CS-LFR-12", "cubedsphere_dual", "cs_dual"); +} + +CASE("cubedsphere_jacobian_test") { + testJacobian("CS-LFR-24", "cubedsphere", "cs_primal"); + testJacobian("CS-LFR-24", "cubedsphere_dual", "cs_dual"); +} + +CASE("test_tiles") { + int resolution(2); + Grid gEA{"CS-EA-L-" + std::to_string(resolution)}; + Grid gLFR{"CS-LFR-L-" + std::to_string(resolution)}; + + using util::Constants; + + util::Config params; + grid::CubedSphereTiles f("cubedsphere_fv3"); + grid::CubedSphereTiles l("cubedsphere_lfric"); + + double cd[2]; + + idx_t jn(0); + + std::array EAOffset{0, + resolution * resolution + 1, + 2 * resolution * resolution + 2, + 3 * resolution * resolution + 2, + 4 * resolution * resolution + 2, + 5 * resolution * resolution + 2, + 6 * resolution * resolution + 2}; + + for (auto crd : gEA.lonlat()) { + atlas::PointLonLat pointLonLat = crd; + cd[LON] = pointLonLat.lon(); + cd[LAT] = pointLonLat.lat(); + + int t = f.indexFromLonLat(cd); + + gEA.projection().lonlat2xy(crd); + cd[LON] = crd.lon(); + cd[LAT] = crd.lat(); + + int t2 = f.indexFromXY(cd); + + for (std::size_t i = 0; i < 6; ++i) { + if (jn >= EAOffset[i] && jn < EAOffset[i + 1]) { + EXPECT(t == static_cast(i)); + EXPECT(t2 == static_cast(i)); + } + } + ++jn; + } + + std::array LFRicOffset{0, + resolution * resolution, + 2 * resolution * resolution, + 3 * resolution * resolution, + 4 * resolution * resolution, + 4 * resolution * resolution + (resolution + 1) * (resolution + 1), + 6 * resolution * resolution + 2}; + + jn = 0; + for (auto crd : gLFR.lonlat()) { + atlas::PointLonLat pointLonLat = crd; + + cd[LON] = pointLonLat.lon(); + cd[LAT] = pointLonLat.lat(); + + int t3 = l.indexFromLonLat(cd); + + gLFR.projection().lonlat2xy(crd); + cd[LON] = crd.lon(); + cd[LAT] = crd.lat(); + + int t4 = l.indexFromXY(cd); + + for (std::size_t i = 0; i < 6; ++i) { + if (jn >= LFRicOffset[i] && jn < LFRicOffset[i + 1]) { + EXPECT(t3 == static_cast(i)); + EXPECT(t4 == static_cast(i)); + } + } + ++jn; + } +} + +CASE("test_projection_cubedsphere_xy_latlon") { + int resolution(12); + std::vector grid_names{"CS-EA-L-" + std::to_string(resolution), + "CS-ED-L-" + std::to_string(resolution)}; + + for (std::string& s : grid_names) { + Grid g{s}; + for (auto crd : g.lonlat()) { + Point2 lonlat{crd}; + g->projection().lonlat2xy(crd); + g->projection().xy2lonlat(crd); + // except for point lonlat (90,82.5) on compiler pgc++ + // we have a maximum error tolerance of 1e-11 + EXPECT_APPROX_EQ(lonlat, crd, 1e-6); + } + for (auto crd : g.xy()) { + Point2 xy{crd}; + g->projection().xy2lonlat(crd); + g->projection().lonlat2xy(crd); + EXPECT_APPROX_EQ(xy, crd, 1e-6); + } + } +} + + +} // namespace test +} // namespace atlas + +int main(int argc, char** argv) { + return atlas::test::run(argc, argv); +} diff --git a/src/tests/projection/test_projection_cubed_sphere.cc b/src/tests/projection/test_projection_cubed_sphere.cc deleted file mode 100644 index f19ea633e..000000000 --- a/src/tests/projection/test_projection_cubed_sphere.cc +++ /dev/null @@ -1,141 +0,0 @@ -/* - * (C) Crown Copyright 2021 MetOffice - * - * 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 -#include -#include "atlas/grid.h" -#include "atlas/grid/Grid.h" -#include "atlas/projection/Projection.h" -#include "atlas/util/Constants.h" -#include "atlas/util/Point.h" - - -#include "atlas/grid/Tiles.h" -#include "atlas/grid/detail/tiles/FV3Tiles.h" -#include "atlas/grid/detail/tiles/LFRicTiles.h" -#include "atlas/grid/detail/tiles/Tiles.h" -#include "tests/AtlasTestEnvironment.h" - -namespace atlas { -namespace test { - - -//----------------------------------------------------------------------------- - -CASE("test_tiles") { - int resolution(2); - Grid gEA{"CS-EA-L-" + std::to_string(resolution)}; - Grid gLFR{"CS-LFR-L-" + std::to_string(resolution)}; - - using util::Constants; - - util::Config params; - grid::CubedSphereTiles f("cubedsphere_fv3"); - grid::CubedSphereTiles l("cubedsphere_lfric"); - - double cd[2]; - - idx_t jn(0); - - std::array EAOffset{0, - resolution * resolution + 1, - 2 * resolution * resolution + 2, - 3 * resolution * resolution + 2, - 4 * resolution * resolution + 2, - 5 * resolution * resolution + 2, - 6 * resolution * resolution + 2}; - - for (auto crd : gEA.lonlat()) { - atlas::PointLonLat pointLonLat = crd; - cd[LON] = pointLonLat.lon(); - cd[LAT] = pointLonLat.lat(); - - int t = f.indexFromLonLat(cd); - - gEA.projection().lonlat2xy(crd); - cd[LON] = crd.lon(); - cd[LAT] = crd.lat(); - - int t2 = f.indexFromXY(cd); - - for (std::size_t i = 0; i < 6; ++i) { - if (jn >= EAOffset[i] && jn < EAOffset[i + 1]) { - EXPECT(t == static_cast(i)); - EXPECT(t2 == static_cast(i)); - } - } - ++jn; - } - - std::array LFRicOffset{0, - resolution * resolution, - 2 * resolution * resolution, - 3 * resolution * resolution, - 4 * resolution * resolution, - 4 * resolution * resolution + (resolution + 1) * (resolution + 1), - 6 * resolution * resolution + 2}; - - jn = 0; - for (auto crd : gLFR.lonlat()) { - atlas::PointLonLat pointLonLat = crd; - - cd[LON] = pointLonLat.lon(); - cd[LAT] = pointLonLat.lat(); - - int t3 = l.indexFromLonLat(cd); - - gLFR.projection().lonlat2xy(crd); - cd[LON] = crd.lon(); - cd[LAT] = crd.lat(); - - int t4 = l.indexFromXY(cd); - - for (std::size_t i = 0; i < 6; ++i) { - if (jn >= LFRicOffset[i] && jn < LFRicOffset[i + 1]) { - EXPECT(t3 == static_cast(i)); - EXPECT(t4 == static_cast(i)); - } - } - ++jn; - } -} - -CASE("test_projection_cubedsphere_xy_latlon") { - int resolution(12); - std::vector grid_names{"CS-EA-L-" + std::to_string(resolution), - "CS-ED-L-" + std::to_string(resolution)}; - - for (std::string& s : grid_names) { - Grid g{s}; - for (auto crd : g.lonlat()) { - Point2 lonlat{crd}; - g->projection().lonlat2xy(crd); - g->projection().xy2lonlat(crd); - // except for point lonlat (90,82.5) on compiler pgc++ - // we have a maximum error tolerance of 1e-11 - EXPECT_APPROX_EQ(lonlat, crd, 1e-6); - } - for (auto crd : g.xy()) { - Point2 xy{crd}; - g->projection().xy2lonlat(crd); - g->projection().lonlat2xy(crd); - EXPECT_APPROX_EQ(xy, crd, 1e-6); - } - } -} - -//----------------------------------------------------------------------------- - -} // namespace test -} // namespace atlas - -int main(int argc, char** argv) { - return atlas::test::run(argc, argv); -} diff --git a/src/tests/trans/test_transgeneral.cc b/src/tests/trans/test_transgeneral.cc index 8cbfb7858..63b523f27 100644 --- a/src/tests/trans/test_transgeneral.cc +++ b/src/tests/trans/test_transgeneral.cc @@ -1582,7 +1582,7 @@ CASE("test_trans_levels") { #endif -#if ATLAS_HAVE_TRANS && defined(TRANS_HAVE_INVTRANS_ADJ) +#if ATLAS_HAVE_TRANS && (ATLAS_HAVE_ECTRANS || defined(TRANS_HAVE_INVTRANS_ADJ)) CASE("test_2level_adjoint_test_with_powerspectrum_convolution") { std::string grid_uid("F64"); // Regular Gaussian F ( 8 N^2) atlas::StructuredGrid g2(grid_uid);