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/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 {