From 8170eb4ba5cd4a847aa80035bc8202255ecca1a0 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Fri, 11 Nov 2022 11:59:09 +0100 Subject: [PATCH 01/25] Propagate field names between global and distributed fields --- src/atlas/functionspace/CellColumns.cc | 6 +++++- src/atlas/functionspace/NodeColumns.cc | 6 ++++-- src/tests/functionspace/fctest_functionspace.F90 | 8 ++++---- 3 files changed, 13 insertions(+), 7 deletions(-) diff --git a/src/atlas/functionspace/CellColumns.cc b/src/atlas/functionspace/CellColumns.cc index 85bb5c062..1c2858886 100644 --- a/src/atlas/functionspace/CellColumns.cc +++ b/src/atlas/functionspace/CellColumns.cc @@ -318,7 +318,7 @@ Field CellColumns::createField(const eckit::Configuration& options) const { } Field CellColumns::createField(const Field& other, const eckit::Configuration& config) const { - return createField(option::datatype(other.datatype()) | option::levels(other.levels()) | + return createField(option::name(other.name()) | option::datatype(other.datatype()) | option::levels(other.levels()) | option::variables(other.variables()) | config); } @@ -472,8 +472,12 @@ void CellColumns::scatter(const FieldSet& global_fieldset, FieldSet& local_field throw_Exception("datatype not supported", Here()); } + auto name = loc.name(); glb.metadata().broadcast(loc.metadata(), root); loc.metadata().set("global", false); + if( !name.empty() ) { + loc.metadata().set("name", name); + } } } void CellColumns::scatter(const Field& global, Field& local) const { diff --git a/src/atlas/functionspace/NodeColumns.cc b/src/atlas/functionspace/NodeColumns.cc index 9565f2409..6847ce0e5 100644 --- a/src/atlas/functionspace/NodeColumns.cc +++ b/src/atlas/functionspace/NodeColumns.cc @@ -333,7 +333,7 @@ Field NodeColumns::createField(const eckit::Configuration& config) const { } Field NodeColumns::createField(const Field& other, const eckit::Configuration& config) const { - return createField(option::datatype(other.datatype()) | option::levels(other.levels()) | + return createField(option::name(other.name()) | option::datatype(other.datatype()) | option::levels(other.levels()) | option::variables(other.variables()) | config); } @@ -540,7 +540,9 @@ void NodeColumns::scatter(const FieldSet& global_fieldset, FieldSet& local_field auto name = loc.name(); glb.metadata().broadcast(loc.metadata(), root); loc.metadata().set("global", false); - loc.metadata().set("name", name); + if( !name.empty() ) { + loc.metadata().set("name", name); + } } } diff --git a/src/tests/functionspace/fctest_functionspace.F90 b/src/tests/functionspace/fctest_functionspace.F90 index a7203cc2b..e46053042 100644 --- a/src/tests/functionspace/fctest_functionspace.F90 +++ b/src/tests/functionspace/fctest_functionspace.F90 @@ -84,7 +84,7 @@ module fcta_FunctionSpace_fxt field = fs%create_field(template) FCTEST_CHECK_EQUAL( field%rank() , 2 ) -FCTEST_CHECK_EQUAL( field%name() , "" ) +FCTEST_CHECK_EQUAL( field%name() , template%name() ) call field%final() field = fs%create_field(template,name="field") @@ -118,7 +118,7 @@ module fcta_FunctionSpace_fxt field = fs%create_field(template,global=.True.) FCTEST_CHECK_EQUAL( field%rank() , 2 ) -FCTEST_CHECK_EQUAL( field%name() , "" ) +FCTEST_CHECK_EQUAL( field%name() , template%name() ) call field%final() field = fs%create_field(template,name="field",global=.True.) @@ -182,7 +182,7 @@ module fcta_FunctionSpace_fxt field = fs%create_field(template) FCTEST_CHECK_EQUAL( field%rank() , 3 ) -FCTEST_CHECK_EQUAL( field%name() , "" ) +FCTEST_CHECK_EQUAL( field%name() , template%name() ) call field%final() field = fs%create_field(template,name="field") @@ -216,7 +216,7 @@ module fcta_FunctionSpace_fxt field = fs%create_field(template,global=.True.) FCTEST_CHECK_EQUAL( field%rank() , 3 ) -FCTEST_CHECK_EQUAL( field%name() , "" ) +FCTEST_CHECK_EQUAL( field%name() , template%name() ) call field%final() field = fs%create_field(template,name="field",global=.True.) From ad53e0e5b0043490aec6f4537f40d69409bb68b8 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Tue, 15 Nov 2022 09:48:56 +0100 Subject: [PATCH 02/25] Fix undefined behaviour. Thanks ubsan/asan :) --- .../detail/CubedSphereStructure.cc | 5 ++++- src/atlas/interpolation/method/Method.cc | 18 ++++++++++-------- src/atlas/interpolation/method/Method.h | 2 +- src/tests/grid/test_largegrid.cc | 4 ++-- 4 files changed, 17 insertions(+), 12 deletions(-) diff --git a/src/atlas/functionspace/detail/CubedSphereStructure.cc b/src/atlas/functionspace/detail/CubedSphereStructure.cc index dc369d1af..cd404d05c 100644 --- a/src/atlas/functionspace/detail/CubedSphereStructure.cc +++ b/src/atlas/functionspace/detail/CubedSphereStructure.cc @@ -61,7 +61,10 @@ CubedSphereStructure::CubedSphereStructure(const Field& tij, const Field& ghost, // Set tijToIdx vectors for (idx_t t = 0; t < 6; ++t) { // Set data array. - const size_t vecSize = static_cast((j_end(t) - j_begin(t)) * (i_end(t) - i_begin(t))); + size_t vecSize = 0; + if( j_end(t) >= j_begin(t) && i_end(t) >= i_begin(t) ) { + vecSize = static_cast((j_end(t) - j_begin(t)) * (i_end(t) - i_begin(t))); + } tijToIdx_.emplace_back(vecSize, invalid_index()); } diff --git a/src/atlas/interpolation/method/Method.cc b/src/atlas/interpolation/method/Method.cc index 70aaeb847..fb7f90221 100644 --- a/src/atlas/interpolation/method/Method.cc +++ b/src/atlas/interpolation/method/Method.cc @@ -383,14 +383,16 @@ void Method::do_execute(const Field& src, Field& tgt, Metadata&) const { haloExchange(src); - if (src.datatype().kind() == array::DataType::KIND_REAL64) { - interpolate_field(src, tgt, *matrix_); - } - else if (src.datatype().kind() == array::DataType::KIND_REAL32) { - interpolate_field(src, tgt, *matrix_); - } - else { - ATLAS_NOTIMPLEMENTED; + if( matrix_ ) { // (matrix == nullptr) when a partition is empty + if (src.datatype().kind() == array::DataType::KIND_REAL64) { + interpolate_field(src, tgt, *matrix_); + } + else if (src.datatype().kind() == array::DataType::KIND_REAL32) { + interpolate_field(src, tgt, *matrix_); + } + else { + ATLAS_NOTIMPLEMENTED; + } } // carry over missing value metadata diff --git a/src/atlas/interpolation/method/Method.h b/src/atlas/interpolation/method/Method.h index 5a0f6eac1..43bd1f4d8 100644 --- a/src/atlas/interpolation/method/Method.h +++ b/src/atlas/interpolation/method/Method.h @@ -155,7 +155,7 @@ class Method : public util::Object { void check_compatibility(const Field& src, const Field& tgt, const Matrix& W) const; private: - const Matrix* matrix_; + const Matrix* matrix_ = nullptr; std::shared_ptr matrix_shared_; interpolation::MatrixCache matrix_cache_; NonLinear nonLinear_; diff --git a/src/tests/grid/test_largegrid.cc b/src/tests/grid/test_largegrid.cc index c273fafb0..82716e7ef 100644 --- a/src/tests/grid/test_largegrid.cc +++ b/src/tests/grid/test_largegrid.cc @@ -48,14 +48,14 @@ CASE("test resources for cropping large grids") { std::vector gridnames{"L40000x20000", "N8000", "O8000"}; for (auto& gridname : gridnames) { - EXPECT(Trace::peakMemory() < 100 * MB); + EXPECT(Trace::peakMemory() < 128 * MB); SECTION(std::string("section") + gridname) { Trace trace(Here(), "Grid{" + gridname + ", GlobalDomain{-180}}"); auto grid = Grid{gridname, GlobalDomain{-180}}; trace.stopAndReport(); EXPECT(trace.elapsed() < 10.); - EXPECT(trace.peakMemory() < 100 * MB); + EXPECT(trace.peakMemory() < 128 * MB); } } } From 9a34035c475b023bb0b7d8deb92e9dce93723bb1 Mon Sep 17 00:00:00 2001 From: Slavko Brdar Date: Tue, 15 Nov 2022 09:48:56 +0100 Subject: [PATCH 03/25] ATLAS-363 Implement BlockStructuredColumns FunctionSpace without halo exchange functionality Co-authored-by: Willem Deconinck --- src/atlas/CMakeLists.txt | 6 + .../functionspace/BlockStructuredColumns.cc | 63 ++ .../functionspace/BlockStructuredColumns.h | 100 +++ .../detail/BlockStructuredColumns.cc | 389 ++++++++++++ .../detail/BlockStructuredColumns.h | 144 +++++ .../detail/BlockStructuredColumnsInterface.cc | 235 +++++++ .../detail/BlockStructuredColumnsInterface.h | 118 ++++ .../functionspace/detail/StructuredColumns.cc | 9 +- .../functionspace/detail/StructuredColumns.h | 8 + src/atlas_f/CMakeLists.txt | 4 + ...ionspace_BlockStructuredColumns_module.F90 | 580 ++++++++++++++++++ src/tests/functionspace/CMakeLists.txt | 17 + .../fctest_blockstructuredcolumns.F90 | 210 +++++++ .../functionspace/fctest_functionspace.F90 | 4 +- .../test_blockstructuredcolumns.cc | 198 ++++++ 15 files changed, 2076 insertions(+), 9 deletions(-) create mode 100644 src/atlas/functionspace/BlockStructuredColumns.cc create mode 100644 src/atlas/functionspace/BlockStructuredColumns.h create mode 100644 src/atlas/functionspace/detail/BlockStructuredColumns.cc create mode 100644 src/atlas/functionspace/detail/BlockStructuredColumns.h create mode 100644 src/atlas/functionspace/detail/BlockStructuredColumnsInterface.cc create mode 100644 src/atlas/functionspace/detail/BlockStructuredColumnsInterface.h create mode 100644 src/atlas_f/functionspace/atlas_functionspace_BlockStructuredColumns_module.F90 create mode 100644 src/tests/functionspace/fctest_blockstructuredcolumns.F90 create mode 100644 src/tests/functionspace/test_blockstructuredcolumns.cc diff --git a/src/atlas/CMakeLists.txt b/src/atlas/CMakeLists.txt index 258c14ec8..fc1d2bec5 100644 --- a/src/atlas/CMakeLists.txt +++ b/src/atlas/CMakeLists.txt @@ -438,6 +438,8 @@ field/detail/MissingValue.h list( APPEND atlas_functionspace_srcs functionspace.h +functionspace/BlockStructuredColumns.h +functionspace/BlockStructuredColumns.cc functionspace/CellColumns.h functionspace/CellColumns.cc functionspace/EdgeColumns.h @@ -454,6 +456,10 @@ functionspace/PointCloud.h functionspace/PointCloud.cc functionspace/CubedSphereColumns.h functionspace/CubedSphereColumns.cc +functionspace/detail/BlockStructuredColumns.h +functionspace/detail/BlockStructuredColumns.cc +functionspace/detail/BlockStructuredColumnsInterface.h +functionspace/detail/BlockStructuredColumnsInterface.cc functionspace/detail/FunctionSpaceImpl.h functionspace/detail/FunctionSpaceImpl.cc functionspace/detail/FunctionSpaceInterface.h diff --git a/src/atlas/functionspace/BlockStructuredColumns.cc b/src/atlas/functionspace/BlockStructuredColumns.cc new file mode 100644 index 000000000..746b3a12f --- /dev/null +++ b/src/atlas/functionspace/BlockStructuredColumns.cc @@ -0,0 +1,63 @@ +/* + * (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. + */ + +#include "atlas/functionspace/BlockStructuredColumns.h" + +namespace atlas { +namespace functionspace { + +// ---------------------------------------------------------------------------- + +BlockStructuredColumns::BlockStructuredColumns(): FunctionSpace(), functionspace_(nullptr) {} + +BlockStructuredColumns::BlockStructuredColumns(const FunctionSpace& functionspace): + FunctionSpace(functionspace), + functionspace_(dynamic_cast(get())) {} + +BlockStructuredColumns::BlockStructuredColumns(const Grid& grid, const eckit::Configuration& config): + FunctionSpace(new detail::BlockStructuredColumns(grid, config)), + functionspace_(dynamic_cast(get())) {} + +BlockStructuredColumns::BlockStructuredColumns(const Grid& grid, const grid::Partitioner& partitioner, + const eckit::Configuration& config): + FunctionSpace(new detail::BlockStructuredColumns(grid, partitioner, config)), + functionspace_(dynamic_cast(get())) {} + +BlockStructuredColumns::BlockStructuredColumns(const Grid& grid, const grid::Distribution& distribution, + const eckit::Configuration& config): + FunctionSpace(new detail::BlockStructuredColumns(grid, distribution, config)), + functionspace_(dynamic_cast(get())) {} + +BlockStructuredColumns::BlockStructuredColumns(const Grid& grid, const Vertical& vertical, const eckit::Configuration& config): + FunctionSpace(new detail::BlockStructuredColumns(grid, vertical, config)), + functionspace_(dynamic_cast(get())) {} + +BlockStructuredColumns::BlockStructuredColumns(const Grid& grid, const Vertical& vertical, const grid::Partitioner& partitioner, + const eckit::Configuration& config): + FunctionSpace(new detail::BlockStructuredColumns(grid, vertical, partitioner, config)), + functionspace_(dynamic_cast(get())) {} + +BlockStructuredColumns::BlockStructuredColumns(const Grid& grid, const grid::Distribution& distribution, const Vertical& vertical, + const eckit::Configuration& config): + FunctionSpace(new detail::BlockStructuredColumns(grid, distribution, vertical, config)), + functionspace_(dynamic_cast(get())) {} + +std::string BlockStructuredColumns::checksum(const FieldSet& fieldset) const { + return functionspace_->checksum(fieldset); +} + +std::string BlockStructuredColumns::checksum(const Field& field) const { + return functionspace_->checksum(field); +} + +// ---------------------------------------------------------------------------- + +} // namespace functionspace +} // namespace atlas diff --git a/src/atlas/functionspace/BlockStructuredColumns.h b/src/atlas/functionspace/BlockStructuredColumns.h new file mode 100644 index 000000000..f4a2f3c93 --- /dev/null +++ b/src/atlas/functionspace/BlockStructuredColumns.h @@ -0,0 +1,100 @@ +/* + * (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 "atlas/functionspace/FunctionSpace.h" +#include "atlas/functionspace/StructuredColumns.h" +#include "atlas/functionspace/detail/BlockStructuredColumns.h" + +namespace atlas { +namespace functionspace { + +// ------------------------------------------------------------------- + +class BlockStructuredColumns : public FunctionSpace { +public: + class Block; +public: + BlockStructuredColumns(); + BlockStructuredColumns(const FunctionSpace&); + BlockStructuredColumns(const Grid&, const eckit::Configuration& = util::NoConfig()); + BlockStructuredColumns(const Grid&, const grid::Partitioner&, const eckit::Configuration& = util::NoConfig()); + BlockStructuredColumns(const Grid&, const grid::Distribution&, const eckit::Configuration& = util::NoConfig()); + BlockStructuredColumns(const Grid&, const Vertical&, const eckit::Configuration& = util::NoConfig()); + BlockStructuredColumns(const Grid&, const Vertical&, const grid::Partitioner&, + const eckit::Configuration& = util::NoConfig()); + BlockStructuredColumns(const Grid&, const grid::Distribution&, const Vertical&, + const eckit::Configuration& = util::NoConfig()); + + static std::string type() { return detail::BlockStructuredColumns::static_type(); } + + operator bool() const { return valid(); } + bool valid() const { return functionspace_; } + + idx_t size() const { return functionspace_->size(); } +// idx_t sizeOwned() const { return functionspace_->sizeOwned(); } +// idx_t sizeHalo() const { return functionspace_->sizeHalo(); } + idx_t levels() const { return functionspace_->levels(); } + + const Vertical& vertical() const { return functionspace_->vertical(); } + + const StructuredGrid& grid() const { return functionspace_->grid(); } + + std::string checksum(const FieldSet&) const; + std::string checksum(const Field&) const; + + idx_t index(idx_t blk, idx_t rof) const { return functionspace_->index(blk, rof); } +// idx_t i_begin(idx_t j) const { return functionspace_->i_begin(j); } +// idx_t i_end(idx_t j) const { return functionspace_->i_end(j); } +// idx_t j_begin() const { return functionspace_->j_begin(); } +// idx_t j_end() const { return functionspace_->j_end(); } + idx_t k_begin() const { return functionspace_->k_begin(); } + idx_t k_end() const { return functionspace_->k_end(); } + idx_t nproma() const { return functionspace_->nproma(); } + idx_t nblks() const { return functionspace_->nblks(); } + + Field xy() const { return functionspace_->xy(); } + Field partition() const { return functionspace_->partition(); } + Field global_index() const { return functionspace_->global_index(); } + Field remote_index() const { return functionspace_->remote_index(); } + Field index_i() const { return functionspace_->index_i(); } + Field index_j() const { return functionspace_->index_j(); } + Field ghost() const { return functionspace_->ghost(); } + + const Block block(idx_t jblk) const { + return Block(functionspace_->block_begin(jblk), functionspace_->block_size(jblk)); + } + + size_t footprint() const { return functionspace_->footprint(); } + + class Block { + public: + Block(idx_t begin, idx_t size) : begin_(begin), size_(size) {} + idx_t index(idx_t j) const { return begin_ + j; }; + idx_t size() const { return size_; } + private: + idx_t begin_; + idx_t size_; + }; + +private: + const detail::BlockStructuredColumns* functionspace_; + void setup(const Grid& grid, const Vertical& vertical, const grid::Distribution& distribution, + const eckit::Configuration& config); +}; + +// ------------------------------------------------------------------- + + +} // namespace functionspace +} // namespace atlas diff --git a/src/atlas/functionspace/detail/BlockStructuredColumns.cc b/src/atlas/functionspace/detail/BlockStructuredColumns.cc new file mode 100644 index 000000000..4f1b84f5e --- /dev/null +++ b/src/atlas/functionspace/detail/BlockStructuredColumns.cc @@ -0,0 +1,389 @@ +/* + * (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. + */ + +#include "atlas/functionspace/BlockStructuredColumns.h" + +#include +#include +#include +#include +#include + +#include "eckit/utils/MD5.h" + +#include "atlas/array/Array.h" +#include "atlas/array/MakeView.h" +#include "atlas/domain.h" +#include "atlas/field/FieldSet.h" +#include "atlas/grid/Distribution.h" +#include "atlas/grid/Partitioner.h" +#include "atlas/grid/StructuredGrid.h" +#include "atlas/grid/StructuredPartitionPolygon.h" +#include "atlas/library/Library.h" +#include "atlas/mesh/Mesh.h" +#include "atlas/parallel/Checksum.h" +#include "atlas/parallel/mpi/mpi.h" +#include "atlas/parallel/omp/fill.h" +#include "atlas/parallel/omp/omp.h" +#include "atlas/runtime/Exception.h" +#include "atlas/runtime/Trace.h" +#include "atlas/util/Checksum.h" +#include "atlas/util/CoordinateEnums.h" +#include "atlas/util/detail/Cache.h" + +namespace { +using namespace atlas; + +template +void block_copy(const Field sloc, Field loc, const functionspace::detail::BlockStructuredColumns& fs) { + if (sloc.variables() and sloc.levels()) { + auto loc_v = array::make_view(loc); + auto sloc_v = array::make_view(sloc); + for (idx_t jblk = 0; jblk < fs.nblks(); ++jblk) { + const idx_t blk_size = fs.block_size(jblk); + const idx_t blk_begin = fs.block_begin(jblk); + for (idx_t jvar = 0; jvar < sloc.shape(2); ++jvar) { + for (idx_t jlev = 0; jlev < sloc.shape(1); ++jlev) { + for (idx_t jrof = 0; jrof < blk_size; ++jrof) { + loc_v(jblk, jvar, jlev, jrof) = sloc_v(blk_begin+jrof, jlev, jvar); + } + } + } + } + } + else if (not sloc.variables() and sloc.levels()) { + auto loc_v = array::make_view(loc); + auto sloc_v = array::make_view(sloc); + for (idx_t jblk = 0; jblk < fs.nblks(); ++jblk) { + const idx_t blk_size = fs.block_size(jblk); + const idx_t blk_begin = fs.block_begin(jblk); + for (idx_t jlev = 0; jlev < sloc.shape(1); ++jlev) { + for (idx_t jrof = 0; jrof < blk_size; ++jrof) { + loc_v(jblk, jlev, jrof) = sloc_v(blk_begin+jrof, jlev); + } + } + } + } + else if (sloc.variables() and not sloc.levels()) { + auto loc_v = array::make_view(loc); + auto sloc_v = array::make_view(sloc); + for (idx_t jblk = 0; jblk < fs.nblks(); ++jblk) { + const idx_t blk_size = fs.block_size(jblk); + const idx_t blk_begin = fs.block_begin(jblk); + for (idx_t jvar = 0; jvar < sloc.shape(1); ++jvar) { + for (idx_t jrof = 0; jrof < blk_size; ++jrof) { + loc_v(jblk, jvar, jrof) = sloc_v(blk_begin+jrof, jvar); + } + } + } + } + else { + auto loc_v = array::make_view(loc); + auto sloc_v = array::make_view(sloc); + for (idx_t jblk = 0; jblk < fs.nblks(); ++jblk) { + const idx_t blk_size = fs.block_size(jblk); + const idx_t blk_begin = fs.block_begin(jblk); + for (idx_t jrof = 0; jrof < blk_size; ++jrof) { + loc_v(jblk, jrof) = sloc_v(blk_begin+jrof); + } + } + } +} + +template +void rev_block_copy(const Field loc, Field sloc, const functionspace::detail::BlockStructuredColumns& fs) { + if (loc.variables() and loc.levels()) { + auto loc_v = array::make_view(loc); + auto sloc_v = array::make_view(sloc); + for (idx_t jblk = 0; jblk < fs.nblks(); ++jblk) { + const idx_t blk_size = fs.block_size(jblk); + const idx_t blk_begin = fs.block_begin(jblk); + for (idx_t jvar = 0; jvar < sloc.shape(2); ++jvar) { + for (idx_t jlev = 0; jlev < sloc.shape(1); ++jlev) { + for (idx_t jrof = 0; jrof < blk_size; ++jrof) { + sloc_v(blk_begin+jrof, jlev, jvar) = loc_v(jblk, jvar, jlev, jrof); + } + } + } + } + } + else if (not loc.variables() and loc.levels()) { + auto loc_v = array::make_view(loc); + auto sloc_v = array::make_view(sloc); + for (idx_t jblk = 0; jblk < fs.nblks(); ++jblk) { + const idx_t blk_size = fs.block_size(jblk); + const idx_t blk_begin = fs.block_begin(jblk); + for (idx_t jlev = 0; jlev < sloc.shape(1); ++jlev) { + for (idx_t jrof = 0; jrof < blk_size; ++jrof) { + sloc_v(blk_begin+jrof, jlev) = loc_v(jblk, jlev, jrof); + } + } + } + } + else if (loc.variables() and not loc.levels()) { + auto loc_v = array::make_view(loc); + auto sloc_v = array::make_view(sloc); + for (idx_t jblk = 0; jblk < fs.nblks(); ++jblk) { + const idx_t blk_size = fs.block_size(jblk); + const idx_t blk_begin = fs.block_begin(jblk); + for (idx_t jvar = 0; jvar < sloc.shape(1); ++jvar) { + for (idx_t jrof = 0; jrof < blk_size; ++jrof) { + sloc_v(blk_begin+jrof, jvar) = loc_v(jblk, jvar, jrof); + } + } + } + } + else { + auto loc_v = array::make_view(loc); + auto sloc_v = array::make_view(sloc); + for (idx_t jblk = 0; jblk < fs.nblks(); ++jblk) { + const idx_t blk_size = fs.block_size(jblk); + const idx_t blk_begin = fs.block_begin(jblk); + for (idx_t jrof = 0; jrof < blk_size; ++jrof) { + sloc_v(blk_begin+jrof) = loc_v(jblk, jrof); + } + } + } +} + +}// namespace + +namespace atlas { +namespace functionspace { +namespace detail { + + +array::ArrayShape BlockStructuredColumns::config_shape(const eckit::Configuration& config) const { + array::ArrayShape shape; + + bool global = false; + config.get("global", global); + if (global) { + return structuredcolumns_->config_shape(config); + } + else { + shape.emplace_back(nblks_); + idx_t variables(0); + config.get("variables", variables); + if (variables > 0) { + shape.emplace_back(variables); + } + idx_t levels(structuredcolumns_->levels()); + config.get("levels", levels); + if (levels > 0) { + shape.emplace_back(levels); + } + shape.emplace_back(nproma_); + } + return shape; +} + +array::ArraySpec BlockStructuredColumns::config_spec(const eckit::Configuration& config) const { + return array::ArraySpec(config_shape(config), structuredcolumns_->config_alignment(config)); +} + +// ---------------------------------------------------------------------------- +// Constructor +// ---------------------------------------------------------------------------- +BlockStructuredColumns::BlockStructuredColumns(const Grid& grid, const eckit::Configuration& config): + BlockStructuredColumns::BlockStructuredColumns(grid, grid::Partitioner(), config) {} + +BlockStructuredColumns::BlockStructuredColumns(const Grid& grid, const grid::Partitioner& p, const eckit::Configuration& config): + BlockStructuredColumns(grid, Vertical(config), p, config) {} + +BlockStructuredColumns::BlockStructuredColumns(const Grid& grid, const grid::Distribution& distribution, + const eckit::Configuration& config): + BlockStructuredColumns(grid, distribution, Vertical(config), config) {} + +BlockStructuredColumns::BlockStructuredColumns(const Grid& grid, const grid::Distribution& distribution, const Vertical& vertical, + const eckit::Configuration& config): + structuredcolumns_(new StructuredColumns(grid, distribution, vertical, config)), + structuredcolumns_handle_(structuredcolumns_){ + setup(config); +} + +BlockStructuredColumns::BlockStructuredColumns(const Grid& grid, const Vertical& vertical, const eckit::Configuration& config): + BlockStructuredColumns(grid, vertical, grid::Partitioner(), config) {} + +BlockStructuredColumns::BlockStructuredColumns(const Grid& grid, const Vertical& vertical, const grid::Partitioner& p, + const eckit::Configuration& config): + structuredcolumns_(new StructuredColumns(grid, vertical, p, config)), + structuredcolumns_handle_(structuredcolumns_){ + setup(config); +} + +// ---------------------------------------------------------------------------- +// Create Field +// ---------------------------------------------------------------------------- +Field BlockStructuredColumns::createField(const eckit::Configuration& options) const { + Field field(structuredcolumns_->config_name(options), structuredcolumns_->config_datatype(options), config_spec(options)); + structuredcolumns_->set_field_metadata(options, field); + field.set_functionspace(this); + return field; +} + +Field BlockStructuredColumns::createField(const Field& other, const eckit::Configuration& config) const { + return createField(option::datatype(other.datatype()) | option::levels(other.levels()) | + option::variables(other.variables()) | + option::type(other.metadata().getString("type", "scalar")) | config); +} +// ---------------------------------------------------------------------------- + +// ---------------------------------------------------------------------------- +// Scatter FieldSet +// ---------------------------------------------------------------------------- +void BlockStructuredColumns::scatter(const FieldSet& global_fieldset, FieldSet& local_fieldset) const { + ATLAS_ASSERT(local_fieldset.size() == global_fieldset.size()); + for (idx_t f = 0; f < local_fieldset.size(); ++f) { + const Field& glb = global_fieldset[f]; + auto config = option::datatype(glb.datatype()) | option::levels(glb.levels()) | option::variables(glb.variables()); + auto sloc = structuredcolumns_->createField(config); + const idx_t nb_fields = 1; + idx_t root(0); + glb.metadata().get("owner", root); + + Field& loc = local_fieldset[f]; + glb.metadata().broadcast(loc.metadata(), root); + loc.metadata().set("global", false); + + if (sloc.datatype().kind() == array::DataType::kind()) { + parallel::Field glb_field(structuredcolumns_->make_leveled_view(glb)); + parallel::Field sloc_field(structuredcolumns_->make_leveled_view(sloc)); + structuredcolumns_->scatter().scatter(&glb_field, &sloc_field, nb_fields, root); + block_copy(sloc, loc, *this); + } + else if (sloc.datatype().kind() == array::DataType::kind()) { + parallel::Field glb_field(structuredcolumns_->make_leveled_view(glb)); + parallel::Field sloc_field(structuredcolumns_->make_leveled_view(sloc)); + structuredcolumns_->scatter().scatter(&glb_field, &sloc_field, nb_fields, root); + block_copy(sloc, loc, *this); + } + else if (sloc.datatype().kind() == array::DataType::kind()) { + parallel::Field glb_field(structuredcolumns_->make_leveled_view(glb)); + parallel::Field sloc_field(structuredcolumns_->make_leveled_view(sloc)); + structuredcolumns_->scatter().scatter(&glb_field, &sloc_field, nb_fields, root); + block_copy(sloc, loc, *this); + } + else if (sloc.datatype().kind() == array::DataType::kind()) { + parallel::Field glb_field(structuredcolumns_->make_leveled_view(glb)); + parallel::Field sloc_field(structuredcolumns_->make_leveled_view(sloc)); + structuredcolumns_->scatter().scatter(&glb_field, &sloc_field, nb_fields, root); + block_copy(sloc, loc, *this); + } + else { + throw_Exception("datatype not supported", Here()); + } + } +} +// ---------------------------------------------------------------------------- + +// ---------------------------------------------------------------------------- +// Scatter Field +// ---------------------------------------------------------------------------- +void BlockStructuredColumns::scatter(const Field& global, Field& local) const { + FieldSet global_fields; + FieldSet local_fields; + global_fields.add(global); + local_fields.add(local); + scatter(global_fields, local_fields); +} +// ---------------------------------------------------------------------------- + +// ---------------------------------------------------------------------------- +// Gather FieldSet +// ---------------------------------------------------------------------------- +void BlockStructuredColumns::gather(const FieldSet& local_fieldset, FieldSet& global_fieldset) const { + ATLAS_ASSERT(local_fieldset.size() == global_fieldset.size()); + for (idx_t f = 0; f < local_fieldset.size(); ++f) { + const Field& loc = local_fieldset[f]; + auto config = option::datatype(loc.datatype()) | option::levels(loc.levels()) | option::variables(loc.variables()); + auto sloc = structuredcolumns_->createField(config); + + + idx_t root(0); + sloc.metadata().set("global", false); + Field& glb = global_fieldset[f]; + glb.metadata().broadcast(sloc.metadata(), root); + const idx_t nb_fields = 1; + glb.metadata().get("owner", root); + + if (sloc.datatype().kind() == array::DataType::kind()) { + rev_block_copy(loc, sloc, *this); + parallel::Field sloc_field(structuredcolumns_->make_leveled_view(sloc)); + parallel::Field glb_field(structuredcolumns_->make_leveled_view(glb)); + structuredcolumns_->gather().gather(&sloc_field, &glb_field, nb_fields, root); + } + else if (sloc.datatype().kind() == array::DataType::kind()) { + rev_block_copy(loc, sloc, *this); + parallel::Field sloc_field(structuredcolumns_->make_leveled_view(sloc)); + parallel::Field glb_field(structuredcolumns_->make_leveled_view(glb)); + structuredcolumns_->gather().gather(&sloc_field, &glb_field, nb_fields, root); + } + else if (sloc.datatype().kind() == array::DataType::kind()) { + rev_block_copy(loc, sloc, *this); + parallel::Field sloc_field(structuredcolumns_->make_leveled_view(sloc)); + parallel::Field glb_field(structuredcolumns_->make_leveled_view(glb)); + structuredcolumns_->gather().gather(&sloc_field, &glb_field, nb_fields, root); + } + else if (sloc.datatype().kind() == array::DataType::kind()) { + rev_block_copy(loc, sloc, *this); + parallel::Field sloc_field(structuredcolumns_->make_leveled_view(sloc)); + parallel::Field glb_field(structuredcolumns_->make_leveled_view(glb)); + structuredcolumns_->gather().gather(&sloc_field, &glb_field, nb_fields, root); + } + else { + throw_Exception("datatype not supported", Here()); + } + } +} +// ---------------------------------------------------------------------------- + +void BlockStructuredColumns::setup(const eckit::Configuration &config) { + nproma_ = 1; + idx_t tmp_nproma; + if (config.get("nproma", tmp_nproma)) { + ATLAS_ASSERT(tmp_nproma > 0); + nproma_ = tmp_nproma; + } + + nblks_ = std::floor( structuredcolumns_->size() / nproma_ ); + endblk_size_ = nproma_; + if (structuredcolumns_->size() % nproma_ > 0) { + endblk_size_ = structuredcolumns_->size() - nblks_ * nproma_; + nblks_++; + } +} + +// ---------------------------------------------------------------------------- +// Gather Field +// ---------------------------------------------------------------------------- +void BlockStructuredColumns::gather(const Field& local, Field& global) const { + FieldSet local_fields; + FieldSet global_fields; + local_fields.add(local); + global_fields.add(global); + gather(local_fields, global_fields); +} + +std::string BlockStructuredColumns::checksum(const Field&) const { + ATLAS_NOTIMPLEMENTED; +} + +std::string BlockStructuredColumns::checksum(const FieldSet&) const { + ATLAS_NOTIMPLEMENTED; +} + + +// ---------------------------------------------------------------------------- + + +} // namespace detail +} // namespace functionspace +} // namespace atlas diff --git a/src/atlas/functionspace/detail/BlockStructuredColumns.h b/src/atlas/functionspace/detail/BlockStructuredColumns.h new file mode 100644 index 000000000..48fb9aee9 --- /dev/null +++ b/src/atlas/functionspace/detail/BlockStructuredColumns.h @@ -0,0 +1,144 @@ +/* + * (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 "atlas/array/DataType.h" +#include "atlas/field/Field.h" +#include "atlas/functionspace/detail/FunctionSpaceImpl.h" +#include "atlas/functionspace/detail/StructuredColumns.h" +#include "atlas/functionspace/StructuredColumns.h" +#include "atlas/grid/StructuredGrid.h" +#include "atlas/grid/Vertical.h" +#include "atlas/library/config.h" +#include "atlas/option.h" +#include "atlas/parallel/mpi/mpi.h" +#include "atlas/parallel/omp/omp.h" +#include "atlas/runtime/Exception.h" +#include "atlas/util/Config.h" +#include "atlas/util/ObjectHandle.h" +#include "atlas/util/Point.h" +#include "atlas/util/Polygon.h" +#include "atlas/util/vector.h" + +namespace atlas { +namespace parallel { +class GatherScatter; +//class HaloExchange; +//class Checksum; +} // namespace parallel +} // namespace atlas + +namespace atlas { +class Field; +//class FieldSet; +class Grid; +class StructuredGrid; +} // namespace atlas + +namespace atlas { +namespace grid { +class Distribution; +class Partitioner; +} // namespace grid +} // namespace atlas + +namespace atlas { +namespace functionspace { +namespace detail { + +// ------------------------------------------------------------------- + +class BlockStructuredColumns : public FunctionSpaceImpl { +public: + BlockStructuredColumns(const Grid&, const eckit::Configuration& = util::NoConfig()); + + BlockStructuredColumns(const Grid&, const grid::Partitioner&, const eckit::Configuration& = util::NoConfig()); + + BlockStructuredColumns(const Grid&, const grid::Distribution&, const eckit::Configuration& = util::NoConfig()); + + BlockStructuredColumns(const Grid&, const grid::Distribution&, const Vertical&, + const eckit::Configuration& = util::NoConfig()); + + BlockStructuredColumns(const Grid&, const Vertical&, const eckit::Configuration& = util::NoConfig()); + + BlockStructuredColumns(const Grid&, const Vertical&, const grid::Partitioner&, + const eckit::Configuration& = util::NoConfig()); + + static std::string static_type() { return "BlockStructuredColumns"; } + std::string type() const override { return static_type(); } + std::string distribution() const override { return structuredcolumns_->distribution(); } + + Field createField(const eckit::Configuration&) const override; + Field createField(const Field&, const eckit::Configuration&) const override; + + void scatter(const FieldSet&, FieldSet&) const override; + void scatter(const Field&, Field&) const override; + void gather(const FieldSet&, FieldSet&) const override; + void gather(const Field&, Field&) const override; + + idx_t size() const override { return structuredcolumns_->size(); } + idx_t index(idx_t jblk, idx_t jrof) const { + return jblk * nproma_ + jrof; // local index; + } + idx_t nproma() const { return nproma_; } + idx_t nblks() const { return nblks_; } + + const Vertical& vertical() const { return structuredcolumns_->vertical(); } + const StructuredGrid& grid() const { return structuredcolumns_->grid(); } + + idx_t levels() const { return structuredcolumns_->levels(); } + Field lonlat() const override { return structuredcolumns_->lonlat(); } + Field xy() const { return structuredcolumns_->xy(); } + Field z() const { return structuredcolumns_->z(); } + Field partition() const { return structuredcolumns_->partition(); } + Field global_index() const override { return structuredcolumns_->global_index(); } + Field remote_index() const override { return structuredcolumns_->remote_index(); } + Field index_i() const { return structuredcolumns_->index_i(); } + Field index_j() const { return structuredcolumns_->index_j(); } + Field ghost() const override { return structuredcolumns_->ghost(); } + size_t footprint() const override { return structuredcolumns_->footprint(); } + idx_t nb_partitions() const override { return structuredcolumns_->nb_partitions(); } + const StructuredColumns& structuredcolumns() const { return *structuredcolumns_; } + idx_t block_begin(idx_t jblk) const { return jblk * nproma_; } + idx_t block_size(idx_t jblk) const { return (jblk != nblks_-1 ? nproma_ : endblk_size_); } + idx_t k_begin() const { return vertical().k_begin(); } + idx_t k_end() const { return vertical().k_end(); } + + std::string checksum(const FieldSet&) const; + std::string checksum(const Field&) const; + +private: // methods + array::ArrayShape config_shape(const eckit::Configuration&) const; + array::ArrayAlignment config_alignment(const eckit::Configuration&) const; + array::ArraySpec config_spec(const eckit::Configuration&) const; + + void create_remote_index() const { structuredcolumns_->create_remote_index(); } + +private: // data + idx_t nproma_; + idx_t endblk_size_; + idx_t nblks_; + + detail::StructuredColumns* structuredcolumns_; + functionspace::StructuredColumns structuredcolumns_handle_; + + void setup(const eckit::Configuration& config); +}; + +// ------------------------------------------------------------------- + +} // namespace detail +} // namespace functionspace +} // namespace atlas diff --git a/src/atlas/functionspace/detail/BlockStructuredColumnsInterface.cc b/src/atlas/functionspace/detail/BlockStructuredColumnsInterface.cc new file mode 100644 index 000000000..988bc1c86 --- /dev/null +++ b/src/atlas/functionspace/detail/BlockStructuredColumnsInterface.cc @@ -0,0 +1,235 @@ +/* + * (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. + */ + +#include + +#include "BlockStructuredColumnsInterface.h" + +#include "atlas/field/FieldSet.h" +#include "atlas/field/detail/FieldImpl.h" +#include "atlas/grid/Distribution.h" +#include "atlas/grid/Grid.h" +#include "atlas/grid/Partitioner.h" +#include "atlas/grid/detail/distribution/DistributionImpl.h" +#include "atlas/runtime/Exception.h" + +namespace atlas { +namespace functionspace { + +// ---------------------------------------------------------------------------- +// Fortran interfaces +// ---------------------------------------------------------------------------- + +namespace detail { +struct BlockStructuredColumnsFortranAccess { + detail::StructuredColumns::Map2to1& ij2gp_; + BlockStructuredColumnsFortranAccess(const detail::BlockStructuredColumns& fs) : + ij2gp_(const_cast(fs.structuredcolumns()).ij2gp_) {}; +}; +} // namespace detail + + +extern "C" { + +const detail::BlockStructuredColumns* atlas__functionspace__BStructuredColumns__new__grid( + const Grid::Implementation* grid, const eckit::Configuration* config) { + return new detail::BlockStructuredColumns(Grid(grid), grid::Partitioner(), *config); +} + +const detail::BlockStructuredColumns* atlas__functionspace__BStructuredColumns__new__grid_dist( + const Grid::Implementation* grid, const grid::DistributionImpl* dist, const eckit::Configuration* config) { + return new detail::BlockStructuredColumns(Grid(grid), grid::Distribution(dist), *config); +} + +const detail::BlockStructuredColumns* atlas__functionspace__BStructuredColumns__new__grid_dist_vert( + const Grid::Implementation* grid, const grid::DistributionImpl* dist, const Vertical* vert, + const eckit::Configuration* config) { + return new detail::BlockStructuredColumns(Grid(grid), grid::Distribution(dist), *vert, *config); +} +const detail::BlockStructuredColumns* atlas__functionspace__BStructuredColumns__new__grid_dist_config( + const Grid::Implementation* grid, const grid::DistributionImpl* dist, const eckit::Configuration* config) { + return new detail::BlockStructuredColumns(Grid(grid), grid::Distribution(dist), *config); +} + +const detail::BlockStructuredColumns* atlas__functionspace__BStructuredColumns__new__grid_part( + const Grid::Implementation* grid, const PartitionerImpl* partitioner, const eckit::Configuration* config) { + return new detail::BlockStructuredColumns(Grid(grid), grid::Partitioner(partitioner), *config); +} + +const detail::BlockStructuredColumns* atlas__functionspace__BStructuredColumns__new__grid_part_vert( + const Grid::Implementation* grid, const PartitionerImpl* partitioner, const Vertical* vert, + const eckit::Configuration* config) { + return new detail::BlockStructuredColumns(Grid(grid), *vert, grid::Partitioner(partitioner), *config); +} + +void atlas__functionspace__BStructuredColumns__gather_field(const detail::BlockStructuredColumns* This, + const field::FieldImpl* local, field::FieldImpl* global) { + ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_functionspace_BlockStructuredColumns"); + ATLAS_ASSERT(global != nullptr, "Cannot access uninitialised atlas_Field"); + ATLAS_ASSERT(local != nullptr, "Cannot access uninitialised atlas_Field"); + const Field l(local); + Field g(global); + This->gather(l, g); +} + +void atlas__functionspace__BStructuredColumns__scatter_field(const detail::BlockStructuredColumns* This, + const field::FieldImpl* global, field::FieldImpl* local) { + ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_functionspace_BlockStructuredColumns"); + ATLAS_ASSERT(global != nullptr, "Cannot access uninitialised atlas_Field"); + ATLAS_ASSERT(local != nullptr, "Cannot access uninitialised atlas_Field"); + const Field g(global); + Field l(local); + This->scatter(g, l); +} + +void atlas__functionspace__BStructuredColumns__gather_fieldset(const detail::BlockStructuredColumns* This, + const field::FieldSetImpl* local, + field::FieldSetImpl* global) { + ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_functionspace_BlockStructuredColumns"); + ATLAS_ASSERT(global != nullptr, "Cannot access uninitialised atlas_FieldSet"); + ATLAS_ASSERT(local != nullptr, "Cannot access uninitialised atlas_FieldSet"); + const FieldSet l(local); + FieldSet g(global); + This->gather(l, g); +} + +void atlas__functionspace__BStructuredColumns__scatter_fieldset(const detail::BlockStructuredColumns* This, + const field::FieldSetImpl* global, + field::FieldSetImpl* local) { + ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_functionspace_BlockStructuredColumns"); + ATLAS_ASSERT(global != nullptr, "Cannot access uninitialised atlas_FieldSet"); + ATLAS_ASSERT(local != nullptr, "Cannot access uninitialised atlas_FieldSet"); + const FieldSet g(global); + FieldSet l(local); + This->scatter(g, l); +} + +void atlas__fs__BStructuredColumns__checksum_fieldset(const detail::BlockStructuredColumns* This, + const field::FieldSetImpl* fieldset, char*& checksum, idx_t& size, + int& allocated) { + ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_functionspace_BlockStructuredColumns"); + ATLAS_ASSERT(fieldset != nullptr, "Cannot access uninitialised atlas_FieldSet"); + std::string checksum_str(This->checksum(fieldset)); + size = static_cast(checksum_str.size()); + checksum = new char[size + 1]; + allocated = true; + std::strncpy(checksum, checksum_str.c_str(), size + 1); +} + +void atlas__fs__BStructuredColumns__checksum_field(const detail::BlockStructuredColumns* This, const field::FieldImpl* field, + char*& checksum, idx_t& size, int& allocated) { + ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_functionspace_BlockStructuredColumns"); + ATLAS_ASSERT(field != nullptr, "Cannot access uninitialised atlas_Field"); + std::string checksum_str(This->checksum(field)); + size = static_cast(checksum_str.size()); + checksum = new char[size + 1]; + allocated = true; + std::strncpy(checksum, checksum_str.c_str(), size + 1); +} + +void atlas__fs__BStructuredColumns__index_host(const detail::BlockStructuredColumns* This, idx_t*& data, idx_t& i_min, + idx_t& i_max, idx_t& j_min, idx_t& j_max) { + ATLAS_ASSERT(This != nullptr, "Cannot access uninitialised atlas_functionspace_BlockStructuredColumns"); + auto _This = detail::BlockStructuredColumnsFortranAccess{*This}; + data = _This.ij2gp_.data_.data(); + i_min = _This.ij2gp_.i_min_ + 1; + i_max = _This.ij2gp_.i_max_ + 1; + j_min = _This.ij2gp_.j_min_ + 1; + j_max = _This.ij2gp_.j_max_ + 1; +} + +idx_t atlas__fs__BStructuredColumns__j_begin(const detail::BlockStructuredColumns* This) { + return This->structuredcolumns().j_begin() + 1; +} +idx_t atlas__fs__BStructuredColumns__j_end(const detail::BlockStructuredColumns* This) { + return This->structuredcolumns().j_end(); +} +idx_t atlas__fs__BStructuredColumns__i_begin(const detail::BlockStructuredColumns* This, idx_t j) { + return This->structuredcolumns().i_begin(j - 1) + 1; +} +idx_t atlas__fs__BStructuredColumns__i_end(const detail::BlockStructuredColumns* This, idx_t j) { + return This->structuredcolumns().i_end(j - 1); +} +idx_t atlas__fs__BStructuredColumns__j_begin_halo(const detail::BlockStructuredColumns* This) { + return This->structuredcolumns().j_begin_halo() + 1; +} +idx_t atlas__fs__BStructuredColumns__j_end_halo(const detail::BlockStructuredColumns* This) { + return This->structuredcolumns().j_end_halo(); +} +idx_t atlas__fs__BStructuredColumns__i_begin_halo(const detail::BlockStructuredColumns* This, idx_t j) { + return This->structuredcolumns().i_begin_halo(j - 1) + 1; +} +idx_t atlas__fs__BStructuredColumns__i_end_halo(const detail::BlockStructuredColumns* This, idx_t j) { + return This->structuredcolumns().i_end_halo(j - 1); +} +idx_t atlas__fs__BStructuredColumns__block_begin(const detail::BlockStructuredColumns* This, idx_t jblk) { + return This->block_begin(jblk); +} +idx_t atlas__fs__BStructuredColumns__block_size(const detail::BlockStructuredColumns* This, idx_t jblk) { + return This->block_size(jblk); +} +idx_t atlas__fs__BStructuredColumns__nblks(const detail::BlockStructuredColumns* This) { + return This->nblks(); +} + +field::FieldImpl* atlas__fs__BStructuredColumns__xy(const detail::BlockStructuredColumns* This) { + return This->xy().get(); +} + +field::FieldImpl* atlas__fs__BStructuredColumns__z(const detail::BlockStructuredColumns* This) { + ATLAS_ASSERT(This != nullptr); + field::FieldImpl* field; + { + Field f = This->z(); + field = f.get(); + field->attach(); + } + field->detach(); + return field; +} + +field::FieldImpl* atlas__fs__BStructuredColumns__partition(const detail::BlockStructuredColumns* This) { + return This->partition().get(); +} + +field::FieldImpl* atlas__fs__BStructuredColumns__global_index(const detail::BlockStructuredColumns* This) { + return This->global_index().get(); +} + +field::FieldImpl* atlas__fs__BStructuredColumns__index_i(const detail::BlockStructuredColumns* This) { + return This->index_i().get(); +} + +field::FieldImpl* atlas__fs__BStructuredColumns__index_j(const detail::BlockStructuredColumns* This) { + return This->index_j().get(); +} + +idx_t atlas__fs__BStructuredColumns__size(const detail::BlockStructuredColumns* This) { + return This->size(); +} + +idx_t atlas__fs__BStructuredColumns__sizeOwned(const detail::BlockStructuredColumns* This) { + return This->structuredcolumns().sizeOwned(); +} + +idx_t atlas__fs__BStructuredColumns__levels(const detail::BlockStructuredColumns* This) { + return This->levels(); +} + +const GridImpl* atlas__fs__BStructuredColumns__grid(const detail::BlockStructuredColumns* This) { + return This->grid().get(); +} +} + + +// ---------------------------------------------------------------------------- + +} // namespace functionspace +} // namespace atlas diff --git a/src/atlas/functionspace/detail/BlockStructuredColumnsInterface.h b/src/atlas/functionspace/detail/BlockStructuredColumnsInterface.h new file mode 100644 index 000000000..7428afdf1 --- /dev/null +++ b/src/atlas/functionspace/detail/BlockStructuredColumnsInterface.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 "atlas/functionspace/BlockStructuredColumns.h" + + +namespace atlas { +namespace field { +class FieldSetImpl; +} +namespace grid { +class DistributionImpl; +} +} // namespace atlas + +namespace atlas { +namespace grid { +namespace detail { +namespace grid { +class Grid; +} // namespace grid +} // namespace detail +} // namespace grid +using GridImpl = grid::detail::grid::Grid; +} // namespace atlas + +namespace atlas { +namespace grid { +namespace detail { +namespace partitioner { +class Partitioner; +} // namespace partitioner +} // namespace detail +} // namespace grid +using PartitionerImpl = grid::detail::partitioner::Partitioner; +} // namespace atlas + + +namespace atlas { +namespace functionspace { + +// ------------------------------------------------------------------- +// C wrapper interfaces to C++ routines +extern "C" { + +const detail::BlockStructuredColumns* atlas__functionspace__BStructuredColumns__new__grid(const GridImpl* grid, + const eckit::Configuration* config); + +const detail::BlockStructuredColumns* atlas__functionspace__BStructuredColumns__new__grid_dist( + const GridImpl* grid, const grid::DistributionImpl* dist, const eckit::Configuration* config); + +const detail::BlockStructuredColumns* atlas__functionspace__BStructuredColumns__new__grid_dist_vert( + const GridImpl* grid, const grid::DistributionImpl* dist, const Vertical* vert, const eckit::Configuration* config); +const detail::BlockStructuredColumns* atlas__functionspace__BStructuredColumns__new__grid_dist_config( + const GridImpl* grid, const grid::DistributionImpl* dist, const eckit::Configuration* config); + +const detail::BlockStructuredColumns* atlas__functionspace__BStructuredColumns__new__grid_part( + const GridImpl* grid, const PartitionerImpl* dist, const eckit::Configuration* config); + +const detail::BlockStructuredColumns* atlas__functionspace__BStructuredColumns__new__grid_part_vert( + const GridImpl* grid, const PartitionerImpl* dist, const Vertical* vert, const eckit::Configuration* config); + +void atlas__functionspace__BStructuredColumns__delete(detail::BlockStructuredColumns* This); +field::FieldImpl* atlas__fs__BStructuredColumns__create_field(const detail::BlockStructuredColumns* This, + const eckit::Configuration* options); + +void atlas__functionspace__BStructuredColumns__gather_field(const detail::BlockStructuredColumns* This, + const field::FieldImpl* local, field::FieldImpl* global); +void atlas__functionspace__BStructuredColumns__scatter_field(const detail::BlockStructuredColumns* This, + const field::FieldImpl* global, field::FieldImpl* local); +void atlas__functionspace__BStructuredColumns__gather_fieldset(const detail::BlockStructuredColumns* This, + const field::FieldSetImpl* local, + field::FieldSetImpl* global); +void atlas__functionspace__BStructuredColumns__scatter_fieldset(const detail::BlockStructuredColumns* This, + const field::FieldSetImpl* global, + field::FieldSetImpl* local); +void atlas__fs__BStructuredColumns__checksum_fieldset(const detail::BlockStructuredColumns* This, + const field::FieldSetImpl* fieldset, char*& checksum, idx_t& size, + int& allocated); +void atlas__fs__BStructuredColumns__checksum_field(const detail::BlockStructuredColumns* This, const field::FieldImpl* field, + char*& checksum, idx_t& size, int& allocated); +void atlas__fs__BStructuredColumns__index_host(const detail::BlockStructuredColumns* This, idx_t*& data, idx_t& i_min, + idx_t& i_max, idx_t& j_min, idx_t& j_max); +idx_t atlas__fs__BStructuredColumns__size(const detail::BlockStructuredColumns* This); +idx_t atlas__fs__BStructuredColumns__sizeOwned(const detail::BlockStructuredColumns* This); +idx_t atlas__fs__BStructuredColumns__j_begin(const detail::BlockStructuredColumns* This); +idx_t atlas__fs__BStructuredColumns__j_end(const detail::BlockStructuredColumns* This); +idx_t atlas__fs__BStructuredColumns__i_begin(const detail::BlockStructuredColumns* This, idx_t j); +idx_t atlas__fs__BStructuredColumns__i_end(const detail::BlockStructuredColumns* This, idx_t j); +idx_t atlas__fs__BStructuredColumns__j_begin_halo(const detail::BlockStructuredColumns* This); +idx_t atlas__fs__BStructuredColumns__j_end_halo(const detail::BlockStructuredColumns* This); +idx_t atlas__fs__BStructuredColumns__i_begin_halo(const detail::BlockStructuredColumns* This, idx_t j); +idx_t atlas__fs__BStructuredColumns__i_end_halo(const detail::BlockStructuredColumns* This, idx_t j); +idx_t atlas__fs__BStructuredColumns__levels(const detail::BlockStructuredColumns* This); +idx_t atlas__fs__BStructuredColumns__block_begin(const detail::BlockStructuredColumns* This, idx_t jblk); +idx_t atlas__fs__BStructuredColumns__block_size(const detail::BlockStructuredColumns* This, idx_t jblk); +idx_t atlas__fs__BStructuredColumns__nblks(const detail::BlockStructuredColumns* This); + +field::FieldImpl* atlas__fs__BStructuredColumns__xy(const detail::BlockStructuredColumns* This); +field::FieldImpl* atlas__fs__BStructuredColumns__z(const detail::BlockStructuredColumns* This); +field::FieldImpl* atlas__fs__BStructuredColumns__partition(const detail::BlockStructuredColumns* This); +field::FieldImpl* atlas__fs__BStructuredColumns__global_index(const detail::BlockStructuredColumns* This); +field::FieldImpl* atlas__fs__BStructuredColumns__index_i(const detail::BlockStructuredColumns* This); +field::FieldImpl* atlas__fs__BStructuredColumns__index_j(const detail::BlockStructuredColumns* This); + +const GridImpl* atlas__fs__BStructuredColumns__grid(const detail::BlockStructuredColumns* This); +} + +} // namespace functionspace +} // namespace atlas diff --git a/src/atlas/functionspace/detail/StructuredColumns.cc b/src/atlas/functionspace/detail/StructuredColumns.cc index 1afc083a7..645c16b7b 100644 --- a/src/atlas/functionspace/detail/StructuredColumns.cc +++ b/src/atlas/functionspace/detail/StructuredColumns.cc @@ -46,11 +46,9 @@ namespace atlas { namespace functionspace { namespace detail { -namespace { - template -array::LocalView make_leveled_view(Field& field) { +array::LocalView StructuredColumns::make_leveled_view(Field& field) const { using namespace array; if (field.levels()) { if (field.variables()) { @@ -71,7 +69,7 @@ array::LocalView make_leveled_view(Field& field) { } template -std::string checksum_3d_field(const parallel::Checksum& checksum, const Field& field) { +std::string StructuredColumns::checksum_3d_field(const parallel::Checksum& checksum, const Field& field) const { auto values = make_leveled_view(field); array::ArrayT surface_field(values.shape(0), values.shape(2)); auto surface = array::make_view(surface_field); @@ -87,9 +85,6 @@ std::string checksum_3d_field(const parallel::Checksum& checksum, const Field& f return checksum.execute(surface.data(), surface_field.stride(0)); } -} // namespace - - class StructuredColumnsHaloExchangeCache : public util::Cache, public grid::detail::grid::GridObserver { private: diff --git a/src/atlas/functionspace/detail/StructuredColumns.h b/src/atlas/functionspace/detail/StructuredColumns.h index c575c318d..210357786 100644 --- a/src/atlas/functionspace/detail/StructuredColumns.h +++ b/src/atlas/functionspace/detail/StructuredColumns.h @@ -298,6 +298,12 @@ class StructuredColumns : public FunctionSpaceImpl { } }; + template + array::LocalView make_leveled_view(Field& field) const; + + template + std::string checksum_3d_field(const parallel::Checksum& checksum, const Field& field) const; + idx_t j_begin_; idx_t j_end_; std::vector i_begin_; @@ -314,8 +320,10 @@ class StructuredColumns : public FunctionSpaceImpl { idx_t nb_partitions_; friend struct StructuredColumnsFortranAccess; + friend struct BlockStructuredColumnsFortranAccess; Map2to1 ij2gp_; + friend class BlockStructuredColumns; void setup(const grid::Distribution& distribution, const eckit::Configuration& config); }; diff --git a/src/atlas_f/CMakeLists.txt b/src/atlas_f/CMakeLists.txt index 7833e5d4c..f015676a5 100644 --- a/src/atlas_f/CMakeLists.txt +++ b/src/atlas_f/CMakeLists.txt @@ -115,6 +115,9 @@ generate_fortran_bindings(FORTRAN_BINDINGS ../atlas/functionspace/detail/Spectra generate_fortran_bindings(FORTRAN_BINDINGS ../atlas/functionspace/detail/StructuredColumnsInterface.h MODULE atlas_functionspace_StructuredColumns_c_binding OUTPUT functionspace_StructuredColumns_c_binding.f90 ) +generate_fortran_bindings(FORTRAN_BINDINGS ../atlas/functionspace/detail/BlockStructuredColumnsInterface.h + MODULE atlas_functionspace_BlockStructuredColumns_c_binding + OUTPUT functionspace_BlockStructuredColumns_c_binding.f90 ) generate_fortran_bindings(FORTRAN_BINDINGS ../atlas/functionspace/detail/NodeColumnsInterface.h MODULE atlas_functionspace_NodeColumns_c_binding OUTPUT functionspace_NodeColumns_c_binding.f90) @@ -169,6 +172,7 @@ ecbuild_add_library( TARGET atlas_f functionspace/atlas_functionspace_EdgeColumns_module.F90 functionspace/atlas_functionspace_NodeColumns_module.fypp functionspace/atlas_functionspace_StructuredColumns_module.F90 + functionspace/atlas_functionspace_BlockStructuredColumns_module.F90 functionspace/atlas_functionspace_Spectral_module.F90 functionspace/atlas_functionspace_PointCloud_module.F90 field/atlas_FieldSet_module.F90 diff --git a/src/atlas_f/functionspace/atlas_functionspace_BlockStructuredColumns_module.F90 b/src/atlas_f/functionspace/atlas_functionspace_BlockStructuredColumns_module.F90 new file mode 100644 index 000000000..f25f76321 --- /dev/null +++ b/src/atlas_f/functionspace/atlas_functionspace_BlockStructuredColumns_module.F90 @@ -0,0 +1,580 @@ +! (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. + +#include "atlas/atlas_f.h" + +module atlas_functionspace_BlockStructuredColumns_module + +use, intrinsic :: iso_c_binding, only : c_ptr, c_double +use fckit_c_interop_module, only : c_str, c_ptr_to_string, c_ptr_free +use atlas_functionspace_module, only : atlas_FunctionSpace +use atlas_Field_module, only: atlas_Field +use atlas_FieldSet_module, only: atlas_FieldSet +use atlas_Grid_module, only: atlas_Grid +use atlas_Config_module, only: atlas_Config +use atlas_kinds_module, only : ATLAS_KIND_IDX +use fckit_owned_object_module, only : fckit_owned_object +use atlas_GridDistribution_module, only : atlas_GridDistribution +use atlas_Partitioner_module, only : atlas_Partitioner +use atlas_Vertical_module, only : atlas_Vertical + +implicit none + +private :: c_ptr, c_double +private :: c_str, c_ptr_to_string, c_ptr_free +private :: atlas_FunctionSpace +private :: atlas_Field +private :: atlas_FieldSet +private :: atlas_Grid +private :: atlas_GridDistribution +private :: atlas_Partitioner +private :: atlas_Vertical +private :: atlas_Config +private :: fckit_owned_object + +public :: atlas_functionspace_BlockStructuredColumns + +private + +!------------------------------------------------------------------------------ +TYPE, extends(atlas_FunctionSpace) :: atlas_functionspace_BlockStructuredColumns + +! Purpose : +! ------- +! *atlas_functionspace_BlockStructuredColumns* : Interpretes spectral fields + +! Methods : +! ------- + +! Author : +! ------ +! August-2015 Willem Deconinck *ECMWF* + +!------------------------------------------------------------------------------ + integer(ATLAS_KIND_IDX), pointer, public :: index(:,:) + +contains + + procedure, public :: assignment_operator_hook + + procedure, private :: gather_fieldset + procedure, private :: gather_field + generic, public :: gather => gather_fieldset, gather_field + + procedure, private :: scatter_fieldset + procedure, private :: scatter_field + generic, public :: scatter => scatter_fieldset, scatter_field + + procedure, private :: checksum_fieldset + procedure, private :: checksum_field + generic, public :: checksum => checksum_fieldset, checksum_field + + procedure :: j_begin + procedure :: j_end + procedure :: i_begin + procedure :: i_end + procedure :: j_begin_halo + procedure :: j_end_halo + procedure :: i_begin_halo + procedure :: i_end_halo + + procedure :: size => get_size + procedure :: size_owned => get_size_owned + + procedure :: levels + procedure :: block_begin + procedure :: block_size + procedure :: nblks + + procedure :: xy + !! Return xy coordinate field + procedure :: z + !! Return z coordinate field + procedure :: partition + !! Return partition field + procedure :: global_index + !! Return global_index field + procedure :: index_i + !! Return index_i field + procedure :: index_j + !! Return index_j field + + procedure :: grid + + procedure, private :: set_index + +#if FCKIT_FINAL_NOT_INHERITING + final :: BStructuredColumns__final_auto +#endif + +END TYPE atlas_functionspace_BlockStructuredColumns + +interface atlas_functionspace_BlockStructuredColumns + module procedure ctor_cptr + module procedure ctor_grid + module procedure ctor_grid_config + module procedure ctor_grid_dist + module procedure ctor_grid_dist_config + module procedure ctor_grid_dist_levels + module procedure ctor_grid_dist_vertical + module procedure ctor_grid_part + module procedure ctor_grid_part_levels + module procedure ctor_grid_part_vertical +end interface + + +!------------------------------------------------------------------------------ + +!======================================================== +contains +!======================================================== + +subroutine assignment_operator_hook(this,other) + class(atlas_functionspace_BlockStructuredColumns) :: this + class(fckit_owned_object) :: other + call this%set_index() + FCKIT_SUPPRESS_UNUSED(other) +end subroutine + +function ctor_cptr(cptr) result(this) + type(atlas_functionspace_BlockStructuredColumns) :: this + type(c_ptr), intent(in) :: cptr + call this%reset_c_ptr( cptr ) + call this%set_index() + call this%return() +end function + +function empty_config() result(config) + type(atlas_Config) :: config + config = atlas_Config() + call config%return() +end function + +function ctor_grid(grid, halo, nproma, levels) result(this) + use atlas_functionspace_BlockStructuredColumns_c_binding + type(atlas_functionspace_BlockStructuredColumns) :: this + class(atlas_Grid), intent(in) :: grid + integer, optional :: halo + integer, optional :: nproma + integer, optional :: levels + type(atlas_Config) :: config + config = empty_config() ! Due to PGI compiler bug, we have to do this instead of "config = atlas_Config()"" + if( present(halo) ) call config%set("halo",halo) + if( present(nproma) ) call config%set("nproma",nproma) + if( present(levels) ) call config%set("levels",levels) + call this%reset_c_ptr( atlas__functionspace__BStructuredColumns__new__grid( grid%CPTR_PGIBUG_A, & + & config%CPTR_PGIBUG_B ) ) + call this%set_index() + call config%final() + call this%return() +end function + +function ctor_grid_config(grid, config) result(this) + use atlas_functionspace_BlockStructuredColumns_c_binding + type(atlas_functionspace_BlockStructuredColumns) :: this + class(atlas_Grid), intent(in) :: grid + type(atlas_Config), intent (in) :: config + call this%reset_c_ptr(atlas__functionspace__BStructuredColumns__new__grid( & + & grid%CPTR_PGIBUG_A, & + & config%CPTR_PGIBUG_B ) ) + call this%set_index() + call this%return() +end function + +function ctor_grid_dist(grid, distribution, halo, levels) result(this) + use atlas_functionspace_BlockStructuredColumns_c_binding + type(atlas_functionspace_BlockStructuredColumns) :: this + class(atlas_Grid), intent(in) :: grid + type(atlas_griddistribution), intent(in) :: distribution + integer, optional :: halo + integer, optional :: levels + type(atlas_Config) :: config + config = empty_config() ! Due to PGI compiler bug, we have to do this instead of "config = atlas_Config()"" + if( present(halo) ) call config%set("halo",halo) + if( present(levels) ) call config%set("levels",levels) + call this%reset_c_ptr( atlas__functionspace__BStructuredColumns__new__grid_dist( & + & grid%CPTR_PGIBUG_A, distribution%CPTR_PGIBUG_A, config%CPTR_PGIBUG_B ) ) + call this%set_index() + call config%final() + call this%return() +end function + +function ctor_grid_dist_levels(grid, distribution, levels, halo) result(this) + use atlas_functionspace_BlockStructuredColumns_c_binding + type(atlas_functionspace_BlockStructuredColumns) :: this + class(atlas_Grid), intent(in) :: grid + type(atlas_griddistribution), intent(in) :: distribution + integer, optional :: halo + real(c_double) :: levels(:) + type(atlas_Config) :: config + type(atlas_Vertical) :: vertical + config = empty_config() ! Due to PGI compiler bug, we have to do this insted of "config = atlas_Config()"" + if( present(halo) ) call config%set("halo",halo) + call config%set("levels",size(levels)) + vertical = atlas_Vertical(levels) + call this%reset_c_ptr( atlas__functionspace__BStructuredColumns__new__grid_dist_vert( & + & grid%CPTR_PGIBUG_A, distribution%CPTR_PGIBUG_A, vertical%CPTR_PGIBUG_B, & + & config%CPTR_PGIBUG_B ) ) + call this%set_index() + call config%final() + call vertical%final() + call this%return() +end function + +function ctor_grid_dist_vertical(grid, distribution, vertical, halo) result(this) + use atlas_functionspace_BlockStructuredColumns_c_binding + type(atlas_functionspace_BlockStructuredColumns) :: this + class(atlas_Grid), intent(in) :: grid + type(atlas_griddistribution), intent(in) :: distribution + integer, optional :: halo + type(atlas_Vertical) :: vertical + type(atlas_Config) :: config + config = empty_config() ! Due to PGI compiler bug, we have to do this insted of "config = atlas_Config()"" + if( present(halo) ) call config%set("halo",halo) + call config%set("levels",vertical%size()) + call this%reset_c_ptr( atlas__functionspace__BStructuredColumns__new__grid_dist_vert( & + & grid%CPTR_PGIBUG_A, distribution%CPTR_PGIBUG_A, vertical%CPTR_PGIBUG_B, & + & config%CPTR_PGIBUG_B ) ) + call this%set_index() + call config%final() + call this%return() +end function + +function ctor_grid_dist_config(grid, distribution, config) result(this) + use atlas_functionspace_BlockStructuredColumns_c_binding + type(atlas_functionspace_BlockStructuredColumns) :: this + class(atlas_Grid), intent(in) :: grid + type(atlas_griddistribution), intent(in) :: distribution + type(atlas_Config), intent (in) :: config + call this%reset_c_ptr(atlas__functionspace__BStructuredColumns__new__grid_dist_config( & + & grid%CPTR_PGIBUG_A, distribution%CPTR_PGIBUG_A, & + & config%CPTR_PGIBUG_B ) ) + call this%set_index() + call this%return() +end function + + +function ctor_grid_part(grid, partitioner, halo, levels) result(this) + use atlas_functionspace_BlockStructuredColumns_c_binding + type(atlas_functionspace_BlockStructuredColumns) :: this + class(atlas_Grid), intent(in) :: grid + type(atlas_Partitioner), intent(in) :: partitioner + integer, optional :: halo + integer, optional :: levels + type(atlas_Config) :: config + config = empty_config() ! Due to PGI compiler bug, we have to do this instead of "config = atlas_Config()"" + if( present(halo) ) call config%set("halo",halo) + if( present(levels) ) call config%set("levels",levels) + call this%reset_c_ptr( atlas__functionspace__BStructuredColumns__new__grid_part( & + & grid%CPTR_PGIBUG_A, partitioner%CPTR_PGIBUG_A, config%CPTR_PGIBUG_B ) ) + call this%set_index() + call config%final() + call this%return() +end function + +function ctor_grid_part_levels(grid, partitioner, levels, halo) result(this) + use atlas_functionspace_BlockStructuredColumns_c_binding + type(atlas_functionspace_BlockStructuredColumns) :: this + class(atlas_Grid), intent(in) :: grid + type(atlas_Partitioner), intent(in) :: partitioner + integer, optional :: halo + real(c_double) :: levels(:) + type(atlas_Config) :: config + type(atlas_Vertical) :: vertical + config = empty_config() ! Due to PGI compiler bug, we have to do this insted of "config = atlas_Config()"" + if( present(halo) ) call config%set("halo",halo) + call config%set("levels",size(levels)) + vertical = atlas_Vertical(levels) + call this%reset_c_ptr( atlas__functionspace__BStructuredColumns__new__grid_part_vert( & + & grid%CPTR_PGIBUG_A, partitioner%CPTR_PGIBUG_A, vertical%CPTR_PGIBUG_B, & + & config%CPTR_PGIBUG_B ) ) + call this%set_index() + call config%final() + call vertical%final() + call this%return() +end function + +function ctor_grid_part_vertical(grid, partitioner, vertical, halo) result(this) + use atlas_functionspace_BlockStructuredColumns_c_binding + type(atlas_functionspace_BlockStructuredColumns) :: this + class(atlas_Grid), intent(in) :: grid + type(atlas_Partitioner), intent(in) :: partitioner + integer, optional :: halo + type(atlas_Vertical) :: vertical + type(atlas_Config) :: config + config = empty_config() ! Due to PGI compiler bug, we have to do this insted of "config = atlas_Config()"" + if( present(halo) ) call config%set("halo",halo) + call config%set("levels",vertical%size()) + call this%reset_c_ptr( atlas__functionspace__BStructuredColumns__new__grid_part_vert( & + & grid%CPTR_PGIBUG_A, partitioner%CPTR_PGIBUG_A, vertical%CPTR_PGIBUG_B, & + & config%CPTR_PGIBUG_B ) ) + call this%set_index() + call config%final() + call this%return() +end function + + +subroutine gather_field(this,local,global) + use atlas_functionspace_BlockStructuredColumns_c_binding + class(atlas_functionspace_BlockStructuredColumns), intent(in) :: this + type(atlas_Field), intent(in) :: local + type(atlas_Field), intent(inout) :: global + call atlas__functionspace__BStructuredColumns__gather_field(this%CPTR_PGIBUG_A,local%CPTR_PGIBUG_A,global%CPTR_PGIBUG_A) +end subroutine + +subroutine gather_fieldset(this,local,global) + use atlas_functionspace_BlockStructuredColumns_c_binding + class(atlas_functionspace_BlockStructuredColumns), intent(in) :: this + type(atlas_FieldSet), intent(in) :: local + type(atlas_FieldSet), intent(inout) :: global + call atlas__functionspace__BStructuredColumns__gather_fieldset(this%CPTR_PGIBUG_A,local%CPTR_PGIBUG_A,global%CPTR_PGIBUG_A) +end subroutine + +subroutine scatter_field(this,global,local) + use atlas_functionspace_BlockStructuredColumns_c_binding + class(atlas_functionspace_BlockStructuredColumns), intent(in) :: this + type(atlas_Field), intent(in) :: global + type(atlas_Field), intent(inout) :: local + call atlas__functionspace__BStructuredColumns__scatter_field(this%CPTR_PGIBUG_A,global%CPTR_PGIBUG_A,local%CPTR_PGIBUG_A) +end subroutine + +subroutine scatter_fieldset(this,global,local) + use atlas_functionspace_BlockStructuredColumns_c_binding + class(atlas_functionspace_BlockStructuredColumns), intent(in) :: this + type(atlas_FieldSet), intent(in) :: global + type(atlas_FieldSet), intent(inout) :: local + call atlas__functionspace__BStructuredColumns__scatter_fieldset(this%CPTR_PGIBUG_A,global%CPTR_PGIBUG_A,local%CPTR_PGIBUG_A) +end subroutine + +function checksum_fieldset(this,fieldset) result(checksum) + use, intrinsic :: iso_c_binding + use atlas_functionspace_BlockStructuredColumns_c_binding + character(len=:), allocatable :: checksum + class(atlas_functionspace_BlockStructuredColumns), intent(in) :: this + type(atlas_FieldSet), intent(in) :: fieldset + type(c_ptr) :: checksum_cptr + integer(ATLAS_KIND_IDX) :: checksum_size + integer(c_int) :: checksum_allocated + call atlas__fs__BStructuredColumns__checksum_fieldset( & + & this%CPTR_PGIBUG_A,fieldset%CPTR_PGIBUG_A,checksum_cptr,checksum_size,checksum_allocated) + allocate(character(len=checksum_size) :: checksum ) + checksum = c_ptr_to_string(checksum_cptr) + if( checksum_allocated == 1 ) call c_ptr_free(checksum_cptr) +end function + + +function checksum_field(this,field) result(checksum) + use, intrinsic :: iso_c_binding + use atlas_functionspace_BlockStructuredColumns_c_binding + character(len=:), allocatable :: checksum + class(atlas_functionspace_BlockStructuredColumns), intent(in) :: this + type(atlas_Field), intent(in) :: field + type(c_ptr) :: checksum_cptr + integer(ATLAS_KIND_IDX) :: checksum_size + integer(c_int) :: checksum_allocated + call atlas__fs__BStructuredColumns__checksum_field( & + & this%CPTR_PGIBUG_A,field%CPTR_PGIBUG_A,checksum_cptr,checksum_size,checksum_allocated) + allocate(character(len=checksum_size) :: checksum ) + checksum = c_ptr_to_string(checksum_cptr) + if( checksum_allocated == 1 ) call c_ptr_free(checksum_cptr) +end function + +subroutine set_index(this) + use, intrinsic :: iso_c_binding, only : c_ptr, c_f_pointer + use atlas_functionspace_BlockStructuredColumns_c_binding + class(atlas_functionspace_BlockStructuredColumns), intent(inout) :: this + type(c_ptr) :: index_cptr + integer(ATLAS_KIND_IDX), pointer :: index_fptr(:) + integer(ATLAS_KIND_IDX) :: i_min, i_max, j_min, j_max + integer(ATLAS_KIND_IDX) :: ni, nj + call atlas__fs__BStructuredColumns__index_host( this%CPTR_PGIBUG_A, index_cptr, i_min, i_max, j_min, j_max ) + ni = i_max-i_min+1; + nj = j_max-j_min+1; + call c_f_pointer( index_cptr, index_fptr, (/ni*nj/) ) + this%index(i_min:i_max,j_min:j_max) => index_fptr(:) +end subroutine + +function j_begin(this) result(j) + use atlas_functionspace_BlockStructuredColumns_c_binding + integer(ATLAS_KIND_IDX) :: j + class(atlas_functionspace_BlockStructuredColumns), intent(in) :: this + j = atlas__fs__BStructuredColumns__j_begin(this%CPTR_PGIBUG_A) +end function + +function j_end(this) result(j) + use atlas_functionspace_BlockStructuredColumns_c_binding + integer(ATLAS_KIND_IDX) :: j + class(atlas_functionspace_BlockStructuredColumns), intent(in) :: this + j = atlas__fs__BStructuredColumns__j_end(this%CPTR_PGIBUG_A) +end function + +function i_begin(this,j) result(i) + use atlas_functionspace_BlockStructuredColumns_c_binding + integer(ATLAS_KIND_IDX) :: i + integer(ATLAS_KIND_IDX), intent(in) :: j + class(atlas_functionspace_BlockStructuredColumns), intent(in) :: this + i = atlas__fs__BStructuredColumns__i_begin(this%CPTR_PGIBUG_A,j) +end function + +function i_end(this,j) result(i) + use atlas_functionspace_BlockStructuredColumns_c_binding + integer(ATLAS_KIND_IDX) :: i + integer(ATLAS_KIND_IDX), intent(in) :: j + class(atlas_functionspace_BlockStructuredColumns), intent(in) :: this + i = atlas__fs__BStructuredColumns__i_end(this%CPTR_PGIBUG_A,j) +end function + + +function j_begin_halo(this) result(j) + use atlas_functionspace_BlockStructuredColumns_c_binding + integer(ATLAS_KIND_IDX) :: j + class(atlas_functionspace_BlockStructuredColumns), intent(in) :: this + j = atlas__fs__BStructuredColumns__j_begin_halo(this%CPTR_PGIBUG_A) +end function + +function j_end_halo(this) result(j) + use atlas_functionspace_BlockStructuredColumns_c_binding + integer(ATLAS_KIND_IDX) :: j + class(atlas_functionspace_BlockStructuredColumns), intent(in) :: this + j = atlas__fs__BStructuredColumns__j_end_halo(this%CPTR_PGIBUG_A) +end function + +function i_begin_halo(this,j) result(i) + use atlas_functionspace_BlockStructuredColumns_c_binding + integer(ATLAS_KIND_IDX) :: i + integer(ATLAS_KIND_IDX), intent(in) :: j + class(atlas_functionspace_BlockStructuredColumns), intent(in) :: this + i = atlas__fs__BStructuredColumns__i_begin_halo(this%CPTR_PGIBUG_A,j) +end function + +function i_end_halo(this,j) result(i) + use atlas_functionspace_BlockStructuredColumns_c_binding + integer(ATLAS_KIND_IDX) :: i + integer(ATLAS_KIND_IDX), intent(in) :: j + class(atlas_functionspace_BlockStructuredColumns), intent(in) :: this + i = atlas__fs__BStructuredColumns__i_end_halo(this%CPTR_PGIBUG_A,j) +end function + +function get_size(this) result(size) + use atlas_functionspace_BlockStructuredColumns_c_binding + integer(ATLAS_KIND_IDX) :: size + class(atlas_functionspace_BlockStructuredColumns), intent(in) :: this + size = atlas__fs__BStructuredColumns__size(this%CPTR_PGIBUG_A) +end function + +function get_size_owned(this) result(size) + use atlas_functionspace_BlockStructuredColumns_c_binding + integer(ATLAS_KIND_IDX) :: size + class(atlas_functionspace_BlockStructuredColumns), intent(in) :: this + size = atlas__fs__BStructuredColumns__sizeOwned(this%CPTR_PGIBUG_A) +end function + +function levels(this) + use atlas_functionspace_BlockStructuredColumns_c_binding + integer(ATLAS_KIND_IDX) :: levels + class(atlas_functionspace_BlockStructuredColumns), intent(in) :: this + levels = atlas__fs__BStructuredColumns__levels(this%CPTR_PGIBUG_A) +end function + +function xy(this) result(field) + use atlas_functionspace_BlockStructuredColumns_c_binding + type(atlas_Field) :: field + class(atlas_functionspace_BlockStructuredColumns), intent(in) :: this + field = atlas_Field( atlas__fs__BStructuredColumns__xy(this%CPTR_PGIBUG_A) ) + call field%return() +end function + +function z(this) result(field) + use atlas_functionspace_BlockStructuredColumns_c_binding + type(atlas_Field) :: field + class(atlas_functionspace_BlockStructuredColumns), intent(in) :: this + field = atlas_Field( atlas__fs__BStructuredColumns__z(this%CPTR_PGIBUG_A) ) + call field%return() +end function + + +function partition(this) result(field) + use atlas_functionspace_BlockStructuredColumns_c_binding + type(atlas_Field) :: field + class(atlas_functionspace_BlockStructuredColumns), intent(in) :: this + field = atlas_Field( atlas__fs__BStructuredColumns__partition(this%CPTR_PGIBUG_A) ) + call field%return() +end function + +function global_index(this) result(field) + use atlas_functionspace_BlockStructuredColumns_c_binding + type(atlas_Field) :: field + class(atlas_functionspace_BlockStructuredColumns), intent(in) :: this + field = atlas_Field( atlas__fs__BStructuredColumns__global_index(this%CPTR_PGIBUG_A) ) + call field%return() +end function + +function index_i(this) result(field) + use atlas_functionspace_BlockStructuredColumns_c_binding + type(atlas_Field) :: field + class(atlas_functionspace_BlockStructuredColumns), intent(in) :: this + field = atlas_Field( atlas__fs__BStructuredColumns__index_i(this%CPTR_PGIBUG_A) ) + call field%return() +end function + +function index_j(this) result(field) + use atlas_functionspace_BlockStructuredColumns_c_binding + type(atlas_Field) :: field + class(atlas_functionspace_BlockStructuredColumns), intent(in) :: this + field = atlas_Field( atlas__fs__BStructuredColumns__index_j(this%CPTR_PGIBUG_A) ) + call field%return() +end function + +function grid(this) + use atlas_functionspace_BlockStructuredColumns_c_binding + type(atlas_Grid) :: grid + class(atlas_functionspace_BlockStructuredColumns), intent(in) :: this + grid = atlas_Grid( atlas__fs__BStructuredColumns__grid(this%CPTR_PGIBUG_A) ) + call grid%return() +end function + +function block_begin(this,j) result(i) + use atlas_functionspace_BlockStructuredColumns_c_binding + integer(ATLAS_KIND_IDX) :: i + integer(ATLAS_KIND_IDX), intent(in) :: j + class(atlas_functionspace_BlockStructuredColumns), intent(in) :: this + i = atlas__fs__BStructuredColumns__block_begin(this%CPTR_PGIBUG_A,j-1) + 1 +end function + +function block_size(this,j) result(i) + use atlas_functionspace_BlockStructuredColumns_c_binding + integer(ATLAS_KIND_IDX) :: i + integer(ATLAS_KIND_IDX), intent(in) :: j + class(atlas_functionspace_BlockStructuredColumns), intent(in) :: this + i = atlas__fs__BStructuredColumns__block_size(this%CPTR_PGIBUG_A,j-1) +end function + +function nblks(this) result(i) + use atlas_functionspace_BlockStructuredColumns_c_binding + integer(ATLAS_KIND_IDX) :: i + class(atlas_functionspace_BlockStructuredColumns), intent(in) :: this + i = atlas__fs__BStructuredColumns__nblks(this%CPTR_PGIBUG_A) +end function + +!------------------------------------------------------------------------------- + +#if FCKIT_FINAL_NOT_INHERITING +ATLAS_FINAL subroutine BStructuredColumns__final_auto(this) + type(atlas_functionspace_BlockStructuredColumns), intent(inout) :: this +#if FCKIT_FINAL_DEBUGGING + write(0,*) "atlas_functionspace_BlockStructuredColumns__final_auto" +#endif +#if FCKIT_FINAL_NOT_PROPAGATING + call this%final() +#endif + FCKIT_SUPPRESS_UNUSED( this ) +end subroutine +#endif + +end module atlas_functionspace_BlockStructuredColumns_module + diff --git a/src/tests/functionspace/CMakeLists.txt b/src/tests/functionspace/CMakeLists.txt index 7d81a7b9d..318a7f285 100644 --- a/src/tests/functionspace/CMakeLists.txt +++ b/src/tests/functionspace/CMakeLists.txt @@ -22,6 +22,15 @@ if( HAVE_FCTEST ) ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} ) + add_fctest( TARGET atlas_fctest_blockstructuredcolumns + MPI 4 + CONDITION eckit_HAVE_MPI + LINKER_LANGUAGE Fortran + SOURCES fctest_blockstructuredcolumns.F90 + LIBS atlas_f + ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} +) + endif() ecbuild_add_test( TARGET atlas_test_functionspace @@ -67,6 +76,14 @@ ecbuild_add_test( TARGET atlas_test_structuredcolumns CONDITION eckit_HAVE_MPI AND MPI_SLOTS GREATER_EQUAL 5 ) +ecbuild_add_test( TARGET atlas_test_blockstructuredcolumns + MPI 3 + SOURCES test_blockstructuredcolumns.cc + LIBS atlas + ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} + CONDITION eckit_HAVE_MPI AND MPI_SLOTS GREATER_EQUAL 3 +) + ecbuild_add_test( TARGET atlas_test_pointcloud SOURCES test_pointcloud.cc LIBS atlas diff --git a/src/tests/functionspace/fctest_blockstructuredcolumns.F90 b/src/tests/functionspace/fctest_blockstructuredcolumns.F90 new file mode 100644 index 000000000..708f6e60e --- /dev/null +++ b/src/tests/functionspace/fctest_blockstructuredcolumns.F90 @@ -0,0 +1,210 @@ +! (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 File contains Unit Tests for testing the +! C++ / Fortran Interfaces to the State Datastructure +! @author Willem Deconinck +! @author Slavko Brdar + +#include "fckit/fctest.h" + +! ----------------------------------------------------------------------------- + +module fcta_BlockStructuredColumns_fxt +use atlas_module +use, intrinsic :: iso_c_binding +implicit none +character(len=1024) :: msg +contains + +end module + +! ----------------------------------------------------------------------------- + +TESTSUITE_WITH_FIXTURE(fcta_BlockStructuredColumns,fcta_BlockStructuredColumns_fxt) + +! ----------------------------------------------------------------------------- + +TESTSUITE_INIT + call atlas_library%initialise() +END_TESTSUITE_INIT + +! ----------------------------------------------------------------------------- + +TESTSUITE_FINALIZE + call atlas_library%finalise() +END_TESTSUITE_FINALIZE + +TEST( test_blockstructuredcolumns ) +use atlas_functionspace_blockstructuredcolumns_module +implicit none +type(atlas_StructuredGrid) :: grid +type(atlas_functionspace_BlockStructuredColumns) :: fs +type(atlas_functionspace) :: fs_base +integer(ATLAS_KIND_IDX) :: i, j, jbegin, jblk, jrof +character(len=256) str + +type(atlas_Field) :: field +type(atlas_Field) :: field_xy +type(atlas_Field) :: field_global_index +type(atlas_Field) :: field_index_j +real(8), pointer :: xy(:,:), x(:,:) +integer(ATLAS_KIND_GIDX), pointer :: global_index(:) +integer(ATLAS_KIND_IDX), pointer :: index_j(:) +integer(ATLAS_KIND_IDX) :: nblks, blksize +integer, parameter :: XX=1 + +grid = atlas_StructuredGrid("O16") +fs = atlas_functionspace_BlockStructuredColumns(grid,halo=0,nproma=12) + +field = fs%create_field(name="field",kind=atlas_real(8)) +FCTEST_CHECK_EQUAL( field%owners(), 1 ) +field_xy = fs%xy() +FCTEST_CHECK_EQUAL( field_xy%owners(), 2 ) +call field%data(x) +call field_xy%data(xy) +field_global_index = fs%global_index() +call field_global_index%data(global_index) + +field_index_j = fs%index_j() +call field_index_j%data(index_j) + +nblks=fs%nblks() +do jblk=1,nblks + write(str,'(A,I4,A)') 'block', jblk, ' : ' + call atlas_log%info(str,newl=.false.) + jbegin = fs%block_begin(jblk) + blksize = fs%block_size(jblk) + do jrof=1,blksize + write(str,'(I6)') global_index( jbegin+jrof-1 ) + call atlas_log%info(str,newl=.false.) + x(jrof,jblk) = xy(XX, jbegin+jrof-1) + enddo + call atlas_log%info("",newl=.true.) +enddo + +!call fs%halo_exchange(field) + +FCTEST_CHECK_EQUAL( field_xy%owners(), 2 ) +fs = atlas_functionspace_BlockStructuredColumns(grid,levels=5) + +field = fs%create_field(atlas_real(c_float)) +FCTEST_CHECK_EQUAL( field%rank() , 3 ) +FCTEST_CHECK_EQUAL( field%name() , "" ) +FCTEST_CHECK_EQUAL( field%kind() , atlas_real(c_float) ) + + +write(str,*) "before: name = ", fs%name(); call atlas_log%info(str) +write(str,*) "before: owners = ", fs%owners(); call atlas_log%info(str) +FCTEST_CHECK_EQUAL( fs%name(), "BlockStructuredColumns" ) +fs_base = field%functionspace(); call atlas_log%info(str) +write(str,*) "after: name = " , fs_base%name(); call atlas_log%info(str) +write(str,*) "after: owners = " , fs_base%owners(); call atlas_log%info(str) +FCTEST_CHECK_EQUAL( fs_base%name(), "BlockStructuredColumns" ) + +FCTEST_CHECK_EQUAL( field%owners(), 1 ) +call field%final() +FCTEST_CHECK_EQUAL( field_xy%owners(), 1 ) +call field_xy%final() +call field_global_index%final() +call fs%final() +call fs_base%final() +call grid%final() +END_TEST + +! ----------------------------------------------------------------------------- + +TEST( test_collectives ) +use fckit_mpi_module +use atlas_functionspace_blockstructuredcolumns_module +implicit none +type(atlas_StructuredGrid) :: grid +type(atlas_functionspace_BlockStructuredColumns) :: fs +type(atlas_Field) :: field, global, scal +type(atlas_Config) :: config +type(atlas_Metadata) :: metadata +type(fckit_mpi_comm) :: mpi +integer(c_int) :: test_broadcast, halo_size, levels, nproma, nvar +call atlas_log%info("test_collectives") +mpi = fckit_mpi_comm() +halo_size = 0 +levels = 10 +nproma = 12 +nvar = 2 +FCTEST_CHECK_EQUAL( halo_size, 0 ) +grid = atlas_StructuredGrid("O16") +config = atlas_Config() +fs = atlas_functionspace_BlockStructuredColumns(grid,nproma=nproma,levels=levels,halo=halo_size) + +! type c_int +field = fs%create_field(kind=atlas_integer(c_int),variables=2) +global = fs%create_field(field,global=.True.) +scal = fs%create_field(kind=atlas_integer(c_int)) +call fs%gather(field,global) +metadata = global%metadata() +if( mpi%rank() == 0 ) then + call metadata%set("test_broadcast",123) +endif +call fs%scatter(global,field) +metadata = field%metadata() +call metadata%get("test_broadcast",test_broadcast) +FCTEST_CHECK_EQUAL( test_broadcast, 123 ) + +! type c_long +field = fs%create_field(kind=c_long,variables=nvar) +global = fs%create_field(field,global=.True.) +scal = fs%create_field(kind=c_long) +call fs%gather(field,global) +metadata = global%metadata() +if( mpi%rank() == 0 ) then + call metadata%set("test_broadcast",123) +endif +call fs%scatter(global,field) +metadata = field%metadata() +call metadata%get("test_broadcast",test_broadcast) +FCTEST_CHECK_EQUAL( test_broadcast, 123 ) + +! type c_float +field = fs%create_field(kind=atlas_real(c_float),variables=nvar) +global = fs%create_field(field,global=.True.) +scal = fs%create_field(kind=atlas_real(c_float)) +write(msg,*) "field: rank",field%rank(), " size ", field%size(), " shape [",field%shape(), "]" +call atlas_log%info(msg) +write(msg,*) "nblk = ",fs%nblks(); call atlas_log%info(msg) +write(msg,*) "global: rank",global%rank()," size ", global%size(), " shape [",global%shape(),"]" +call atlas_log%info(msg) +call fs%gather(field,global) +!call fs%halo_exchange(field) +metadata = global%metadata() +if( mpi%rank() == 0 ) then + call metadata%set("test_broadcast",123) + FCTEST_CHECK_EQUAL(global%shape(3), grid%size() ) + FCTEST_CHECK_EQUAL(global%shape(2), levels) + FCTEST_CHECK_EQUAL(global%shape(1), nvar ) +endif +call fs%scatter(global,field) +metadata = field%metadata() +call metadata%get("test_broadcast",test_broadcast) +FCTEST_CHECK_EQUAL( test_broadcast, 123 ) + +! type c_double +field = fs%create_field(kind=atlas_real(c_double),variables=2) +global = fs%create_field(field,global=.True.) +scal = fs%create_field(kind=atlas_real(c_double)) +call fs%gather(field,global) +metadata = global%metadata() +if( mpi%rank() == 0 ) then + call metadata%set("test_broadcast",123) +endif +call fs%scatter(global,field) +metadata = field%metadata() +call metadata%get("test_broadcast",test_broadcast) +FCTEST_CHECK_EQUAL( test_broadcast, 123 ) +END_TEST + +END_TESTSUITE + diff --git a/src/tests/functionspace/fctest_functionspace.F90 b/src/tests/functionspace/fctest_functionspace.F90 index e46053042..a705790a6 100644 --- a/src/tests/functionspace/fctest_functionspace.F90 +++ b/src/tests/functionspace/fctest_functionspace.F90 @@ -569,8 +569,8 @@ module fcta_FunctionSpace_fxt write(0,*) "before: name = ", fs%name() write(0,*) "before: owners = ", fs%owners() fs_base = field%functionspace() -write(0,*) "after: name = " , fs%name() -write(0,*) "after: owners = " , fs%owners() +write(0,*) "after: name = " , fs_base%name() +write(0,*) "after: owners = " , fs_base%owners() FCTEST_CHECK_EQUAL( field%owners(), 1 ) call field%final() diff --git a/src/tests/functionspace/test_blockstructuredcolumns.cc b/src/tests/functionspace/test_blockstructuredcolumns.cc new file mode 100644 index 000000000..953dd9bd2 --- /dev/null +++ b/src/tests/functionspace/test_blockstructuredcolumns.cc @@ -0,0 +1,198 @@ +/* + * (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. + */ + +#include "atlas/array/ArrayView.h" +#include "atlas/array/MakeView.h" +#include "atlas/field/Field.h" +#include "atlas/functionspace/BlockStructuredColumns.h" +#include "atlas/grid/Partitioner.h" +#include "atlas/grid/StructuredGrid.h" +#include "atlas/parallel/mpi/mpi.h" + +#include "tests/AtlasTestEnvironment.h" + +using namespace eckit; +using namespace atlas; +using namespace atlas::functionspace; +using namespace atlas::test; +using namespace atlas::util; + +namespace { +template +void run_scatter_gather(const Grid& grid, const BlockStructuredColumns& fs, int nlev, int nvar) { + auto glb_field = fs.createField( option::global() | option::levels(nlev) | option::variables(nvar) ); + // Only allocated on rank 0 + if( atlas::mpi::comm().rank() != 0 ) { + EXPECT_EQ( glb_field.shape(0), 0); + } + else { + EXPECT_EQ( glb_field.shape(0), grid.size()); + } + + + auto glb_field_v = array::make_view(glb_field); + + for (gidx_t p = 0; p < glb_field.shape(0); p++) { + for (idx_t jlev = 0; jlev < nlev; ++jlev) { + for (idx_t jvar = 0; jvar < nvar; ++jvar) { + glb_field_v(p, jlev, jvar) = p*nvar*nlev + jlev*nvar + jvar; + } + } + } + + auto field = fs.createField(glb_field); + + fs.scatter(glb_field, field); + + auto glb_field_2 = fs.createField(glb_field , option::global(atlas::mpi::comm().size()-1)); + // Only allocated on rank MPI_SIZE-1 + if( atlas::mpi::comm().rank() != atlas::mpi::comm().size() - 1 ) { + EXPECT_EQ( glb_field_2.shape(0), 0); + } + else { + EXPECT_EQ( glb_field_2.shape(0), grid.size()); + } + + fs.gather(field, glb_field_2); + + auto glb_field_2_v = array::make_view(glb_field_2); + for (gidx_t p = 0; p < glb_field_2.shape(0); p++) { + for (idx_t jlev = 0; jlev < nlev; ++jlev) { + for (idx_t jvar = 0; jvar < nvar; ++jvar) { + EXPECT_EQ(glb_field_2_v(p, jlev, jvar), p*nvar*nlev + jlev*nvar + jvar); + } + } + } +} + +} // namespace + +namespace atlas { +namespace test { + +//----------------------------------------------------------------------------- + + +CASE("test_BlockStructuredColumns") { + std::string gridname = eckit::Resource("--grid", "O8"); + idx_t halo = eckit::Resource("--halo", 0); + idx_t nproma = eckit::Resource("--nproma", 12); + idx_t nvar = 3; + idx_t nlev = 4; + + auto grid = StructuredGrid(gridname); + + util::Config config; + config.set("halo", halo); + config.set("nproma", nproma); + config.set("levels", nlev); + auto fs = functionspace::BlockStructuredColumns(grid, config); + + SECTION("blocking indices parallel without halo") { + ATLAS_DEBUG_VAR(mpi::comm().size()); + ATLAS_DEBUG_VAR(halo); + ATLAS_DEBUG_VAR(fs.size()); + ATLAS_DEBUG_VAR(fs.nproma()); + ATLAS_DEBUG_VAR(fs.nblks()); + ATLAS_DEBUG_VAR(fs.block(fs.nblks()-1).size()); + ATLAS_DEBUG_VAR(fs.levels()); + ATLAS_DEBUG_VAR(nvar); + + idx_t ngptot = fs.size(); + idx_t nblk = ngptot / nproma; + idx_t nrof_last = nproma; + if (ngptot % nproma > 0) { + // Correction when ngptot is not divisible by nproma + nrof_last = ngptot - nblk * nproma; + ++nblk; + } + EXPECT_EQ(fs.nproma(), nproma); + EXPECT_EQ(fs.levels(), nlev); + EXPECT_EQ(fs.nblks(), nblk); + EXPECT_EQ(fs.block(nblk-1).size(), nrof_last); + + ATLAS_DEBUG("Test cover full iteration space"); + idx_t jpoint = 0; + for (idx_t jblk = 0; jblk < nblk; ++jblk) { + auto blk = fs.block(jblk); + for (idx_t jrof = 0; jrof < blk.size(); ++jrof, ++jpoint) { + idx_t h_idx = fs.index(jblk, jrof); + idx_t hh_idx = blk.index(jrof); + EXPECT_EQ(h_idx, jpoint); + EXPECT_EQ(h_idx, hh_idx); + } + } + EXPECT_EQ(jpoint, ngptot); + + Field field = fs.createField(option::name("field") | option::variables(nvar)); + EXPECT_EQ(field.shape(0), nblk ); + EXPECT_EQ(field.shape(1), nvar ); + EXPECT_EQ(field.shape(2), nlev ); + EXPECT_EQ(field.shape(3), nproma ); + + // When above is failing, we need to really stop here. + REQUIRE( not current_test().failed() ); + + auto value = array::make_view(field); + auto g = array::make_view(fs.global_index()); + + ATLAS_DEBUG("Set"); + for (idx_t jblk = 0; jblk < nblk; ++jblk) { + auto blk = fs.block(jblk); + for (idx_t jvar = 0; jvar < nvar; ++jvar) { + for (idx_t jlev = 0; jlev < nlev; ++jlev) { + for (idx_t jrof = 0; jrof < blk.size(); ++jrof) { + idx_t h_idx = blk.index(jrof); + value(jblk,jvar,jlev,jrof) = g(h_idx); + } + } + } + } + + ATLAS_DEBUG("Check"); + std::vector counters; + counters.resize(nvar * nlev); + std::fill(counters.begin(), counters.end(), 0); + for (idx_t jblk = 0; jblk < nblk; ++jblk) { + auto blk = fs.block(jblk); + for (idx_t jvar = 0; jvar < nvar; ++jvar) { + for (idx_t jlev = 0; jlev < nlev; ++jlev) { + for (idx_t jrof = 0; jrof < blk.size(); ++jrof) { + idx_t& icheck = counters[jvar*nlev + jlev]; + EXPECT_EQ(value(jblk,jvar,jlev,jrof), g(icheck)); + icheck++; + } + } + } + } + + auto xy = array::make_view(fs.xy()); + auto r = array::make_view(fs.remote_index()); + auto p = array::make_view(fs.partition()); + } + + SECTION("test_BlockStructuredColumns scatter/gather") { + auto fs = functionspace::StructuredColumns(grid, config); + run_scatter_gather(grid, fs, nlev, nvar); + run_scatter_gather(grid, fs, nlev, nvar); + run_scatter_gather(grid, fs, nlev, nvar); + run_scatter_gather(grid, fs, nlev, nvar); + } +} + +//----------------------------------------------------------------------------- + + +} // namespace test +} // namespace atlas + +int main(int argc, char** argv) { + return atlas::test::run(argc, argv); +} From bf5a663441c1f20ad7f20a9119a850db5358c3e2 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Tue, 15 Nov 2022 13:44:06 +0100 Subject: [PATCH 04/25] ATLAS-363 Fix compilation and tidy --- .../detail/BlockStructuredColumns.cc | 148 ++++++++---------- .../detail/BlockStructuredColumns.h | 2 - .../functionspace/detail/StructuredColumns.cc | 7 +- .../functionspace/detail/StructuredColumns.h | 9 +- 4 files changed, 69 insertions(+), 97 deletions(-) diff --git a/src/atlas/functionspace/detail/BlockStructuredColumns.cc b/src/atlas/functionspace/detail/BlockStructuredColumns.cc index 4f1b84f5e..a291e0294 100644 --- a/src/atlas/functionspace/detail/BlockStructuredColumns.cc +++ b/src/atlas/functionspace/detail/BlockStructuredColumns.cc @@ -38,8 +38,11 @@ #include "atlas/util/CoordinateEnums.h" #include "atlas/util/detail/Cache.h" +namespace atlas { +namespace functionspace { +namespace detail { + namespace { -using namespace atlas; template void block_copy(const Field sloc, Field loc, const functionspace::detail::BlockStructuredColumns& fs) { @@ -153,13 +156,48 @@ void rev_block_copy(const Field loc, Field sloc, const functionspace::detail::Bl } } -}// namespace -namespace atlas { -namespace functionspace { -namespace detail { +void transpose_nonblocked_to_blocked(const Field& nonblocked, Field& blocked, const functionspace::detail::BlockStructuredColumns& fs) { + auto kind = nonblocked.datatype().kind(); + if (kind == array::DataType::kind()) { + block_copy(nonblocked, blocked, fs); + } + else if (kind == array::DataType::kind()) { + block_copy(nonblocked, blocked, fs); + } + else if (kind == array::DataType::kind()) { + block_copy(nonblocked, blocked, fs); + } + else if (kind == array::DataType::kind()) { + block_copy(nonblocked, blocked, fs); + } + else { + throw_Exception("datatype not supported", Here()); + } +} + +void transpose_blocked_to_nonblocked(const Field& blocked, Field& nonblocked, const functionspace::detail::BlockStructuredColumns& fs) { + auto kind = blocked.datatype().kind(); + if (kind == array::DataType::kind()) { + rev_block_copy(blocked, nonblocked, fs); + } + else if (kind == array::DataType::kind()) { + rev_block_copy(blocked, nonblocked, fs); + } + else if (kind == array::DataType::kind()) { + rev_block_copy(blocked, nonblocked, fs); + } + else if (kind == array::DataType::kind()) { + rev_block_copy(blocked, nonblocked, fs); + } + else { + throw_Exception("datatype not supported", Here()); + } +} +}// namespace + array::ArrayShape BlockStructuredColumns::config_shape(const eckit::Configuration& config) const { array::ArrayShape shape; @@ -230,8 +268,8 @@ Field BlockStructuredColumns::createField(const eckit::Configuration& options) c } Field BlockStructuredColumns::createField(const Field& other, const eckit::Configuration& config) const { - return createField(option::datatype(other.datatype()) | option::levels(other.levels()) | - option::variables(other.variables()) | + return createField(option::name(other.name()) | option::datatype(other.datatype()) | + option::levels(other.levels()) |option::variables(other.variables()) | option::type(other.metadata().getString("type", "scalar")) | config); } // ---------------------------------------------------------------------------- @@ -242,47 +280,14 @@ Field BlockStructuredColumns::createField(const Field& other, const eckit::Confi void BlockStructuredColumns::scatter(const FieldSet& global_fieldset, FieldSet& local_fieldset) const { ATLAS_ASSERT(local_fieldset.size() == global_fieldset.size()); for (idx_t f = 0; f < local_fieldset.size(); ++f) { - const Field& glb = global_fieldset[f]; - auto config = option::datatype(glb.datatype()) | option::levels(glb.levels()) | option::variables(glb.variables()); - auto sloc = structuredcolumns_->createField(config); - const idx_t nb_fields = 1; - idx_t root(0); - glb.metadata().get("owner", root); - - Field& loc = local_fieldset[f]; - glb.metadata().broadcast(loc.metadata(), root); - loc.metadata().set("global", false); - - if (sloc.datatype().kind() == array::DataType::kind()) { - parallel::Field glb_field(structuredcolumns_->make_leveled_view(glb)); - parallel::Field sloc_field(structuredcolumns_->make_leveled_view(sloc)); - structuredcolumns_->scatter().scatter(&glb_field, &sloc_field, nb_fields, root); - block_copy(sloc, loc, *this); - } - else if (sloc.datatype().kind() == array::DataType::kind()) { - parallel::Field glb_field(structuredcolumns_->make_leveled_view(glb)); - parallel::Field sloc_field(structuredcolumns_->make_leveled_view(sloc)); - structuredcolumns_->scatter().scatter(&glb_field, &sloc_field, nb_fields, root); - block_copy(sloc, loc, *this); - } - else if (sloc.datatype().kind() == array::DataType::kind()) { - parallel::Field glb_field(structuredcolumns_->make_leveled_view(glb)); - parallel::Field sloc_field(structuredcolumns_->make_leveled_view(sloc)); - structuredcolumns_->scatter().scatter(&glb_field, &sloc_field, nb_fields, root); - block_copy(sloc, loc, *this); - } - else if (sloc.datatype().kind() == array::DataType::kind()) { - parallel::Field glb_field(structuredcolumns_->make_leveled_view(glb)); - parallel::Field sloc_field(structuredcolumns_->make_leveled_view(sloc)); - structuredcolumns_->scatter().scatter(&glb_field, &sloc_field, nb_fields, root); - block_copy(sloc, loc, *this); - } - else { - throw_Exception("datatype not supported", Here()); - } + const Field& glb = global_fieldset[f]; + Field& loc = local_fieldset[f]; + auto sloc = structuredcolumns_->createField(glb, util::Config("global",false)); + structuredcolumns_->scatter(glb, sloc); + loc.metadata() = sloc.metadata(); + transpose_nonblocked_to_blocked(sloc, loc, *this); } } -// ---------------------------------------------------------------------------- // ---------------------------------------------------------------------------- // Scatter Field @@ -294,7 +299,6 @@ void BlockStructuredColumns::scatter(const Field& global, Field& local) const { local_fields.add(local); scatter(global_fields, local_fields); } -// ---------------------------------------------------------------------------- // ---------------------------------------------------------------------------- // Gather FieldSet @@ -302,45 +306,13 @@ void BlockStructuredColumns::scatter(const Field& global, Field& local) const { void BlockStructuredColumns::gather(const FieldSet& local_fieldset, FieldSet& global_fieldset) const { ATLAS_ASSERT(local_fieldset.size() == global_fieldset.size()); for (idx_t f = 0; f < local_fieldset.size(); ++f) { - const Field& loc = local_fieldset[f]; - auto config = option::datatype(loc.datatype()) | option::levels(loc.levels()) | option::variables(loc.variables()); - auto sloc = structuredcolumns_->createField(config); - - - idx_t root(0); - sloc.metadata().set("global", false); - Field& glb = global_fieldset[f]; - glb.metadata().broadcast(sloc.metadata(), root); - const idx_t nb_fields = 1; - glb.metadata().get("owner", root); - - if (sloc.datatype().kind() == array::DataType::kind()) { - rev_block_copy(loc, sloc, *this); - parallel::Field sloc_field(structuredcolumns_->make_leveled_view(sloc)); - parallel::Field glb_field(structuredcolumns_->make_leveled_view(glb)); - structuredcolumns_->gather().gather(&sloc_field, &glb_field, nb_fields, root); - } - else if (sloc.datatype().kind() == array::DataType::kind()) { - rev_block_copy(loc, sloc, *this); - parallel::Field sloc_field(structuredcolumns_->make_leveled_view(sloc)); - parallel::Field glb_field(structuredcolumns_->make_leveled_view(glb)); - structuredcolumns_->gather().gather(&sloc_field, &glb_field, nb_fields, root); - } - else if (sloc.datatype().kind() == array::DataType::kind()) { - rev_block_copy(loc, sloc, *this); - parallel::Field sloc_field(structuredcolumns_->make_leveled_view(sloc)); - parallel::Field glb_field(structuredcolumns_->make_leveled_view(glb)); - structuredcolumns_->gather().gather(&sloc_field, &glb_field, nb_fields, root); - } - else if (sloc.datatype().kind() == array::DataType::kind()) { - rev_block_copy(loc, sloc, *this); - parallel::Field sloc_field(structuredcolumns_->make_leveled_view(sloc)); - parallel::Field glb_field(structuredcolumns_->make_leveled_view(glb)); - structuredcolumns_->gather().gather(&sloc_field, &glb_field, nb_fields, root); - } - else { - throw_Exception("datatype not supported", Here()); - } + const Field& loc = local_fieldset[f]; + Field& glb = global_fieldset[f]; + auto sloc = structuredcolumns_->createField(loc, util::Config("global",false)); + transpose_blocked_to_nonblocked(loc, sloc, *this); + structuredcolumns_->gather(sloc, glb); + glb.metadata() = loc.metadata(); + glb.metadata().set("global", false); } } // ---------------------------------------------------------------------------- @@ -372,6 +344,10 @@ void BlockStructuredColumns::gather(const Field& local, Field& global) const { gather(local_fields, global_fields); } +// ---------------------------------------------------------------------------- +// Checksum Field +// ---------------------------------------------------------------------------- + std::string BlockStructuredColumns::checksum(const Field&) const { ATLAS_NOTIMPLEMENTED; } diff --git a/src/atlas/functionspace/detail/BlockStructuredColumns.h b/src/atlas/functionspace/detail/BlockStructuredColumns.h index 48fb9aee9..c34a0aaa7 100644 --- a/src/atlas/functionspace/detail/BlockStructuredColumns.h +++ b/src/atlas/functionspace/detail/BlockStructuredColumns.h @@ -124,8 +124,6 @@ class BlockStructuredColumns : public FunctionSpaceImpl { array::ArrayAlignment config_alignment(const eckit::Configuration&) const; array::ArraySpec config_spec(const eckit::Configuration&) const; - void create_remote_index() const { structuredcolumns_->create_remote_index(); } - private: // data idx_t nproma_; idx_t endblk_size_; diff --git a/src/atlas/functionspace/detail/StructuredColumns.cc b/src/atlas/functionspace/detail/StructuredColumns.cc index 645c16b7b..140e47c0f 100644 --- a/src/atlas/functionspace/detail/StructuredColumns.cc +++ b/src/atlas/functionspace/detail/StructuredColumns.cc @@ -46,9 +46,10 @@ namespace atlas { namespace functionspace { namespace detail { +namespace { template -array::LocalView StructuredColumns::make_leveled_view(Field& field) const { +array::LocalView make_leveled_view(Field& field) { using namespace array; if (field.levels()) { if (field.variables()) { @@ -69,7 +70,7 @@ array::LocalView StructuredColumns::make_leveled_view(Field& field) const } template -std::string StructuredColumns::checksum_3d_field(const parallel::Checksum& checksum, const Field& field) const { +std::string checksum_3d_field(const parallel::Checksum& checksum, const Field& field) { auto values = make_leveled_view(field); array::ArrayT surface_field(values.shape(0), values.shape(2)); auto surface = array::make_view(surface_field); @@ -85,6 +86,8 @@ std::string StructuredColumns::checksum_3d_field(const parallel::Checksum& check return checksum.execute(surface.data(), surface_field.stride(0)); } +} // namespace + class StructuredColumnsHaloExchangeCache : public util::Cache, public grid::detail::grid::GridObserver { private: diff --git a/src/atlas/functionspace/detail/StructuredColumns.h b/src/atlas/functionspace/detail/StructuredColumns.h index 210357786..f3d3d6850 100644 --- a/src/atlas/functionspace/detail/StructuredColumns.h +++ b/src/atlas/functionspace/detail/StructuredColumns.h @@ -298,12 +298,6 @@ class StructuredColumns : public FunctionSpaceImpl { } }; - template - array::LocalView make_leveled_view(Field& field) const; - - template - std::string checksum_3d_field(const parallel::Checksum& checksum, const Field& field) const; - idx_t j_begin_; idx_t j_end_; std::vector i_begin_; @@ -323,8 +317,9 @@ class StructuredColumns : public FunctionSpaceImpl { friend struct BlockStructuredColumnsFortranAccess; Map2to1 ij2gp_; - friend class BlockStructuredColumns; void setup(const grid::Distribution& distribution, const eckit::Configuration& config); + + friend class BlockStructuredColumns; }; // ------------------------------------------------------------------- From bedc88d33778d755663bb85d9e86a8ba9d89f1cd Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Fri, 18 Nov 2022 09:26:21 +0100 Subject: [PATCH 05/25] Fix Mesh::operator bool() to reflect empty mesh --- src/atlas/mesh/Mesh.cc | 3 +++ src/atlas/mesh/Mesh.h | 2 ++ 2 files changed, 5 insertions(+) diff --git a/src/atlas/mesh/Mesh.cc b/src/atlas/mesh/Mesh.cc index 0e02781de..e7e782e32 100644 --- a/src/atlas/mesh/Mesh.cc +++ b/src/atlas/mesh/Mesh.cc @@ -9,6 +9,7 @@ */ #include "atlas/mesh/Mesh.h" +#include "atlas/mesh/Nodes.h" #include "atlas/grid/Grid.h" #include "atlas/grid/Partitioner.h" #include "atlas/meshgenerator/MeshGenerator.h" @@ -42,6 +43,8 @@ Mesh::Mesh(const Grid& grid, const grid::Partitioner& partitioner): Mesh::Mesh(eckit::Stream& stream): Handle(new Implementation(stream)) {} +Mesh::operator bool() const { return get()->nodes().size() > 0; } + //---------------------------------------------------------------------------------------------------------------------- } // namespace atlas diff --git a/src/atlas/mesh/Mesh.h b/src/atlas/mesh/Mesh.h index 27ffc59af..3c4e09af7 100644 --- a/src/atlas/mesh/Mesh.h +++ b/src/atlas/mesh/Mesh.h @@ -72,6 +72,8 @@ class Mesh : DOXYGEN_HIDE(public util::ObjectHandle) { using Handle::Handle; Mesh(); + operator bool() const; + /// @brief Generate a mesh from a Grid with recommended mesh generator and partitioner strategy Mesh(const Grid&); From 10866b6bfedf22e949088ecf4994e8ef83f04d7d Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Fri, 18 Nov 2022 09:27:17 +0100 Subject: [PATCH 06/25] Add missing CellColumns::halo() function --- src/atlas/functionspace/CellColumns.cc | 4 ++++ src/atlas/functionspace/CellColumns.h | 4 ++++ 2 files changed, 8 insertions(+) diff --git a/src/atlas/functionspace/CellColumns.cc b/src/atlas/functionspace/CellColumns.cc index 1c2858886..51c542281 100644 --- a/src/atlas/functionspace/CellColumns.cc +++ b/src/atlas/functionspace/CellColumns.cc @@ -839,5 +839,9 @@ const parallel::Checksum& CellColumns::checksum() const { return functionspace_->checksum(); } +const mesh::Halo& CellColumns::halo() const { + return functionspace_->halo(); +} + } // namespace functionspace } // namespace atlas diff --git a/src/atlas/functionspace/CellColumns.h b/src/atlas/functionspace/CellColumns.h index ff620d9c9..407d47439 100644 --- a/src/atlas/functionspace/CellColumns.h +++ b/src/atlas/functionspace/CellColumns.h @@ -69,6 +69,8 @@ class CellColumns : public functionspace::FunctionSpaceImpl { // -- Parallelisation aware methods + const mesh::Halo& halo() const { return halo_; } + void haloExchange(const FieldSet&, bool on_device = false) const override; void haloExchange(const Field&, bool on_device = false) const override; const parallel::HaloExchange& halo_exchange() const; @@ -181,6 +183,8 @@ class CellColumns : public FunctionSpace { const mesh::HybridElements& cells() const; // -- Parallelisation aware methods + const mesh::Halo& halo() const; + const parallel::HaloExchange& halo_exchange() const; std::string checksum(const FieldSet&) const; From 5f871d1c67306de984a0bb6df0f98c4c9885cc94 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Fri, 18 Nov 2022 09:30:46 +0100 Subject: [PATCH 07/25] Fix interpolation method UnstructuredBilinearLonLat for failing triangular elements --- .../UnstructuredBilinearLonLat.cc | 26 ++++++++++++++++--- 1 file changed, 23 insertions(+), 3 deletions(-) diff --git a/src/atlas/interpolation/method/unstructured/UnstructuredBilinearLonLat.cc b/src/atlas/interpolation/method/unstructured/UnstructuredBilinearLonLat.cc index 9814560f2..bd8c0c579 100644 --- a/src/atlas/interpolation/method/unstructured/UnstructuredBilinearLonLat.cc +++ b/src/atlas/interpolation/method/unstructured/UnstructuredBilinearLonLat.cc @@ -326,7 +326,7 @@ Method::Triplets UnstructuredBilinearLonLat::projectPointToElements(size_t ip, c } PointLonLat o_loc{o_lon, (*olonlat_)(ip, LAT)}; // lookup point - auto inv_dist_weight = [](element::Quad2D& q, const PointXY& loc, std::array& w) { + auto inv_dist_weight_quad = [](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); @@ -342,6 +342,22 @@ Method::Triplets UnstructuredBilinearLonLat::projectPointToElements(size_t ip, c w[i] *= suminv; } }; + auto inv_dist_weight_triag = [](element::Triag2D& q, const PointXY& loc, std::array& w) { + double d[3]; + 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); + w[0] = d[1] * d[2]; + w[1] = d[0] * d[2]; + w[2] = d[1] * d[0]; + w[3] = 0.; + + double suminv = 1. / (w[0] + w[1] + w[2]); + for (size_t i = 0; i < 3; ++i) { + w[i] *= suminv; + } + }; + for (ElemIndex3::NodeList::const_iterator itc = elems.begin(); itc != elems.end(); ++itc) { const idx_t elem_id = idx_t((*itc).value().payload()); @@ -361,6 +377,10 @@ Method::Triplets UnstructuredBilinearLonLat::projectPointToElements(size_t ip, c PointLonLat{(*ilonlat_)(idx[1], LON), (*ilonlat_)(idx[1], LAT)}, PointLonLat{(*ilonlat_)(idx[2], LON), (*ilonlat_)(idx[2], LAT)}); + if (itc == elems.begin()) { + inv_dist_weight_triag(triag, 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(triag.area()); @@ -410,7 +430,7 @@ Method::Triplets UnstructuredBilinearLonLat::projectPointToElements(size_t ip, c 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); + inv_dist_weight_quad(quad, o_loc, inv_dist_w); } // pick an epsilon based on a characteristic length (sqrt(area)) @@ -444,7 +464,7 @@ Method::Triplets UnstructuredBilinearLonLat::projectPointToElements(size_t ip, c // 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) { + for (size_t i = 0; i < connectivity_->cols(elem_id); ++i) { idx[i] = (*connectivity_)(elem_id, i); triplets.emplace_back(ip, idx[i], inv_dist_w[i]); } From d7dcaa7aa5446a71985308f6cf18f0a10de22edb Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Fri, 18 Nov 2022 09:32:05 +0100 Subject: [PATCH 08/25] Rename StructuredMeshGenerator option 3d -> three_dimensional, but still backward compatible --- .../meshgenerator/detail/StructuredMeshGenerator.cc | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/atlas/meshgenerator/detail/StructuredMeshGenerator.cc b/src/atlas/meshgenerator/detail/StructuredMeshGenerator.cc index 4a4b773fe..87420300f 100644 --- a/src/atlas/meshgenerator/detail/StructuredMeshGenerator.cc +++ b/src/atlas/meshgenerator/detail/StructuredMeshGenerator.cc @@ -98,7 +98,7 @@ StructuredMeshGenerator::StructuredMeshGenerator(const eckit::Parametrisation& p bool three_dimensional; if (p.get("three_dimensional", three_dimensional) || p.get("3d", three_dimensional)) { - options.set("3d", three_dimensional); + options.set("three_dimensional", three_dimensional); } size_t nb_parts; @@ -151,14 +151,14 @@ void StructuredMeshGenerator::configure_defaults() { // connects elements // to the first node only. Note this option will only be looked at in case // other option - // "3d"==true + // "three_dimensional"==true options.set("unique_pole", true); // This option creates elements that connect east to west at greenwich // meridian // when true, instead of creating periodic ghost-points at east boundary when // false - options.set("3d", false); + options.set("three_dimensional", false); // This option sets number of parts the mesh will be split in options.set("nb_parts", mpi::size()); @@ -266,7 +266,7 @@ void StructuredMeshGenerator::generate_region(const StructuredGrid& rg, const gr double max_angle = options.getDouble("angle"); bool triangulate_quads = options.getBool("triangulate"); - bool three_dimensional = options.getBool("3d"); + bool three_dimensional = options.getBool("three_dimensional"); bool has_north_pole = eckit::types::is_approximately_equal(rg.y().front(), 90.); bool has_south_pole = eckit::types::is_approximately_equal(rg.y().back(), -90.); bool unique_pole = options.getBool("unique_pole") && three_dimensional && has_north_pole && has_south_pole; @@ -874,7 +874,7 @@ void StructuredMeshGenerator::generate_mesh(const StructuredGrid& rg, const grid bool mypart_at_north = std::any_of(part_north.begin(), part_north.end(), [&](int p) { return p == mypart; }); bool mypart_at_south = std::any_of(part_south.begin(), part_south.end(), [&](int p) { return p == mypart; }); - bool three_dimensional = options.getBool("3d"); + bool three_dimensional = options.getBool("three_dimensional"); bool periodic_east_west = rg.periodic(); bool include_periodic_ghost_points = periodic_east_west && !three_dimensional; bool remove_periodic_ghost_points = periodic_east_west && three_dimensional; From 5cd6a445fcfddf0f3f5827c1f70c8656f2526083 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Fri, 18 Nov 2022 09:33:53 +0100 Subject: [PATCH 09/25] Support more functionspaces for interpolation with FiniteElement, UnstructuredBilinearLonLat, GridBoxMethod, (K)NearestNeighbour --- .../interpolation/method/knn/GridBoxMethod.cc | 27 ++++++++-- .../method/knn/KNearestNeighbours.cc | 20 +++----- .../method/knn/KNearestNeighboursBase.cc | 51 +++++++++++++++++-- .../method/knn/NearestNeighbour.cc | 16 ++---- .../method/unstructured/FiniteElement.cc | 31 ++++++----- .../UnstructuredBilinearLonLat.cc | 29 +++++------ 6 files changed, 110 insertions(+), 64 deletions(-) diff --git a/src/atlas/interpolation/method/knn/GridBoxMethod.cc b/src/atlas/interpolation/method/knn/GridBoxMethod.cc index 24357f60b..729146e20 100644 --- a/src/atlas/interpolation/method/knn/GridBoxMethod.cc +++ b/src/atlas/interpolation/method/knn/GridBoxMethod.cc @@ -45,14 +45,31 @@ void GridBoxMethod::print(std::ostream& out) const { out << "GridBoxMethod[]"; } +namespace { +StructuredGrid extract_structured_grid(const FunctionSpace& fs) { + if (functionspace::StructuredColumns{fs}) { + return functionspace::StructuredColumns(fs).grid(); + } + if (functionspace::NodeColumns{fs}) { + return functionspace::NodeColumns(fs).mesh().grid(); + } + return Grid(); +} +} void GridBoxMethod::do_setup(const FunctionSpace& source, const FunctionSpace& target) { - if (functionspace::StructuredColumns(source) && functionspace::StructuredColumns(target)) { - do_setup(functionspace::StructuredColumns(source).grid(), functionspace::StructuredColumns(target).grid(), - Cache()); - return; + if( mpi::size() > 1 ) { + ATLAS_THROW_EXCEPTION("Cannot use GridBoxMethod in parallel yet."); + } + auto src_grid = extract_structured_grid(source); + auto tgt_grid = extract_structured_grid(target); + if( not src_grid ) { + ATLAS_THROW_EXCEPTION("Could not extract StructuredGrid from source function space " << source.type() ); + } + if( not tgt_grid ) { + ATLAS_THROW_EXCEPTION("Could not extract StructuredGrid from target function space " << target.type() ); } - ATLAS_NOTIMPLEMENTED; + do_setup(src_grid,tgt_grid,Cache()); } bool GridBoxMethod::intersect(size_t i, const GridBox& box, const util::IndexKDTree::ValueList& closest, diff --git a/src/atlas/interpolation/method/knn/KNearestNeighbours.cc b/src/atlas/interpolation/method/knn/KNearestNeighbours.cc index d39ffbf88..07df3faea 100644 --- a/src/atlas/interpolation/method/knn/KNearestNeighbours.cc +++ b/src/atlas/interpolation/method/knn/KNearestNeighbours.cc @@ -65,22 +65,14 @@ void KNearestNeighbours::do_setup(const Grid& source, const Grid& target, const void KNearestNeighbours::do_setup(const FunctionSpace& source, const FunctionSpace& target) { source_ = source; target_ = target; - functionspace::NodeColumns src = source; - functionspace::NodeColumns tgt = target; - ATLAS_ASSERT(src); - ATLAS_ASSERT(tgt); - - Mesh meshSource = src.mesh(); - Mesh meshTarget = tgt.mesh(); // build point-search tree - buildPointSearchTree(meshSource, src.halo()); + buildPointSearchTree(source); - array::ArrayView lonlat = array::make_view(meshTarget.nodes().lonlat()); + array::ArrayView lonlat = array::make_view(target.lonlat()); - size_t inp_npts = meshSource.nodes().size(); - meshSource.metadata().get("nb_nodes_including_halo[" + std::to_string(src.halo().size()) + "]", inp_npts); - size_t out_npts = meshTarget.nodes().size(); + size_t inp_npts = source.size(); + size_t out_npts = target.size(); // fill the sparse matrix std::vector weights_triplets; @@ -93,10 +85,12 @@ void KNearestNeighbours::do_setup(const FunctionSpace& source, const FunctionSpa Log::debug() << "Computing interpolation weights for " << out_npts << " points." << std::endl; for (size_t ip = 0; ip < out_npts; ++ip) { if (ip && (ip % 1000 == 0)) { + timer.stop(); auto elapsed = timer.elapsed(); + timer.start(); auto rate = eckit::types::is_approximately_equal(elapsed, 0.) ? std::numeric_limits::infinity() : (ip / elapsed); - Log::debug() << eckit::BigNum(ip) << " (at " << rate << " points/s)..." << std::endl; + Log::debug() << eckit::BigNum(ip) << " (at " << size_t(rate) << " points/s)... after " << elapsed << " s" << std::endl; } // find the closest input points to the output point diff --git a/src/atlas/interpolation/method/knn/KNearestNeighboursBase.cc b/src/atlas/interpolation/method/knn/KNearestNeighboursBase.cc index 6b1fd1908..f002589d6 100644 --- a/src/atlas/interpolation/method/knn/KNearestNeighboursBase.cc +++ b/src/atlas/interpolation/method/knn/KNearestNeighboursBase.cc @@ -13,6 +13,9 @@ #include "atlas/array.h" #include "atlas/functionspace/PointCloud.h" +#include "atlas/functionspace/NodeColumns.h" +#include "atlas/functionspace/CellColumns.h" +#include "atlas/functionspace/StructuredColumns.h" #include "atlas/interpolation/method/knn/KNearestNeighboursBase.h" #include "atlas/library/Library.h" #include "atlas/mesh/Nodes.h" @@ -46,8 +49,8 @@ void KNearestNeighboursBase::buildPointSearchTree(Mesh& meshSource, const mesh:: pTree_.build(); - // generate 3D point coordinates - mesh::actions::BuildXYZField("xyz")(meshSource); +// // generate 3D point coordinates +// mesh::actions::BuildXYZField("xyz")(meshSource); } namespace { @@ -60,20 +63,62 @@ void insert_tree(util::IndexKDTree& tree, const FunctionSpace_type& functionspac } } +void insert_tree(util::IndexKDTree& tree, const functionspace::NodeColumns& functionspace) { + auto lonlat = array::make_view(functionspace.lonlat()); + auto halo = array::make_view(functionspace.nodes().halo()); + int h = functionspace.halo().size(); + + for (idx_t ip = 0; ip < lonlat.shape(0); ++ip) { + if (halo(ip) <= h) { + tree.insert(PointLonLat(lonlat(ip, LON), lonlat(ip, LAT)), ip); + } + } +} + +void insert_tree(util::IndexKDTree& tree, const functionspace::CellColumns& functionspace) { + auto lonlat = array::make_view(functionspace.lonlat()); + auto halo = array::make_view(functionspace.cells().halo()); + int h = functionspace.halo().size(); + + for (idx_t ip = 0; ip < lonlat.shape(0); ++ip) { + if (halo(ip) <= h) { + tree.insert(PointLonLat(lonlat(ip, LON), lonlat(ip, LAT)), ip); + } + } +} + +void insert_tree(util::IndexKDTree& tree, const functionspace::StructuredColumns& functionspace) { + auto lonlat = array::make_view(functionspace.lonlat()); + + for (idx_t ip = 0; ip < lonlat.shape(0); ++ip) { + tree.insert(PointLonLat(lonlat(ip, LON), lonlat(ip, LAT)), ip); + } +} + } // namespace void KNearestNeighboursBase::buildPointSearchTree(const FunctionSpace& functionspace) { ATLAS_TRACE(); - eckit::TraceTimer tim("KNearestNeighboursBase::buildPointSearchTree()"); + eckit::TraceTimer timer("KNearestNeighboursBase::buildPointSearchTree()"); static bool fastBuildKDTrees = eckit::Resource("$ATLAS_FAST_BUILD_KDTREES", true); if (fastBuildKDTrees) { pTree_.reserve(functionspace.size()); } + if (functionspace::PointCloud fs = functionspace) { insert_tree(pTree_, fs); } + else if (functionspace::NodeColumns fs = functionspace) { + insert_tree(pTree_, fs); + } + else if (functionspace::CellColumns fs = functionspace) { + insert_tree(pTree_, fs); + } + else if (functionspace::StructuredColumns fs = functionspace) { + insert_tree(pTree_, fs); + } else { ATLAS_NOTIMPLEMENTED; } diff --git a/src/atlas/interpolation/method/knn/NearestNeighbour.cc b/src/atlas/interpolation/method/knn/NearestNeighbour.cc index fb500ba83..4e2811582 100644 --- a/src/atlas/interpolation/method/knn/NearestNeighbour.cc +++ b/src/atlas/interpolation/method/knn/NearestNeighbour.cc @@ -59,22 +59,14 @@ void NearestNeighbour::do_setup(const Grid& source, const Grid& target, const Ca void NearestNeighbour::do_setup(const FunctionSpace& source, const FunctionSpace& target) { source_ = source; target_ = target; - functionspace::NodeColumns src = source; - functionspace::NodeColumns tgt = target; - ATLAS_ASSERT(src); - ATLAS_ASSERT(tgt); - - Mesh meshSource = src.mesh(); - Mesh meshTarget = tgt.mesh(); // build point-search tree - buildPointSearchTree(meshSource, src.halo()); + buildPointSearchTree(source); - array::ArrayView lonlat = array::make_view(meshTarget.nodes().lonlat()); + array::ArrayView lonlat = array::make_view(target.lonlat()); - size_t inp_npts = meshSource.nodes().size(); - meshSource.metadata().get("nb_nodes_including_halo[" + std::to_string(src.halo().size()) + "]", inp_npts); - size_t out_npts = meshTarget.nodes().size(); + size_t inp_npts = source.size(); + size_t out_npts = target.size(); // fill the sparse matrix std::vector weights_triplets; diff --git a/src/atlas/interpolation/method/unstructured/FiniteElement.cc b/src/atlas/interpolation/method/unstructured/FiniteElement.cc index c949e3858..8584ae1a1 100644 --- a/src/atlas/interpolation/method/unstructured/FiniteElement.cc +++ b/src/atlas/interpolation/method/unstructured/FiniteElement.cc @@ -89,32 +89,31 @@ void FiniteElement::do_setup(const FunctionSpace& source, const FunctionSpace& t target_ = target; ATLAS_TRACE_SCOPE("Setup target") { - if (functionspace::NodeColumns tgt = target) { - Mesh meshTarget = tgt.mesh(); - // generate 3D point coordinates - 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_); + + auto create_xyz = [](Field lonlat_field) { + auto xyz_field = Field("xyz", array::make_datatype(), array::make_shape(lonlat_field.shape(0), 3)); + auto lonlat = array::make_view(lonlat_field); + auto xyz = array::make_view(xyz_field); PointXYZ p2; - for (idx_t n = 0; n < N; ++n) { + for (idx_t n = 0; n < lonlat.shape(0); ++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(); } + return xyz_field; + }; + + target_ghost_ = target.ghost(); + target_lonlat_ = target.lonlat(); + if (functionspace::NodeColumns tgt = target) { + auto meshTarget = tgt.mesh(); + target_xyz_ = mesh::actions::BuildXYZField("xyz")(meshTarget); } else { - ATLAS_NOTIMPLEMENTED; + target_xyz_ = create_xyz(target_lonlat_); } } diff --git a/src/atlas/interpolation/method/unstructured/UnstructuredBilinearLonLat.cc b/src/atlas/interpolation/method/unstructured/UnstructuredBilinearLonLat.cc index bd8c0c579..62276db31 100644 --- a/src/atlas/interpolation/method/unstructured/UnstructuredBilinearLonLat.cc +++ b/src/atlas/interpolation/method/unstructured/UnstructuredBilinearLonLat.cc @@ -86,31 +86,30 @@ void UnstructuredBilinearLonLat::do_setup(const FunctionSpace& source, const Fun target_ = target; ATLAS_TRACE_SCOPE("Setup target") { - 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_); + auto create_xyz = [](Field lonlat_field) { + auto xyz_field = Field("xyz", array::make_datatype(), array::make_shape(lonlat_field.shape(0), 3)); + auto lonlat = array::make_view(lonlat_field); + auto xyz = array::make_view(xyz_field); PointXYZ p2; - for (idx_t n = 0; n < N; ++n) { + for (idx_t n = 0; n < lonlat.shape(0); ++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(); } + return xyz_field; + }; + + target_ghost_ = target.ghost(); + target_lonlat_ = target.lonlat(); + if (functionspace::NodeColumns tgt = target) { + auto meshTarget = tgt.mesh(); + target_xyz_ = mesh::actions::BuildXYZField("xyz")(meshTarget); } else { - ATLAS_NOTIMPLEMENTED; + target_xyz_ = create_xyz(target_lonlat_); } } From 3b3fc3d6ad35ffe2248212ec73ecd37efed18711 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Sat, 19 Nov 2022 22:29:28 +0100 Subject: [PATCH 10/25] Fixup previous commit, fixing segmentation fault --- src/atlas/interpolation/method/knn/KNearestNeighbours.cc | 6 +++--- src/atlas/interpolation/method/knn/NearestNeighbour.cc | 4 +++- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/src/atlas/interpolation/method/knn/KNearestNeighbours.cc b/src/atlas/interpolation/method/knn/KNearestNeighbours.cc index 07df3faea..d17839489 100644 --- a/src/atlas/interpolation/method/knn/KNearestNeighbours.cc +++ b/src/atlas/interpolation/method/knn/KNearestNeighbours.cc @@ -78,16 +78,16 @@ void KNearestNeighbours::do_setup(const FunctionSpace& source, const FunctionSpa std::vector weights_triplets; weights_triplets.reserve(out_npts * k_); { - Trace timer(Here(), "atlas::interpolation::method::NearestNeighbour::do_setup()"); + Trace timer(Here(), "atlas::interpolation::method::KNearestNeighbour::do_setup()"); std::vector weights; Log::debug() << "Computing interpolation weights for " << out_npts << " points." << std::endl; for (size_t ip = 0; ip < out_npts; ++ip) { if (ip && (ip % 1000 == 0)) { - timer.stop(); + timer.pause(); auto elapsed = timer.elapsed(); - timer.start(); + timer.resume(); auto rate = eckit::types::is_approximately_equal(elapsed, 0.) ? std::numeric_limits::infinity() : (ip / elapsed); Log::debug() << eckit::BigNum(ip) << " (at " << size_t(rate) << " points/s)... after " << elapsed << " s" << std::endl; diff --git a/src/atlas/interpolation/method/knn/NearestNeighbour.cc b/src/atlas/interpolation/method/knn/NearestNeighbour.cc index 4e2811582..515f9249c 100644 --- a/src/atlas/interpolation/method/knn/NearestNeighbour.cc +++ b/src/atlas/interpolation/method/knn/NearestNeighbour.cc @@ -75,10 +75,12 @@ void NearestNeighbour::do_setup(const FunctionSpace& source, const FunctionSpace Trace timer(Here(), "atlas::interpolation::method::NearestNeighbour::do_setup()"); for (size_t ip = 0; ip < out_npts; ++ip) { if (ip && (ip % 1000 == 0)) { + timer.pause(); auto elapsed = timer.elapsed(); + timer.resume(); auto rate = eckit::types::is_approximately_equal(elapsed, 0.) ? std::numeric_limits::infinity() : (ip / elapsed); - Log::debug() << eckit::BigNum(ip) << " (at " << rate << " points/s)..." << std::endl; + Log::debug() << eckit::BigNum(ip) << " (at " << rate << " points/s)... after " << elapsed << " s" << std::endl; } // find the closest input point to the output point From 2a87caeb70ab98ab03b64844d89629951b1c8c2f Mon Sep 17 00:00:00 2001 From: Paul Cresswell Date: Tue, 22 Nov 2022 13:15:51 +0000 Subject: [PATCH 11/25] ATLAS-372: Add missing header file --- src/atlas/grid/detail/distribution/BandsDistribution.cc | 1 + 1 file changed, 1 insertion(+) diff --git a/src/atlas/grid/detail/distribution/BandsDistribution.cc b/src/atlas/grid/detail/distribution/BandsDistribution.cc index 237f382fa..1926fe8f1 100644 --- a/src/atlas/grid/detail/distribution/BandsDistribution.cc +++ b/src/atlas/grid/detail/distribution/BandsDistribution.cc @@ -12,6 +12,7 @@ #include #include +#include #include "atlas/grid/Grid.h" #include "atlas/runtime/Exception.h" From 1f8182ee140a0b1d9ce9c4af26002a6070f31696 Mon Sep 17 00:00:00 2001 From: Slavko Brdar Date: Thu, 15 Dec 2022 14:54:17 +0100 Subject: [PATCH 12/25] Add Fortran constructors for ShiftedLonLat, ShiftedLon, ShiftedLat grids commit 7e4adacb7938b733dc909adcac7dda55e9f6ad16 Author: Slavko Brdar Date: Thu Dec 15 14:35:08 2022 +0100 code simplification (tnx Willem) commit 4e64d0873b86a2ff5441b7896cfdcefadfd85ec3 Author: Slavko Brdar Date: Thu Dec 15 14:21:53 2022 +0100 fix run time error for shifted grids - correct usage of configuration parameters commit 9b2ebd1ffa1e523d9a6b30874e1760819906e8f9 Author: Slavko Brdar Date: Thu Dec 15 13:37:13 2022 +0100 add Shifted* grids with a unit test to Fortran interface commit 4f053e61f7a76b24886c6629f54b1e696187a986 Author: Slavko Brdar Date: Wed Dec 14 16:49:26 2022 +0100 1) replace NOTIMPLEMENTED in fortran api for RegularLonLat and Shifted* grids; 2) add unit test --- src/atlas/grid/detail/grid/Structured.cc | 27 ++-- src/atlas/grid/detail/grid/Structured.h | 12 +- src/atlas_f/atlas_module.F90 | 3 + src/atlas_f/grid/atlas_Grid_module.F90 | 168 +++++++++++++++++++++++ src/tests/grid/fctest_grids.F90 | 6 + 5 files changed, 199 insertions(+), 17 deletions(-) diff --git a/src/atlas/grid/detail/grid/Structured.cc b/src/atlas/grid/detail/grid/Structured.cc index e79d0b56b..02b26384e 100644 --- a/src/atlas/grid/detail/grid/Structured.cc +++ b/src/atlas/grid/detail/grid/Structured.cc @@ -852,7 +852,7 @@ int atlas__grid__Structured__reduced(Structured* This) { return This->reduced(); } -const Structured* atlas__grid__Structured(char* identifier) { +const Structured* atlas__grid__Structured(const char* identifier) { const Structured* grid = dynamic_cast(Grid::create(std::string(identifier))); ATLAS_ASSERT(grid != nullptr); return grid; @@ -870,20 +870,25 @@ void atlas__grid__Structured__delete(Structured* This) { delete This; } -Structured* atlas__grid__regular__RegularGaussian(long N) { - ATLAS_NOTIMPLEMENTED; +const Structured* atlas__grid__regular__RegularGaussian(long N) { + std::string gridname = "F"+std::to_string(N); + return atlas__grid__Structured(gridname.c_str()); } -Structured* atlas__grid__regular__RegularLonLat(long nlon, long nlat) { - ATLAS_NOTIMPLEMENTED; +const Structured* atlas__grid__regular__RegularLonLat(long nlon, long nlat) { + std::string gridname = "L"+std::to_string(nlon)+"x"+std::to_string(nlat); + return atlas__grid__Structured(gridname.c_str()); } -Structured* atlas__grid__regular__ShiftedLonLat(long nlon, long nlat) { - ATLAS_NOTIMPLEMENTED; +const Structured* atlas__grid__regular__ShiftedLonLat(long nlon, long nlat) { + std::string gridname = "S"+std::to_string(nlon)+"x"+std::to_string(nlat); + return atlas__grid__Structured(gridname.c_str()); } -Structured* atlas__grid__regular__ShiftedLon(long nlon, long nlat) { - ATLAS_NOTIMPLEMENTED; +const Structured* atlas__grid__regular__ShiftedLon(long nlon, long nlat) { + std::string gridname = "Slon"+std::to_string(nlon)+"x"+std::to_string(nlat); + return atlas__grid__Structured(gridname.c_str()); } -Structured* atlas__grid__regular__ShiftedLat(long nlon, long nlat) { - ATLAS_NOTIMPLEMENTED; +const Structured* atlas__grid__regular__ShiftedLat(long nlon, long nlat) { + std::string gridname = "Slat"+std::to_string(nlon)+"x"+std::to_string(nlat); + return atlas__grid__Structured(gridname.c_str()); } idx_t atlas__grid__Gaussian__N(Structured* This) { diff --git a/src/atlas/grid/detail/grid/Structured.h b/src/atlas/grid/detail/grid/Structured.h index 23a3d590c..2783c022a 100644 --- a/src/atlas/grid/detail/grid/Structured.h +++ b/src/atlas/grid/detail/grid/Structured.h @@ -426,19 +426,19 @@ class Structured : public Grid { extern "C" { void atlas__grid__Structured__delete(Structured* This); -const Structured* atlas__grid__Structured(char* identifier); +const Structured* atlas__grid__Structured(const char* identifier); const Structured* atlas__grid__Structured__config(util::Config* conf); -Structured* atlas__grid__regular__RegularGaussian(long N); +const Structured* atlas__grid__regular__RegularGaussian(long N); Structured* atlas__grid__reduced__ReducedGaussian_int(int nx[], long ny); Structured* atlas__grid__reduced__ReducedGaussian_long(long nx[], long ny); Structured* atlas__grid__reduced__ReducedGaussian_int_projection(int nx[], long ny, const Projection::Implementation* projection); Structured* atlas__grid__reduced__ReducedGaussian_long_projection(long nx[], long ny, const Projection::Implementation* projection); -Structured* atlas__grid__regular__RegularLonLat(long nx, long ny); -Structured* atlas__grid__regular__ShiftedLonLat(long nx, long ny); -Structured* atlas__grid__regular__ShiftedLon(long nx, long ny); -Structured* atlas__grid__regular__ShiftedLat(long nx, long ny); +const Structured* atlas__grid__regular__RegularLonLat(long nx, long ny); +const Structured* atlas__grid__regular__ShiftedLonLat(long nx, long ny); +const Structured* atlas__grid__regular__ShiftedLon(long nx, long ny); +const Structured* atlas__grid__regular__ShiftedLat(long nx, long ny); void atlas__grid__Structured__nx_array(Structured* This, const idx_t*& nx, idx_t& size); idx_t atlas__grid__Structured__nx(Structured* This, idx_t j); diff --git a/src/atlas_f/atlas_module.F90 b/src/atlas_f/atlas_module.F90 index 58d1b1540..8931791df 100644 --- a/src/atlas_f/atlas_module.F90 +++ b/src/atlas_f/atlas_module.F90 @@ -81,6 +81,9 @@ module atlas_module & atlas_ReducedGaussianGrid, & & atlas_RegularGaussianGrid, & & atlas_RegularLonLatGrid, & + & atlas_ShiftedLonLatGrid, & + & atlas_ShiftedLonGrid, & + & atlas_ShiftedLatGrid, & & atlas_RegionalGrid use atlas_Vertical_module, only :& & atlas_Vertical diff --git a/src/atlas_f/grid/atlas_Grid_module.F90 b/src/atlas_f/grid/atlas_Grid_module.F90 index aae3e384d..42dd9e458 100644 --- a/src/atlas_f/grid/atlas_Grid_module.F90 +++ b/src/atlas_f/grid/atlas_Grid_module.F90 @@ -31,6 +31,9 @@ module atlas_Grid_module public :: atlas_ReducedGaussianGrid public :: atlas_RegularGaussianGrid public :: atlas_RegularLonLatGrid +public :: atlas_ShiftedLonLatGrid +public :: atlas_ShiftedLonGrid +public :: atlas_ShiftedLatGrid public :: atlas_RegionalGrid private @@ -276,6 +279,87 @@ module atlas_Grid_module !------------------------------------------------------------------------------ +TYPE, extends(atlas_StructuredGrid) :: atlas_ShiftedLonLatGrid + +! Purpose : +! ------- +! *atlas_ShiftedLonLatGrid* : Object Grid specifications for LonLat Grids + +! Methods : +! ------- + +! Author : +! ------ +! 9-Oct-2014 Willem Deconinck *ECMWF* + +!------------------------------------------------------------------------------ +contains +#if FCKIT_FINAL_NOT_INHERITING + final :: atlas_ShiftedLonLatGrid__final_auto +#endif +END TYPE atlas_ShiftedLonLatGrid + +interface atlas_ShiftedLonLatGrid + module procedure atlas_grid_ShiftedLonLat__ctor_int32 + module procedure atlas_grid_ShiftedLonLat__ctor_int64 +end interface + +!------------------------------------------------------------------------------ + +TYPE, extends(atlas_StructuredGrid) :: atlas_ShiftedLonGrid + +! Purpose : +! ------- +! *atlas_ShiftedLonGrid* : Object Grid specifications for LonLat Grids + +! Methods : +! ------- + +! Author : +! ------ +! 9-Oct-2014 Willem Deconinck *ECMWF* + +!------------------------------------------------------------------------------ +contains +#if FCKIT_FINAL_NOT_INHERITING + final :: atlas_ShiftedLonGrid__final_auto +#endif +END TYPE atlas_ShiftedLonGrid + +interface atlas_ShiftedLonGrid + module procedure atlas_grid_ShiftedLon__ctor_int32 + module procedure atlas_grid_ShiftedLon__ctor_int64 +end interface + +!------------------------------------------------------------------------------ + +TYPE, extends(atlas_StructuredGrid) :: atlas_ShiftedLatGrid + +! Purpose : +! ------- +! *atlas_ShiftedLatGrid* : Object Grid specifications for LonLat Grids + +! Methods : +! ------- + +! Author : +! ------ +! 9-Oct-2014 Willem Deconinck *ECMWF* + +!------------------------------------------------------------------------------ +contains +#if FCKIT_FINAL_NOT_INHERITING + final :: atlas_ShiftedLatGrid__final_auto +#endif +END TYPE atlas_ShiftedLatGrid + +interface atlas_ShiftedLatGrid + module procedure atlas_grid_ShiftedLat__ctor_int32 + module procedure atlas_grid_ShiftedLat__ctor_int64 +end interface + +!------------------------------------------------------------------------------ + interface atlas_RegionalGrid module procedure atlas_RegionalGrid_ctor_int32 module procedure atlas_RegionalGrid_ctor_int64 @@ -382,6 +466,30 @@ ATLAS_FINAL subroutine atlas_RegularLonLatGrid__final_auto(this) FCKIT_SUPPRESS_UNUSED( this ) end subroutine +ATLAS_FINAL subroutine atlas_ShiftedLonLatGrid__final_auto(this) + type(atlas_ShiftedLonLatGrid), intent(inout) :: this +#if FCKIT_FINAL_NOT_PROPAGATING + call this%final() +#endif + FCKIT_SUPPRESS_UNUSED( this ) +end subroutine + +ATLAS_FINAL subroutine atlas_ShiftedLonGrid__final_auto(this) + type(atlas_ShiftedLonGrid), intent(inout) :: this +#if FCKIT_FINAL_NOT_PROPAGATING + call this%final() +#endif + FCKIT_SUPPRESS_UNUSED( this ) +end subroutine + +ATLAS_FINAL subroutine atlas_ShiftedLatGrid__final_auto(this) + type(atlas_ShiftedLatGrid), intent(inout) :: this +#if FCKIT_FINAL_NOT_PROPAGATING + call this%final() +#endif + FCKIT_SUPPRESS_UNUSED( this ) +end subroutine + ATLAS_FINAL subroutine atlas_RegularGaussianGrid__final_auto(this) type(atlas_RegularGaussianGrid), intent(inout) :: this #if FCKIT_FINAL_NOT_PROPAGATING @@ -580,6 +688,66 @@ function atlas_grid_RegularLonLat__ctor_int64(nlon,nlat) result(this) call this%return() end function +!----------------------------------------------------------------------------- + +function atlas_grid_ShiftedLonLat__ctor_int32(nlon,nlat) result(this) + use, intrinsic :: iso_c_binding, only: c_int, c_long + use atlas_grid_Structured_c_binding + type(atlas_ShiftedLonLatGrid) :: this + integer(c_int), intent(in) :: nlon, nlat + call this%reset_c_ptr( atlas__grid__regular__ShiftedLonLat(int(nlon,c_long),int(nlat,c_long)) ) + call this%return() +end function + +function atlas_grid_ShiftedLonLat__ctor_int64(nlon,nlat) result(this) + use, intrinsic :: iso_c_binding, only: c_long + use atlas_grid_Structured_c_binding + type(atlas_ShiftedLonLatGrid) :: this + integer(c_long), intent(in) :: nlon, nlat + call this%reset_c_ptr( atlas__grid__regular__ShiftedLonLat( nlon, nlat ) ) + call this%return() +end function + +!----------------------------------------------------------------------------- + +function atlas_grid_ShiftedLon__ctor_int32(nlon,nlat) result(this) + use, intrinsic :: iso_c_binding, only: c_int, c_long + use atlas_grid_Structured_c_binding + type(atlas_ShiftedLonGrid) :: this + integer(c_int), intent(in) :: nlon, nlat + call this%reset_c_ptr( atlas__grid__regular__ShiftedLon(int(nlon,c_long),int(nlat,c_long)) ) + call this%return() +end function + +function atlas_grid_ShiftedLon__ctor_int64(nlon,nlat) result(this) + use, intrinsic :: iso_c_binding, only: c_long + use atlas_grid_Structured_c_binding + type(atlas_ShiftedLonGrid) :: this + integer(c_long), intent(in) :: nlon, nlat + call this%reset_c_ptr( atlas__grid__regular__ShiftedLon( nlon, nlat ) ) + call this%return() +end function + +!----------------------------------------------------------------------------- + +function atlas_grid_ShiftedLat__ctor_int32(nlon,nlat) result(this) + use, intrinsic :: iso_c_binding, only: c_int, c_long + use atlas_grid_Structured_c_binding + type(atlas_ShiftedLatGrid) :: this + integer(c_int), intent(in) :: nlon, nlat + call this%reset_c_ptr( atlas__grid__regular__ShiftedLat(int(nlon,c_long),int(nlat,c_long)) ) + call this%return() +end function + +function atlas_grid_ShiftedLat__ctor_int64(nlon,nlat) result(this) + use, intrinsic :: iso_c_binding, only: c_long + use atlas_grid_Structured_c_binding + type(atlas_ShiftedLatGrid) :: this + integer(c_long), intent(in) :: nlon, nlat + call this%reset_c_ptr( atlas__grid__regular__ShiftedLat( nlon, nlat ) ) + call this%return() +end function + ! ----------------------------------------------------------------------------- ! Structured members diff --git a/src/tests/grid/fctest_grids.F90 b/src/tests/grid/fctest_grids.F90 index 96997123a..0e298438c 100644 --- a/src/tests/grid/fctest_grids.F90 +++ b/src/tests/grid/fctest_grids.F90 @@ -67,6 +67,12 @@ json = spec%json() FCTEST_CHECK( json == json_sorted .or. json == json_ordered ) + grid = atlas_RegularGaussianGrid(8) + grid = atlas_RegularLonLatGrid(8,8) + grid = atlas_ShiftedLonLatGrid(8,8) + grid = atlas_ShiftedLonGrid(8,8) + grid = atlas_ShiftedLatGrid(8,8) + call spec%final() call grid%final() From f9d326e8ac814cef95afca8afa6c40500f1f3854 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Mon, 9 Jan 2023 16:07:13 +0100 Subject: [PATCH 13/25] Cosmetic change --- src/atlas/grid/detail/grid/Healpix.cc | 22 +++++++++++++++++----- 1 file changed, 17 insertions(+), 5 deletions(-) diff --git a/src/atlas/grid/detail/grid/Healpix.cc b/src/atlas/grid/detail/grid/Healpix.cc index d4e435db3..08baa68b6 100644 --- a/src/atlas/grid/detail/grid/Healpix.cc +++ b/src/atlas/grid/detail/grid/Healpix.cc @@ -65,18 +65,23 @@ namespace grid { Healpix::XSpace healpix_xspace(long N) { std::vector vp(4 * N - 1); + + // Polar caps for (int r = 1; r < N; r++) { - double start = 180. * M_1_PI * M_PI_4 / r; + double start = 45./r; vp[r - 1] = LinearSpacing(start, start + 360., 4 * r, false); - vp[4 * N - r - 1] = LinearSpacing(start, start + 360., 4 * r, false); + vp[4 * N - r - 1] = vp[r - 1]; } - double start = 180. * M_1_PI * M_PI_4 / N; + // Equatorial belt + const double start = 45. / N; for (int r = N; r < 2 * N; r++) { double r_start = start * (2. - (r - N + 1) % 2); vp[r - 1] = LinearSpacing(r_start, r_start + 360., 4 * N, false); - vp[4 * N - r - 1] = LinearSpacing(r_start, r_start + 360., 4 * N, false); + vp[4 * N - r - 1] = vp[r - 1]; } + + // Equator double r_start = start * (1 - (N % 2 ? 1 : 0)); vp[2 * N - 1] = LinearSpacing(r_start, r_start + 360., 4 * N, false); @@ -86,15 +91,22 @@ Healpix::XSpace healpix_xspace(long N) { Healpix::YSpace healpix_yspace(long N) { constexpr double rad2deg = util::Constants::radiansToDegrees(); std::vector y(4 * N - 1); - y[2 * N - 1] = 0.; + + // Polar caps for (int r = 1; r < N; r++) { y[r - 1] = 90. - rad2deg * std::acos(1. - r * r / (3. * N * N)); y[4 * N - 1 - r] = -y[r - 1]; } + + // Equatorial belt for (int r = N; r < 2 * N; r++) { y[r - 1] = 90. - rad2deg * std::acos((4. * N - 2. * r) / (3. * N)); y[4 * N - 1 - r] = -y[r - 1]; } + + // Equator + y[2 * N - 1] = 0.; + return new spacing::CustomSpacing(y.size(), y.data()); } From 5680f67ec14eafd06713c62a525d2926b3d6360d Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Mon, 9 Jan 2023 17:21:30 +0100 Subject: [PATCH 14/25] ATLAS-373 Enable gaussian latitudes computation for very high resolution, improving memory requirement with orders of magnitude --- src/apps/atlas-gaussian-latitudes.cc | 3 + .../grid/detail/spacing/gaussian/Latitudes.cc | 127 ++++++++---------- 2 files changed, 56 insertions(+), 74 deletions(-) diff --git a/src/apps/atlas-gaussian-latitudes.cc b/src/apps/atlas-gaussian-latitudes.cc index c365d5ee4..7641d13c4 100644 --- a/src/apps/atlas-gaussian-latitudes.cc +++ b/src/apps/atlas-gaussian-latitudes.cc @@ -39,6 +39,8 @@ class AtlasGaussianLatitudes : public eckit::Tool { public: AtlasGaussianLatitudes(int argc, char** argv): eckit::Tool(argc, argv) { + atlas::initialise(argc,argv); + do_run = false; bool help = Resource("--help", false); @@ -132,6 +134,7 @@ void AtlasGaussianLatitudes::run() { } std::cout << "};" << std::endl; } + atlas::finalise(); } //------------------------------------------------------------------------------ diff --git a/src/atlas/grid/detail/spacing/gaussian/Latitudes.cc b/src/atlas/grid/detail/spacing/gaussian/Latitudes.cc index 26dffdd55..0f9d003c0 100644 --- a/src/atlas/grid/detail/spacing/gaussian/Latitudes.cc +++ b/src/atlas/grid/detail/spacing/gaussian/Latitudes.cc @@ -97,7 +97,7 @@ namespace { // Anonymous namespace //----------------------------------------------------------------------------- -void legpol_newton_iteration(int kn, const double pfn[], double px, double& pxn, double& pxmod) { +void legpol_newton_iteration(size_t kn, const double pfn[], double px, double& pxn, double& pxmod) { // Routine to perform a single Newton iteration step to find // the zero of the ordinary Legendre polynomial of degree N @@ -112,8 +112,8 @@ void legpol_newton_iteration(int kn, const double pfn[], double px, double& pxn, // PXMOD : PXN-PX (out) double zdlx, zdlk, zdlldn, zdlxn, zdlmod; - int ik; - int kodd = kn % 2; // mod(kn,2) + size_t ik; + size_t kodd = kn % 2; // mod(kn,2) zdlx = px; zdlk = 0.; @@ -126,7 +126,7 @@ void legpol_newton_iteration(int kn, const double pfn[], double px, double& pxn, // Newton interation step // ---------------------- - for (int jn = 2 - kodd; jn <= kn; jn += 2) { + for (size_t jn = 2 - kodd; jn <= kn; jn += 2) { // normalised ordinary Legendre polynomial == \overbar{P_n}^0 zdlk += pfn[ik] * std::cos(static_cast(jn) * zdlx); // normalised derivative == d/d\theta(\overbar{P_n}^0) @@ -140,7 +140,7 @@ void legpol_newton_iteration(int kn, const double pfn[], double px, double& pxn, pxmod = zdlmod; } -void legpol_weight(const int kn, const double pfn[], const double px, double& pw) { +void legpol_weight(const size_t kn, const double pfn[], const double px, double& pw) { // Routine to compute the quadrature weight of the legendre polynomial // Explicit arguments : @@ -155,8 +155,8 @@ void legpol_weight(const int kn, const double pfn[], const double px, double& pw // PXMOD : PXN-PX (out) double zdlx, zdlldn; - int ik; - int kodd = kn % 2; + size_t ik; + size_t kodd = kn % 2; zdlx = px; zdlldn = 0.; @@ -164,7 +164,7 @@ void legpol_weight(const int kn, const double pfn[], const double px, double& pw // Compute weights // --------------- - for (int jn = 2 - kodd; jn <= kn; jn += 2) { + for (size_t jn = 2 - kodd; jn <= kn; jn += 2) { zdlldn -= pfn[ik] * static_cast(jn) * std::sin(static_cast(jn) * zdlx); ++ik; } @@ -173,7 +173,7 @@ void legpol_weight(const int kn, const double pfn[], const double px, double& pw //----------------------------------------------------------------------------- -void legpol_quadrature(const int kn, const double pfn[], double& pl, double& pw, int& kiter, double& pmod) { +void legpol_quadrature(const size_t kn, const double pfn[], double& pl, double& pw, size_t& kiter, double& pmod) { //**** *GAWL * - Routine to perform the Newton loop // Purpose. @@ -193,40 +193,29 @@ void legpol_quadrature(const int kn, const double pfn[], double& pl, double& pw, // KITER Number of iterations (out) // PMOD Last modification (inout) - int iflag, itemax; + constexpr size_t itemax = 20; + constexpr double ztol = std::numeric_limits::epsilon() * 1000.; - double zx = 0; - double zw = 0; - double zxn = 0; + double zx = pl; + double zxn; + double zw; + bool tol_reached = false; - //* 1. Initialization. - // --------------- - - itemax = 20; - zx = pl; - iflag = 0; - pw = 0.; - - //* 2. Newton iteration. - // ----------------- - - const double zeps = std::numeric_limits::epsilon(); - - for (int jter = 1; jter <= itemax + 1; ++jter) { + for (size_t jter = 1; jter <= itemax + 1; ++jter) { kiter = jter; legpol_newton_iteration(kn, pfn, zx, zxn, pmod); zx = zxn; - if (iflag == 1) { + if (tol_reached) { legpol_weight(kn, pfn, zx, zw); break; } - if (std::abs(pmod) <= zeps * 1000.) { - iflag = 1; + if (std::abs(pmod) <= ztol) { + tol_reached = true; } } - if (iflag != 1) { + if (not tol_reached) { std::stringstream s; - s << "Could not converge gaussian latitude to accuracy [" << zeps * 1000 << "]\n"; + s << "Could not converge gaussian latitude to accuracy [" << ztol << "]\n"; s << "after " << itemax << " iterations. Consequently also failed to compute quadrature weight."; throw_Exception(s.str(), Here()); } @@ -242,59 +231,49 @@ void legpol_quadrature(const int kn, const double pfn[], double& pl, double& pw, //----------------------------------------------------------------------------- void compute_gaussian_quadrature_npole_equator(const size_t N, double lats[], double weights[]) { - Log::debug() << "Atlas computing Gaussian latitudes for N " << N << " which requires temporary memory of " - << eckit::Bytes(sizeof(double) * (2 * N + 1) * (2 * N + 1)) << std::endl; + Log::debug() << "Atlas computing Gaussian latitudes for N " << N << std::endl; ATLAS_TRACE(); - // Compute first guess for colatitudes in radians - double z; - for (size_t i = 0; i < N; ++i) { - z = (4. * (i + 1.) - 1.) * M_PI / (4. * 2. * N + 2.); - lats[i] = (z + 1. / (tan(z) * (8. * (2. * N) * (2. * N)))); - } - - int kdgl = 2 * N; - array::ArrayT zfn_(kdgl + 1, kdgl + 1); - // WARNING: potentially HUGE allocation ( N=16000 --> 7.6 Gbytes ) - - ArrayView zfn = make_view(zfn_); - - int iodd; - - // Belousov, Swarztrauber use zfn(0,0)=std::sqrt(2.) - // IFS normalisation chosen to be 0.5*Integral(Pnm**2) = 1 - zfn(0, 0) = 2.; - for (int jn = 1; jn <= kdgl; ++jn) { - double zfnn = zfn(0, 0); - for (int jgl = 1; jgl <= jn; ++jgl) { - zfnn *= std::sqrt(1. - 0.25 / (static_cast(jgl * jgl))); - } - iodd = jn % 2; - zfn(jn, jn) = zfnn; - for (int jgl = 2; jgl <= jn - iodd; jgl += 2) { - zfn(jn, jn - jgl) = zfn(jn, jn - jgl + 2) * static_cast((jgl - 1) * (2 * jn - jgl + 2)) / - static_cast(jgl * (2 * jn - jgl + 1)); - } - } - - iodd = kdgl % 2; - int ik = iodd; + size_t kdgl = 2 * N; + // Computation of Fourier coefficients of series expansion for the ordinary Legendre polynomials std::vector zzfn(N + 1); - std::vector zmod(kdgl); - std::vector iter(kdgl); + { + // Belousov, Swarztrauber use zfn(0)=std::sqrt(2.) + // IFS normalisation chosen to be 0.5*Integral(Pnm**2) = 1 + constexpr double zfn_0 = 2.; + std::vector zfn(kdgl+1); + zfn[0] = zfn_0; + for (size_t jn = 1; jn <= kdgl; ++jn) { + zfn[jn] = zfn_0; + for (size_t jgl = 1; jgl <= jn; ++jgl) { + zfn[jn] *= std::sqrt(1. - 0.25 / (static_cast(jgl * jgl))); + } + size_t iodd = jn % 2; + for (size_t jgl = 2; jgl <= jn - iodd; jgl += 2) { + zfn[jn - jgl] = zfn[jn - jgl + 2] * static_cast((jgl - 1) * (2 * jn - jgl + 2)) / + static_cast(jgl * (2 * jn - jgl + 1)); + } + } - for (int jgl = iodd; jgl <= kdgl; jgl += 2) { - zzfn[ik] = zfn(kdgl, jgl); - ++ik; + size_t iodd = kdgl % 2; + for (size_t jgl = iodd, ik = iodd; jgl <= kdgl; jgl += 2, ++ik) { + zzfn[ik] = zfn[jgl]; + } } - const double pole = 90.; for (size_t jgl = 0; jgl < N; ++jgl) { + // Compute first guess for colatitudes in radians + double z = (4. * (jgl + 1.) - 1.) * M_PI / (4. * 2. * N + 2.); + lats[jgl] = (z + 1. / (std::tan(z) * (8. * (2. * N) * (2. * N)))); + // refine colat first guess here via Newton's method - legpol_quadrature(kdgl, zzfn.data(), lats[jgl], weights[jgl], iter[jgl], zmod[jgl]); + size_t iter; + double zmod; + legpol_quadrature(kdgl, zzfn.data(), lats[jgl], weights[jgl], iter, zmod); // Convert colat to lat, in degrees + constexpr double pole = 90.; lats[jgl] = pole - lats[jgl] * util::Constants::radiansToDegrees(); } } From 256786711015e70c95ec71d1fbdcb9b744906efb Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Tue, 10 Jan 2023 10:02:00 +0100 Subject: [PATCH 15/25] Github Actions: Update macos-11 to macos-12 (latest available) --- .github/workflows/build.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index b37ee6934..601fc29f6 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -85,8 +85,8 @@ jobs: - name: macos # Xcode compiler requires empty environment variables, so we pass null (~) here - os: macos-11 - compiler: clang-12 + os: macos-12 + compiler: clang-14 compiler_cc: ~ compiler_cxx: ~ compiler_fc: gfortran-11 From b8aba2b189d0a04013b0c97b915346d492568dc6 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Tue, 10 Jan 2023 10:23:44 +0100 Subject: [PATCH 16/25] Github Actions: Prevent homebrew from auto-updating --- .github/workflows/build.yml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index 601fc29f6..ef6315dde 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -110,6 +110,9 @@ jobs: echo "FC=${{ matrix.compiler_fc }}" >> $GITHUB_ENV if [[ "${{ matrix.os }}" =~ macos ]]; then + export HOMEBREW_NO_INSTALLED_DEPENDENTS_CHECK=1 + export HOMEBREW_NO_AUTO_UPDATE=1 + export HOMEBREW_NO_INSTALL_CLEANUP=1 brew install ninja brew install libomp else From a63cff3cca830d554fba181d5cf85aaaac6e0293 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Tue, 10 Jan 2023 10:40:03 +0100 Subject: [PATCH 17/25] Github Actions: Prevent homebrew from auto-updating --- .github/workflows/build.yml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index ef6315dde..f61efce08 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -113,6 +113,9 @@ jobs: export HOMEBREW_NO_INSTALLED_DEPENDENTS_CHECK=1 export HOMEBREW_NO_AUTO_UPDATE=1 export HOMEBREW_NO_INSTALL_CLEANUP=1 + echo "HOMEBREW_NO_INSTALLED_DEPENDENTS_CHECK=1" >> $GITHUB_ENV + echo "HOMEBREW_NO_AUTO_UPDATE=1" >> $GITHUB_ENV + echo "HOMEBREW_NO_INSTALL_CLEANUP=1" >> $GITHUB_ENV brew install ninja brew install libomp else From f8ac720e6fd9de390ec65100487edcfe521bfb8b Mon Sep 17 00:00:00 2001 From: Toby Searle <14909402+twsearle@users.noreply.github.com> Date: Wed, 11 Jan 2023 09:01:26 +0000 Subject: [PATCH 18/25] Wrap-around the longitudes in the polygon locator (#108) * Wrap-around the longitudes in the polygon locator * correct error message --- src/atlas/mesh/detail/PartitionGraph.cc | 1 + src/atlas/util/PolygonLocator.h | 24 +++++++++++++++++++++++- 2 files changed, 24 insertions(+), 1 deletion(-) diff --git a/src/atlas/mesh/detail/PartitionGraph.cc b/src/atlas/mesh/detail/PartitionGraph.cc index c498b8aac..ac3bb31cb 100644 --- a/src/atlas/mesh/detail/PartitionGraph.cc +++ b/src/atlas/mesh/detail/PartitionGraph.cc @@ -29,6 +29,7 @@ namespace detail { //---------------------------------------------------------------------------------------------------------------------- PartitionGraph* build_partition_graph(const MeshImpl& mesh) { + ATLAS_TRACE("build_partition_graph..."); const eckit::mpi::Comm& comm = mpi::comm(); const int mpi_size = int(comm.size()); diff --git a/src/atlas/util/PolygonLocator.h b/src/atlas/util/PolygonLocator.h index b6e75720d..a119ce303 100644 --- a/src/atlas/util/PolygonLocator.h +++ b/src/atlas/util/PolygonLocator.h @@ -95,9 +95,31 @@ class PolygonLocator { Log::info() << "NOT_FOUND" << std::endl; #endif } + + auto try_shifted_point = [&](double x) { + const Point2 shifted{x, point.y()}; + const auto shifted_found = kdtree_.closestPoints(shifted, k_); + for (size_t i = 0; i < shifted_found.size(); ++i) { + idx_t ii = shifted_found[i].payload(); + if (polygons_[ii].contains(lonlat2xy(shifted))) { + return ii; + } + } + return -1; + }; + + if (partition < 0) { + partition = try_shifted_point(point.x() + 360.0); + } + if (partition < 0) { + partition = try_shifted_point(point.x() - 360.0); + } if (partition < 0) { std::stringstream out; - out << "Could not find find point {lon,lat} = " << point << " in `k=" << k_ << "` \"nearest\" polygons ["; + out << "Could not find find point {lon,lat} = " << point + << " or " << Point2{point.x() + 360.0, point.y()} + << " or " << Point2{point.x() - 360.0, point.y()} + << " in `k=" << k_ << "` \"nearest\" polygons ["; for (size_t i = 0; i < found.size(); ++i) { if (i > 0) { out << ", "; From 6b194dff56fc40b068d9b2775d3d96db195fd644 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Wed, 11 Jan 2023 10:04:14 +0100 Subject: [PATCH 19/25] Remove assertion checking for empty mesh partition in NodeColumns (see #108) --- src/atlas/functionspace/NodeColumns.cc | 1 - 1 file changed, 1 deletion(-) diff --git a/src/atlas/functionspace/NodeColumns.cc b/src/atlas/functionspace/NodeColumns.cc index 6847ce0e5..8286f1eed 100644 --- a/src/atlas/functionspace/NodeColumns.cc +++ b/src/atlas/functionspace/NodeColumns.cc @@ -206,7 +206,6 @@ NodeColumns::NodeColumns(Mesh mesh, const eckit::Configuration& config): else { halo_ = mesh::Halo(mesh); } - ATLAS_ASSERT(mesh_); mesh::actions::build_nodes_parallel_fields(mesh_.nodes()); mesh::actions::build_periodic_boundaries(mesh_); From 1585d7a1cdc1d6b8306966d2c8e8fc5b18bf9daa Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Mon, 16 Jan 2023 10:44:06 +0000 Subject: [PATCH 20/25] Allow interpolation to partitions without points or elements --- .../method/cubedsphere/CubedSphereBilinear.cc | 7 +++++++ src/atlas/interpolation/method/knn/KNearestNeighbours.cc | 7 +++++++ src/atlas/interpolation/method/knn/NearestNeighbour.cc | 7 +++++++ .../method/structured/StructuredInterpolation2D.tcc | 6 ++++++ .../interpolation/method/unstructured/FiniteElement.cc | 7 +++++++ 5 files changed, 34 insertions(+) diff --git a/src/atlas/interpolation/method/cubedsphere/CubedSphereBilinear.cc b/src/atlas/interpolation/method/cubedsphere/CubedSphereBilinear.cc index 0fb1a7793..ac0d76353 100644 --- a/src/atlas/interpolation/method/cubedsphere/CubedSphereBilinear.cc +++ b/src/atlas/interpolation/method/cubedsphere/CubedSphereBilinear.cc @@ -33,6 +33,13 @@ void CubedSphereBilinear::do_setup(const FunctionSpace& source, const FunctionSp ATLAS_ASSERT(ncSource); ATLAS_ASSERT(target_); + // return early if no output points on this partition reserve is called on + // the triplets but also during the sparseMatrix constructor. This won't + // work for empty matrices + if (target_.size() == 0) { + return; + } + const auto finder = cubedsphere::CellFinder(ncSource.mesh(), util::Config("halo", halo_)); // Numeric tolerance should scale with N. diff --git a/src/atlas/interpolation/method/knn/KNearestNeighbours.cc b/src/atlas/interpolation/method/knn/KNearestNeighbours.cc index d17839489..4ac350962 100644 --- a/src/atlas/interpolation/method/knn/KNearestNeighbours.cc +++ b/src/atlas/interpolation/method/knn/KNearestNeighbours.cc @@ -74,6 +74,13 @@ void KNearestNeighbours::do_setup(const FunctionSpace& source, const FunctionSpa size_t inp_npts = source.size(); size_t out_npts = target.size(); + // return early if no output points on this partition reserve is called on + // the triplets but also during the sparseMatrix constructor. This won't + // work for empty matrices + if (out_npts == 0) { + return; + } + // fill the sparse matrix std::vector weights_triplets; weights_triplets.reserve(out_npts * k_); diff --git a/src/atlas/interpolation/method/knn/NearestNeighbour.cc b/src/atlas/interpolation/method/knn/NearestNeighbour.cc index 515f9249c..65f36a46d 100644 --- a/src/atlas/interpolation/method/knn/NearestNeighbour.cc +++ b/src/atlas/interpolation/method/knn/NearestNeighbour.cc @@ -68,6 +68,13 @@ void NearestNeighbour::do_setup(const FunctionSpace& source, const FunctionSpace size_t inp_npts = source.size(); size_t out_npts = target.size(); + // return early if no output points on this partition reserve is called on + // the triplets but also during the sparseMatrix constructor. This won't + // work for empty matrices + if (out_npts == 0) { + return; + } + // fill the sparse matrix std::vector weights_triplets; weights_triplets.reserve(out_npts); diff --git a/src/atlas/interpolation/method/structured/StructuredInterpolation2D.tcc b/src/atlas/interpolation/method/structured/StructuredInterpolation2D.tcc index 2ee6cc2c0..4375ba4e0 100644 --- a/src/atlas/interpolation/method/structured/StructuredInterpolation2D.tcc +++ b/src/atlas/interpolation/method/structured/StructuredInterpolation2D.tcc @@ -152,6 +152,12 @@ void StructuredInterpolation2D::setup( const FunctionSpace& source ) { idx_t inp_npts = source.size(); idx_t out_npts = target_lonlat_.shape( 0 ); + // return early if no output points on this partition reserve is called on + // the triplets but also during the sparseMatrix constructor. This won't + // work for empty matrices + if (out_npts == 0) { + return; + } auto lonlat = array::make_view( target_lonlat_ ); diff --git a/src/atlas/interpolation/method/unstructured/FiniteElement.cc b/src/atlas/interpolation/method/unstructured/FiniteElement.cc index 8584ae1a1..74c3e63c4 100644 --- a/src/atlas/interpolation/method/unstructured/FiniteElement.cc +++ b/src/atlas/interpolation/method/unstructured/FiniteElement.cc @@ -227,6 +227,13 @@ void FiniteElement::setup(const FunctionSpace& source) { idx_t inp_npts = i_nodes.size(); idx_t out_npts = ocoords_->shape(0); + // return early if no output points on this partition reserve is called on + // the triplets but also during the sparseMatrix constructor. This won't + // work for empty matrices + if (out_npts == 0) { + return; + } + array::ArrayView out_ghosts = array::make_view(target_ghost_); array::ArrayView out_lonlat = array::make_view(target_lonlat_); From e507bf2299237ecf7cadd4b03b30cc7f931d52d6 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Mon, 16 Jan 2023 10:47:03 +0000 Subject: [PATCH 21/25] Adapt StructuredMeshGenerator to allow creation of meshes where some mpi tasks have no elements or nodes --- .../detail/StructuredMeshGenerator.cc | 32 +++++++++++++------ 1 file changed, 22 insertions(+), 10 deletions(-) diff --git a/src/atlas/meshgenerator/detail/StructuredMeshGenerator.cc b/src/atlas/meshgenerator/detail/StructuredMeshGenerator.cc index 87420300f..8547f59f1 100644 --- a/src/atlas/meshgenerator/detail/StructuredMeshGenerator.cc +++ b/src/atlas/meshgenerator/detail/StructuredMeshGenerator.cc @@ -54,12 +54,12 @@ static double to_deg = 180. * M_1_PI; } // namespace struct Region { - int north; - int south; + int north {-1}; + int south {-1}; std::unique_ptr elems; - int ntriags; - int nquads; - int nnodes; + int ntriags {0}; + int nquads {0}; + int nnodes {0}; std::vector lat_begin; std::vector lat_end; std::vector nb_lat_elems; @@ -272,6 +272,9 @@ void StructuredMeshGenerator::generate_region(const StructuredGrid& rg, const gr bool unique_pole = options.getBool("unique_pole") && three_dimensional && has_north_pole && has_south_pole; bool periodic_east_west = rg.periodic(); + region.lat_begin.resize(rg.ny(), -1); + region.lat_end.resize(rg.ny(), -1); + region.nb_lat_elems.resize(rg.ny(), 0); idx_t lat_north = -1; idx_t lat_south = -1; @@ -351,6 +354,13 @@ Find min and max latitudes used by this part. } } + if (lat_north == -1 && lat_south == -1 ) { + return; + } + + ATLAS_ASSERT(lat_north >= 0); + ATLAS_ASSERT(lat_south < rg.ny() ); + std::vector offset(rg.ny(), 0); int n = 0; @@ -368,9 +378,7 @@ We need to connect to next region if (idx_t(lat_south + 1) <= rg.ny() - 1 && rg.nx(lat_south + 1) > 0) { ++lat_south; } - region.lat_begin.resize(rg.ny(), -1); - region.lat_end.resize(rg.ny(), -1); - region.nb_lat_elems.resize(rg.ny(), 0); + region.north = lat_north; region.south = lat_south; @@ -993,8 +1001,11 @@ void StructuredMeshGenerator::generate_mesh(const StructuredGrid& rg, const grid auto flags = array::make_view(nodes.flags()); auto halo = array::make_view(nodes.halo()); + idx_t jnorth = -1; + idx_t jsouth = -1; std::vector node_numbering(node_numbering_size, -1); + if (nnodes > 0) { if (options.getBool("ghost_at_end")) { std::vector ghost_nodes; ghost_nodes.reserve(nnodes); @@ -1165,7 +1176,6 @@ void StructuredMeshGenerator::generate_mesh(const StructuredGrid& rg, const grid } } - idx_t jnorth = -1; if (include_north_pole) { idx_t inode = node_numbering.at(jnode); jnorth = jnode; @@ -1189,7 +1199,6 @@ void StructuredMeshGenerator::generate_mesh(const StructuredGrid& rg, const grid ++jnode; } - idx_t jsouth = -1; if (include_south_pole) { idx_t inode = node_numbering.at(jnode); jsouth = jnode; @@ -1212,6 +1221,7 @@ void StructuredMeshGenerator::generate_mesh(const StructuredGrid& rg, const grid Topology::set(flags(inode), Topology::SOUTH); ++jnode; } + } mesh.metadata().set("nb_nodes_including_halo[0]", nodes.size()); nodes.metadata().set("NbRealPts", size_t(nnodes - nnewnodes)); @@ -1268,6 +1278,7 @@ void StructuredMeshGenerator::generate_mesh(const StructuredGrid& rg, const grid regular_cells_glb_idx = false; } + if ((region.nquads + region.ntriags) > 0) { for (idx_t jlat = region.north; jlat < region.south; ++jlat) { idx_t ilat = jlat - region.north; idx_t jlatN = jlat; @@ -1568,6 +1579,7 @@ void StructuredMeshGenerator::generate_mesh(const StructuredGrid& rg, const grid } } } + } if (not regular_cells_glb_idx) { generateGlobalElementNumbering(mesh); } From d9f284832c5c20f4aeacdadcceaec6ad770095d0 Mon Sep 17 00:00:00 2001 From: Toby Searle <14909402+twsearle@users.noreply.github.com> Date: Mon, 16 Jan 2023 12:20:33 +0000 Subject: [PATCH 22/25] Serial distribution to use entire mesh on every partition (#114) * Use entire mesh on every partition Co-authored-by: Willem Deconinck --- src/atlas/grid/detail/distribution/SerialDistribution.cc | 2 ++ src/atlas/grid/detail/distribution/SerialDistribution.h | 6 ++++-- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/src/atlas/grid/detail/distribution/SerialDistribution.cc b/src/atlas/grid/detail/distribution/SerialDistribution.cc index 71e915a32..a9b1b16a9 100644 --- a/src/atlas/grid/detail/distribution/SerialDistribution.cc +++ b/src/atlas/grid/detail/distribution/SerialDistribution.cc @@ -14,6 +14,7 @@ #include #include "atlas/grid/Grid.h" +#include "atlas/parallel/mpi/mpi.h" namespace atlas { namespace grid { @@ -23,6 +24,7 @@ namespace distribution { SerialDistribution::SerialDistribution(const Grid& grid): DistributionFunctionT(grid) { type_ = "serial"; nb_partitions_ = 1; + rank_ = mpi::rank(); size_ = grid.size(); nb_pts_.resize(nb_partitions_, grid.size()); max_pts_ = *std::max_element(nb_pts_.begin(), nb_pts_.end()); diff --git a/src/atlas/grid/detail/distribution/SerialDistribution.h b/src/atlas/grid/detail/distribution/SerialDistribution.h index 6732777f3..090d1a54c 100644 --- a/src/atlas/grid/detail/distribution/SerialDistribution.h +++ b/src/atlas/grid/detail/distribution/SerialDistribution.h @@ -10,7 +10,6 @@ #pragma once - #include "atlas/grid/detail/distribution/DistributionFunction.h" namespace atlas { @@ -22,7 +21,10 @@ class SerialDistribution : public DistributionFunctionT { public: SerialDistribution(const Grid& grid); - ATLAS_ALWAYS_INLINE int function(gidx_t gidx) const { return 0; } + ATLAS_ALWAYS_INLINE int function(gidx_t gidx) const { return rank_; } + +private: + int rank_{0}; }; } // namespace distribution From 9e0b26d026208916b8d9bd79c461865cdfb004d6 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Mon, 16 Jan 2023 16:09:22 +0100 Subject: [PATCH 23/25] Add function SolidBodyRotation --- src/atlas/CMakeLists.txt | 2 + src/atlas/util/function/SolidBodyRotation.cc | 149 ++++++++++++++++++ src/atlas/util/function/SolidBodyRotation.h | 54 +++++++ .../atlas-conservative-interpolation.cc | 10 +- .../numerics/test_fvm_nabla_validation.cc | 80 +--------- 5 files changed, 221 insertions(+), 74 deletions(-) create mode 100644 src/atlas/util/function/SolidBodyRotation.cc create mode 100644 src/atlas/util/function/SolidBodyRotation.h diff --git a/src/atlas/CMakeLists.txt b/src/atlas/CMakeLists.txt index fc1d2bec5..938ba1ab2 100644 --- a/src/atlas/CMakeLists.txt +++ b/src/atlas/CMakeLists.txt @@ -757,6 +757,8 @@ util/detail/BlackMagic.h util/detail/Cache.h util/detail/Debug.h util/detail/KDTree.h +util/function/SolidBodyRotation.h +util/function/SolidBodyRotation.cc util/function/SphericalHarmonic.h util/function/SphericalHarmonic.cc util/function/VortexRollup.h diff --git a/src/atlas/util/function/SolidBodyRotation.cc b/src/atlas/util/function/SolidBodyRotation.cc new file mode 100644 index 000000000..721399953 --- /dev/null +++ b/src/atlas/util/function/SolidBodyRotation.cc @@ -0,0 +1,149 @@ +/* + * (C) Copyright 2023 ECMWF. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + * In applying this licence, ECMWF does not waive the privileges and immunities + * granted to it by virtue of its status as an intergovernmental organisation + * nor does it submit to any jurisdiction. + */ + +#include + +#include "atlas/util/Constants.h" + +#include "atlas/util/function/SolidBodyRotation.h" + +namespace atlas { + +namespace util { + +namespace function { + +SolidBodyRotation::SolidBodyRotation(const double beta, const double radius) { + sin_beta_ = std::sin(beta * Constants::degreesToRadians()); + cos_beta_ = std::cos(beta * Constants::degreesToRadians()); + radius_ = radius; +} + +void SolidBodyRotation::wind(const double lon, const double lat, double& u, double& v) const { + double x = lon * Constants::degreesToRadians(); + double y = lat * Constants::degreesToRadians(); + double cos_x = std::cos(x); + double cos_y = std::cos(y); + double sin_x = std::sin(x); + double sin_y = std::sin(y); + u = cos_y * cos_beta_ + cos_x * sin_y * sin_beta_; + v = -sin_x * sin_beta_; +} + +double SolidBodyRotation::u(const double lon, const double lat) const { + double x = lon * Constants::degreesToRadians(); + double y = lat * Constants::degreesToRadians(); + double cos_x = std::cos(x); + double cos_y = std::cos(y); + double sin_y = std::sin(y); + return cos_y * cos_beta_ + cos_x * sin_y * sin_beta_; +} + +double SolidBodyRotation::v(const double lon, const double lat) const { + double x = lon * Constants::degreesToRadians(); + double sin_x = std::sin(x); + return -sin_x * sin_beta_; +} + +void SolidBodyRotation::vordiv(const double lon, const double lat, double& vor, double& div) const { + double x = lon * Constants::degreesToRadians(); + double y = lat * Constants::degreesToRadians(); + + double cos_x = std::cos(x); + double cos_y = std::cos(y); + double sin_x = std::sin(x); + double sin_y = std::sin(y); + + // Divergence = 1./(R*cos(y)) * ( d/dx( u ) + d/dy( v * cos(y) ) ) + // Vorticity = 1./(R*cos(y)) * ( d/dx( v ) - d/dy( u * cos(y) ) ) + double ddx_u = -sin_x * sin_y * sin_beta_; + double ddy_cosy_v = (-sin_x * sin_beta_) * (-sin_y); + double ddx_v = -cos_x * sin_beta_; + double ddy_cosy_u = + 2 * cos_y * (-sin_y) * cos_beta_ + (-sin_y) * cos_x * sin_y * sin_beta_ + cos_y * cos_x * cos_y * sin_beta_; + + double metric = 1. / (radius_ * cos_y); + + div = metric * (ddx_u + ddy_cosy_v); + vor = metric * (ddx_v - ddy_cosy_u); +} + +double SolidBodyRotation::vorticity(const double lon, const double lat) const { + double x = lon * Constants::degreesToRadians(); + double y = lat * Constants::degreesToRadians(); + + double cos_x = std::cos(x); + double cos_y = std::cos(y); + double sin_y = std::sin(y); + + // Vorticity = 1./(R*cos(y)) * ( d/dx( v ) - d/dy( u * cos(y) ) ) + double ddx_v = -cos_x * sin_beta_; + double ddy_cosy_u = + 2 * cos_y * (-sin_y) * cos_beta_ + (-sin_y) * cos_x * sin_y * sin_beta_ + cos_y * cos_x * cos_y * sin_beta_; + + double metric = 1. / (radius_ * cos_y); + + return metric * (ddx_v - ddy_cosy_u); +} + +double SolidBodyRotation::divergence(const double lon, const double lat) const { + double x = lon * Constants::degreesToRadians(); + double y = lat * Constants::degreesToRadians(); + + double cos_y = std::cos(y); + double sin_x = std::sin(x); + double sin_y = std::sin(y); + + // Divergence = 1./(R*cos(y)) * ( d/dx( u ) + d/dy( v * cos(y) ) ) + double ddx_u = -sin_x * sin_y * sin_beta_; + double ddy_cosy_v = (-sin_x * sin_beta_) * (-sin_y); + + double metric = 1. / (radius_ * cos_y); + + return metric * (ddx_u + ddy_cosy_v); +} + +double SolidBodyRotation::windMagnitude(const double lon, const double lat) const { + return std::sqrt(windMagnitudeSquared(lon, lat)); +} + +double SolidBodyRotation::windMagnitudeSquared(const double lon, const double lat) const { + double u, v; + wind(lon, lat, u, v); + return u * u + v * v; +} + +void SolidBodyRotation::windMagnitudeSquaredGradient(const double lon, const double lat, double& dfdx, double& dfdy) const { + double x = lon * Constants::degreesToRadians(); + double y = lat * Constants::degreesToRadians(); + + double cos_x = std::cos(x); + double cos_y = std::cos(y); + double sin_x = std::sin(x); + double sin_y = std::sin(y); + + double metric_y = 1. / radius_; + double metric_x = metric_y / cos_y; + + double u = cos_y * cos_beta_ + cos_x * sin_y * sin_beta_; + double v = -sin_x * sin_beta_; + double dudx = metric_x * (-sin_x * sin_y * sin_beta_); + double dudy = metric_y * (-sin_y * cos_beta_ + cos_x * cos_y * sin_beta_); + double dvdx = metric_x * (-cos_x * sin_beta_); + double dvdy = metric_y * (0.); + dfdx = 2 * u * dudx + 2 * v * dvdx; + dfdy = 2 * u * dudy + 2 * v * dvdy; +} + +} // namespace function + +} // namespace util + +} // namespace atlas diff --git a/src/atlas/util/function/SolidBodyRotation.h b/src/atlas/util/function/SolidBodyRotation.h new file mode 100644 index 000000000..e705e1dd1 --- /dev/null +++ b/src/atlas/util/function/SolidBodyRotation.h @@ -0,0 +1,54 @@ +/* + * (C) Copyright 2023 ECMWF. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + * In applying this licence, ECMWF does not waive the privileges and immunities + * granted to it by virtue of its status as an intergovernmental organisation nor + * does it submit to any jurisdiction. + */ + + +#pragma once + +namespace atlas { + +namespace util { + +namespace function { + +/// \brief An analytic function that provides solid body rotation winds on a sphere. +/// +/// All angles must be provided in degrees. +/// +class SolidBodyRotation { +public: + + SolidBodyRotation() : SolidBodyRotation(0., 1.) {} + + + SolidBodyRotation(const double beta) : SolidBodyRotation(beta, 1.) {} + SolidBodyRotation(const double beta, const double radius); + + void wind(const double lon, const double lat, double& u, double& v) const; + void vordiv(const double lon, const double lat, double& vor, double& div) const; + double windMagnitude(const double lon, const double lat) const; + double u(const double lon, const double lat) const; + double v(const double lon, const double lat) const; + double vorticity(const double lon, const double lat) const; + double divergence(const double lon, const double lat) const; + + double windMagnitudeSquared(const double lon, const double lat) const; + void windMagnitudeSquaredGradient(const double lon, const double lat, double& dfdx, double& dfdy) const; + +private: + double sin_beta_; + double cos_beta_; + double radius_; +}; + +} // namespace function + +} // namespace util + +} // namespace atlas diff --git a/src/sandbox/interpolation/atlas-conservative-interpolation.cc b/src/sandbox/interpolation/atlas-conservative-interpolation.cc index 36f8b67b4..afe570c92 100644 --- a/src/sandbox/interpolation/atlas-conservative-interpolation.cc +++ b/src/sandbox/interpolation/atlas-conservative-interpolation.cc @@ -45,6 +45,7 @@ #include "atlas/output/Gmsh.h" #include "atlas/runtime/AtlasTool.h" #include "atlas/util/Config.h" +#include "atlas/util/function/SolidBodyRotation.h" #include "atlas/util/function/SphericalHarmonic.h" #include "atlas/util/function/VortexRollup.h" @@ -107,7 +108,8 @@ class AtlasParallelInterpolation : public AtlasTool { add_option(new eckit::option::Separator("Initial condition options")); add_option(new SimpleOption( - "init", "Setup initial source field [ constant, spherical_harmonic, vortex_rollup (default) ]")); + "init", "Setup initial source field [ constant, spherical_harmonic, vortex_rollup (default), solid_body_rotation_wind_magnitude ]")); + add_option(new SimpleOption("solid_body_rotation.angle", "Angle of solid body rotation (default = 0.)")); add_option(new SimpleOption("vortex_rollup.t", "Value that controls vortex rollup (default = 0.5)")); add_option(new SimpleOption("constant.value", "Value that is assigned in case init==constant)")); add_option(new SimpleOption("spherical_harmonic.n", "total wave number 'n' of a spherical harmonic")); @@ -147,6 +149,12 @@ std::function get_init(const AtlasTool::Args& args) args.get("constant.value", value = 1.); return [value](const PointLonLat&) { return value; }; } + else if (init == "solid_body_rotation_wind_magnitude") { + double beta; + args.get("solid_body_rotation.angle", beta = 0.); + util::function::SolidBodyRotation sbr(beta); + return [sbr](const PointLonLat& p) { return sbr.windMagnitude(p.lon(), p.lat()); }; + } else { if (args.has("init")) { Log::error() << "Bad value for \"init\": \"" << init << "\" not recognised." << std::endl; diff --git a/src/tests/numerics/test_fvm_nabla_validation.cc b/src/tests/numerics/test_fvm_nabla_validation.cc index 389740b63..0b9bc351c 100644 --- a/src/tests/numerics/test_fvm_nabla_validation.cc +++ b/src/tests/numerics/test_fvm_nabla_validation.cc @@ -32,6 +32,7 @@ #include "atlas/util/Constants.h" #include "atlas/util/CoordinateEnums.h" #include "atlas/util/Earth.h" +#include "atlas/util/function/SolidBodyRotation.h" #include "tests/AtlasTestEnvironment.h" @@ -79,77 +80,10 @@ static int metric_approach() { return 0; } -class SolidBodyRotation { - double radius; - double beta; - double sin_beta; - double cos_beta; - -public: - SolidBodyRotation(const double _radius, const double _beta): radius{_radius}, beta{_beta} { - sin_beta = std::sin(beta); - cos_beta = std::cos(beta); - } - void wind(const double x, const double y, double& u, double& v) { - double cos_x = std::cos(x); - double cos_y = std::cos(y); - double sin_x = std::sin(x); - double sin_y = std::sin(y); - u = cos_y * cos_beta + cos_x * sin_y * sin_beta; - v = -sin_x * sin_beta; - } - - void vordiv(const double x, const double y, double& vor, double& div) { - double cos_x = std::cos(x); - double cos_y = std::cos(y); - double sin_x = std::sin(x); - double sin_y = std::sin(y); - - // Divergence = 1./(R*cos(y)) * ( d/dx( u ) + d/dy( v * cos(y) ) ) - // Vorticity = 1./(R*cos(y)) * ( d/dx( v ) - d/dy( u * cos(y) ) ) - double ddx_u = -sin_x * sin_y * sin_beta; - double ddy_cosy_v = (-sin_x * sin_beta) * (-sin_y); - double ddx_v = -cos_x * sin_beta; - double ddy_cosy_u = - 2 * cos_y * (-sin_y) * cos_beta + (-sin_y) * cos_x * sin_y * sin_beta + cos_y * cos_x * cos_y * sin_beta; - - double metric = 1. / (radius * cos_y); - - div = metric * (ddx_u + ddy_cosy_v); - vor = metric * (ddx_v - ddy_cosy_u); - } - - void wind_magnitude_squared(const double x, const double y, double& f) { - double u, v; - wind(x, y, u, v); - f = u * u + v * v; - } - - void wind_magnitude_squared_gradient(const double x, const double y, double& dfdx, double& dfdy) { - double cos_x = std::cos(x); - double cos_y = std::cos(y); - double sin_x = std::sin(x); - double sin_y = std::sin(y); - - double metric_y = 1. / radius; - double metric_x = metric_y / cos_y; - - double u = cos_y * cos_beta + cos_x * sin_y * sin_beta; - double v = -sin_x * sin_beta; - double dudx = metric_x * (-sin_x * sin_y * sin_beta); - double dudy = metric_y * (-sin_y * cos_beta + cos_x * cos_y * sin_beta); - double dvdx = metric_x * (-cos_x * sin_beta); - double dvdy = metric_y * (0.); - dfdx = 2 * u * dudx + 2 * v * dvdx; - dfdy = 2 * u * dudy + 2 * v * dvdy; - } -}; - FieldSet analytical_fields(const fvm::Method& fvm) { - constexpr double deg2rad = M_PI / 180.; const double radius = fvm.radius(); - auto lonlat_deg = array::make_view(fvm.mesh().nodes().lonlat()); + auto lonlat = array::make_view(fvm.mesh().nodes().lonlat()); FieldSet fields; auto add_scalar_field = [&](const std::string& name) { @@ -169,20 +103,20 @@ FieldSet analytical_fields(const fvm::Method& fvm) { auto div = add_scalar_field("ref_div"); auto vor = add_scalar_field("ref_vor"); - auto flow = SolidBodyRotation{radius, beta_in_degrees() * deg2rad}; + auto flow = atlas::util::function::SolidBodyRotation{beta_in_degrees(), radius}; auto is_ghost = array::make_view(fvm.mesh().nodes().ghost()); const idx_t nnodes = fvm.mesh().nodes().size(); for (idx_t jnode = 0; jnode < nnodes; ++jnode) { if (is_ghost(jnode)) { continue; } - double x = lonlat_deg(jnode, LON) * deg2rad; - double y = lonlat_deg(jnode, LAT) * deg2rad; + double x = lonlat(jnode, LON); + double y = lonlat(jnode, LAT); flow.wind(x, y, u(jnode), v(jnode)); flow.vordiv(x, y, vor(jnode), div(jnode)); - flow.wind_magnitude_squared(x, y, f(jnode)); - flow.wind_magnitude_squared_gradient(x, y, dfdx(jnode), dfdy(jnode)); + f(jnode) = flow.windMagnitudeSquared(x, y); + flow.windMagnitudeSquaredGradient(x, y, dfdx(jnode), dfdy(jnode)); uv(jnode, XX) = u(jnode); uv(jnode, YY) = v(jnode); From aac67c9d37cd35defa9d2ef7a17ed2ec800eaa2c Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Mon, 23 Jan 2023 11:52:06 +0100 Subject: [PATCH 24/25] Fix compilation with ATLAS_BITS_LOCAL=64 --- src/atlas/util/PolygonLocator.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/atlas/util/PolygonLocator.h b/src/atlas/util/PolygonLocator.h index a119ce303..b64749038 100644 --- a/src/atlas/util/PolygonLocator.h +++ b/src/atlas/util/PolygonLocator.h @@ -105,7 +105,7 @@ class PolygonLocator { return ii; } } - return -1; + return idx_t{-1}; }; if (partition < 0) { From df95b8e9a0ee2f582ad04f98ace799c5a838b217 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Mon, 16 Jan 2023 17:43:51 +0100 Subject: [PATCH 25/25] Version 0.32.0 --- CHANGELOG.md | 18 ++++++++++++++++++ VERSION | 2 +- 2 files changed, 19 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 401372d81..6e19fee71 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,23 @@ This project adheres to [Semantic Versioning](http://semver.org/spec/v2.0.0.html ## [Unreleased] +## [0.32.0] - 2023-01-23 +### Added +- Added BlockStructuredColumns FunctionSpace +- Added function SolidBodyRotation +- Added convenience Fortran constructors for ShiftedLonLat, SHiftedLon, ShiftedLat +- Support for more FunctionSpaces in various interpolation methods + +### Changed +- PolygonLocator now wraps around longitudes for non-matching domains +- StructuredMeshGenerator can generate meshes with partitions without elements +- SerialDistribution makes every MPI task have the entire mesh + +### Fixed +- Remove assertion for checking of empty mesh in NodeColumns FunctionSpace +- Gaussian latitudes can now be computed for very high resolution, without running out of memory. +- Interpolation to functionspaces with empty partitions + ## [0.31.1] - 2022-11-11 ### Fixed - Fix bug introduced in StructuredMeshGenerator global numbering with release 0.31.0 in commit a63fc62a2 @@ -411,6 +428,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.32.0]: https://github.com/ecmwf/atlas/compare/0.31.1...0.32.0 [0.31.1]: https://github.com/ecmwf/atlas/compare/0.31.0...0.31.1 [0.31.0]: https://github.com/ecmwf/atlas/compare/0.30.0...0.31.0 [0.30.0]: https://github.com/ecmwf/atlas/compare/0.29.0...0.30.0 diff --git a/VERSION b/VERSION index f176c9441..9eb2aa3f1 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -0.31.1 +0.32.0