From f1ffac22c8919bac4c847803b0f79e6699bfdecb Mon Sep 17 00:00:00 2001 From: Giacomo Serafini Date: Tue, 17 Mar 2020 15:57:17 +0100 Subject: [PATCH 1/2] update generated and interface --- .gitignore | 1 + README.md | 6 +- scripts/convergence_plot_dec.m | 4 +- stencils/CMakeLists.txt | 5 +- stencils/atlasIconLaplaceDriver.cpp | 6 +- .../{interfaces => driver-includes}/defs.hpp | 0 stencils/driver-includes/extent.hpp | 58 +++++++ .../unstructured_interface.hpp | 1 + stencils/generated_iconLaplace.hpp | 124 +++++++++------ stencils/interfaces/atlas_interface.hpp | 149 ++++++------------ stencils/interfaces/mylib_interface.hpp | 2 +- 11 files changed, 196 insertions(+), 160 deletions(-) create mode 100644 .gitignore rename stencils/{interfaces => driver-includes}/defs.hpp (100%) create mode 100644 stencils/driver-includes/extent.hpp rename stencils/{interfaces => driver-includes}/unstructured_interface.hpp (98%) diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..30d388a --- /dev/null +++ b/.gitignore @@ -0,0 +1 @@ +build* \ No newline at end of file diff --git a/README.md b/README.md index aedb929..b90df91 100644 --- a/README.md +++ b/README.md @@ -6,7 +6,7 @@ This repoistory contains some utilities to transform Atlas meshes as well as the Requirements: -* libnetcdc-c++4 +* libnetcdf-c++4-dev * Atlas * eckit @@ -60,8 +60,8 @@ python3 convergence_plot.py ./(mylib|atlas)IconLaplaceDriver Various error norms are collected for divergence, curl and (normal) vector Laplacian. -Various Octave scripts are provided to assist in visualizing debug and output data. For example, to get a convergence plot: +Various Octave scripts (inside `scripts` directory) are provided to assist in visualizing debug and output data. For example, to get a convergence plot: ``` -octave:1> convergence_plot_dec(' convergence_plot_dec('path/to/conv.csv') ``` diff --git a/scripts/convergence_plot_dec.m b/scripts/convergence_plot_dec.m index edc194a..9b3d922 100644 --- a/scripts/convergence_plot_dec.m +++ b/scripts/convergence_plot_dec.m @@ -12,7 +12,7 @@ function convergence_plot_dec(fname) hcb = legend(['L_{\infty}';'L_1';'L_2';'lin. conv';'quad. conv']); xlabel('log(\Delta x)[°]') ylabel('log(L) [-]') - set(gca,'FontSize',50) - set(hcb,'FontSize',50) + set(gca,'FontSize',20) + set(hcb,'FontSize',20) ylim([-6,0]); endfunction diff --git a/stencils/CMakeLists.txt b/stencils/CMakeLists.txt index 28f701b..321e2b2 100644 --- a/stencils/CMakeLists.txt +++ b/stencils/CMakeLists.txt @@ -1,8 +1,11 @@ add_executable(atlasIconLaplaceDriver atlasIconLaplaceDriver.cpp) target_link_libraries(atlasIconLaplaceDriver atlas eckit atlasUtilsLib) +target_include_directories(atlasIconLaplaceDriver PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}) add_executable(atlasShallowWater shallowWater.cpp) target_link_libraries(atlasShallowWater atlas eckit atlasUtilsLib) +target_include_directories(atlasShallowWater PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}) add_executable(mylibIconLaplaceDriver mylibIconLaplaceDriver.cpp) -target_link_libraries(mylibIconLaplaceDriver atlasUtilsLib myLib) \ No newline at end of file +target_link_libraries(mylibIconLaplaceDriver atlasUtilsLib myLib) +target_include_directories(mylibIconLaplaceDriver PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}) \ No newline at end of file diff --git a/stencils/atlasIconLaplaceDriver.cpp b/stencils/atlasIconLaplaceDriver.cpp index b86d0d3..4418fed 100644 --- a/stencils/atlasIconLaplaceDriver.cpp +++ b/stencils/atlasIconLaplaceDriver.cpp @@ -115,7 +115,7 @@ int main(int argc, char const* argv[]) { // dump a whole bunch of debug output (meant to be visualized using Octave, but gnuplot and the // like will certainly work too) - const bool dbg_out = false; + const bool dbg_out = true; const bool readMeshFromDisk = false; atlas::Mesh mesh; @@ -490,8 +490,8 @@ int main(int argc, char const* argv[]) { primal_edge_length, dual_edge_length, tangent_orientation, geofac_rot, geofac_div) .run(); clock_t end = clock(); - std::cout << "run time Laplacian at resolution " << w << " " - << (end - start) / (double CLOCKS_PER_SEC) << "\n"; + // std::cout << "run time Laplacian at resolution " << w << " " + // << (end - start) / (double CLOCKS_PER_SEC) << "\n"; if(dbg_out) { dumpEdgeField("laplICONatlas_nabla2t1.txt", mesh, wrapper, nabla2t1_vec, level, diff --git a/stencils/interfaces/defs.hpp b/stencils/driver-includes/defs.hpp similarity index 100% rename from stencils/interfaces/defs.hpp rename to stencils/driver-includes/defs.hpp diff --git a/stencils/driver-includes/extent.hpp b/stencils/driver-includes/extent.hpp new file mode 100644 index 0000000..dafe3c0 --- /dev/null +++ b/stencils/driver-includes/extent.hpp @@ -0,0 +1,58 @@ +//===--------------------------------------------------------------------------------*- C++ -*-===// +// _ +// | | +// __| | __ ___ ___ ___ +// / _` |/ _` \ \ /\ / / '_ | +// | (_| | (_| |\ V V /| | | | +// \__,_|\__,_| \_/\_/ |_| |_| - Compiler Toolchain +// +// +// This file is distributed under the MIT License (MIT). +// See LICENSE.txt for details. +// +//===------------------------------------------------------------------------------------------===// + +#pragma once + +#include + +namespace dawn { +namespace driver { + +// Access dimensions with `auto dim_extent = get();`, where D = 0,1,2 (=i,j,k) +// and minus/plus component with `get(dim_extent)`, where C = 0 (minus), 1 (plus) +// (Note, need `using std::get` for ADL to kick-in if std < c++20) +using cartesian_extent = std::array, 3>; + +// Access dimensions with get(), where D = 0,1 (=horizontal, vertical) +// for D == 0 it is either true (have extent) or false (no extent) +// for D == 1, see Cartesian case. +struct unstructured_extent { + bool horizontal; + std::array vertical; +}; + +namespace detail { +template +struct get_impl; + +template <> +struct get_impl<0> { + using type = bool; + constexpr type operator()(unstructured_extent const& extent) { return extent.horizontal; } +}; +template <> +struct get_impl<1> { + using type = std::array; + constexpr type operator()(unstructured_extent const& extent) { return extent.vertical; } +}; + +} // namespace detail + +template +constexpr auto get(unstructured_extent const& extent) -> decltype(detail::get_impl{}(extent)) { + return detail::get_impl{}(extent); +} + +} // namespace driver +} // namespace dawn diff --git a/stencils/interfaces/unstructured_interface.hpp b/stencils/driver-includes/unstructured_interface.hpp similarity index 98% rename from stencils/interfaces/unstructured_interface.hpp rename to stencils/driver-includes/unstructured_interface.hpp index c0911f1..945567e 100644 --- a/stencils/interfaces/unstructured_interface.hpp +++ b/stencils/driver-includes/unstructured_interface.hpp @@ -15,6 +15,7 @@ #pragma once #include "defs.hpp" +#include "extent.hpp" namespace dawn { diff --git a/stencils/generated_iconLaplace.hpp b/stencils/generated_iconLaplace.hpp index 9a434eb..fbbdcf3 100644 --- a/stencils/generated_iconLaplace.hpp +++ b/stencils/generated_iconLaplace.hpp @@ -1,12 +1,18 @@ +//---- Preprocessor defines ---- #define DAWN_GENERATED 1 +#undef DAWN_BACKEND_T #define DAWN_BACKEND_T CXXNAIVEICO -#include "interfaces/unstructured_interface.hpp" +#include + +//---- Globals ---- + +//---- Stencils ---- namespace dawn_generated { namespace cxxnaiveico { template class icon { private: - struct stencil_399 { + struct stencil_64 { dawn::mesh_t const& m_mesh; int m_k_size; dawn::edge_field_t& m_vec; @@ -22,83 +28,101 @@ class icon { dawn::sparse_cell_field_t& m_geofac_div; public: - stencil_399(dawn::mesh_t const& mesh, int k_size, - dawn::edge_field_t& vec, - dawn::cell_field_t& div_vec, - dawn::vertex_field_t& rot_vec, - dawn::edge_field_t& nabla2t1_vec, - dawn::edge_field_t& nabla2t2_vec, - dawn::edge_field_t& nabla2_vec, - dawn::edge_field_t& primal_edge_length, - dawn::edge_field_t& dual_edge_length, - dawn::edge_field_t& tangent_orientation, - dawn::sparse_vertex_field_t& geofac_rot, - dawn::sparse_cell_field_t& geofac_div) + stencil_64(dawn::mesh_t const& mesh, int k_size, + dawn::edge_field_t& vec, dawn::cell_field_t& div_vec, + dawn::vertex_field_t& rot_vec, + dawn::edge_field_t& nabla2t1_vec, + dawn::edge_field_t& nabla2t2_vec, + dawn::edge_field_t& nabla2_vec, + dawn::edge_field_t& primal_edge_length, + dawn::edge_field_t& dual_edge_length, + dawn::edge_field_t& tangent_orientation, + dawn::sparse_vertex_field_t& geofac_rot, + dawn::sparse_cell_field_t& geofac_div) : m_mesh(mesh), m_k_size(k_size), m_vec(vec), m_div_vec(div_vec), m_rot_vec(rot_vec), m_nabla2t1_vec(nabla2t1_vec), m_nabla2t2_vec(nabla2t2_vec), m_nabla2_vec(nabla2_vec), m_primal_edge_length(primal_edge_length), m_dual_edge_length(dual_edge_length), m_tangent_orientation(tangent_orientation), m_geofac_rot(geofac_rot), m_geofac_div(geofac_div) {} - ~stencil_399() {} + ~stencil_64() {} void sync_storages() {} + static constexpr dawn::driver::unstructured_extent vec_extent = {false, 0, 0}; + static constexpr dawn::driver::unstructured_extent div_vec_extent = {false, 0, 0}; + static constexpr dawn::driver::unstructured_extent rot_vec_extent = {false, 0, 0}; + static constexpr dawn::driver::unstructured_extent nabla2t1_vec_extent = {false, 0, 0}; + static constexpr dawn::driver::unstructured_extent nabla2t2_vec_extent = {false, 0, 0}; + static constexpr dawn::driver::unstructured_extent nabla2_vec_extent = {false, 0, 0}; + static constexpr dawn::driver::unstructured_extent primal_edge_length_extent = {false, 0, 0}; + static constexpr dawn::driver::unstructured_extent dual_edge_length_extent = {false, 0, 0}; + static constexpr dawn::driver::unstructured_extent tangent_orientation_extent = {false, 0, 0}; + static constexpr dawn::driver::unstructured_extent geofac_rot_extent = {false, 0, 0}; + static constexpr dawn::driver::unstructured_extent geofac_div_extent = {false, 0, 0}; void run() { using dawn::deref; - ; { for(int k = 0 + 0; k <= (m_k_size == 0 ? 0 : (m_k_size - 1)) + 0 + 0; ++k) { - for(auto const& loc : getVertices(LibTag{}, m_mesh)) { + for(auto const& loc : getCells(LibTag{}, m_mesh)) { int m_sparse_dimension_idx = 0; - m_rot_vec(deref(LibTag{}, loc), k + 0) = reduceEdgeToVertex( - LibTag{}, m_mesh, loc, (::dawn::float_type)0.000000, + m_div_vec(deref(LibTag{}, loc), k + 0) = reduceEdgeToCell( + LibTag{}, m_mesh, loc, (::dawn::float_type)0.0, [&](auto& lhs, auto const& red_loc) { lhs += (m_vec(deref(LibTag{}, red_loc), k + 0) * - m_geofac_rot(deref(LibTag{}, loc), m_sparse_dimension_idx, k + 0)); + m_geofac_div(deref(LibTag{}, loc), m_sparse_dimension_idx, k + 0)); m_sparse_dimension_idx++; return lhs; }); } - for(auto const& loc : getCells(LibTag{}, m_mesh)) { + for(auto const& loc : getEdges(LibTag{}, m_mesh)) { int m_sparse_dimension_idx = 0; - m_div_vec(deref(LibTag{}, loc), k + 0) = reduceEdgeToCell( - LibTag{}, m_mesh, loc, (::dawn::float_type)0.000000, + m_nabla2t2_vec(deref(LibTag{}, loc), k + 0) = + reduceCellToEdge(LibTag{}, m_mesh, loc, (::dawn::float_type)0.0, + [&](auto& lhs, auto const& red_loc, auto const& weight) { + lhs += weight * m_div_vec(deref(LibTag{}, red_loc), k + 0); + m_sparse_dimension_idx++; + return lhs; + }, + std::vector({-1.000000, 1.000000})); + } + for(auto const& loc : getEdges(LibTag{}, m_mesh)) { + m_nabla2t2_vec(deref(LibTag{}, loc), k + 0) = + (m_nabla2t2_vec(deref(LibTag{}, loc), k + 0) / + m_dual_edge_length(deref(LibTag{}, loc), k + 0)); + } + for(auto const& loc : getVertices(LibTag{}, m_mesh)) { + int m_sparse_dimension_idx = 0; + m_rot_vec(deref(LibTag{}, loc), k + 0) = reduceEdgeToVertex( + LibTag{}, m_mesh, loc, (::dawn::float_type)0.0, [&](auto& lhs, auto const& red_loc) { lhs += (m_vec(deref(LibTag{}, red_loc), k + 0) * - m_geofac_div(deref(LibTag{}, loc), m_sparse_dimension_idx, k + 0)); + m_geofac_rot(deref(LibTag{}, loc), m_sparse_dimension_idx, k + 0)); m_sparse_dimension_idx++; return lhs; }); } for(auto const& loc : getEdges(LibTag{}, m_mesh)) { int m_sparse_dimension_idx = 0; - m_nabla2t1_vec(deref(LibTag{}, loc), k + 0) = reduceVertexToEdge( - LibTag{}, m_mesh, loc, (::dawn::float_type)0.000000, - [&](auto& lhs, auto const& red_loc, auto const& weight) { - lhs += weight * m_rot_vec(deref(LibTag{}, red_loc), k + 0); - m_sparse_dimension_idx++; - return lhs; - }, - std::vector({-1.000000, 1.000000})); + m_nabla2t1_vec(deref(LibTag{}, loc), k + 0) = + reduceVertexToEdge(LibTag{}, m_mesh, loc, (::dawn::float_type)0.0, + [&](auto& lhs, auto const& red_loc, auto const& weight) { + lhs += weight * m_rot_vec(deref(LibTag{}, red_loc), k + 0); + m_sparse_dimension_idx++; + return lhs; + }, + std::vector({-1.000000, 1.000000})); + } + for(auto const& loc : getEdges(LibTag{}, m_mesh)) { m_nabla2t1_vec(deref(LibTag{}, loc), k + 0) = ((m_tangent_orientation(deref(LibTag{}, loc), k + 0) * m_nabla2t1_vec(deref(LibTag{}, loc), k + 0)) / m_primal_edge_length(deref(LibTag{}, loc), k + 0)); - m_nabla2t2_vec(deref(LibTag{}, loc), k + 0) = reduceCellToEdge( - LibTag{}, m_mesh, loc, (::dawn::float_type)0.000000, - [&](auto& lhs, auto const& red_loc, auto const& weight) { - lhs += weight * m_div_vec(deref(LibTag{}, red_loc), k + 0); - m_sparse_dimension_idx++; - return lhs; - }, - std::vector({-1.000000, 1.000000})); - m_nabla2t2_vec(deref(LibTag{}, loc), k + 0) = - (m_nabla2t2_vec(deref(LibTag{}, loc), k + 0) / - m_dual_edge_length(deref(LibTag{}, loc), k + 0)); + } + for(auto const& loc : getEdges(LibTag{}, m_mesh)) { m_nabla2_vec(deref(LibTag{}, loc), k + 0) = - (-m_nabla2t1_vec(deref(LibTag{}, loc), k + 0) + - m_nabla2t2_vec(deref(LibTag{}, loc), k + 0)); + (m_nabla2t2_vec(deref(LibTag{}, loc), k + 0) - + m_nabla2t1_vec(deref(LibTag{}, loc), k + 0)); } } } @@ -106,7 +130,7 @@ class icon { } }; static constexpr const char* s_name = "icon"; - stencil_399 m_stencil_399; + stencil_64 m_stencil_64; public: icon(const icon&) = delete; @@ -123,12 +147,12 @@ class icon { dawn::edge_field_t& tangent_orientation, dawn::sparse_vertex_field_t& geofac_rot, dawn::sparse_cell_field_t& geofac_div) - : m_stencil_399(mesh, k_size, vec, div_vec, rot_vec, nabla2t1_vec, nabla2t2_vec, nabla2_vec, - primal_edge_length, dual_edge_length, tangent_orientation, geofac_rot, - geofac_div) {} + : m_stencil_64(mesh, k_size, vec, div_vec, rot_vec, nabla2t1_vec, nabla2t2_vec, nabla2_vec, + primal_edge_length, dual_edge_length, tangent_orientation, geofac_rot, + geofac_div) {} void run() { - m_stencil_399.run(); + m_stencil_64.run(); ; } }; diff --git a/stencils/interfaces/atlas_interface.hpp b/stencils/interfaces/atlas_interface.hpp index d592eab..e1d4bd8 100644 --- a/stencils/interfaces/atlas_interface.hpp +++ b/stencils/interfaces/atlas_interface.hpp @@ -68,17 +68,14 @@ struct atlasTag {}; template class Field { -private: - atlas::array::ArrayView atlas_field_; - public: T const& operator()(int f, int k) const { return atlas_field_(f, k); } T& operator()(int f, int k) { return atlas_field_(f, k); } Field(atlas::array::ArrayView const& atlas_field) : atlas_field_(atlas_field) {} - T* data() { return atlas_field_.data(); } - const T* data() const { return atlas_field_.data(); } - int numElements() const { return atlas_field_.shape(0) * atlas_field_.shape(1); } + +private: + atlas::array::ArrayView atlas_field_; }; template @@ -90,9 +87,6 @@ Field vertexFieldType(atlasTag); template class SparseDimension { -private: - atlas::array::ArrayView sparse_dimension_; - public: T const& operator()(int elem_idx, int sparse_dim_idx, int level) const { return sparse_dimension_(elem_idx, level, sparse_dim_idx); @@ -101,14 +95,11 @@ class SparseDimension { return sparse_dimension_(elem_idx, level, sparse_dim_idx); } - T* data() { return sparse_dimension_.data(); } - const T* data() const { return sparse_dimension_.data(); } - int numElements() const { - return sparse_dimension_.shape(0) * sparse_dimension_.shape(1) * sparse_dimension_.shape(2); - } - SparseDimension(atlas::array::ArrayView const& sparse_dimension) : sparse_dimension_(sparse_dimension) {} + +private: + atlas::array::ArrayView sparse_dimension_; }; template @@ -120,17 +111,11 @@ SparseDimension sparseVertexFieldType(atlasTag); atlas::Mesh meshType(atlasTag); -inline auto getCells(atlasTag, atlas::Mesh const& m) { - return utility::irange(0, m.cells().size()); -} -inline auto getEdges(atlasTag, atlas::Mesh const& m) { - return utility::irange(0, m.edges().size()); -} -inline auto getVertices(atlasTag, atlas::Mesh const& m) { - return utility::irange(0, m.nodes().size()); -} +auto getCells(atlasTag, atlas::Mesh const& m) { return utility::irange(0, m.cells().size()); } +auto getEdges(atlasTag, atlas::Mesh const& m) { return utility::irange(0, m.edges().size()); } +auto getVertices(atlasTag, atlas::Mesh const& m) { return utility::irange(0, m.nodes().size()); } -inline std::vector getNeighs(const atlas::Mesh::HybridElements::Connectivity& conn, int idx) { +std::vector getNeighs(const atlas::Mesh::HybridElements::Connectivity& conn, int idx) { std::vector neighs; for(int n = 0; n < conn.cols(idx); ++n) { neighs.emplace_back(conn(idx, n)); @@ -138,7 +123,7 @@ inline std::vector getNeighs(const atlas::Mesh::HybridElements::Connectivit return neighs; } -inline std::vector getNeighs(const atlas::mesh::Nodes::Connectivity& conn, int idx) { +std::vector getNeighs(const atlas::mesh::Nodes::Connectivity& conn, int idx) { std::vector neighs; for(int n = 0; n < conn.cols(idx); ++n) { neighs.emplace_back(conn(idx, n)); @@ -146,101 +131,65 @@ inline std::vector getNeighs(const atlas::mesh::Nodes::Connectivity& conn, return neighs; } -inline std::vector const& cellNeighboursOfCell(atlas::Mesh const& m, int const& idx) { - // note this is only a workaround and does only work as long as we have only one mesh - static std::map> neighs; - if(neighs.count(idx) == 0) { - const auto& conn = m.cells().edge_connectivity(); - neighs[idx] = std::vector{}; - for(int n = 0; n < conn.cols(idx); ++n) { - int initialEdge = conn(idx, n); - for(int c1 = 0; c1 < m.cells().size(); ++c1) { - for(int n1 = 0; n1 < conn.cols(c1); ++n1) { - int compareEdge = conn(c1, n1); - if(initialEdge == compareEdge && c1 != idx) { - neighs[idx].emplace_back(c1); - } +std::vector const cellNeighboursOfCell(atlas::Mesh const& m, int const& idx) { + const auto& conn = m.cells().edge_connectivity(); + auto neighs = std::vector{}; + for(int n = 0; n < conn.cols(idx); ++n) { + int initialEdge = conn(idx, n); + for(int c1 = 0; c1 < m.cells().size(); ++c1) { + for(int n1 = 0; n1 < conn.cols(c1); ++n1) { + int compareEdge = conn(c1, n1); + if(initialEdge == compareEdge && c1 != idx) { + neighs.emplace_back(c1); } } } } - return neighs[idx]; + return neighs; } -inline std::vector const& edgeNeighboursOfCell(atlas::Mesh const& m, int const& idx) { - // note this is only a workaround and does only work as long as we have only one mesh - static std::map> neighs; - if(neighs.count(idx) == 0) { - neighs[idx] = getNeighs(m.cells().edge_connectivity(), idx); - } - return neighs[idx]; +std::vector const edgeNeighboursOfCell(atlas::Mesh const& m, int const& idx) { + return getNeighs(m.cells().edge_connectivity(), idx); } -inline std::vector const& nodeNeighboursOfCell(atlas::Mesh const& m, int const& idx) { - // note this is only a workaround and does only work as long as we have only one mesh - static std::map> neighs; - if(neighs.count(idx) == 0) { - neighs[idx] = getNeighs(m.cells().node_connectivity(), idx); - } - return neighs[idx]; +std::vector const nodeNeighboursOfCell(atlas::Mesh const& m, int const& idx) { + return getNeighs(m.cells().node_connectivity(), idx); } -inline std::vector cellNeighboursOfEdge(atlas::Mesh const& m, int const& idx) { - // note this is only a workaround and does only work as long as we have only one mesh - static std::map> neighs; - if(neighs.count(idx) == 0) { - neighs[idx] = getNeighs(m.edges().cell_connectivity(), idx); - } - assert(neighs[idx].size() == 2); - return neighs[idx]; +std::vector const cellNeighboursOfEdge(atlas::Mesh const& m, int const& idx) { + auto neighs = getNeighs(m.edges().cell_connectivity(), idx); + assert(neighs.size() == 2); + return neighs; } -inline std::vector nodeNeighboursOfEdge(atlas::Mesh const& m, int const& idx) { - // note this is only a workaround and does only work as long as we have only one mesh - static std::map> neighs; - if(neighs.count(idx) == 0) { - neighs[idx] = getNeighs(m.edges().node_connectivity(), idx); - } - assert(neighs[idx].size() == 2); - return neighs[idx]; +std::vector const nodeNeighboursOfEdge(atlas::Mesh const& m, int const& idx) { + auto neighs = getNeighs(m.edges().node_connectivity(), idx); + assert(neighs.size() == 2); + return neighs; } -inline std::vector cellNeighboursOfNode(atlas::Mesh const& m, int const& idx) { - // note this is only a workaround and does only work as long as we have only one mesh - static std::map> neighs; - if(neighs.count(idx) == 0) { - neighs[idx] = getNeighs(m.nodes().cell_connectivity(), idx); - } - return neighs[idx]; +std::vector const cellNeighboursOfNode(atlas::Mesh const& m, int const& idx) { + return getNeighs(m.nodes().cell_connectivity(), idx); } -inline std::vector edgeNeighboursOfNode(atlas::Mesh const& m, int const& idx) { - // note this is only a workaround and does only work as long as we have only one mesh - static std::map> neighs; - if(neighs.count(idx) == 0) { - neighs[idx] = getNeighs(m.nodes().edge_connectivity(), idx); - } - return neighs[idx]; +std::vector const edgeNeighboursOfNode(atlas::Mesh const& m, int const& idx) { + return getNeighs(m.cells().edge_connectivity(), idx); } -inline std::vector nodeNeighboursOfNode(atlas::Mesh const& m, int const& idx) { - // note this is only a workaround and does only work as long as we have only one mesh - static std::map> neighs; - if(neighs.count(idx) == 0) { - const auto& conn_nodes_to_edge = m.nodes().edge_connectivity(); - neighs[idx] = std::vector{}; - for(int ne = 0; ne < conn_nodes_to_edge.cols(idx); ++ne) { - int nbh_edge_idx = conn_nodes_to_edge(idx, ne); - const auto& conn_edge_to_nodes = m.edges().node_connectivity(); - for(int nn = 0; nn < conn_edge_to_nodes.cols(nbh_edge_idx); ++nn) { - int nbhNode = conn_edge_to_nodes(idx, nn); - if(nbhNode != idx) { - neighs[idx].emplace_back(nbhNode); - } +std::vector const nodeNeighboursOfNode(atlas::Mesh const& m, int const& idx) { + const auto& conn_nodes_to_edge = m.nodes().edge_connectivity(); + auto neighs = std::vector{}; + for(int ne = 0; ne < conn_nodes_to_edge.cols(idx); ++ne) { + int nbh_edge_idx = conn_nodes_to_edge(idx, ne); + const auto& conn_edge_to_nodes = m.edges().node_connectivity(); + for(int nn = 0; nn < conn_edge_to_nodes.cols(nbh_edge_idx); ++nn) { + int nbhNode = conn_edge_to_nodes(idx, nn); + if(nbhNode != idx) { + neighs.emplace_back(nbhNode); } } } - return neighs[idx]; + return neighs; } //===------------------------------------------------------------------------------------------===// diff --git a/stencils/interfaces/mylib_interface.hpp b/stencils/interfaces/mylib_interface.hpp index dcbb5e1..d1c15fa 100644 --- a/stencils/interfaces/mylib_interface.hpp +++ b/stencils/interfaces/mylib_interface.hpp @@ -1,8 +1,8 @@ #ifndef DAWN_PROTOTYPE_MYLIB_INTERFACE_H_ #define DAWN_PROTOTYPE_MYLIB_INTERFACE_H_ +#include "driver-includes/unstructured_interface.hpp" #include "mylib.hpp" -#include "unstructured_interface.hpp" #include namespace mylibInterface { From 4b776a945829b0ffc4a065de73ebf11562ba7c2b Mon Sep 17 00:00:00 2001 From: Giacomo Serafini Date: Tue, 17 Mar 2020 17:17:27 +0100 Subject: [PATCH 2/2] Revert atlas_interface.hpp --- stencils/interfaces/atlas_interface.hpp | 149 ++++++++++++++++-------- 1 file changed, 100 insertions(+), 49 deletions(-) diff --git a/stencils/interfaces/atlas_interface.hpp b/stencils/interfaces/atlas_interface.hpp index e1d4bd8..d592eab 100644 --- a/stencils/interfaces/atlas_interface.hpp +++ b/stencils/interfaces/atlas_interface.hpp @@ -68,14 +68,17 @@ struct atlasTag {}; template class Field { +private: + atlas::array::ArrayView atlas_field_; + public: T const& operator()(int f, int k) const { return atlas_field_(f, k); } T& operator()(int f, int k) { return atlas_field_(f, k); } Field(atlas::array::ArrayView const& atlas_field) : atlas_field_(atlas_field) {} - -private: - atlas::array::ArrayView atlas_field_; + T* data() { return atlas_field_.data(); } + const T* data() const { return atlas_field_.data(); } + int numElements() const { return atlas_field_.shape(0) * atlas_field_.shape(1); } }; template @@ -87,6 +90,9 @@ Field vertexFieldType(atlasTag); template class SparseDimension { +private: + atlas::array::ArrayView sparse_dimension_; + public: T const& operator()(int elem_idx, int sparse_dim_idx, int level) const { return sparse_dimension_(elem_idx, level, sparse_dim_idx); @@ -95,11 +101,14 @@ class SparseDimension { return sparse_dimension_(elem_idx, level, sparse_dim_idx); } + T* data() { return sparse_dimension_.data(); } + const T* data() const { return sparse_dimension_.data(); } + int numElements() const { + return sparse_dimension_.shape(0) * sparse_dimension_.shape(1) * sparse_dimension_.shape(2); + } + SparseDimension(atlas::array::ArrayView const& sparse_dimension) : sparse_dimension_(sparse_dimension) {} - -private: - atlas::array::ArrayView sparse_dimension_; }; template @@ -111,11 +120,17 @@ SparseDimension sparseVertexFieldType(atlasTag); atlas::Mesh meshType(atlasTag); -auto getCells(atlasTag, atlas::Mesh const& m) { return utility::irange(0, m.cells().size()); } -auto getEdges(atlasTag, atlas::Mesh const& m) { return utility::irange(0, m.edges().size()); } -auto getVertices(atlasTag, atlas::Mesh const& m) { return utility::irange(0, m.nodes().size()); } +inline auto getCells(atlasTag, atlas::Mesh const& m) { + return utility::irange(0, m.cells().size()); +} +inline auto getEdges(atlasTag, atlas::Mesh const& m) { + return utility::irange(0, m.edges().size()); +} +inline auto getVertices(atlasTag, atlas::Mesh const& m) { + return utility::irange(0, m.nodes().size()); +} -std::vector getNeighs(const atlas::Mesh::HybridElements::Connectivity& conn, int idx) { +inline std::vector getNeighs(const atlas::Mesh::HybridElements::Connectivity& conn, int idx) { std::vector neighs; for(int n = 0; n < conn.cols(idx); ++n) { neighs.emplace_back(conn(idx, n)); @@ -123,7 +138,7 @@ std::vector getNeighs(const atlas::Mesh::HybridElements::Connectivity& conn return neighs; } -std::vector getNeighs(const atlas::mesh::Nodes::Connectivity& conn, int idx) { +inline std::vector getNeighs(const atlas::mesh::Nodes::Connectivity& conn, int idx) { std::vector neighs; for(int n = 0; n < conn.cols(idx); ++n) { neighs.emplace_back(conn(idx, n)); @@ -131,65 +146,101 @@ std::vector getNeighs(const atlas::mesh::Nodes::Connectivity& conn, int idx return neighs; } -std::vector const cellNeighboursOfCell(atlas::Mesh const& m, int const& idx) { - const auto& conn = m.cells().edge_connectivity(); - auto neighs = std::vector{}; - for(int n = 0; n < conn.cols(idx); ++n) { - int initialEdge = conn(idx, n); - for(int c1 = 0; c1 < m.cells().size(); ++c1) { - for(int n1 = 0; n1 < conn.cols(c1); ++n1) { - int compareEdge = conn(c1, n1); - if(initialEdge == compareEdge && c1 != idx) { - neighs.emplace_back(c1); +inline std::vector const& cellNeighboursOfCell(atlas::Mesh const& m, int const& idx) { + // note this is only a workaround and does only work as long as we have only one mesh + static std::map> neighs; + if(neighs.count(idx) == 0) { + const auto& conn = m.cells().edge_connectivity(); + neighs[idx] = std::vector{}; + for(int n = 0; n < conn.cols(idx); ++n) { + int initialEdge = conn(idx, n); + for(int c1 = 0; c1 < m.cells().size(); ++c1) { + for(int n1 = 0; n1 < conn.cols(c1); ++n1) { + int compareEdge = conn(c1, n1); + if(initialEdge == compareEdge && c1 != idx) { + neighs[idx].emplace_back(c1); + } } } } } - return neighs; + return neighs[idx]; } -std::vector const edgeNeighboursOfCell(atlas::Mesh const& m, int const& idx) { - return getNeighs(m.cells().edge_connectivity(), idx); +inline std::vector const& edgeNeighboursOfCell(atlas::Mesh const& m, int const& idx) { + // note this is only a workaround and does only work as long as we have only one mesh + static std::map> neighs; + if(neighs.count(idx) == 0) { + neighs[idx] = getNeighs(m.cells().edge_connectivity(), idx); + } + return neighs[idx]; } -std::vector const nodeNeighboursOfCell(atlas::Mesh const& m, int const& idx) { - return getNeighs(m.cells().node_connectivity(), idx); +inline std::vector const& nodeNeighboursOfCell(atlas::Mesh const& m, int const& idx) { + // note this is only a workaround and does only work as long as we have only one mesh + static std::map> neighs; + if(neighs.count(idx) == 0) { + neighs[idx] = getNeighs(m.cells().node_connectivity(), idx); + } + return neighs[idx]; } -std::vector const cellNeighboursOfEdge(atlas::Mesh const& m, int const& idx) { - auto neighs = getNeighs(m.edges().cell_connectivity(), idx); - assert(neighs.size() == 2); - return neighs; +inline std::vector cellNeighboursOfEdge(atlas::Mesh const& m, int const& idx) { + // note this is only a workaround and does only work as long as we have only one mesh + static std::map> neighs; + if(neighs.count(idx) == 0) { + neighs[idx] = getNeighs(m.edges().cell_connectivity(), idx); + } + assert(neighs[idx].size() == 2); + return neighs[idx]; } -std::vector const nodeNeighboursOfEdge(atlas::Mesh const& m, int const& idx) { - auto neighs = getNeighs(m.edges().node_connectivity(), idx); - assert(neighs.size() == 2); - return neighs; +inline std::vector nodeNeighboursOfEdge(atlas::Mesh const& m, int const& idx) { + // note this is only a workaround and does only work as long as we have only one mesh + static std::map> neighs; + if(neighs.count(idx) == 0) { + neighs[idx] = getNeighs(m.edges().node_connectivity(), idx); + } + assert(neighs[idx].size() == 2); + return neighs[idx]; } -std::vector const cellNeighboursOfNode(atlas::Mesh const& m, int const& idx) { - return getNeighs(m.nodes().cell_connectivity(), idx); +inline std::vector cellNeighboursOfNode(atlas::Mesh const& m, int const& idx) { + // note this is only a workaround and does only work as long as we have only one mesh + static std::map> neighs; + if(neighs.count(idx) == 0) { + neighs[idx] = getNeighs(m.nodes().cell_connectivity(), idx); + } + return neighs[idx]; } -std::vector const edgeNeighboursOfNode(atlas::Mesh const& m, int const& idx) { - return getNeighs(m.cells().edge_connectivity(), idx); +inline std::vector edgeNeighboursOfNode(atlas::Mesh const& m, int const& idx) { + // note this is only a workaround and does only work as long as we have only one mesh + static std::map> neighs; + if(neighs.count(idx) == 0) { + neighs[idx] = getNeighs(m.nodes().edge_connectivity(), idx); + } + return neighs[idx]; } -std::vector const nodeNeighboursOfNode(atlas::Mesh const& m, int const& idx) { - const auto& conn_nodes_to_edge = m.nodes().edge_connectivity(); - auto neighs = std::vector{}; - for(int ne = 0; ne < conn_nodes_to_edge.cols(idx); ++ne) { - int nbh_edge_idx = conn_nodes_to_edge(idx, ne); - const auto& conn_edge_to_nodes = m.edges().node_connectivity(); - for(int nn = 0; nn < conn_edge_to_nodes.cols(nbh_edge_idx); ++nn) { - int nbhNode = conn_edge_to_nodes(idx, nn); - if(nbhNode != idx) { - neighs.emplace_back(nbhNode); +inline std::vector nodeNeighboursOfNode(atlas::Mesh const& m, int const& idx) { + // note this is only a workaround and does only work as long as we have only one mesh + static std::map> neighs; + if(neighs.count(idx) == 0) { + const auto& conn_nodes_to_edge = m.nodes().edge_connectivity(); + neighs[idx] = std::vector{}; + for(int ne = 0; ne < conn_nodes_to_edge.cols(idx); ++ne) { + int nbh_edge_idx = conn_nodes_to_edge(idx, ne); + const auto& conn_edge_to_nodes = m.edges().node_connectivity(); + for(int nn = 0; nn < conn_edge_to_nodes.cols(nbh_edge_idx); ++nn) { + int nbhNode = conn_edge_to_nodes(idx, nn); + if(nbhNode != idx) { + neighs[idx].emplace_back(nbhNode); + } } } } - return neighs; + return neighs[idx]; } //===------------------------------------------------------------------------------------------===//