From 8318f2c4cf992a2d5c21485070548b26514dc0b9 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Sun, 20 Oct 2024 22:39:01 +0200 Subject: [PATCH 01/44] [WIP] Transfer matrix operation --- cpp/dolfinx/CMakeLists.txt | 1 + cpp/dolfinx/multigrid/CMakeLists.txt | 4 + cpp/dolfinx/multigrid/transfer_matrix.h | 259 ++++++++++++++ cpp/test/CMakeLists.txt | 1 + cpp/test/multigrid/transfer_matrix.cpp | 386 +++++++++++++++++++++ python/CMakeLists.txt | 1 + python/dolfinx/transfer.py | 29 ++ python/dolfinx/wrappers/dolfinx.cpp | 6 + python/dolfinx/wrappers/multigrid.cpp | 67 ++++ python/test/unit/transfer/test_transfer.py | 63 ++++ 10 files changed, 817 insertions(+) create mode 100644 cpp/dolfinx/multigrid/CMakeLists.txt create mode 100644 cpp/dolfinx/multigrid/transfer_matrix.h create mode 100644 cpp/test/multigrid/transfer_matrix.cpp create mode 100644 python/dolfinx/transfer.py create mode 100644 python/dolfinx/wrappers/multigrid.cpp create mode 100644 python/test/unit/transfer/test_transfer.py diff --git a/cpp/dolfinx/CMakeLists.txt b/cpp/dolfinx/CMakeLists.txt index 7adaaeb73cf..001cb2116be 100644 --- a/cpp/dolfinx/CMakeLists.txt +++ b/cpp/dolfinx/CMakeLists.txt @@ -17,6 +17,7 @@ set(DOLFINX_DIRS io la mesh + multigrid nls refinement ) diff --git a/cpp/dolfinx/multigrid/CMakeLists.txt b/cpp/dolfinx/multigrid/CMakeLists.txt new file mode 100644 index 00000000000..36d60438135 --- /dev/null +++ b/cpp/dolfinx/multigrid/CMakeLists.txt @@ -0,0 +1,4 @@ +set(HEADERS_multigrid + ${CMAKE_CURRENT_SOURCE_DIR}/transfer_matrix.h + PARENT_SCOPE +) diff --git a/cpp/dolfinx/multigrid/transfer_matrix.h b/cpp/dolfinx/multigrid/transfer_matrix.h new file mode 100644 index 00000000000..bb68fc49fa3 --- /dev/null +++ b/cpp/dolfinx/multigrid/transfer_matrix.h @@ -0,0 +1,259 @@ +// Copyright (C) 2024 Paul T. Kühner +// +// This file is part of DOLFINX (https://www.fenicsproject.org) +// +// SPDX-License-Identifier: LGPL-3.0-or-later + +#pragma once + +#include +#include +#include +#include +#include +#include + +#include + +#include "dolfinx/common/IndexMap.h" +#include "dolfinx/fem/FunctionSpace.h" +#include "dolfinx/la/SparsityPattern.h" +#include "dolfinx/la/utils.h" +#include "dolfinx/mesh/Mesh.h" + +namespace dolfinx::multigrid +{ + +template +std::vector +inclusion_mapping(const dolfinx::mesh::Mesh& mesh_from, + const dolfinx::mesh::Mesh& mesh_to) +{ + + const common::IndexMap& im_from = *mesh_from.topology()->index_map(0); + const common::IndexMap& im_to = *mesh_to.topology()->index_map(0); + + std::vector map(im_from.size_global(), -1); + + std::span x_from = mesh_from.geometry().x(); + std::span x_to = mesh_to.geometry().x(); + + auto build_global_to_local = [&](const auto& im) + { + return [&](std::int32_t idx) + { + std::array tmp; + im.local_to_global(std::vector{idx}, tmp); + return tmp[0]; + }; + }; + + auto to_global_to = build_global_to_local(im_to); + auto to_global_from = build_global_to_local(im_from); + + for (std::int32_t i = 0; i < im_from.size_local(); i++) + { + std::ranges::subrange vertex_from(std::next(x_from.begin(), 3 * i), + std::next(x_from.begin(), 3 * (i + 1))); + for (std::int64_t j = 0; j < im_to.size_local() + im_to.num_ghosts(); j++) + { + std::ranges::subrange vertex_to(std::next(x_to.begin(), 3 * j), + std::next(x_to.begin(), 3 * (j + 1))); + + if (std::ranges::equal( + vertex_from, vertex_to, [](T a, T b) + { return std::abs(a - b) <= std::numeric_limits::epsilon(); })) + { + map[to_global_from(i)] = to_global_to(j); + break; + } + } + } + + if (dolfinx::MPI::size(mesh_to.comm()) == 1) + { + // no communication required + assert(std::ranges::all_of(map, [](auto e) { return e >= 0; })); + return map; + } + + // map holds at this point for every original local index the corresponding + // mapped global index. All other entries are still -1, but are available on + // other processes. + std::vector result(map.size(), -1); + MPI_Allreduce(map.data(), result.data(), map.size(), + dolfinx::MPI::mpi_t, MPI_MAX, mesh_from.comm()); + + assert(std::ranges::all_of(result, [](auto e) { return e >= 0; })); + return result; +} + +template +la::SparsityPattern +create_sparsity_pattern(const dolfinx::fem::FunctionSpace& V_from, + const dolfinx::fem::FunctionSpace& V_to, + const std::vector& inclusion_map) +{ + if (*V_from.element() != *V_to.element()) + throw std::runtime_error( + "Transfer between different element types not supported"); + + if (V_from.element()->basix_element().family() != basix::element::family::P) + throw std::runtime_error("Only Lagrange elements supported"); + + if (V_from.element()->basix_element().degree() != 1) + throw std::runtime_error("Only piecewise linear elements supported"); + + // TODO: mixed elements? value shapes? DG? + + auto mesh_from = V_from.mesh(); + auto mesh_to = V_to.mesh(); + assert(mesh_from); + assert(mesh_to); + + MPI_Comm comm = mesh_from->comm(); + { + // Check comms equal + int result; + MPI_Comm_compare(comm, mesh_to->comm(), &result); + assert(result == MPI_CONGRUENT); + } + assert(mesh_from->topology()->dim() == mesh_to->topology()->dim()); + + auto to_v_to_f = mesh_to->topology()->connectivity(0, 1); + auto to_f_to_v = mesh_to->topology()->connectivity(1, 0); + assert(to_v_to_f); + assert(to_f_to_v); + + auto dofmap_from = V_from.dofmap(); + auto dofmap_to = V_to.dofmap(); + assert(dofmap_from); + assert(dofmap_to); + + assert(mesh_to->topology()->index_map(0)); + assert(mesh_from->topology()->index_map(0)); + const common::IndexMap& im_to = *mesh_to->topology()->index_map(0); + const common::IndexMap& im_from = *mesh_from->topology()->index_map(0); + + dolfinx::la::SparsityPattern sp( + comm, {dofmap_from->index_map, dofmap_to->index_map}, + {dofmap_from->index_map_bs(), dofmap_to->index_map_bs()}); + + assert(inclusion_map.size() == im_from.size_global()); + for (int dof_from_global = 0; dof_from_global < im_from.size_global(); + dof_from_global++) + { + std::int64_t dof_to_global = inclusion_map[dof_from_global]; + + std::vector local_dof_to_v{0}; + im_to.global_to_local(std::vector{dof_to_global}, + local_dof_to_v); + + auto local_dof_to = local_dof_to_v[0]; + + bool is_remote = (local_dof_to == -1); + bool is_ghost = local_dof_to >= im_to.size_local(); + if (is_remote || is_ghost) + continue; + + std::vector dof_from_v{0}; + im_from.global_to_local(std::vector{dof_from_global}, + dof_from_v); + + std::ranges::for_each( + to_v_to_f->links(local_dof_to), + [&](auto e) + { + sp.insert( + std::vector(to_f_to_v->links(e).size(), dof_from_v[0]), + to_f_to_v->links(e)); + }); + } + sp.finalize(); + return sp; +} + +template +void assemble_transfer_matrix(la::MatSet auto mat_set, + const dolfinx::fem::FunctionSpace& V_from, + const dolfinx::fem::FunctionSpace& V_to, + const std::vector& inclusion_map, + std::function weight) +{ + if (*V_from.element() != *V_to.element()) + throw std::runtime_error( + "Transfer between different element types not supported"); + + if (V_from.element()->basix_element().family() != basix::element::family::P) + throw std::runtime_error("Only Lagrange elements supported"); + + if (V_from.element()->basix_element().degree() != 1) + throw std::runtime_error("Only piecewise linear elements supported"); + + // TODO: mixed elements? value shapes? DG? + + auto mesh_from = V_from.mesh(); + auto mesh_to = V_to.mesh(); + assert(mesh_from); + + MPI_Comm comm = mesh_from->comm(); + { + // Check comms equal + int result; + MPI_Comm_compare(comm, mesh_to->comm(), &result); + assert(result == MPI_CONGRUENT); + } + assert(mesh_from->topology()->dim() == mesh_to->topology()->dim()); + + auto to_v_to_f = mesh_to->topology()->connectivity(0, 1); + auto to_f_to_v = mesh_to->topology()->connectivity(1, 0); + assert(to_v_to_f); + assert(to_f_to_v); + + auto dofmap_from = V_from.dofmap(); + auto dofmap_to = V_to.dofmap(); + assert(dofmap_from); + assert(dofmap_to); + + assert(mesh_to->topology()->index_map(0)); + assert(mesh_from->topology()->index_map(0)); + const common::IndexMap& im_to = *mesh_to->topology()->index_map(0); + const common::IndexMap& im_from = *mesh_from->topology()->index_map(0); + + assert(inclusion_map.size() == im_from.size_global()); + + for (int dof_from_global = 0; dof_from_global < im_from.size_global(); + dof_from_global++) + { + std::int64_t dof_to_global = inclusion_map[dof_from_global]; + + std::vector local_dof_to_v{0}; + im_to.global_to_local(std::vector{dof_to_global}, + local_dof_to_v); + + auto local_dof_to = local_dof_to_v[0]; + + bool is_remote = (local_dof_to == -1); + bool is_ghost = local_dof_to >= im_to.size_local(); + if (is_remote || is_ghost) + continue; + + std::vector dof_from_v{0}; + im_from.global_to_local(std::vector{dof_from_global}, + dof_from_v); + + for (auto e : to_v_to_f->links(local_dof_to)) + { + for (auto n : to_f_to_v->links(e)) + { + // For now we only support distance <= 1 -> this should be somehow + // asserted. + std::int32_t distance = (n == local_dof_to) ? 0 : 1; + mat_set(std::vector{dof_from_v[0]}, std::vector{n}, + std::vector{weight(distance)}); + } + } + } +} + +} // namespace dolfinx::transfer \ No newline at end of file diff --git a/cpp/test/CMakeLists.txt b/cpp/test/CMakeLists.txt index ada545792a1..e956f3fda46 100644 --- a/cpp/test/CMakeLists.txt +++ b/cpp/test/CMakeLists.txt @@ -48,6 +48,7 @@ add_executable( mesh/refinement/interval.cpp mesh/refinement/option.cpp mesh/refinement/rectangle.cpp + multigrid/transfer_matrix.cpp ${CMAKE_CURRENT_BINARY_DIR}/poisson.c ) target_link_libraries(unittests PRIVATE Catch2::Catch2WithMain dolfinx) diff --git a/cpp/test/multigrid/transfer_matrix.cpp b/cpp/test/multigrid/transfer_matrix.cpp new file mode 100644 index 00000000000..60fa6dd3b04 --- /dev/null +++ b/cpp/test/multigrid/transfer_matrix.cpp @@ -0,0 +1,386 @@ +// Copyright (C) 2024 Paul T. Kühner +// +// This file is part of DOLFINX (https://www.fenicsproject.org) +// +// SPDX-License-Identifier: LGPL-3.0-or-later + +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace dolfinx; +using namespace Catch::Matchers; + +template +constexpr auto weight + = [](std::int32_t d) -> T { return d == 0 ? T(1) : T(.5); }; + +template +std::shared_ptr> +P1_element(basix::cell::type cell_type) +{ + auto element + = basix::create_element(basix::element::family::P, cell_type, 1, + basix::element::lagrange_variant::unset, + basix::element::dpc_variant::unset, false); + return std::make_shared>(element); +} + +TEMPLATE_TEST_CASE("Transfer Matrix 1D", "[transfer_matrix]", double, float) +{ + using T = TestType; + + auto mesh_coarse + = std::make_shared>(dolfinx::mesh::create_interval( + MPI_COMM_SELF, 2, {0.0, 1.0}, mesh::GhostMode::none)); + + auto [mesh_fine, parent_cell, parent_facet] + = refinement::refine(*mesh_coarse, std::nullopt); + + auto element = P1_element(basix::cell::type::interval); + auto V_coarse = std::make_shared>( + fem::create_functionspace(mesh_coarse, element)); + auto V_fine + = std::make_shared>(fem::create_functionspace( + std::make_shared>(mesh_fine), element)); + + mesh_fine.topology()->create_connectivity(1, 0); + mesh_fine.topology()->create_connectivity(0, 1); + + std::vector inclusion_map + = multigrid::inclusion_mapping(*mesh_coarse, mesh_fine); + + CHECK_THAT(inclusion_map, RangeEquals(std::vector{0, 2, 4})); + + la::SparsityPattern sp + = multigrid::create_sparsity_pattern(*V_coarse, *V_fine, inclusion_map); + + la::MatrixCSR transfer_matrix(std::move(sp), la::BlockMode::compact); + multigrid::assemble_transfer_matrix(transfer_matrix.mat_set_values(), + *V_coarse, *V_fine, inclusion_map, + weight); + + // clang-format off + std::vector expected{1.0, 0.5, 0.0, 0.0, 0.0, + 0.0, 0.5, 1.0, 0.5, 0.0, + 0.0, 0.0, 0.0, 0.5, 1.0}; + // clang-format on + CHECK_THAT(transfer_matrix.to_dense(), + RangeEquals(expected, + [](auto a, auto b) { + return std::abs(a - b) + <= std::numeric_limits::epsilon(); + })); +} + +// TEMPLATE_TEST_CASE("Transfer Matrix 1D (parallel)", "[transfer_matrix]", +// double) // TODO: float +// { +// using T = TestType; + +// int comm_size = dolfinx::MPI::size(MPI_COMM_WORLD); +// int rank = dolfinx::MPI::rank(MPI_COMM_WORLD); + +// // TODO: see https://github.com/FEniCS/dolfinx/issues/3358 +// // auto part = mesh::create_cell_partitioner(mesh::GhostMode::shared_facet); +// mesh::CellPartitionFunction part +// = [&](MPI_Comm /* comm */, int /* nparts */, +// const std::vector& /* cell_types */, +// const std::vector>& /* cells */) +// { +// std::vector> data; +// if (comm_size == 1) +// data = {{0}, {0}, {0}, {0}}; +// else if (comm_size == 2) +// data = {{0}, {0}, {0}, {0}, {0, 1}, {1, 0}, {1}, {1}, {1}}; +// else if (comm_size == 3) +// data = {{1}, {1}, {1}, {1}, {1, 2}, {2, 1}, {2}, +// {2}, {2, 0}, {0, 2}, {0}, {0}, {0}, {0}}; +// else +// FAIL("Test only supports <= 3 processes"); +// return graph::AdjacencyList(std::move(data)); +// }; + +// auto mesh_coarse = std::make_shared>( +// mesh::create_interval(MPI_COMM_WORLD, 5 * comm_size - 1, {0., 1.}, +// mesh::GhostMode::none, part)); + +// // TODO: see https://github.com/FEniCS/dolfinx/issues/3358 +// // auto part_refine +// // = mesh::create_cell_partitioner(mesh::GhostMode::shared_facet); +// mesh::CellPartitionFunction part_refine +// = [&](MPI_Comm comm, int /* nparts */, +// const std::vector& /* cell_types */, +// const std::vector>& /* cells */) +// { +// std::vector> data; +// if (dolfinx::MPI::size(comm) == 1) +// data = {{0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}}; +// else if (dolfinx::MPI::size(comm) == 2) +// { +// if (rank == 0) +// data = {{0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0, 1}, {1, 0}}; +// else +// data = {{1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}}; +// } +// else if (dolfinx::MPI::size(comm) == 3) +// { +// if (rank == 0) +// data = {{2, 0}, {0, 2}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}}; +// else if (rank == 1) +// data = {{1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1, 2}}; +// else +// data = {{2, 1}, {2}, {2}, {2}, {2}, {2}, {2}, {2}}; +// } +// else +// FAIL("Test only supports <= 3 processes"); +// return graph::AdjacencyList(std::move(data)); +// }; + +// auto [mesh_fine, parent_cell, parent_facet] = refinement::refine( +// *mesh_coarse, std::nullopt, part_refine, refinement::Option::none); + +// auto element = P1_element(basix::cell::type::interval); + +// auto V_coarse = std::make_shared>( +// fem::create_functionspace(mesh_coarse, element)); +// auto V_fine +// = std::make_shared>(fem::create_functionspace( +// std::make_shared>(mesh_fine), element)); + +// mesh_fine.topology()->create_connectivity(1, 0); +// mesh_fine.topology()->create_connectivity(0, 1); + +// std::vector inclusion_map; +// switch (comm_size) +// { +// case 1: +// inclusion_map = {0, 2, 4, 6, 8}; +// break; +// case 2: +// inclusion_map = {0, 2, 4, 6, 8, 10, 12, 14, 16, 18}; +// break; +// case 3: +// inclusion_map = {0, 2, 4, 6, 8, 9, 11, 13, 15, 17, 19, 27, 25, 23, 20}; +// break; +// } + +// CHECK_THAT(multigrid::inclusion_mapping(*mesh_coarse, mesh_fine), +// RangeEquals(inclusion_map)); + +// la::SparsityPattern sp +// = multigrid::create_sparsity_pattern(*V_coarse, *V_fine, inclusion_map); + +// la::MatrixCSR transfer_matrix(std::move(sp), la::BlockMode::compact); +// multigrid::assemble_transfer_matrix(transfer_matrix.mat_set_values(), +// *V_coarse, *V_fine, inclusion_map, +// weight); + +// std::array, 3> expected; +// if (comm_size == 1) +// { +// // clang-format off +// expected[0] = { +// /* row_0 */ 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* row_1 */ 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* row_2 */ 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, +// /* row_3 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, +// /* row_4 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0 +// }; +// // clang-format on +// } +// else if (comm_size == 2) +// { +// // clang-format off +// expected[0] = { +// /* row_0 */ 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* row_1 */ 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* row_2 */ 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* row_3 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* row_4 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* row_5 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* row_6 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 +// }; +// expected[1] = { +// /* row_0 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* row_1 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, +// /* row_2 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, +// /* row_3 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, +// /* row_4 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* row_5 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 +// }; +// // clang-format on +// } +// else if (comm_size == 3) +// { +// // clang-format off +// expected[0] = { +// /* row_0 */ 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* row_1 */ 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* row_2 */ 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* row_3 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* row_4 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* row_5 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* row_6 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// }; +// expected[1] = { +// /* row_0 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* row_1 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* row_2 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* row_3 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* row_4 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* row_5 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, +// /* row_6 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// }; +// expected[2] = { +// /* row_0 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, +// /* row_1 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, +// /* row_2 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, +// /* row_3 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* row_4 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* row_5 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* row_6 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 +// }; +// // clang-format on +// } +// else +// { +// CHECK(false); +// } + +// CHECK_THAT(transfer_matrix.to_dense(), +// RangeEquals(expected[rank], +// [](auto a, auto b) { +// return std::abs(a - b) +// <= std::numeric_limits::epsilon(); +// })); +// } + +TEMPLATE_TEST_CASE("Transfer Matrix 2D", "[transfer_matrix]", double) +{ + using T = TestType; + + auto mesh_coarse + = std::make_shared>(dolfinx::mesh::create_rectangle( + MPI_COMM_SELF, {{{0, 0}, {1, 1}}}, {1, 1}, mesh::CellType::triangle)); + mesh_coarse->topology()->create_entities(1); + + auto [mesh_fine, parent_cell, parent_facet] + = refinement::refine(*mesh_coarse, std::nullopt); + + auto element = P1_element(basix::cell::type::triangle); + + auto V_coarse = std::make_shared>( + fem::create_functionspace(mesh_coarse, element)); + auto V_fine + = std::make_shared>(fem::create_functionspace( + std::make_shared>(mesh_fine), element)); + + mesh_fine.topology()->create_connectivity(1, 0); + mesh_fine.topology()->create_connectivity(0, 1); + + std::vector inclusion_map + = multigrid::inclusion_mapping(*mesh_coarse, mesh_fine); + CHECK_THAT(inclusion_map, RangeEquals(std::vector{4, 1, 5, 8})); + + la::SparsityPattern sp + = multigrid::create_sparsity_pattern(*V_coarse, *V_fine, inclusion_map); + la::MatrixCSR transfer_matrix(std::move(sp), la::BlockMode::compact); + multigrid::assemble_transfer_matrix(transfer_matrix.mat_set_values(), + *V_coarse, *V_fine, inclusion_map, + weight); + transfer_matrix.scatter_rev(); + std::vector expected{0.5, 0.0, 0.5, 0.0, 1.0, 0.0, 0.5, 0.0, 0.0, + 0.5, 1.0, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.5, 0.0, 0.0, 0.5, 0.0, 1.0, 0.0, 0.5, 0.0, + 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.5, 1.0}; + CHECK_THAT(transfer_matrix.to_dense(), + RangeEquals(expected, + [](auto a, auto b) { + return std::abs(a - b) + <= std::numeric_limits::epsilon(); + })); +} + +TEMPLATE_TEST_CASE("Transfer Matrix 3D", "[transfer_matrix]", double) +{ + using T = TestType; + + auto mesh_coarse = std::make_shared>( + dolfinx::mesh::create_box(MPI_COMM_SELF, {{{0, 0, 0}, {1, 1, 1}}}, + {1, 1, 1}, mesh::CellType::tetrahedron)); + mesh_coarse->topology()->create_entities(1); + + auto [mesh_fine, parent_cell, parent_facet] + = refinement::refine(*mesh_coarse, std::nullopt); + + auto element = P1_element(basix::cell::type::tetrahedron); + + auto V_coarse = std::make_shared>( + fem::create_functionspace(mesh_coarse, element)); + auto V_fine + = std::make_shared>(fem::create_functionspace( + std::make_shared>(mesh_fine), element)); + + mesh_fine.topology()->create_connectivity(1, 0); + mesh_fine.topology()->create_connectivity(0, 1); + + std::vector inclusion_map + = multigrid::inclusion_mapping(*mesh_coarse, mesh_fine); + CHECK_THAT(inclusion_map, RangeEquals(std::vector{ + 0, 6, 15, 25, 17, 9, 11, 22})); + + la::SparsityPattern sp + = multigrid::create_sparsity_pattern(*V_coarse, *V_fine, inclusion_map); + la::MatrixCSR transfer_matrix(std::move(sp), la::BlockMode::compact); + multigrid::assemble_transfer_matrix(transfer_matrix.mat_set_values(), + *V_coarse, *V_fine, inclusion_map, + weight); + + std::vector expected{ + 1.0, 0.5, 0.5, 0.5, 0.5, 0.5, 0.0, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.5, + 0.5, 0.5, 0.0, 1.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.5, 0.0, 0.5, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.5, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.5, 0.5, 1.0, 0.0, 0.0, 0.0, 0.5, 0.0, + 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.5, 0.0, 0.5, + 0.5, 1.0, 0.5, 0.0, 0.5, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.5, 0.0, 0.5, 0.0, 0.5, 1.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.5, 0.0, 0.0, + 0.0, 0.5, 0.5, 0.0, 0.0, 0.5, 0.0, 0.0, 0.5, 1.0, 0.0, 0.0, 0.0, 0.5, 0.0, + 0.0, 0.0, 0.0, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, + 0.0, 0.5, 0.0, 0.0, 0.5, 0.5, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, + 0.0, 0.0, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.5, 0.5, + 0.5, 1.0, 0.0, 0.0, 0.0, 0.5}; + CHECK_THAT(transfer_matrix.to_dense(), + RangeEquals(expected, + [](auto a, auto b) { + return std::abs(a - b) + <= std::numeric_limits::epsilon(); + })); +} diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt index 0f440260214..5425d10896f 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -44,6 +44,7 @@ nanobind_add_module( dolfinx/wrappers/la.cpp dolfinx/wrappers/log.cpp dolfinx/wrappers/mesh.cpp + dolfinx/wrappers/multigrid.cpp dolfinx/wrappers/petsc.cpp dolfinx/wrappers/refinement.cpp ) diff --git a/python/dolfinx/transfer.py b/python/dolfinx/transfer.py new file mode 100644 index 00000000000..7b9e8b59e7c --- /dev/null +++ b/python/dolfinx/transfer.py @@ -0,0 +1,29 @@ +# Copyright (C) 2024 Paul T. Kühner +# +# This file is part of DOLFINx (https://www.fenicsproject.org) +# +# SPDX-License-Identifier: LGPL-3.0-or-later + +import numpy as np +from numpy.typing import NDArray + +from dolfinx.cpp.la import SparsityPattern +from dolfinx.cpp.transfer import ( + assemble_transfer_matrix, +) +from dolfinx.cpp.transfer import create_sparsity_pattern as _create_sparsity_pattern +from dolfinx.cpp.transfer import inclusion_mapping as _inclusion_mapping +from dolfinx.fem import FunctionSpace +from dolfinx.mesh import Mesh + +__all__ = ["assemble_transfer_matrix", "create_sparsity_pattern", "inclusion_mapping"] + + +def inclusion_mapping(mesh_from: Mesh, mesh_to: Mesh) -> NDArray[np.int64]: + return _inclusion_mapping(mesh_from._cpp_object, mesh_to._cpp_object) + + +def create_sparsity_pattern( + V_from: FunctionSpace, V_to: FunctionSpace, inclusion_map: NDArray[np.int64] +) -> SparsityPattern: + return _create_sparsity_pattern(V_from._cpp_object, V_to._cpp_object, inclusion_map) diff --git a/python/dolfinx/wrappers/dolfinx.cpp b/python/dolfinx/wrappers/dolfinx.cpp index 9c0864f9ef1..9d2ca3e716d 100644 --- a/python/dolfinx/wrappers/dolfinx.cpp +++ b/python/dolfinx/wrappers/dolfinx.cpp @@ -25,6 +25,8 @@ void la(nb::module_& m); void mesh(nb::module_& m); void nls(nb::module_& m); void refinement(nb::module_& m); +void transfer(nb::module_& m); + } // namespace dolfinx_wrappers NB_MODULE(cpp, m) @@ -74,6 +76,10 @@ NB_MODULE(cpp, m) nb::module_ refinement = m.def_submodule("refinement", "Refinement module"); dolfinx_wrappers::refinement(refinement); + // Create transfer submodule + nb::module_ transfer = m.def_submodule("transfer", "Transfer module"); + dolfinx_wrappers::transfer(transfer); + #if defined(HAS_PETSC) && defined(HAS_PETSC4PY) // PETSc-specific wrappers nb::module_ nls = m.def_submodule("nls", "Nonlinear solver module"); diff --git a/python/dolfinx/wrappers/multigrid.cpp b/python/dolfinx/wrappers/multigrid.cpp new file mode 100644 index 00000000000..cabc8e53726 --- /dev/null +++ b/python/dolfinx/wrappers/multigrid.cpp @@ -0,0 +1,67 @@ +// Copyright (C) 2024 Paul T. Kühner +// +// This file is part of DOLFINX (https://www.fenicsproject.org) +// +// SPDX-License-Identifier: LGPL-3.0-or-later + +#include +#include + +#include +#include + +#include +#include + +#include "array.h" + +namespace nb = nanobind; + +namespace dolfinx_wrappers +{ + +void transfer(nb::module_& m) +{ + m.def( + "inclusion_mapping", + [](const dolfinx::mesh::Mesh& mesh_from, + const dolfinx::mesh::Mesh& mesh_to) + { + std::vector map + = dolfinx::multigrid::inclusion_mapping(mesh_from, mesh_to); + return dolfinx_wrappers::as_nbarray(std::move(map)); + }, + nb::arg("mesh_from"), nb::arg("mesh_to"), + "Computes inclusion mapping between two meshes"); + + m.def( + "create_sparsity_pattern", + [](const dolfinx::fem::FunctionSpace& V_from, + const dolfinx::fem::FunctionSpace& V_to, + nb::ndarray& inclusion_map) + { + // TODO: cahnge to accepting range; + auto vec = std::vector(inclusion_map.data(), + inclusion_map.data() + inclusion_map.size()); + return dolfinx::multigrid::create_sparsity_pattern(V_from, V_to, + vec); + }, + nb::arg("V_from"), nb::arg("V_to"), nb::arg("inclusion_map")); + + m.def( + "assemble_transfer_matrix", + [](dolfinx::la::MatrixCSR& A, + const dolfinx::fem::FunctionSpace& V_from, + const dolfinx::fem::FunctionSpace& V_to, + const std::vector& inclusion_map, + std::function weight) + { + dolfinx::multigrid::assemble_transfer_matrix( + A.mat_set_values(), V_from, V_to, inclusion_map, weight); + }, + nb::arg("A"), nb::arg("V_from"), nb::arg("V_to"), + nb::arg("inclusion_map"), nb::arg("weight"), + "Assembles a transfer matrix between two function spaces."); +} + +} // namespace dolfinx_wrappers diff --git a/python/test/unit/transfer/test_transfer.py b/python/test/unit/transfer/test_transfer.py new file mode 100644 index 00000000000..8d60dfb12e2 --- /dev/null +++ b/python/test/unit/transfer/test_transfer.py @@ -0,0 +1,63 @@ +# Copyright (C) 2024 Paul T. Kühner +# +# This file is part of DOLFINx (https://www.fenicsproject.org) +# +# SPDX-License-Identifier: LGPL-3.0-or-later + +from mpi4py import MPI + +import pytest + +from dolfinx.fem import functionspace +from dolfinx.mesh import GhostMode, create_interval, refine +from dolfinx.transfer import create_sparsity_pattern, inclusion_mapping + + +@pytest.mark.parametrize( + "ghost_mode", [GhostMode.none, GhostMode.shared_vertex, GhostMode.shared_facet] +) +def test_1d(ghost_mode): + mesh = create_interval(MPI.COMM_WORLD, 10, (0, 1), ghost_mode=ghost_mode) + mesh_fine, _, _ = refine(mesh) + inclusion_map = inclusion_mapping(mesh, mesh_fine) + V = functionspace(mesh, ("Lagrange", 1, (1,))) + V_fine = functionspace(mesh_fine, ("Lagrange", 1, (1,))) + assert V.element == V_fine.element + V_fine.mesh.topology.create_connectivity(1, 0) + V_fine.mesh.topology.create_connectivity(0, 1) + create_sparsity_pattern(V, V_fine, inclusion_map) + # T = matrix_csr(sp) + # assemble_transfer_matrix( + # T._cpp_object, + # V._cpp_object, + # V_fine._cpp_object, + # inclusion_map, + # lambda i: 1.0 if i == 0 else 0.5, + # ) + # continue with assembly of matrix + + +@pytest.mark.parametrize("degree", [2, 3, 4]) +@pytest.mark.parametrize( + "ghost_mode", [GhostMode.none, GhostMode.shared_vertex, GhostMode.shared_facet] +) +def test_1d_higher_order(degree, ghost_mode): + mesh = create_interval(MPI.COMM_WORLD, 10, (0, 1), ghost_mode=ghost_mode) + mesh_fine, _, _ = refine(mesh) + + # this is a strictly geometric operation, and thus should pass + inclusion_map = inclusion_mapping(mesh, mesh_fine) + + V = functionspace(mesh, ("Lagrange", degree, (1,))) + V_fine = functionspace(mesh_fine, ("Lagrange", degree, (1,))) + + V_fine.mesh.topology.create_connectivity(1, 0) + V_fine.mesh.topology.create_connectivity(0, 1) + + # Check not supported throws + with pytest.raises(Exception): + create_sparsity_pattern(V, V_fine, inclusion_map) + + +if __name__ == "__main__": + pytest.main([__file__]) From 20a374c30efd57c9a68ab4d684fa0c44861e1c03 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Thu, 2 Jan 2025 15:41:15 +0100 Subject: [PATCH 02/44] Fix: interval refinement assert ignores optional type --- cpp/dolfinx/refinement/interval.h | 2 +- cpp/test/mesh/refinement/interval.cpp | 10 ++++++++++ 2 files changed, 11 insertions(+), 1 deletion(-) diff --git a/cpp/dolfinx/refinement/interval.h b/cpp/dolfinx/refinement/interval.h index 0b1505acd43..dc718470ad2 100644 --- a/cpp/dolfinx/refinement/interval.h +++ b/cpp/dolfinx/refinement/interval.h @@ -176,7 +176,7 @@ compute_refinement_data(const mesh::Mesh& mesh, } assert(cell_topology.size() == 2 * refined_cell_count); - assert(parent_cell->size() == (compute_parent_cell ? refined_cell_count : 0)); + assert((!compute_parent_cell) || (parent_cell->size() == refined_cell_count)); std::vector offsets(refined_cell_count + 1); std::ranges::generate(offsets, [i = 0]() mutable { return 2 * i++; }); diff --git a/cpp/test/mesh/refinement/interval.cpp b/cpp/test/mesh/refinement/interval.cpp index 1253fe00b04..d356aeb683e 100644 --- a/cpp/test/mesh/refinement/interval.cpp +++ b/cpp/test/mesh/refinement/interval.cpp @@ -239,3 +239,13 @@ TEMPLATE_TEST_CASE("Interval Refinement (parallel)", != v_to_e->links((center_index + 2) % 3)[0]); } } + +TEMPLATE_TEST_CASE("Interval uniform refinement", "[refinement][interva]", + double, float) +{ + auto interval + = dolfinx::mesh::create_interval(MPI_COMM_WORLD, 2, {0.0, 1.0}); + auto [refined, parent_edge, parent_facet] + = dolfinx::refinement::refine(interval, std::nullopt); + CHECK(refined.topology()->index_map(0)->size_global() == 5); +} From 2b72911b82c6666c8b375bf1c01d19f93d2d368e Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Sun, 20 Oct 2024 22:39:01 +0200 Subject: [PATCH 03/44] Introduce multigrid inclusion map computation --- cpp/dolfinx/CMakeLists.txt | 1 + cpp/dolfinx/multigrid/CMakeLists.txt | 4 + cpp/dolfinx/multigrid/inclusion.h | 152 ++++++++++++++++++++ cpp/test/CMakeLists.txt | 1 + cpp/test/multigrid/inclusion.cpp | 115 +++++++++++++++ python/CMakeLists.txt | 1 + python/dolfinx/transfer.py | 17 +++ python/dolfinx/wrappers/dolfinx.cpp | 6 + python/dolfinx/wrappers/multigrid.cpp | 38 +++++ python/test/unit/transfer/test_inclusion.py | 26 ++++ 10 files changed, 361 insertions(+) create mode 100644 cpp/dolfinx/multigrid/CMakeLists.txt create mode 100644 cpp/dolfinx/multigrid/inclusion.h create mode 100644 cpp/test/multigrid/inclusion.cpp create mode 100644 python/dolfinx/transfer.py create mode 100644 python/dolfinx/wrappers/multigrid.cpp create mode 100644 python/test/unit/transfer/test_inclusion.py diff --git a/cpp/dolfinx/CMakeLists.txt b/cpp/dolfinx/CMakeLists.txt index 7adaaeb73cf..001cb2116be 100644 --- a/cpp/dolfinx/CMakeLists.txt +++ b/cpp/dolfinx/CMakeLists.txt @@ -17,6 +17,7 @@ set(DOLFINX_DIRS io la mesh + multigrid nls refinement ) diff --git a/cpp/dolfinx/multigrid/CMakeLists.txt b/cpp/dolfinx/multigrid/CMakeLists.txt new file mode 100644 index 00000000000..fa5d718cb74 --- /dev/null +++ b/cpp/dolfinx/multigrid/CMakeLists.txt @@ -0,0 +1,4 @@ +set(HEADERS_multigrid + ${CMAKE_CURRENT_SOURCE_DIR}/inclusion.h + PARENT_SCOPE +) diff --git a/cpp/dolfinx/multigrid/inclusion.h b/cpp/dolfinx/multigrid/inclusion.h new file mode 100644 index 00000000000..0325ee26a74 --- /dev/null +++ b/cpp/dolfinx/multigrid/inclusion.h @@ -0,0 +1,152 @@ +// Copyright (C) 2024 Paul T. Kühner +// +// This file is part of DOLFINX (https://www.fenicsproject.org) +// +// SPDX-License-Identifier: LGPL-3.0-or-later + +#pragma once + +#include +#include +#include +#include +#include +#include + +#include + +#include "dolfinx/common/IndexMap.h" +#include "dolfinx/la/SparsityPattern.h" +#include "dolfinx/mesh/Mesh.h" + +namespace dolfinx::multigrid +{ + +template +std::vector +inclusion_mapping(const dolfinx::mesh::Mesh& mesh_from, + const dolfinx::mesh::Mesh& mesh_to) +{ + { + // Check comms equal + int result; + MPI_Comm_compare(mesh_from.comm(), mesh_to.comm(), &result); + assert(result == MPI_CONGRUENT); + } + + const common::IndexMap& im_from = *mesh_from.topology()->index_map(0); + const common::IndexMap& im_to = *mesh_to.topology()->index_map(0); + + std::vector map(im_from.size_global(), -1); + + std::span x_from = mesh_from.geometry().x(); + std::span x_to = mesh_to.geometry().x(); + + auto build_global_to_local = [&](const auto& im) + { + return [&](std::int32_t idx) + { + std::array tmp; + im.local_to_global(std::vector{idx}, tmp); + return tmp[0]; + }; + }; + + auto to_global_to = build_global_to_local(im_to); + auto to_global_from = build_global_to_local(im_from); + + for (std::int32_t i = 0; i < im_from.size_local(); i++) + { + std::ranges::subrange vertex_from(std::next(x_from.begin(), 3 * i), + std::next(x_from.begin(), 3 * (i + 1))); + for (std::int64_t j = 0; j < im_to.size_local() + im_to.num_ghosts(); j++) + { + std::ranges::subrange vertex_to(std::next(x_to.begin(), 3 * j), + std::next(x_to.begin(), 3 * (j + 1))); + + if (std::ranges::equal( + vertex_from, vertex_to, [](T a, T b) + { return std::abs(a - b) <= std::numeric_limits::epsilon(); })) + { + map[to_global_from(i)] = to_global_to(j); + break; + } + } + } + + if (dolfinx::MPI::size(mesh_to.comm()) == 1) + { + // no communication required + assert(std::ranges::all_of(map, [](auto e) { return e >= 0; })); + return map; + } + + // map holds at this point for every original local index the corresponding + // mapped global index (if it was available on the same process on the to + // mesh). + std::vector result(map.size(), -1); + MPI_Allreduce(map.data(), result.data(), map.size(), + dolfinx::MPI::mpi_t, MPI_MAX, mesh_from.comm()); + + if (std::ranges::all_of(result, [](auto e) { return e >= 0; })) + // All vertices indentified + return result; + + // Build global to vertex list + + // 1) exchange local sizes + std::vector local_sizes(dolfinx::MPI::size(mesh_from.comm())); + { + std::array tmp{im_to.size_local() * 3}; + MPI_Allgather(&tmp, 1, MPI_INT32_T, local_sizes.data(), 1, MPI_INT32_T, + mesh_from.comm()); + } + + // for (auto ls : local_sizes) + // std::cout << ls << ", "; + + // 2) compute displacement vector + std::vector displacements(local_sizes.size() + 1, 0); + std::partial_sum(local_sizes.begin(), local_sizes.end(), + displacements.begin() + 1); + + // for (auto ls : displacements) + // std::cout << ls << ", "; + + // 3) Allgather global x vector + std::vector global_x_to(im_to.size_global() * 3); + MPI_Allgatherv(mesh_to.geometry().x().data(), im_to.size_local() * 3, + dolfinx::MPI::mpi_t, global_x_to.data(), local_sizes.data(), + displacements.data(), dolfinx::MPI::mpi_t, + mesh_from.comm()); + + // Recheck indices on global data structure + for (std::int32_t i = 0; i < im_from.size_local(); i++) + { + std::ranges::subrange vertex_from(std::next(x_from.begin(), 3 * i), + std::next(x_from.begin(), 3 * (i + 1))); + for (std::int64_t j = 0; j < im_to.size_global(); j++) + { + std::ranges::subrange vertex_to( + std::next(global_x_to.begin(), 3 * j), + std::next(global_x_to.begin(), 3 * (j + 1))); + + if (std::ranges::equal( + vertex_from, vertex_to, [](T a, T b) + { return std::abs(a - b) <= std::numeric_limits::epsilon(); })) + { + map[to_global_from(i)] = j; + break; + } + } + } + + MPI_Allreduce(map.data(), result.data(), map.size(), + dolfinx::MPI::mpi_t, MPI_MAX, mesh_from.comm()); + + assert(std::ranges::all_of(result, [&](auto e) + { return e >= 0 && e < im_to.size_global(); })); + return result; +} + +} // namespace dolfinx::multigrid diff --git a/cpp/test/CMakeLists.txt b/cpp/test/CMakeLists.txt index ada545792a1..8d900d02dc9 100644 --- a/cpp/test/CMakeLists.txt +++ b/cpp/test/CMakeLists.txt @@ -48,6 +48,7 @@ add_executable( mesh/refinement/interval.cpp mesh/refinement/option.cpp mesh/refinement/rectangle.cpp + multigrid/inclusion.cpp ${CMAKE_CURRENT_BINARY_DIR}/poisson.c ) target_link_libraries(unittests PRIVATE Catch2::Catch2WithMain dolfinx) diff --git a/cpp/test/multigrid/inclusion.cpp b/cpp/test/multigrid/inclusion.cpp new file mode 100644 index 00000000000..7b0b533735a --- /dev/null +++ b/cpp/test/multigrid/inclusion.cpp @@ -0,0 +1,115 @@ +// Copyright (C) 2024 Paul T. Kühner +// +// This file is part of DOLFINX (https://www.fenicsproject.org) +// +// SPDX-License-Identifier: LGPL-3.0-or-later + +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include + +#include + +#include +#include +#include +#include + +using namespace dolfinx; +using namespace Catch::Matchers; + +template +void CHECK_inclusion_map(const dolfinx::mesh::Mesh& from, + const dolfinx::mesh::Mesh& to, + const std::vector& map) +{ + const common::IndexMap& im_from = *from.topology()->index_map(0); + const common::IndexMap& im_to = *to.topology()->index_map(0); + + // 1) exchange local sizes + std::vector local_sizes(dolfinx::MPI::size(from.comm())); + { + std::array tmp{im_to.size_local() * 3}; + MPI_Allgather(&tmp, 1, MPI_INT32_T, local_sizes.data(), 1, MPI_INT32_T, + from.comm()); + } + + // 2) compute displacement vector + std::vector displacements(local_sizes.size() + 1, 0); + std::partial_sum(local_sizes.begin(), local_sizes.end(), + displacements.begin() + 1); + + // 3) Allgather global x vector + std::vector global_x_to(im_to.size_global() * 3); + MPI_Allgatherv(to.geometry().x().data(), im_to.size_local() * 3, + dolfinx::MPI::mpi_t, global_x_to.data(), local_sizes.data(), + displacements.data(), dolfinx::MPI::mpi_t, from.comm()); + + REQUIRE(static_cast(map.size()) == im_from.size_global()); + for (std::int64_t i = 0; i < map.size(); i++) + { + std::array local{-1}; + im_from.global_to_local(std::array{i}, local); + if (local[0] == -1) + continue; + + CHECK(std::abs(from.geometry().x()[3 * local[0]] - global_x_to[3 * map[i]]) + < std::numeric_limits::epsilon()); + CHECK(std::abs(from.geometry().x()[3 * local[0] + 1] + - global_x_to[3 * map[i] + 1]) + < std::numeric_limits::epsilon()); + CHECK(std::abs(from.geometry().x()[3 * local[0] + 2] + - global_x_to[3 * map[i] + 2]) + < std::numeric_limits::epsilon()); + } +} + +/// Performs one uniform refinement and checks the inclusion map between coarse +/// and fine mesh against the (provided) list. +template +void TEST_inclusion(dolfinx::mesh::Mesh&& mesh_coarse) +{ + mesh_coarse.topology()->create_entities(1); + + auto [mesh_fine, parent_cell, parent_facet] + = refinement::refine(mesh_coarse, std::nullopt); + mesh_fine.topology()->create_connectivity(1, 0); + mesh_fine.topology()->create_connectivity(0, 1); + std::vector inclusion_map + = multigrid::inclusion_mapping(mesh_coarse, mesh_fine); + + CHECK_inclusion_map(mesh_coarse, mesh_fine, inclusion_map); +} + +TEMPLATE_TEST_CASE("Inclusion (interval)", "[multigrid][inclusion]", double, + float) +{ + TEST_inclusion( + dolfinx::mesh::create_interval(MPI_COMM_WORLD, 2, {0.0, 1.0})); +} + +TEMPLATE_TEST_CASE("Inclusion (triangle)", "[multigrid][inclusion]", double, + float) +{ + TEST_inclusion(dolfinx::mesh::create_rectangle( + MPI_COMM_WORLD, {{{0, 0}, {1, 1}}}, {1, 1}, mesh::CellType::triangle)); +} + +// TODO: fix! +// TEMPLATE_TEST_CASE("Inclusion (tetrahedron)", "[multigrid][inclusion]", +// double, +// float) +// { +// TEST_inclusion(dolfinx::mesh::create_box( +// MPI_COMM_WORLD, {{{0, 0, 0}, {1, 1, 1}}}, {1, 1, 1}, +// mesh::CellType::tetrahedron)); +// } diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt index 0f440260214..5425d10896f 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -44,6 +44,7 @@ nanobind_add_module( dolfinx/wrappers/la.cpp dolfinx/wrappers/log.cpp dolfinx/wrappers/mesh.cpp + dolfinx/wrappers/multigrid.cpp dolfinx/wrappers/petsc.cpp dolfinx/wrappers/refinement.cpp ) diff --git a/python/dolfinx/transfer.py b/python/dolfinx/transfer.py new file mode 100644 index 00000000000..c7fc59c0aba --- /dev/null +++ b/python/dolfinx/transfer.py @@ -0,0 +1,17 @@ +# Copyright (C) 2024 Paul T. Kühner +# +# This file is part of DOLFINx (https://www.fenicsproject.org) +# +# SPDX-License-Identifier: LGPL-3.0-or-later + +import numpy as np +from numpy.typing import NDArray + +from dolfinx.cpp.transfer import inclusion_mapping as _inclusion_mapping +from dolfinx.mesh import Mesh + +__all__ = ["inclusion_mapping"] + + +def inclusion_mapping(mesh_from: Mesh, mesh_to: Mesh) -> NDArray[np.int64]: + return _inclusion_mapping(mesh_from._cpp_object, mesh_to._cpp_object) diff --git a/python/dolfinx/wrappers/dolfinx.cpp b/python/dolfinx/wrappers/dolfinx.cpp index 9c0864f9ef1..9d2ca3e716d 100644 --- a/python/dolfinx/wrappers/dolfinx.cpp +++ b/python/dolfinx/wrappers/dolfinx.cpp @@ -25,6 +25,8 @@ void la(nb::module_& m); void mesh(nb::module_& m); void nls(nb::module_& m); void refinement(nb::module_& m); +void transfer(nb::module_& m); + } // namespace dolfinx_wrappers NB_MODULE(cpp, m) @@ -74,6 +76,10 @@ NB_MODULE(cpp, m) nb::module_ refinement = m.def_submodule("refinement", "Refinement module"); dolfinx_wrappers::refinement(refinement); + // Create transfer submodule + nb::module_ transfer = m.def_submodule("transfer", "Transfer module"); + dolfinx_wrappers::transfer(transfer); + #if defined(HAS_PETSC) && defined(HAS_PETSC4PY) // PETSc-specific wrappers nb::module_ nls = m.def_submodule("nls", "Nonlinear solver module"); diff --git a/python/dolfinx/wrappers/multigrid.cpp b/python/dolfinx/wrappers/multigrid.cpp new file mode 100644 index 00000000000..73f09cbd150 --- /dev/null +++ b/python/dolfinx/wrappers/multigrid.cpp @@ -0,0 +1,38 @@ +// Copyright (C) 2024 Paul T. Kühner +// +// This file is part of DOLFINX (https://www.fenicsproject.org) +// +// SPDX-License-Identifier: LGPL-3.0-or-later + +#include +#include + +#include +#include + +#include +#include + +#include "array.h" + +namespace nb = nanobind; + +namespace dolfinx_wrappers +{ + +void transfer(nb::module_& m) +{ + m.def( + "inclusion_mapping", + [](const dolfinx::mesh::Mesh& mesh_from, + const dolfinx::mesh::Mesh& mesh_to) + { + std::vector map + = dolfinx::multigrid::inclusion_mapping(mesh_from, mesh_to); + return dolfinx_wrappers::as_nbarray(std::move(map)); + }, + nb::arg("mesh_from"), nb::arg("mesh_to"), + "Computes inclusion mapping between two meshes"); +} + +} // namespace dolfinx_wrappers diff --git a/python/test/unit/transfer/test_inclusion.py b/python/test/unit/transfer/test_inclusion.py new file mode 100644 index 00000000000..9074d7ed41e --- /dev/null +++ b/python/test/unit/transfer/test_inclusion.py @@ -0,0 +1,26 @@ +# Copyright (C) 2024 Paul T. Kühner +# +# This file is part of DOLFINx (https://www.fenicsproject.org) +# +# SPDX-License-Identifier: LGPL-3.0-or-later + +from mpi4py import MPI + +import pytest + +from dolfinx.mesh import GhostMode, create_interval, refine +from dolfinx.transfer import inclusion_mapping + + +@pytest.mark.parametrize( + "ghost_mode", [GhostMode.none, GhostMode.shared_vertex, GhostMode.shared_facet] +) +def test_1d(ghost_mode): + mesh = create_interval(MPI.COMM_WORLD, 10, (0, 1), ghost_mode=ghost_mode) + mesh_fine, _, _ = refine(mesh) + inclusion_mapping(mesh, mesh_fine) + # TODO: further testing + + +if __name__ == "__main__": + pytest.main([__file__]) From 042b7dbb218a90697a75b4a2f522c025f3ba6e0f Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Thu, 2 Jan 2025 15:00:47 +0100 Subject: [PATCH 04/44] Add cast --- cpp/test/multigrid/inclusion.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/test/multigrid/inclusion.cpp b/cpp/test/multigrid/inclusion.cpp index 7b0b533735a..55665bd9002 100644 --- a/cpp/test/multigrid/inclusion.cpp +++ b/cpp/test/multigrid/inclusion.cpp @@ -55,7 +55,7 @@ void CHECK_inclusion_map(const dolfinx::mesh::Mesh& from, displacements.data(), dolfinx::MPI::mpi_t, from.comm()); REQUIRE(static_cast(map.size()) == im_from.size_global()); - for (std::int64_t i = 0; i < map.size(); i++) + for (std::int64_t i = 0; i < static_cast(map.size()); i++) { std::array local{-1}; im_from.global_to_local(std::array{i}, local); From b301dc02ab999e438a669675abd1d6f4a119bda7 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Thu, 2 Jan 2025 15:02:22 +0100 Subject: [PATCH 05/44] Rename to multigrid --- python/dolfinx/{transfer.py => multigrid.py} | 2 +- python/dolfinx/wrappers/dolfinx.cpp | 8 ++++---- python/dolfinx/wrappers/multigrid.cpp | 2 +- .../test/unit/{transfer => multigrid}/test_inclusion.py | 2 +- 4 files changed, 7 insertions(+), 7 deletions(-) rename python/dolfinx/{transfer.py => multigrid.py} (85%) rename python/test/unit/{transfer => multigrid}/test_inclusion.py (92%) diff --git a/python/dolfinx/transfer.py b/python/dolfinx/multigrid.py similarity index 85% rename from python/dolfinx/transfer.py rename to python/dolfinx/multigrid.py index c7fc59c0aba..ae84ab3992c 100644 --- a/python/dolfinx/transfer.py +++ b/python/dolfinx/multigrid.py @@ -7,7 +7,7 @@ import numpy as np from numpy.typing import NDArray -from dolfinx.cpp.transfer import inclusion_mapping as _inclusion_mapping +from dolfinx.cpp.multigrid import inclusion_mapping as _inclusion_mapping from dolfinx.mesh import Mesh __all__ = ["inclusion_mapping"] diff --git a/python/dolfinx/wrappers/dolfinx.cpp b/python/dolfinx/wrappers/dolfinx.cpp index 9d2ca3e716d..21bcec93f4f 100644 --- a/python/dolfinx/wrappers/dolfinx.cpp +++ b/python/dolfinx/wrappers/dolfinx.cpp @@ -25,7 +25,7 @@ void la(nb::module_& m); void mesh(nb::module_& m); void nls(nb::module_& m); void refinement(nb::module_& m); -void transfer(nb::module_& m); +void multigrid(nb::module_& m); } // namespace dolfinx_wrappers @@ -76,9 +76,9 @@ NB_MODULE(cpp, m) nb::module_ refinement = m.def_submodule("refinement", "Refinement module"); dolfinx_wrappers::refinement(refinement); - // Create transfer submodule - nb::module_ transfer = m.def_submodule("transfer", "Transfer module"); - dolfinx_wrappers::transfer(transfer); + // Create multigrid submodule + nb::module_ multigrid = m.def_submodule("multigrid", "Multigrid module"); + dolfinx_wrappers::multigrid(multigrid); #if defined(HAS_PETSC) && defined(HAS_PETSC4PY) // PETSc-specific wrappers diff --git a/python/dolfinx/wrappers/multigrid.cpp b/python/dolfinx/wrappers/multigrid.cpp index 73f09cbd150..34687f2edbf 100644 --- a/python/dolfinx/wrappers/multigrid.cpp +++ b/python/dolfinx/wrappers/multigrid.cpp @@ -20,7 +20,7 @@ namespace nb = nanobind; namespace dolfinx_wrappers { -void transfer(nb::module_& m) +void multigrid(nb::module_& m) { m.def( "inclusion_mapping", diff --git a/python/test/unit/transfer/test_inclusion.py b/python/test/unit/multigrid/test_inclusion.py similarity index 92% rename from python/test/unit/transfer/test_inclusion.py rename to python/test/unit/multigrid/test_inclusion.py index 9074d7ed41e..e60167e885d 100644 --- a/python/test/unit/transfer/test_inclusion.py +++ b/python/test/unit/multigrid/test_inclusion.py @@ -9,7 +9,7 @@ import pytest from dolfinx.mesh import GhostMode, create_interval, refine -from dolfinx.transfer import inclusion_mapping +from dolfinx.multigrid import inclusion_mapping @pytest.mark.parametrize( From 93efeecfe5397b0d978c4b1bf396f7abc0775884 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Thu, 2 Jan 2025 15:30:29 +0100 Subject: [PATCH 06/44] Partial fix --- cpp/dolfinx/multigrid/inclusion.h | 2 +- cpp/test/multigrid/inclusion.cpp | 26 +++++++++++++++----------- 2 files changed, 16 insertions(+), 12 deletions(-) diff --git a/cpp/dolfinx/multigrid/inclusion.h b/cpp/dolfinx/multigrid/inclusion.h index 0325ee26a74..ca7867bfcf9 100644 --- a/cpp/dolfinx/multigrid/inclusion.h +++ b/cpp/dolfinx/multigrid/inclusion.h @@ -59,7 +59,7 @@ inclusion_mapping(const dolfinx::mesh::Mesh& mesh_from, { std::ranges::subrange vertex_from(std::next(x_from.begin(), 3 * i), std::next(x_from.begin(), 3 * (i + 1))); - for (std::int64_t j = 0; j < im_to.size_local() + im_to.num_ghosts(); j++) + for (std::int64_t j = 0; j < im_to.size_local(); j++) // + im_to.num_ghosts() TODO { std::ranges::subrange vertex_to(std::next(x_to.begin(), 3 * j), std::next(x_to.begin(), 3 * (j + 1))); diff --git a/cpp/test/multigrid/inclusion.cpp b/cpp/test/multigrid/inclusion.cpp index 55665bd9002..fd7e767a22c 100644 --- a/cpp/test/multigrid/inclusion.cpp +++ b/cpp/test/multigrid/inclusion.cpp @@ -62,6 +62,11 @@ void CHECK_inclusion_map(const dolfinx::mesh::Mesh& from, if (local[0] == -1) continue; + std::cout << "(" << from.geometry().x()[3 * local[0]] << ", " + << from.geometry().x()[3 * local[0] + 1] << ", " + << from.geometry().x()[3 * local[0] + 2] << ") == (" + << global_x_to[3 * map[i]] << ", " << global_x_to[3 * map[i] + 1] + << ", " << global_x_to[3 * map[i] + 2] << ")" << std::endl; CHECK(std::abs(from.geometry().x()[3 * local[0]] - global_x_to[3 * map[i]]) < std::numeric_limits::epsilon()); CHECK(std::abs(from.geometry().x()[3 * local[0] + 1] @@ -74,7 +79,7 @@ void CHECK_inclusion_map(const dolfinx::mesh::Mesh& from, } /// Performs one uniform refinement and checks the inclusion map between coarse -/// and fine mesh against the (provided) list. +/// and fine mesh. template void TEST_inclusion(dolfinx::mesh::Mesh&& mesh_coarse) { @@ -101,15 +106,14 @@ TEMPLATE_TEST_CASE("Inclusion (triangle)", "[multigrid][inclusion]", double, float) { TEST_inclusion(dolfinx::mesh::create_rectangle( - MPI_COMM_WORLD, {{{0, 0}, {1, 1}}}, {1, 1}, mesh::CellType::triangle)); + MPI_COMM_WORLD, {{{0, 0}, {1, 1}}}, {2, 2}, mesh::CellType::triangle)); } -// TODO: fix! -// TEMPLATE_TEST_CASE("Inclusion (tetrahedron)", "[multigrid][inclusion]", -// double, -// float) -// { -// TEST_inclusion(dolfinx::mesh::create_box( -// MPI_COMM_WORLD, {{{0, 0, 0}, {1, 1, 1}}}, {1, 1, 1}, -// mesh::CellType::tetrahedron)); -// } +TEMPLATE_TEST_CASE("Inclusion (tetrahedron)", "[multigrid][inclusion]", double, + float) +{ + // TODO: n = {2, 2, 2} fails + TEST_inclusion(dolfinx::mesh::create_box( + MPI_COMM_WORLD, {{{0, 0, 0}, {1, 1, 1}}}, {1, 1, 1}, + mesh::CellType::tetrahedron)); +} From be1987be5cb39eb4b34e4488333131a31e6925d5 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Sat, 4 Jan 2025 08:06:25 +0100 Subject: [PATCH 07/44] Local check run also over ghosts --- cpp/dolfinx/multigrid/inclusion.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cpp/dolfinx/multigrid/inclusion.h b/cpp/dolfinx/multigrid/inclusion.h index ca7867bfcf9..06c217e7d5e 100644 --- a/cpp/dolfinx/multigrid/inclusion.h +++ b/cpp/dolfinx/multigrid/inclusion.h @@ -55,11 +55,11 @@ inclusion_mapping(const dolfinx::mesh::Mesh& mesh_from, auto to_global_to = build_global_to_local(im_to); auto to_global_from = build_global_to_local(im_from); - for (std::int32_t i = 0; i < im_from.size_local(); i++) + for (std::int32_t i = 0; i < im_from.size_local() + im_from.num_ghosts(); i++) { std::ranges::subrange vertex_from(std::next(x_from.begin(), 3 * i), std::next(x_from.begin(), 3 * (i + 1))); - for (std::int64_t j = 0; j < im_to.size_local(); j++) // + im_to.num_ghosts() TODO + for (std::int64_t j = 0; j < im_to.size_local() + im_to.num_ghosts(); j++) { std::ranges::subrange vertex_to(std::next(x_to.begin(), 3 * j), std::next(x_to.begin(), 3 * (j + 1))); From c0f8e3a07b646203ed03adf83773004f67d16a93 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Sat, 4 Jan 2025 08:34:53 +0100 Subject: [PATCH 08/44] Switch to central definition of gather_global and unit test --- cpp/dolfinx/multigrid/inclusion.h | 56 ++++++++++++++------------ cpp/test/multigrid/inclusion.cpp | 67 +++++++++++++++++++++---------- 2 files changed, 75 insertions(+), 48 deletions(-) diff --git a/cpp/dolfinx/multigrid/inclusion.h b/cpp/dolfinx/multigrid/inclusion.h index 06c217e7d5e..04a3fd237e9 100644 --- a/cpp/dolfinx/multigrid/inclusion.h +++ b/cpp/dolfinx/multigrid/inclusion.h @@ -22,6 +22,33 @@ namespace dolfinx::multigrid { +template +std::vector gather_global(std::span local, std::int64_t global_size, + MPI_Comm comm) +{ + // 1) exchange local sizes + std::vector local_sizes(dolfinx::MPI::size(comm)); + { + std::array tmp{local.size()}; + MPI_Allgather(&tmp, 1, MPI_INT32_T, local_sizes.data(), 1, MPI_INT32_T, + comm); + } + + // 2) compute displacement vector + std::vector displacements(local_sizes.size() + 1, 0); + std::partial_sum(local_sizes.begin(), local_sizes.end(), + displacements.begin() + 1); + + // 3) Allgather global x vector + std::vector global_x_to(global_size); + MPI_Allgatherv(local.data(), local.size(), dolfinx::MPI::mpi_t, + global_x_to.data(), local_sizes.data(), displacements.data(), + dolfinx::MPI::mpi_t, comm); + + return global_x_to; +} + + template std::vector inclusion_mapping(const dolfinx::mesh::Mesh& mesh_from, @@ -93,32 +120,9 @@ inclusion_mapping(const dolfinx::mesh::Mesh& mesh_from, return result; // Build global to vertex list - - // 1) exchange local sizes - std::vector local_sizes(dolfinx::MPI::size(mesh_from.comm())); - { - std::array tmp{im_to.size_local() * 3}; - MPI_Allgather(&tmp, 1, MPI_INT32_T, local_sizes.data(), 1, MPI_INT32_T, - mesh_from.comm()); - } - - // for (auto ls : local_sizes) - // std::cout << ls << ", "; - - // 2) compute displacement vector - std::vector displacements(local_sizes.size() + 1, 0); - std::partial_sum(local_sizes.begin(), local_sizes.end(), - displacements.begin() + 1); - - // for (auto ls : displacements) - // std::cout << ls << ", "; - - // 3) Allgather global x vector - std::vector global_x_to(im_to.size_global() * 3); - MPI_Allgatherv(mesh_to.geometry().x().data(), im_to.size_local() * 3, - dolfinx::MPI::mpi_t, global_x_to.data(), local_sizes.data(), - displacements.data(), dolfinx::MPI::mpi_t, - mesh_from.comm()); + auto global_x_to + = gather_global(mesh_to.geometry().x().subspan(0, im_to.size_local() * 3), + im_to.size_global() * 3, mesh_to.comm()); // Recheck indices on global data structure for (std::int32_t i = 0; i < im_from.size_local(); i++) diff --git a/cpp/test/multigrid/inclusion.cpp b/cpp/test/multigrid/inclusion.cpp index fd7e767a22c..5f8f5294ad2 100644 --- a/cpp/test/multigrid/inclusion.cpp +++ b/cpp/test/multigrid/inclusion.cpp @@ -6,6 +6,7 @@ #include #include +#include #include #include #include @@ -13,6 +14,7 @@ #include +#include #include #include #include @@ -27,6 +29,25 @@ using namespace dolfinx; using namespace Catch::Matchers; +TEMPLATE_TEST_CASE("Gather global", "[multigrid]", double, float) +{ + using T = TestType; + MPI_Comm comm = MPI_COMM_WORLD; + auto comm_size = dolfinx::MPI::size(comm); + auto rank = dolfinx::MPI::rank(comm); + std::vector local(static_cast(rank), static_cast(rank + 1)); + + std::vector global = multigrid::gather_global( + std::span(local.data(), local.size()), comm_size * 2, comm); + + CHECK(global.size() == 2 * comm_size); + for (std::size_t i = 0; i < comm_size; i++) + { + CHECK(global[2 * i] == Catch::Approx(i)); + CHECK(global[2 * i + 1] == Catch::Approx(i + 1)); + } +} + template void CHECK_inclusion_map(const dolfinx::mesh::Mesh& from, const dolfinx::mesh::Mesh& to, @@ -95,25 +116,27 @@ void TEST_inclusion(dolfinx::mesh::Mesh&& mesh_coarse) CHECK_inclusion_map(mesh_coarse, mesh_fine, inclusion_map); } -TEMPLATE_TEST_CASE("Inclusion (interval)", "[multigrid][inclusion]", double, - float) -{ - TEST_inclusion( - dolfinx::mesh::create_interval(MPI_COMM_WORLD, 2, {0.0, 1.0})); -} - -TEMPLATE_TEST_CASE("Inclusion (triangle)", "[multigrid][inclusion]", double, - float) -{ - TEST_inclusion(dolfinx::mesh::create_rectangle( - MPI_COMM_WORLD, {{{0, 0}, {1, 1}}}, {2, 2}, mesh::CellType::triangle)); -} - -TEMPLATE_TEST_CASE("Inclusion (tetrahedron)", "[multigrid][inclusion]", double, - float) -{ - // TODO: n = {2, 2, 2} fails - TEST_inclusion(dolfinx::mesh::create_box( - MPI_COMM_WORLD, {{{0, 0, 0}, {1, 1, 1}}}, {1, 1, 1}, - mesh::CellType::tetrahedron)); -} +// TEMPLATE_TEST_CASE("Inclusion (interval)", "[multigrid][inclusion]", double, +// float) +// { +// TEST_inclusion( +// dolfinx::mesh::create_interval(MPI_COMM_WORLD, 2, +// {0.0, 1.0})); +// } + +// TEMPLATE_TEST_CASE("Inclusion (triangle)", "[multigrid][inclusion]", double, +// float) +// { +// TEST_inclusion(dolfinx::mesh::create_rectangle( +// MPI_COMM_WORLD, {{{0, 0}, {1, 1}}}, {3, 3}, mesh::CellType::triangle)); +// } + +// // TEMPLATE_TEST_CASE("Inclusion (tetrahedron)", "[multigrid][inclusion]", +// // double, +// // float) +// // { +// // // TODO: n = {2, 2, 2} fails +// // TEST_inclusion(dolfinx::mesh::create_box( +// // MPI_COMM_WORLD, {{{0, 0, 0}, {1, 1, 1}}}, {1, 1, 1}, +// // mesh::CellType::tetrahedron)); +// // } From 98b32262eb82342a39f0f4ea773f8aa84f81f85a Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Sat, 4 Jan 2025 08:52:44 +0100 Subject: [PATCH 09/44] Tidy gather global --- cpp/dolfinx/multigrid/inclusion.h | 13 ++++----- cpp/test/multigrid/inclusion.cpp | 47 +++++++++++-------------------- 2 files changed, 22 insertions(+), 38 deletions(-) diff --git a/cpp/dolfinx/multigrid/inclusion.h b/cpp/dolfinx/multigrid/inclusion.h index 04a3fd237e9..224b84b3ad2 100644 --- a/cpp/dolfinx/multigrid/inclusion.h +++ b/cpp/dolfinx/multigrid/inclusion.h @@ -23,7 +23,7 @@ namespace dolfinx::multigrid { template -std::vector gather_global(std::span local, std::int64_t global_size, +std::vector gather_global(std::span local, std::int64_t global_size, MPI_Comm comm) { // 1) exchange local sizes @@ -39,16 +39,15 @@ std::vector gather_global(std::span local, std::int64_t global_size, std::partial_sum(local_sizes.begin(), local_sizes.end(), displacements.begin() + 1); - // 3) Allgather global x vector - std::vector global_x_to(global_size); + // 3) Allgather global vector + std::vector global(global_size); MPI_Allgatherv(local.data(), local.size(), dolfinx::MPI::mpi_t, - global_x_to.data(), local_sizes.data(), displacements.data(), + global.data(), local_sizes.data(), displacements.data(), dolfinx::MPI::mpi_t, comm); - return global_x_to; + return global; } - template std::vector inclusion_mapping(const dolfinx::mesh::Mesh& mesh_from, @@ -120,7 +119,7 @@ inclusion_mapping(const dolfinx::mesh::Mesh& mesh_from, return result; // Build global to vertex list - auto global_x_to + std::vector global_x_to = gather_global(mesh_to.geometry().x().subspan(0, im_to.size_local() * 3), im_to.size_global() * 3, mesh_to.comm()); diff --git a/cpp/test/multigrid/inclusion.cpp b/cpp/test/multigrid/inclusion.cpp index 5f8f5294ad2..ae16d8468b8 100644 --- a/cpp/test/multigrid/inclusion.cpp +++ b/cpp/test/multigrid/inclusion.cpp @@ -35,13 +35,13 @@ TEMPLATE_TEST_CASE("Gather global", "[multigrid]", double, float) MPI_Comm comm = MPI_COMM_WORLD; auto comm_size = dolfinx::MPI::size(comm); auto rank = dolfinx::MPI::rank(comm); - std::vector local(static_cast(rank), static_cast(rank + 1)); + std::vector local{static_cast(rank), static_cast(rank + 1)}; - std::vector global = multigrid::gather_global( + std::vector global = multigrid::gather_global( std::span(local.data(), local.size()), comm_size * 2, comm); - CHECK(global.size() == 2 * comm_size); - for (std::size_t i = 0; i < comm_size; i++) + CHECK(global.size() == static_cast(2 * comm_size)); + for (int i = 0; i < comm_size; i++) { CHECK(global[2 * i] == Catch::Approx(i)); CHECK(global[2 * i + 1] == Catch::Approx(i + 1)); @@ -56,24 +56,9 @@ void CHECK_inclusion_map(const dolfinx::mesh::Mesh& from, const common::IndexMap& im_from = *from.topology()->index_map(0); const common::IndexMap& im_to = *to.topology()->index_map(0); - // 1) exchange local sizes - std::vector local_sizes(dolfinx::MPI::size(from.comm())); - { - std::array tmp{im_to.size_local() * 3}; - MPI_Allgather(&tmp, 1, MPI_INT32_T, local_sizes.data(), 1, MPI_INT32_T, - from.comm()); - } - - // 2) compute displacement vector - std::vector displacements(local_sizes.size() + 1, 0); - std::partial_sum(local_sizes.begin(), local_sizes.end(), - displacements.begin() + 1); - - // 3) Allgather global x vector - std::vector global_x_to(im_to.size_global() * 3); - MPI_Allgatherv(to.geometry().x().data(), im_to.size_local() * 3, - dolfinx::MPI::mpi_t, global_x_to.data(), local_sizes.data(), - displacements.data(), dolfinx::MPI::mpi_t, from.comm()); + std::vector global_x_to = multigrid::gather_global( + to.geometry().x().subspan(0, im_to.size_local() * 3), + im_to.size_global() * 3, to.comm()); REQUIRE(static_cast(map.size()) == im_from.size_global()); for (std::int64_t i = 0; i < static_cast(map.size()); i++) @@ -131,12 +116,12 @@ void TEST_inclusion(dolfinx::mesh::Mesh&& mesh_coarse) // MPI_COMM_WORLD, {{{0, 0}, {1, 1}}}, {3, 3}, mesh::CellType::triangle)); // } -// // TEMPLATE_TEST_CASE("Inclusion (tetrahedron)", "[multigrid][inclusion]", -// // double, -// // float) -// // { -// // // TODO: n = {2, 2, 2} fails -// // TEST_inclusion(dolfinx::mesh::create_box( -// // MPI_COMM_WORLD, {{{0, 0, 0}, {1, 1, 1}}}, {1, 1, 1}, -// // mesh::CellType::tetrahedron)); -// // } +// TEMPLATE_TEST_CASE("Inclusion (tetrahedron)", "[multigrid][inclusion]", +// double, +// float) +// { +// // TODO: n = {2, 2, 2} fails +// TEST_inclusion(dolfinx::mesh::create_box( +// MPI_COMM_WORLD, {{{0, 0, 0}, {1, 1, 1}}}, {1, 1, 1}, +// mesh::CellType::tetrahedron)); +// } From 751a242a3344d4bf1ae3e44e6c44953caf9952a6 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Sat, 4 Jan 2025 09:04:14 +0100 Subject: [PATCH 10/44] Extended sequential tests passing --- cpp/test/multigrid/inclusion.cpp | 54 ++++++++++++++++++-------------- 1 file changed, 30 insertions(+), 24 deletions(-) diff --git a/cpp/test/multigrid/inclusion.cpp b/cpp/test/multigrid/inclusion.cpp index ae16d8468b8..cd5ba43daba 100644 --- a/cpp/test/multigrid/inclusion.cpp +++ b/cpp/test/multigrid/inclusion.cpp @@ -101,27 +101,33 @@ void TEST_inclusion(dolfinx::mesh::Mesh&& mesh_coarse) CHECK_inclusion_map(mesh_coarse, mesh_fine, inclusion_map); } -// TEMPLATE_TEST_CASE("Inclusion (interval)", "[multigrid][inclusion]", double, -// float) -// { -// TEST_inclusion( -// dolfinx::mesh::create_interval(MPI_COMM_WORLD, 2, -// {0.0, 1.0})); -// } - -// TEMPLATE_TEST_CASE("Inclusion (triangle)", "[multigrid][inclusion]", double, -// float) -// { -// TEST_inclusion(dolfinx::mesh::create_rectangle( -// MPI_COMM_WORLD, {{{0, 0}, {1, 1}}}, {3, 3}, mesh::CellType::triangle)); -// } - -// TEMPLATE_TEST_CASE("Inclusion (tetrahedron)", "[multigrid][inclusion]", -// double, -// float) -// { -// // TODO: n = {2, 2, 2} fails -// TEST_inclusion(dolfinx::mesh::create_box( -// MPI_COMM_WORLD, {{{0, 0, 0}, {1, 1, 1}}}, {1, 1, 1}, -// mesh::CellType::tetrahedron)); -// } +TEMPLATE_TEST_CASE("Inclusion (interval)", "[multigrid][inclusion]", double, + float) +{ + for (auto n : {1, 2, 3}) + { + TEST_inclusion(dolfinx::mesh::create_interval(MPI_COMM_WORLD, 2, + {0.0, 1.0})); + } +} + +TEMPLATE_TEST_CASE("Inclusion (triangle)", "[multigrid][inclusion]", double, + float) +{ + for (auto n : {1, 2, 3}) + { + TEST_inclusion(dolfinx::mesh::create_rectangle( + MPI_COMM_WORLD, {{{0, 0}, {1, 1}}}, {n, n}, mesh::CellType::triangle)); + } +} + +TEMPLATE_TEST_CASE("Inclusion (tetrahedron)", "[multigrid][inclusion]", double, + float) +{ + for (auto n : {1, 2, 3}) + { + TEST_inclusion(dolfinx::mesh::create_box( + MPI_COMM_WORLD, {{{0, 0, 0}, {1, 1, 1}}}, {n, n, n}, + mesh::CellType::tetrahedron)); + } +} From d0cee96f131fe16caf56184ef1209876eec96bfa Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Sat, 4 Jan 2025 09:13:28 +0100 Subject: [PATCH 11/44] Switch to inplace MPI communication --- cpp/dolfinx/multigrid/inclusion.h | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/cpp/dolfinx/multigrid/inclusion.h b/cpp/dolfinx/multigrid/inclusion.h index 224b84b3ad2..4ec47b70928 100644 --- a/cpp/dolfinx/multigrid/inclusion.h +++ b/cpp/dolfinx/multigrid/inclusion.h @@ -110,13 +110,12 @@ inclusion_mapping(const dolfinx::mesh::Mesh& mesh_from, // map holds at this point for every original local index the corresponding // mapped global index (if it was available on the same process on the to // mesh). - std::vector result(map.size(), -1); - MPI_Allreduce(map.data(), result.data(), map.size(), + MPI_Allreduce(MPI_IN_PLACE, map.data(), map.size(), dolfinx::MPI::mpi_t, MPI_MAX, mesh_from.comm()); - if (std::ranges::all_of(result, [](auto e) { return e >= 0; })) + if (std::ranges::all_of(map, [](auto e) { return e >= 0; })) // All vertices indentified - return result; + return map; // Build global to vertex list std::vector global_x_to @@ -144,12 +143,12 @@ inclusion_mapping(const dolfinx::mesh::Mesh& mesh_from, } } - MPI_Allreduce(map.data(), result.data(), map.size(), + MPI_Allreduce(MPI_IN_PLACE, map.data(), map.size(), dolfinx::MPI::mpi_t, MPI_MAX, mesh_from.comm()); - assert(std::ranges::all_of(result, [&](auto e) + assert(std::ranges::all_of(map, [&](auto e) { return e >= 0 && e < im_to.size_global(); })); - return result; + return map; } } // namespace dolfinx::multigrid From 6fe584a85900f314f41775f81d2178bcab3956c2 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Sat, 4 Jan 2025 09:28:53 +0100 Subject: [PATCH 12/44] Switch to local computation --- cpp/dolfinx/multigrid/inclusion.h | 9 +++++---- cpp/test/multigrid/inclusion.cpp | 31 ++++++++++++------------------- 2 files changed, 17 insertions(+), 23 deletions(-) diff --git a/cpp/dolfinx/multigrid/inclusion.h b/cpp/dolfinx/multigrid/inclusion.h index 4ec47b70928..d6e685a811c 100644 --- a/cpp/dolfinx/multigrid/inclusion.h +++ b/cpp/dolfinx/multigrid/inclusion.h @@ -63,7 +63,8 @@ inclusion_mapping(const dolfinx::mesh::Mesh& mesh_from, const common::IndexMap& im_from = *mesh_from.topology()->index_map(0); const common::IndexMap& im_to = *mesh_to.topology()->index_map(0); - std::vector map(im_from.size_global(), -1); + std::vector map(im_from.size_local() + im_from.num_ghosts(), + -1); std::span x_from = mesh_from.geometry().x(); std::span x_to = mesh_to.geometry().x(); @@ -94,7 +95,7 @@ inclusion_mapping(const dolfinx::mesh::Mesh& mesh_from, vertex_from, vertex_to, [](T a, T b) { return std::abs(a - b) <= std::numeric_limits::epsilon(); })) { - map[to_global_from(i)] = to_global_to(j); + map[i] = to_global_to(j); break; } } @@ -123,7 +124,7 @@ inclusion_mapping(const dolfinx::mesh::Mesh& mesh_from, im_to.size_global() * 3, mesh_to.comm()); // Recheck indices on global data structure - for (std::int32_t i = 0; i < im_from.size_local(); i++) + for (std::int32_t i = 0; i < im_from.size_local() + im_from.num_ghosts(); i++) { std::ranges::subrange vertex_from(std::next(x_from.begin(), 3 * i), std::next(x_from.begin(), 3 * (i + 1))); @@ -137,7 +138,7 @@ inclusion_mapping(const dolfinx::mesh::Mesh& mesh_from, vertex_from, vertex_to, [](T a, T b) { return std::abs(a - b) <= std::numeric_limits::epsilon(); })) { - map[to_global_from(i)] = j; + map[i] = j; break; } } diff --git a/cpp/test/multigrid/inclusion.cpp b/cpp/test/multigrid/inclusion.cpp index cd5ba43daba..3ae804f3dc6 100644 --- a/cpp/test/multigrid/inclusion.cpp +++ b/cpp/test/multigrid/inclusion.cpp @@ -4,7 +4,6 @@ // // SPDX-License-Identifier: LGPL-3.0-or-later -#include #include #include #include @@ -60,26 +59,20 @@ void CHECK_inclusion_map(const dolfinx::mesh::Mesh& from, to.geometry().x().subspan(0, im_to.size_local() * 3), im_to.size_global() * 3, to.comm()); - REQUIRE(static_cast(map.size()) == im_from.size_global()); + REQUIRE(static_cast(map.size()) + == im_from.size_local() + im_from.num_ghosts()); for (std::int64_t i = 0; i < static_cast(map.size()); i++) { - std::array local{-1}; - im_from.global_to_local(std::array{i}, local); - if (local[0] == -1) - continue; - - std::cout << "(" << from.geometry().x()[3 * local[0]] << ", " - << from.geometry().x()[3 * local[0] + 1] << ", " - << from.geometry().x()[3 * local[0] + 2] << ") == (" + std::cout << "(" << from.geometry().x()[3 * i] << ", " + << from.geometry().x()[3 * i + 1] << ", " + << from.geometry().x()[3 * i + 2] << ") == (" << global_x_to[3 * map[i]] << ", " << global_x_to[3 * map[i] + 1] << ", " << global_x_to[3 * map[i] + 2] << ")" << std::endl; - CHECK(std::abs(from.geometry().x()[3 * local[0]] - global_x_to[3 * map[i]]) + CHECK(std::abs(from.geometry().x()[3 * i] - global_x_to[3 * map[i]]) < std::numeric_limits::epsilon()); - CHECK(std::abs(from.geometry().x()[3 * local[0] + 1] - - global_x_to[3 * map[i] + 1]) + CHECK(std::abs(from.geometry().x()[3 * i + 1] - global_x_to[3 * map[i] + 1]) < std::numeric_limits::epsilon()); - CHECK(std::abs(from.geometry().x()[3 * local[0] + 2] - - global_x_to[3 * map[i] + 2]) + CHECK(std::abs(from.geometry().x()[3 * i + 2] - global_x_to[3 * map[i] + 2]) < std::numeric_limits::epsilon()); } } @@ -104,9 +97,9 @@ void TEST_inclusion(dolfinx::mesh::Mesh&& mesh_coarse) TEMPLATE_TEST_CASE("Inclusion (interval)", "[multigrid][inclusion]", double, float) { - for (auto n : {1, 2, 3}) + for (auto n : {10}) { - TEST_inclusion(dolfinx::mesh::create_interval(MPI_COMM_WORLD, 2, + TEST_inclusion(dolfinx::mesh::create_interval(MPI_COMM_WORLD, n, {0.0, 1.0})); } } @@ -114,7 +107,7 @@ TEMPLATE_TEST_CASE("Inclusion (interval)", "[multigrid][inclusion]", double, TEMPLATE_TEST_CASE("Inclusion (triangle)", "[multigrid][inclusion]", double, float) { - for (auto n : {1, 2, 3}) + for (auto n : {5}) { TEST_inclusion(dolfinx::mesh::create_rectangle( MPI_COMM_WORLD, {{{0, 0}, {1, 1}}}, {n, n}, mesh::CellType::triangle)); @@ -124,7 +117,7 @@ TEMPLATE_TEST_CASE("Inclusion (triangle)", "[multigrid][inclusion]", double, TEMPLATE_TEST_CASE("Inclusion (tetrahedron)", "[multigrid][inclusion]", double, float) { - for (auto n : {1, 2, 3}) + for (auto n : {5}) { TEST_inclusion(dolfinx::mesh::create_box( MPI_COMM_WORLD, {{{0, 0, 0}, {1, 1, 1}}}, {n, n, n}, From 1d5e1cde37acf14636a4f578519fb92ee2d3de81 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Sat, 4 Jan 2025 09:51:28 +0100 Subject: [PATCH 13/44] Fix logic --- cpp/dolfinx/multigrid/inclusion.h | 17 ++++++----------- cpp/test/multigrid/inclusion.cpp | 5 ----- 2 files changed, 6 insertions(+), 16 deletions(-) diff --git a/cpp/dolfinx/multigrid/inclusion.h b/cpp/dolfinx/multigrid/inclusion.h index d6e685a811c..61cd40c7924 100644 --- a/cpp/dolfinx/multigrid/inclusion.h +++ b/cpp/dolfinx/multigrid/inclusion.h @@ -108,14 +108,12 @@ inclusion_mapping(const dolfinx::mesh::Mesh& mesh_from, return map; } - // map holds at this point for every original local index the corresponding - // mapped global index (if it was available on the same process on the to - // mesh). - MPI_Allreduce(MPI_IN_PLACE, map.data(), map.size(), - dolfinx::MPI::mpi_t, MPI_MAX, mesh_from.comm()); - - if (std::ranges::all_of(map, [](auto e) { return e >= 0; })) - // All vertices indentified + bool locally_complete + = std::ranges::all_of(map, [](auto e) { return e >= 0; }); + bool globally_complete = false; + MPI_Allreduce(&locally_complete, &locally_complete, 1, MPI_CXX_BOOL, MPI_LAND, + mesh_to.comm()); + if (globally_complete) return map; // Build global to vertex list @@ -144,9 +142,6 @@ inclusion_mapping(const dolfinx::mesh::Mesh& mesh_from, } } - MPI_Allreduce(MPI_IN_PLACE, map.data(), map.size(), - dolfinx::MPI::mpi_t, MPI_MAX, mesh_from.comm()); - assert(std::ranges::all_of(map, [&](auto e) { return e >= 0 && e < im_to.size_global(); })); return map; diff --git a/cpp/test/multigrid/inclusion.cpp b/cpp/test/multigrid/inclusion.cpp index 3ae804f3dc6..66ca7b0bcbb 100644 --- a/cpp/test/multigrid/inclusion.cpp +++ b/cpp/test/multigrid/inclusion.cpp @@ -63,11 +63,6 @@ void CHECK_inclusion_map(const dolfinx::mesh::Mesh& from, == im_from.size_local() + im_from.num_ghosts()); for (std::int64_t i = 0; i < static_cast(map.size()); i++) { - std::cout << "(" << from.geometry().x()[3 * i] << ", " - << from.geometry().x()[3 * i + 1] << ", " - << from.geometry().x()[3 * i + 2] << ") == (" - << global_x_to[3 * map[i]] << ", " << global_x_to[3 * map[i] + 1] - << ", " << global_x_to[3 * map[i] + 2] << ")" << std::endl; CHECK(std::abs(from.geometry().x()[3 * i] - global_x_to[3 * map[i]]) < std::numeric_limits::epsilon()); CHECK(std::abs(from.geometry().x()[3 * i + 1] - global_x_to[3 * map[i] + 1]) From f5d5e1ac77370c495d3f3b9ea1ebc9a14d3247d7 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Sat, 4 Jan 2025 09:58:17 +0100 Subject: [PATCH 14/44] Add all to all control flag --- cpp/dolfinx/multigrid/inclusion.h | 11 ++++++++++- cpp/test/multigrid/inclusion.cpp | 2 +- 2 files changed, 11 insertions(+), 2 deletions(-) diff --git a/cpp/dolfinx/multigrid/inclusion.h b/cpp/dolfinx/multigrid/inclusion.h index 61cd40c7924..a371f9b6031 100644 --- a/cpp/dolfinx/multigrid/inclusion.h +++ b/cpp/dolfinx/multigrid/inclusion.h @@ -11,6 +11,7 @@ #include #include #include +#include #include #include @@ -51,7 +52,8 @@ std::vector gather_global(std::span local, std::int64_t global_size, template std::vector inclusion_mapping(const dolfinx::mesh::Mesh& mesh_from, - const dolfinx::mesh::Mesh& mesh_to) + const dolfinx::mesh::Mesh& mesh_to, + bool allow_all_to_all = false) { { // Check comms equal @@ -116,6 +118,13 @@ inclusion_mapping(const dolfinx::mesh::Mesh& mesh_from, if (globally_complete) return map; + if (!allow_all_to_all) + { + throw std::runtime_error( + "Parallelization of mesh requires all to all communication to compute " + "inclusion map, but allow_all_to_all is not set."); + } + // Build global to vertex list std::vector global_x_to = gather_global(mesh_to.geometry().x().subspan(0, im_to.size_local() * 3), diff --git a/cpp/test/multigrid/inclusion.cpp b/cpp/test/multigrid/inclusion.cpp index 66ca7b0bcbb..d77d14990e3 100644 --- a/cpp/test/multigrid/inclusion.cpp +++ b/cpp/test/multigrid/inclusion.cpp @@ -84,7 +84,7 @@ void TEST_inclusion(dolfinx::mesh::Mesh&& mesh_coarse) mesh_fine.topology()->create_connectivity(1, 0); mesh_fine.topology()->create_connectivity(0, 1); std::vector inclusion_map - = multigrid::inclusion_mapping(mesh_coarse, mesh_fine); + = multigrid::inclusion_mapping(mesh_coarse, mesh_fine, true); CHECK_inclusion_map(mesh_coarse, mesh_fine, inclusion_map); } From 7450830c719fb3ca6c3f304b2419f8d45ad9476d Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Sat, 4 Jan 2025 10:04:05 +0100 Subject: [PATCH 15/44] Add TODO --- cpp/dolfinx/multigrid/inclusion.h | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/cpp/dolfinx/multigrid/inclusion.h b/cpp/dolfinx/multigrid/inclusion.h index a371f9b6031..f4b17f8d06a 100644 --- a/cpp/dolfinx/multigrid/inclusion.h +++ b/cpp/dolfinx/multigrid/inclusion.h @@ -133,6 +133,10 @@ inclusion_mapping(const dolfinx::mesh::Mesh& mesh_from, // Recheck indices on global data structure for (std::int32_t i = 0; i < im_from.size_local() + im_from.num_ghosts(); i++) { + // TODO: + // if (map[i] >= 0) + // continue; + std::ranges::subrange vertex_from(std::next(x_from.begin(), 3 * i), std::next(x_from.begin(), 3 * (i + 1))); for (std::int64_t j = 0; j < im_to.size_global(); j++) From dffd614fe215435bf09831e6a6a874442d5ac0d4 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Sat, 4 Jan 2025 10:22:59 +0100 Subject: [PATCH 16/44] Finish python export and tests --- python/dolfinx/multigrid.py | 6 +++-- python/dolfinx/wrappers/multigrid.cpp | 7 +++--- python/test/unit/multigrid/test_inclusion.py | 26 ++++++++++++++++++-- 3 files changed, 32 insertions(+), 7 deletions(-) diff --git a/python/dolfinx/multigrid.py b/python/dolfinx/multigrid.py index ae84ab3992c..2fde4813f72 100644 --- a/python/dolfinx/multigrid.py +++ b/python/dolfinx/multigrid.py @@ -13,5 +13,7 @@ __all__ = ["inclusion_mapping"] -def inclusion_mapping(mesh_from: Mesh, mesh_to: Mesh) -> NDArray[np.int64]: - return _inclusion_mapping(mesh_from._cpp_object, mesh_to._cpp_object) +def inclusion_mapping( + mesh_from: Mesh, mesh_to: Mesh, allow_all_to_all: bool = False +) -> NDArray[np.int64]: + return _inclusion_mapping(mesh_from._cpp_object, mesh_to._cpp_object, allow_all_to_all) diff --git a/python/dolfinx/wrappers/multigrid.cpp b/python/dolfinx/wrappers/multigrid.cpp index 34687f2edbf..619d8fec848 100644 --- a/python/dolfinx/wrappers/multigrid.cpp +++ b/python/dolfinx/wrappers/multigrid.cpp @@ -25,13 +25,14 @@ void multigrid(nb::module_& m) m.def( "inclusion_mapping", [](const dolfinx::mesh::Mesh& mesh_from, - const dolfinx::mesh::Mesh& mesh_to) + const dolfinx::mesh::Mesh& mesh_to, bool allow_all_to_all) { std::vector map - = dolfinx::multigrid::inclusion_mapping(mesh_from, mesh_to); + = dolfinx::multigrid::inclusion_mapping(mesh_from, mesh_to, + allow_all_to_all); return dolfinx_wrappers::as_nbarray(std::move(map)); }, - nb::arg("mesh_from"), nb::arg("mesh_to"), + nb::arg("mesh_from"), nb::arg("mesh_to"), nb::arg("allow_all_to_all"), "Computes inclusion mapping between two meshes"); } diff --git a/python/test/unit/multigrid/test_inclusion.py b/python/test/unit/multigrid/test_inclusion.py index e60167e885d..92c995de67e 100644 --- a/python/test/unit/multigrid/test_inclusion.py +++ b/python/test/unit/multigrid/test_inclusion.py @@ -8,7 +8,7 @@ import pytest -from dolfinx.mesh import GhostMode, create_interval, refine +from dolfinx.mesh import GhostMode, create_interval, create_unit_cube, create_unit_square, refine from dolfinx.multigrid import inclusion_mapping @@ -19,7 +19,29 @@ def test_1d(ghost_mode): mesh = create_interval(MPI.COMM_WORLD, 10, (0, 1), ghost_mode=ghost_mode) mesh_fine, _, _ = refine(mesh) inclusion_mapping(mesh, mesh_fine) - # TODO: further testing + # TODO: extend with future operations on inclusion mappings + + +@pytest.mark.parametrize( + "ghost_mode", [GhostMode.none, GhostMode.shared_vertex, GhostMode.shared_facet] +) +def test_2d(ghost_mode): + mesh = create_unit_square(MPI.COMM_WORLD, 5, 5, ghost_mode=ghost_mode) + mesh.topology.create_entities(1) + mesh_fine, _, _ = refine(mesh) + inclusion_mapping(mesh, mesh_fine) + # TODO: extend with future operations on inclusion mappings + + +@pytest.mark.parametrize( + "ghost_mode", [GhostMode.none, GhostMode.shared_vertex, GhostMode.shared_facet] +) +def test_3d(ghost_mode): + mesh = create_unit_cube(MPI.COMM_WORLD, 5, 5, 5, ghost_mode=ghost_mode) + mesh.topology.create_entities(1) + mesh_fine, _, _ = refine(mesh) + inclusion_mapping(mesh, mesh_fine) + # TODO: extend with future operations on inclusion mappings if __name__ == "__main__": From 6d5ffa409b99fcedebfaf81dfbf3e5081e2f024d Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Sat, 4 Jan 2025 10:35:28 +0100 Subject: [PATCH 17/44] Start doc --- cpp/doc/source/api.rst | 1 + cpp/doc/source/multigrid.rst | 5 +++++ 2 files changed, 6 insertions(+) create mode 100644 cpp/doc/source/multigrid.rst diff --git a/cpp/doc/source/api.rst b/cpp/doc/source/api.rst index f9b914c732b..8464594f7a2 100644 --- a/cpp/doc/source/api.rst +++ b/cpp/doc/source/api.rst @@ -14,6 +14,7 @@ generated documentation is `here `_. io la mesh + multigrid refinement * :ref:`genindex` diff --git a/cpp/doc/source/multigrid.rst b/cpp/doc/source/multigrid.rst new file mode 100644 index 00000000000..f0fe0e6ecb1 --- /dev/null +++ b/cpp/doc/source/multigrid.rst @@ -0,0 +1,5 @@ +Multigrid (``dolfinx::multigrid``) +==================================== + +.. doxygennamespace:: dolfinx::multigrid + :project: DOLFINx From 5732a1cbb868d5a4461932a22f0b6c1fe0943637 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Sat, 4 Jan 2025 10:56:24 +0100 Subject: [PATCH 18/44] More doc --- cpp/dolfinx/multigrid/dolfinx_multigrid.h | 11 ++++++++++ cpp/dolfinx/multigrid/inclusion.h | 26 ++++++++++++++++++++++- 2 files changed, 36 insertions(+), 1 deletion(-) create mode 100644 cpp/dolfinx/multigrid/dolfinx_multigrid.h diff --git a/cpp/dolfinx/multigrid/dolfinx_multigrid.h b/cpp/dolfinx/multigrid/dolfinx_multigrid.h new file mode 100644 index 00000000000..fb454826c89 --- /dev/null +++ b/cpp/dolfinx/multigrid/dolfinx_multigrid.h @@ -0,0 +1,11 @@ +#pragma once + +/// @brief Multigrid algorithms. +/// +namespace dolfinx::multigrid +{ +} + +// DOLFINx refinement interface + +#include diff --git a/cpp/dolfinx/multigrid/inclusion.h b/cpp/dolfinx/multigrid/inclusion.h index f4b17f8d06a..22f60499013 100644 --- a/cpp/dolfinx/multigrid/inclusion.h +++ b/cpp/dolfinx/multigrid/inclusion.h @@ -23,6 +23,15 @@ namespace dolfinx::multigrid { +/** + * @brief Gathers a global vector from combination of local data. + * @note Performs an all-to-all communication. + * + * @param local local data + * @param global_size number of global data entries + * @param comm MPI communicator + * @return std::vector on communicator gathered global data. + */ template std::vector gather_global(std::span local, std::int64_t global_size, MPI_Comm comm) @@ -49,6 +58,21 @@ std::vector gather_global(std::span local, std::int64_t global_size, return global; } +/** + * @brief Computes an inclusion map, i.e. local list of global vertex indices of + * another mesh, between to meshes. + * + * + * @param mesh_from Coarser mesh (domain of the inclusion map) + * @param mesh_to Finer mesh (range of the inclusion map) + * @param allow_all_to_all If the vertices of `mesh_from` are not equally + * spatially parallelized as `mesh_to` an all-to-all gathering of all vertices + * in `mesh_to` is performed. If true, performs all-to-all gathering, otherwise + * throws an exception if this becomes necessary. + * @return std::vector Map from local vertex index in `mesh_from` + * to global vertex index in `mesh_to`, i.e. `mesh_from.geometry.x()[i:i+3] == + * mesh_to.geometry.x()[map[i]:map[i]+3]`. + */ template std::vector inclusion_mapping(const dolfinx::mesh::Mesh& mesh_from, @@ -124,7 +148,7 @@ inclusion_mapping(const dolfinx::mesh::Mesh& mesh_from, "Parallelization of mesh requires all to all communication to compute " "inclusion map, but allow_all_to_all is not set."); } - + // Build global to vertex list std::vector global_x_to = gather_global(mesh_to.geometry().x().subspan(0, im_to.size_local() * 3), From c36c5d012df3a08328ed315c18b5567b8f1e6081 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Sat, 4 Jan 2025 10:58:28 +0100 Subject: [PATCH 19/44] Copy Typo --- cpp/dolfinx/multigrid/dolfinx_multigrid.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/dolfinx/multigrid/dolfinx_multigrid.h b/cpp/dolfinx/multigrid/dolfinx_multigrid.h index fb454826c89..562acb169ab 100644 --- a/cpp/dolfinx/multigrid/dolfinx_multigrid.h +++ b/cpp/dolfinx/multigrid/dolfinx_multigrid.h @@ -6,6 +6,6 @@ namespace dolfinx::multigrid { } -// DOLFINx refinement interface +// DOLFINx multigrid interface #include From 814a39bf938fb4eb0f55ba24bcadd95f04224fde Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Sat, 4 Jan 2025 11:00:12 +0100 Subject: [PATCH 20/44] Allow all to all in python tests --- python/test/unit/multigrid/test_inclusion.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/python/test/unit/multigrid/test_inclusion.py b/python/test/unit/multigrid/test_inclusion.py index 92c995de67e..06fee848579 100644 --- a/python/test/unit/multigrid/test_inclusion.py +++ b/python/test/unit/multigrid/test_inclusion.py @@ -18,7 +18,7 @@ def test_1d(ghost_mode): mesh = create_interval(MPI.COMM_WORLD, 10, (0, 1), ghost_mode=ghost_mode) mesh_fine, _, _ = refine(mesh) - inclusion_mapping(mesh, mesh_fine) + inclusion_mapping(mesh, mesh_fine, True) # TODO: extend with future operations on inclusion mappings @@ -29,7 +29,7 @@ def test_2d(ghost_mode): mesh = create_unit_square(MPI.COMM_WORLD, 5, 5, ghost_mode=ghost_mode) mesh.topology.create_entities(1) mesh_fine, _, _ = refine(mesh) - inclusion_mapping(mesh, mesh_fine) + inclusion_mapping(mesh, mesh_fine, True) # TODO: extend with future operations on inclusion mappings @@ -40,7 +40,7 @@ def test_3d(ghost_mode): mesh = create_unit_cube(MPI.COMM_WORLD, 5, 5, 5, ghost_mode=ghost_mode) mesh.topology.create_entities(1) mesh_fine, _, _ = refine(mesh) - inclusion_mapping(mesh, mesh_fine) + inclusion_mapping(mesh, mesh_fine, True) # TODO: extend with future operations on inclusion mappings From 48e8f83fddd51b71b8f58f72cc0490d6f2d1db06 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Sat, 4 Jan 2025 11:28:41 +0100 Subject: [PATCH 21/44] Tidy --- cpp/dolfinx/multigrid/inclusion.h | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/cpp/dolfinx/multigrid/inclusion.h b/cpp/dolfinx/multigrid/inclusion.h index 22f60499013..69f543e8e95 100644 --- a/cpp/dolfinx/multigrid/inclusion.h +++ b/cpp/dolfinx/multigrid/inclusion.h @@ -121,6 +121,7 @@ inclusion_mapping(const dolfinx::mesh::Mesh& mesh_from, vertex_from, vertex_to, [](T a, T b) { return std::abs(a - b) <= std::numeric_limits::epsilon(); })) { + assert(map[i] == -1); map[i] = to_global_to(j); break; } @@ -134,12 +135,10 @@ inclusion_mapping(const dolfinx::mesh::Mesh& mesh_from, return map; } - bool locally_complete - = std::ranges::all_of(map, [](auto e) { return e >= 0; }); - bool globally_complete = false; - MPI_Allreduce(&locally_complete, &locally_complete, 1, MPI_CXX_BOOL, MPI_LAND, + bool all_found = std::ranges::all_of(map, [](auto e) { return e >= 0; }); + MPI_Allreduce(MPI_IN_PLACE, &all_found, 1, MPI_CXX_BOOL, MPI_LAND, mesh_to.comm()); - if (globally_complete) + if (all_found) return map; if (!allow_all_to_all) From 6e4376b10bfa76089ff6512d32f146ab232cd07d Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Mon, 10 Mar 2025 18:06:21 +0100 Subject: [PATCH 22/44] Adapt to main --- python/dolfinx/wrappers/multigrid.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/dolfinx/wrappers/multigrid.cpp b/python/dolfinx/wrappers/multigrid.cpp index 619d8fec848..ee098d3b935 100644 --- a/python/dolfinx/wrappers/multigrid.cpp +++ b/python/dolfinx/wrappers/multigrid.cpp @@ -13,7 +13,7 @@ #include #include -#include "array.h" +#include "dolfinx_wrappers/array.h" namespace nb = nanobind; From d2adf740b41995c7958c63d8173d8f83d3e42054 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Mon, 10 Mar 2025 18:32:20 +0100 Subject: [PATCH 23/44] Addapt --- cpp/dolfinx/multigrid/CMakeLists.txt | 1 + cpp/dolfinx/multigrid/dolfinx_multigrid.h | 1 + 2 files changed, 2 insertions(+) diff --git a/cpp/dolfinx/multigrid/CMakeLists.txt b/cpp/dolfinx/multigrid/CMakeLists.txt index fa5d718cb74..55516dc1ef0 100644 --- a/cpp/dolfinx/multigrid/CMakeLists.txt +++ b/cpp/dolfinx/multigrid/CMakeLists.txt @@ -1,4 +1,5 @@ set(HEADERS_multigrid ${CMAKE_CURRENT_SOURCE_DIR}/inclusion.h + ${CMAKE_CURRENT_SOURCE_DIR}/transfer_matrix.h PARENT_SCOPE ) diff --git a/cpp/dolfinx/multigrid/dolfinx_multigrid.h b/cpp/dolfinx/multigrid/dolfinx_multigrid.h index 562acb169ab..5dff5c3502f 100644 --- a/cpp/dolfinx/multigrid/dolfinx_multigrid.h +++ b/cpp/dolfinx/multigrid/dolfinx_multigrid.h @@ -9,3 +9,4 @@ namespace dolfinx::multigrid // DOLFINx multigrid interface #include +#include From 2532d3b98057baa8f3af1df3ec8985001d2a5dc9 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Mon, 10 Mar 2025 19:38:14 +0100 Subject: [PATCH 24/44] Fix python module --- python/dolfinx/wrappers/multigrid.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/python/dolfinx/wrappers/multigrid.cpp b/python/dolfinx/wrappers/multigrid.cpp index a9f1b872a9c..cd655f394b1 100644 --- a/python/dolfinx/wrappers/multigrid.cpp +++ b/python/dolfinx/wrappers/multigrid.cpp @@ -10,8 +10,10 @@ #include #include +#include #include #include +#include #include "dolfinx_wrappers/array.h" @@ -32,7 +34,7 @@ void multigrid(nb::module_& m) allow_all_to_all); return dolfinx_wrappers::as_nbarray(std::move(map)); }, - nb::arg("mesh_from"), nb::arg("mesh_to"), + nb::arg("mesh_from"), nb::arg("mesh_to"), nb::arg("allow_all_to_all"), "Computes inclusion mapping between two meshes"); m.def( From a4f14b617a2c29fd536ea610f455db3621a5e17b Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Mon, 10 Mar 2025 20:18:06 +0100 Subject: [PATCH 25/44] Finish rename in python layer --- python/dolfinx/multigrid.py | 12 +++++++- python/dolfinx/transfer.py | 29 ------------------- .../{transfer => multigrid}/test_transfer.py | 2 +- 3 files changed, 12 insertions(+), 31 deletions(-) delete mode 100644 python/dolfinx/transfer.py rename python/test/unit/{transfer => multigrid}/test_transfer.py (96%) diff --git a/python/dolfinx/multigrid.py b/python/dolfinx/multigrid.py index 2fde4813f72..c50a3abf981 100644 --- a/python/dolfinx/multigrid.py +++ b/python/dolfinx/multigrid.py @@ -7,13 +7,23 @@ import numpy as np from numpy.typing import NDArray +from dolfinx.cpp.la import SparsityPattern +from dolfinx.cpp.multigrid import assemble_transfer_matrix +from dolfinx.cpp.multigrid import create_sparsity_pattern as _create_sparsity_pattern from dolfinx.cpp.multigrid import inclusion_mapping as _inclusion_mapping +from dolfinx.fem import FunctionSpace from dolfinx.mesh import Mesh -__all__ = ["inclusion_mapping"] +__all__ = ["assemble_transfer_matrix", "create_sparsity_pattern", "inclusion_mapping"] def inclusion_mapping( mesh_from: Mesh, mesh_to: Mesh, allow_all_to_all: bool = False ) -> NDArray[np.int64]: return _inclusion_mapping(mesh_from._cpp_object, mesh_to._cpp_object, allow_all_to_all) + + +def create_sparsity_pattern( + V_from: FunctionSpace, V_to: FunctionSpace, inclusion_map: NDArray[np.int64] +) -> SparsityPattern: + return _create_sparsity_pattern(V_from._cpp_object, V_to._cpp_object, inclusion_map) diff --git a/python/dolfinx/transfer.py b/python/dolfinx/transfer.py deleted file mode 100644 index 7b9e8b59e7c..00000000000 --- a/python/dolfinx/transfer.py +++ /dev/null @@ -1,29 +0,0 @@ -# Copyright (C) 2024 Paul T. Kühner -# -# This file is part of DOLFINx (https://www.fenicsproject.org) -# -# SPDX-License-Identifier: LGPL-3.0-or-later - -import numpy as np -from numpy.typing import NDArray - -from dolfinx.cpp.la import SparsityPattern -from dolfinx.cpp.transfer import ( - assemble_transfer_matrix, -) -from dolfinx.cpp.transfer import create_sparsity_pattern as _create_sparsity_pattern -from dolfinx.cpp.transfer import inclusion_mapping as _inclusion_mapping -from dolfinx.fem import FunctionSpace -from dolfinx.mesh import Mesh - -__all__ = ["assemble_transfer_matrix", "create_sparsity_pattern", "inclusion_mapping"] - - -def inclusion_mapping(mesh_from: Mesh, mesh_to: Mesh) -> NDArray[np.int64]: - return _inclusion_mapping(mesh_from._cpp_object, mesh_to._cpp_object) - - -def create_sparsity_pattern( - V_from: FunctionSpace, V_to: FunctionSpace, inclusion_map: NDArray[np.int64] -) -> SparsityPattern: - return _create_sparsity_pattern(V_from._cpp_object, V_to._cpp_object, inclusion_map) diff --git a/python/test/unit/transfer/test_transfer.py b/python/test/unit/multigrid/test_transfer.py similarity index 96% rename from python/test/unit/transfer/test_transfer.py rename to python/test/unit/multigrid/test_transfer.py index 8d60dfb12e2..a1caf7d9736 100644 --- a/python/test/unit/transfer/test_transfer.py +++ b/python/test/unit/multigrid/test_transfer.py @@ -10,7 +10,7 @@ from dolfinx.fem import functionspace from dolfinx.mesh import GhostMode, create_interval, refine -from dolfinx.transfer import create_sparsity_pattern, inclusion_mapping +from dolfinx.multigrid import create_sparsity_pattern, inclusion_mapping @pytest.mark.parametrize( From 77623d6e877b7e429feeabbc1100975d501a93c5 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Mon, 10 Mar 2025 20:33:27 +0100 Subject: [PATCH 26/44] Passing python test with assemble_transfer_matrix --- python/dolfinx/wrappers/multigrid.cpp | 1 + python/test/unit/multigrid/test_transfer.py | 21 +++++++++++---------- 2 files changed, 12 insertions(+), 10 deletions(-) diff --git a/python/dolfinx/wrappers/multigrid.cpp b/python/dolfinx/wrappers/multigrid.cpp index cd655f394b1..e0b22fa4118 100644 --- a/python/dolfinx/wrappers/multigrid.cpp +++ b/python/dolfinx/wrappers/multigrid.cpp @@ -8,6 +8,7 @@ #include #include +#include #include #include diff --git a/python/test/unit/multigrid/test_transfer.py b/python/test/unit/multigrid/test_transfer.py index a1caf7d9736..455e13944e1 100644 --- a/python/test/unit/multigrid/test_transfer.py +++ b/python/test/unit/multigrid/test_transfer.py @@ -9,8 +9,9 @@ import pytest from dolfinx.fem import functionspace +from dolfinx.la import matrix_csr from dolfinx.mesh import GhostMode, create_interval, refine -from dolfinx.multigrid import create_sparsity_pattern, inclusion_mapping +from dolfinx.multigrid import assemble_transfer_matrix, create_sparsity_pattern, inclusion_mapping @pytest.mark.parametrize( @@ -25,15 +26,15 @@ def test_1d(ghost_mode): assert V.element == V_fine.element V_fine.mesh.topology.create_connectivity(1, 0) V_fine.mesh.topology.create_connectivity(0, 1) - create_sparsity_pattern(V, V_fine, inclusion_map) - # T = matrix_csr(sp) - # assemble_transfer_matrix( - # T._cpp_object, - # V._cpp_object, - # V_fine._cpp_object, - # inclusion_map, - # lambda i: 1.0 if i == 0 else 0.5, - # ) + sp = create_sparsity_pattern(V, V_fine, inclusion_map) + T = matrix_csr(sp) + assemble_transfer_matrix( + T._cpp_object, + V._cpp_object, + V_fine._cpp_object, + inclusion_map, + lambda i: 1.0 if i == 0 else 0.5, + ) # continue with assembly of matrix From 6b791cabc3571ae681daf0a1fca4bda0b99dd515 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Tue, 11 Mar 2025 09:03:13 +0100 Subject: [PATCH 27/44] Python wrapping of assemble_transfer_matrix --- python/dolfinx/multigrid.py | 17 ++++++++++++++++- python/test/unit/multigrid/test_transfer.py | 12 +++--------- 2 files changed, 19 insertions(+), 10 deletions(-) diff --git a/python/dolfinx/multigrid.py b/python/dolfinx/multigrid.py index c50a3abf981..407999c3c10 100644 --- a/python/dolfinx/multigrid.py +++ b/python/dolfinx/multigrid.py @@ -4,14 +4,17 @@ # # SPDX-License-Identifier: LGPL-3.0-or-later +from typing import Callable + import numpy as np from numpy.typing import NDArray from dolfinx.cpp.la import SparsityPattern -from dolfinx.cpp.multigrid import assemble_transfer_matrix +from dolfinx.cpp.multigrid import assemble_transfer_matrix as _assemble_transfer_matrix from dolfinx.cpp.multigrid import create_sparsity_pattern as _create_sparsity_pattern from dolfinx.cpp.multigrid import inclusion_mapping as _inclusion_mapping from dolfinx.fem import FunctionSpace +from dolfinx.la import MatrixCSR from dolfinx.mesh import Mesh __all__ = ["assemble_transfer_matrix", "create_sparsity_pattern", "inclusion_mapping"] @@ -23,6 +26,18 @@ def inclusion_mapping( return _inclusion_mapping(mesh_from._cpp_object, mesh_to._cpp_object, allow_all_to_all) +def assemble_transfer_matrix( + T: MatrixCSR, + V_from: FunctionSpace, + V_to: FunctionSpace, + inclusion_map: NDArray[np.int64], + weight: Callable[[int], float], +): + return _assemble_transfer_matrix( + T._cpp_object, V_from._cpp_object, V_to._cpp_object, inclusion_map, weight + ) + + def create_sparsity_pattern( V_from: FunctionSpace, V_to: FunctionSpace, inclusion_map: NDArray[np.int64] ) -> SparsityPattern: diff --git a/python/test/unit/multigrid/test_transfer.py b/python/test/unit/multigrid/test_transfer.py index 455e13944e1..3bfbfcc66be 100644 --- a/python/test/unit/multigrid/test_transfer.py +++ b/python/test/unit/multigrid/test_transfer.py @@ -20,7 +20,7 @@ def test_1d(ghost_mode): mesh = create_interval(MPI.COMM_WORLD, 10, (0, 1), ghost_mode=ghost_mode) mesh_fine, _, _ = refine(mesh) - inclusion_map = inclusion_mapping(mesh, mesh_fine) + inclusion_map = inclusion_mapping(mesh, mesh_fine, True) V = functionspace(mesh, ("Lagrange", 1, (1,))) V_fine = functionspace(mesh_fine, ("Lagrange", 1, (1,))) assert V.element == V_fine.element @@ -28,13 +28,7 @@ def test_1d(ghost_mode): V_fine.mesh.topology.create_connectivity(0, 1) sp = create_sparsity_pattern(V, V_fine, inclusion_map) T = matrix_csr(sp) - assemble_transfer_matrix( - T._cpp_object, - V._cpp_object, - V_fine._cpp_object, - inclusion_map, - lambda i: 1.0 if i == 0 else 0.5, - ) + assemble_transfer_matrix(T, V, V_fine, inclusion_map, lambda i: 1.0 if i == 0 else 0.5) # continue with assembly of matrix @@ -47,7 +41,7 @@ def test_1d_higher_order(degree, ghost_mode): mesh_fine, _, _ = refine(mesh) # this is a strictly geometric operation, and thus should pass - inclusion_map = inclusion_mapping(mesh, mesh_fine) + inclusion_map = inclusion_mapping(mesh, mesh_fine, True) V = functionspace(mesh, ("Lagrange", degree, (1,))) V_fine = functionspace(mesh_fine, ("Lagrange", degree, (1,))) From e13cdcec3d8bb7f2aba2c2033e40eeaa238df4ef Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Tue, 11 Mar 2025 10:13:21 +0100 Subject: [PATCH 28/44] Remove return value --- python/dolfinx/multigrid.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/dolfinx/multigrid.py b/python/dolfinx/multigrid.py index 407999c3c10..b045466e427 100644 --- a/python/dolfinx/multigrid.py +++ b/python/dolfinx/multigrid.py @@ -33,7 +33,7 @@ def assemble_transfer_matrix( inclusion_map: NDArray[np.int64], weight: Callable[[int], float], ): - return _assemble_transfer_matrix( + _assemble_transfer_matrix( T._cpp_object, V_from._cpp_object, V_to._cpp_object, inclusion_map, weight ) From 22f747152393081fb3e6a44b4914ac579225c221 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Tue, 11 Mar 2025 13:34:09 +0100 Subject: [PATCH 29/44] Sketch out strategy and copy maintain_coarse_partitioner (for non ghosted case) --- cpp/dolfinx/refinement/refine.h | 30 ++++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h index 430bcc495bc..618204022b9 100644 --- a/cpp/dolfinx/refinement/refine.h +++ b/cpp/dolfinx/refinement/refine.h @@ -84,4 +84,34 @@ refine(const mesh::Mesh& mesh, return {std::move(mesh1), std::move(parent_cell), std::move(parent_facet)}; } + +template +mesh::CellPartitionFunction +create_identity_partitioner(const mesh::Mesh& parent_mesh, + std::vector parent_cell) +{ + // for every non ghost index maintain_coarse_partitioner below should be doing + // the correct thing. ghost cells need to given the rank, which previously the + // COARSE cell was assigned. These ranks can be computed for parent_mesh with + // compute_destination_ranks and the parent_cell[i] is the cell idx in + // parent_mesh to access it. + + return + [&](MPI_Comm comm, int /*nparts*/, std::vector cell_types, + std::vector> cells) + -> graph::AdjacencyList + { + int num_cell_vertices = mesh::num_cell_vertices(cell_types.front()); + std::int32_t num_cells = cells.front().size() / num_cell_vertices; + std::vector destinations(num_cells, dolfinx::MPI::rank(comm)); + // TODO: iterate over previous ghost cells here and replace destination + // ranks with previous coarse ranks, i.e. do not change ownership of ghosts + // cells to local process + std::vector dest_offsets(num_cells + 1); + std::iota(dest_offsets.begin(), dest_offsets.end(), 0); + return graph::AdjacencyList(std::move(destinations), + std::move(dest_offsets)); + }; +} + } // namespace dolfinx::refinement From 6160a96a61231f32acb1d9922d6740c531237ce5 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Tue, 11 Mar 2025 14:53:34 +0100 Subject: [PATCH 30/44] A first version --- cpp/dolfinx/refinement/refine.h | 229 +++++++++++++++++++++++++++++--- 1 file changed, 214 insertions(+), 15 deletions(-) diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h index 618204022b9..bab24157164 100644 --- a/cpp/dolfinx/refinement/refine.h +++ b/cpp/dolfinx/refinement/refine.h @@ -6,12 +6,14 @@ #pragma once +#include "common/MPI.h" #include "dolfinx/graph/AdjacencyList.h" #include "dolfinx/mesh/Mesh.h" #include "dolfinx/mesh/Topology.h" #include "dolfinx/mesh/cell_types.h" #include "dolfinx/mesh/utils.h" #include "interval.h" +#include "mesh/graphbuild.h" #include "plaza.h" #include #include @@ -19,6 +21,176 @@ #include #include +// TODO: Remove once works +namespace +{ +/// @todo Is it un-documented that the owning rank must come first in +/// reach list of edges? +/// +/// @param[in] comm The communicator +/// @param[in] graph Graph, using global indices for graph edges +/// @param[in] node_disp The distribution of graph nodes across MPI +/// ranks. The global index `gidx` of local index `lidx` is `lidx + +/// node_disp[my_rank]`. +/// @param[in] part The destination rank for owned nodes, i.e. `dest[i]` +/// is the destination of the node with local index `i`. +/// @return Destination ranks for each local node. +template +graph::AdjacencyList compute_destination_ranks( + MPI_Comm comm, const graph::AdjacencyList& graph, + const std::vector& node_disp, const std::vector& part) +{ + common::Timer timer("Extend graph destination ranks for halo"); + + const int rank = dolfinx::MPI::rank(comm); + const std::int64_t range0 = node_disp[rank]; + const std::int64_t range1 = node_disp[rank + 1]; + assert(static_cast(range1 - range0) == graph.num_nodes()); + + // Wherever an owned 'node' goes, so must the nodes connected to it by + // an edge ('node1'). Task is to let the owner of node1 know the extra + // ranks that it needs to send node1 to. + std::vector> node_to_dest; + for (int node0 = 0; node0 < graph.num_nodes(); ++node0) + { + // Wherever 'node' goes to, so must the attached 'node1' + for (auto node1 : graph.links(node0)) + { + if (node1 < range0 or node1 >= range1) + { + auto it = std::ranges::upper_bound(node_disp, node1); + int remote_rank = std::distance(node_disp.begin(), it) - 1; + node_to_dest.push_back( + {remote_rank, node1, static_cast(part[node0])}); + } + else + node_to_dest.push_back( + {rank, node1, static_cast(part[node0])}); + } + } + + std::ranges::sort(node_to_dest); + auto [unique_end, range_end] = std::ranges::unique(node_to_dest); + node_to_dest.erase(unique_end, range_end); + + // Build send data and buffer + std::vector dest, send_sizes; + std::vector send_buffer; + { + auto it = node_to_dest.begin(); + while (it != node_to_dest.end()) + { + // Current destination rank + dest.push_back(it->front()); + + // Find iterator to next destination rank and pack send data + auto it1 + = std::find_if(it, node_to_dest.end(), [r0 = dest.back()](auto& idx) + { return idx[0] != r0; }); + send_sizes.push_back(2 * std::distance(it, it1)); + for (auto itx = it; itx != it1; ++itx) + { + send_buffer.push_back(itx->at(1)); + send_buffer.push_back(itx->at(2)); + } + + it = it1; + } + } + + // Prepare send displacements + std::vector send_disp(send_sizes.size() + 1, 0); + std::partial_sum(send_sizes.begin(), send_sizes.end(), + std::next(send_disp.begin())); + + // Discover src ranks. ParMETIS/KaHIP are not scalable (holding an + // array of size equal to the comm size), so no extra harm in using + // non-scalable neighbourhood detection (which might be faster for + // small rank counts). + const std::vector src + = dolfinx::MPI::compute_graph_edges_pcx(comm, dest); + + // Create neighbourhood communicator + MPI_Comm neigh_comm; + MPI_Dist_graph_create_adjacent(comm, src.size(), src.data(), MPI_UNWEIGHTED, + dest.size(), dest.data(), MPI_UNWEIGHTED, + MPI_INFO_NULL, false, &neigh_comm); + + // Determine receives sizes + std::vector recv_sizes(dest.size()); + send_sizes.reserve(1); + recv_sizes.reserve(1); + MPI_Neighbor_alltoall(send_sizes.data(), 1, MPI_INT, recv_sizes.data(), 1, + MPI_INT, neigh_comm); + + // Prepare receive displacements + std::vector recv_disp(recv_sizes.size() + 1, 0); + std::partial_sum(recv_sizes.begin(), recv_sizes.end(), + std::next(recv_disp.begin())); + + // Send/receive data + std::vector recv_buffer(recv_disp.back()); + MPI_Neighbor_alltoallv(send_buffer.data(), send_sizes.data(), + send_disp.data(), MPI_INT64_T, recv_buffer.data(), + recv_sizes.data(), recv_disp.data(), MPI_INT64_T, + neigh_comm); + MPI_Comm_free(&neigh_comm); + + // Prepare (local node index, destination rank) array. Add local data, + // then add the received data, and the make unique. + std::vector> local_node_to_dest; + for (auto d : part) + { + local_node_to_dest.push_back( + {static_cast(local_node_to_dest.size()), static_cast(d)}); + } + for (std::size_t i = 0; i < recv_buffer.size(); i += 2) + { + std::int64_t idx = recv_buffer[i]; + int d = recv_buffer[i + 1]; + assert(idx >= range0 and idx < range1); + std::int32_t idx_local = idx - range0; + local_node_to_dest.push_back({idx_local, d}); + } + + { + std::ranges::sort(local_node_to_dest); + auto [unique_end, range_end] = std::ranges::unique(local_node_to_dest); + local_node_to_dest.erase(unique_end, range_end); + } + // Compute offsets + std::vector offsets(graph.num_nodes() + 1, 0); + { + std::vector num_dests(graph.num_nodes(), 0); + for (auto x : local_node_to_dest) + ++num_dests[x[0]]; + std::partial_sum(num_dests.begin(), num_dests.end(), + std::next(offsets.begin())); + } + + // Fill data array + std::vector data(offsets.back()); + { + std::vector pos = offsets; + for (auto [x0, x1] : local_node_to_dest) + data[pos[x0]++] = x1; + } + + graph::AdjacencyList g(std::move(data), std::move(offsets)); + + // Make sure the owning rank comes first for each node + for (std::int32_t i = 0; i < g.num_nodes(); ++i) + { + auto d = g.links(i); + auto it = std::find(d.begin(), d.end(), part[i]); + assert(it != d.end()); + std::iter_swap(d.begin(), it); + } + + return g; +} +} // namespace + namespace dolfinx::refinement { /// @brief Refine a mesh with markers. @@ -56,8 +228,7 @@ std::tuple, std::optional>, std::optional>> refine(const mesh::Mesh& mesh, std::optional> edges, - const mesh::CellPartitionFunction& partitioner - = mesh::create_cell_partitioner(mesh::GhostMode::none), + const mesh::CellPartitionFunction& partitioner = nullptr, Option option = Option::none) { auto topology = mesh.topology(); @@ -70,6 +241,11 @@ refine(const mesh::Mesh& mesh, ? interval::compute_refinement_data(mesh, edges, option) : plaza::compute_refinement_data(mesh, edges, option); + if (!partitioner) + { + partitioner = create_identity_partitioner(mesh, parent_cell); + } + mesh::Mesh mesh1 = mesh::create_mesh( mesh.comm(), mesh.comm(), cell_adj.array(), mesh.geometry().cmap(), mesh.comm(), new_vertex_coords, xshape, partitioner); @@ -88,29 +264,52 @@ refine(const mesh::Mesh& mesh, template mesh::CellPartitionFunction create_identity_partitioner(const mesh::Mesh& parent_mesh, - std::vector parent_cell) + std::span parent_cell) { - // for every non ghost index maintain_coarse_partitioner below should be doing - // the correct thing. ghost cells need to given the rank, which previously the - // COARSE cell was assigned. These ranks can be computed for parent_mesh with - // compute_destination_ranks and the parent_cell[i] is the cell idx in - // parent_mesh to access it. + // TODO: optimize for non ghosted mesh? return [&](MPI_Comm comm, int /*nparts*/, std::vector cell_types, std::vector> cells) -> graph::AdjacencyList { + auto parent_top = parent_mesh.topology(); + auto parent_cell_im = parent_top.index_map(parent_top.dim()); + int num_cell_vertices = mesh::num_cell_vertices(cell_types.front()); std::int32_t num_cells = cells.front().size() / num_cell_vertices; - std::vector destinations(num_cells, dolfinx::MPI::rank(comm)); - // TODO: iterate over previous ghost cells here and replace destination - // ranks with previous coarse ranks, i.e. do not change ownership of ghosts - // cells to local process + std::vector destinations(num_cells); + std::vector dest_offsets(num_cells + 1); - std::iota(dest_offsets.begin(), dest_offsets.end(), 0); - return graph::AdjacencyList(std::move(destinations), - std::move(dest_offsets)); + int rank = dolfinx::MPI::rank(comm); + for (std::int32_t i = 0; i < destinations.size(); i++) + { + bool ghost_parent_cell = parent_cell[i] > parent_cell_im.size_local(); + if (ghost_parent_cell) + { + destinations[i] + = parent_cell_im + .owners()[parent_cell[i] - parent_cell_im.size_local()]; + } + else + { + destinations[i] = rank; + } + } + + auto dual_graph = mesh::build_dual_graph(comm, cell_types, cells); + std::vector node_disp; + node_disp = std::vector(MPI::size(comm) + 1, 0); + MPI_Allgather(&cells.front().size(), 1, dolfinx::MPI::mpi_t, + node_disp.data() + 1, 1, dolfinx::MPI::mpi_t, + comm); + std::partial_sum(node_disp.begin(), node_disp.end(), node_disp.begin()); + + return compute_destination_ranks(comm, dual_graph, node_disp, destinations); + + // std::iota(dest_offsets.begin(), dest_offsets.end(), 0); + // return graph::AdjacencyList(std::move(destinations), + // std::move(dest_offsets)); }; } From c7f6e321f05490014c4977df64631db813e601c1 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Tue, 11 Mar 2025 16:49:00 +0100 Subject: [PATCH 31/44] Compiling... --- cpp/dolfinx/refinement/refine.h | 116 +++++++++++++++++--------------- 1 file changed, 60 insertions(+), 56 deletions(-) diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h index bab24157164..f982c76726b 100644 --- a/cpp/dolfinx/refinement/refine.h +++ b/cpp/dolfinx/refinement/refine.h @@ -6,14 +6,14 @@ #pragma once -#include "common/MPI.h" +#include "dolfinx/common/MPI.h" #include "dolfinx/graph/AdjacencyList.h" #include "dolfinx/mesh/Mesh.h" #include "dolfinx/mesh/Topology.h" #include "dolfinx/mesh/cell_types.h" +#include "dolfinx/mesh/graphbuild.h" #include "dolfinx/mesh/utils.h" #include "interval.h" -#include "mesh/graphbuild.h" #include "plaza.h" #include #include @@ -22,6 +22,7 @@ #include // TODO: Remove once works +namespace dolfinx { namespace { /// @todo Is it un-documented that the owning rank must come first in @@ -190,9 +191,63 @@ graph::AdjacencyList compute_destination_ranks( return g; } } // namespace +} namespace dolfinx::refinement { + +template +mesh::CellPartitionFunction +create_identity_partitioner(const mesh::Mesh& parent_mesh, + std::span parent_cell) +{ + // TODO: optimize for non ghosted mesh? + + return + [&](MPI_Comm comm, int /*nparts*/, std::vector cell_types, + std::vector> cells) + -> graph::AdjacencyList + { + auto parent_top = parent_mesh.topology(); + auto parent_cell_im = parent_top->index_map(parent_top->dim()); + + int num_cell_vertices = mesh::num_cell_vertices(cell_types.front()); + std::int32_t num_cells = cells.front().size() / num_cell_vertices; + std::vector destinations(num_cells); + + std::vector dest_offsets(num_cells + 1); + int rank = dolfinx::MPI::rank(comm); + for (std::int32_t i = 0; i < destinations.size(); i++) + { + bool ghost_parent_cell = parent_cell[i] > parent_cell_im->size_local(); + if (ghost_parent_cell) + { + destinations[i] + = parent_cell_im->owners()[parent_cell[i] - parent_cell_im->size_local()]; + } + else + { + destinations[i] = rank; + } + } + + auto dual_graph = mesh::build_dual_graph(comm, cell_types, cells); + std::vector node_disp; + node_disp = std::vector(MPI::size(comm) + 1, 0); + std::int32_t local_size = cells.front().size(); + MPI_Allgather(&local_size, 1, dolfinx::MPI::mpi_t, + node_disp.data() + 1, 1, dolfinx::MPI::mpi_t, + comm); + std::partial_sum(node_disp.begin(), node_disp.end(), node_disp.begin()); + + return compute_destination_ranks(comm, dual_graph, node_disp, destinations); + + // std::iota(dest_offsets.begin(), dest_offsets.end(), 0); + // return graph::AdjacencyList(std::move(destinations), + // std::move(dest_offsets)); + }; +} + /// @brief Refine a mesh with markers. /// /// The refined mesh can be optionally re-partitioned across processes. @@ -228,7 +283,7 @@ std::tuple, std::optional>, std::optional>> refine(const mesh::Mesh& mesh, std::optional> edges, - const mesh::CellPartitionFunction& partitioner = nullptr, + mesh::CellPartitionFunction partitioner = nullptr, Option option = Option::none) { auto topology = mesh.topology(); @@ -243,7 +298,8 @@ refine(const mesh::Mesh& mesh, if (!partitioner) { - partitioner = create_identity_partitioner(mesh, parent_cell); + assert(parent_cell); + partitioner = create_identity_partitioner(mesh, parent_cell.value()); } mesh::Mesh mesh1 = mesh::create_mesh( @@ -261,56 +317,4 @@ refine(const mesh::Mesh& mesh, return {std::move(mesh1), std::move(parent_cell), std::move(parent_facet)}; } -template -mesh::CellPartitionFunction -create_identity_partitioner(const mesh::Mesh& parent_mesh, - std::span parent_cell) -{ - // TODO: optimize for non ghosted mesh? - - return - [&](MPI_Comm comm, int /*nparts*/, std::vector cell_types, - std::vector> cells) - -> graph::AdjacencyList - { - auto parent_top = parent_mesh.topology(); - auto parent_cell_im = parent_top.index_map(parent_top.dim()); - - int num_cell_vertices = mesh::num_cell_vertices(cell_types.front()); - std::int32_t num_cells = cells.front().size() / num_cell_vertices; - std::vector destinations(num_cells); - - std::vector dest_offsets(num_cells + 1); - int rank = dolfinx::MPI::rank(comm); - for (std::int32_t i = 0; i < destinations.size(); i++) - { - bool ghost_parent_cell = parent_cell[i] > parent_cell_im.size_local(); - if (ghost_parent_cell) - { - destinations[i] - = parent_cell_im - .owners()[parent_cell[i] - parent_cell_im.size_local()]; - } - else - { - destinations[i] = rank; - } - } - - auto dual_graph = mesh::build_dual_graph(comm, cell_types, cells); - std::vector node_disp; - node_disp = std::vector(MPI::size(comm) + 1, 0); - MPI_Allgather(&cells.front().size(), 1, dolfinx::MPI::mpi_t, - node_disp.data() + 1, 1, dolfinx::MPI::mpi_t, - comm); - std::partial_sum(node_disp.begin(), node_disp.end(), node_disp.begin()); - - return compute_destination_ranks(comm, dual_graph, node_disp, destinations); - - // std::iota(dest_offsets.begin(), dest_offsets.end(), 0); - // return graph::AdjacencyList(std::move(destinations), - // std::move(dest_offsets)); - }; -} - } // namespace dolfinx::refinement From c4afbcb059c3adb8349dcf0562037f7f5bbef6da Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Tue, 11 Mar 2025 21:58:53 +0100 Subject: [PATCH 32/44] Fix: python layer default value of partitioner does not align with cpp layer --- python/dolfinx/mesh.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py index 46385f0d345..cc11abd7db5 100644 --- a/python/dolfinx/mesh.py +++ b/python/dolfinx/mesh.py @@ -518,7 +518,7 @@ def transfer_meshtag( def refine( msh: Mesh, edges: typing.Optional[np.ndarray] = None, - partitioner: typing.Optional[typing.Callable] = create_cell_partitioner(GhostMode.none), + partitioner: typing.Optional[typing.Callable] = None, option: RefinementOption = RefinementOption.none, ) -> tuple[Mesh, npt.NDArray[np.int32], npt.NDArray[np.int8]]: """Refine a mesh. From d8650569577f49de782e933824bc2943b425fe01 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Tue, 11 Mar 2025 22:10:26 +0100 Subject: [PATCH 33/44] Debug --- cpp/dolfinx/refinement/refine.h | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h index f982c76726b..22e9ed49c5d 100644 --- a/cpp/dolfinx/refinement/refine.h +++ b/cpp/dolfinx/refinement/refine.h @@ -22,7 +22,8 @@ #include // TODO: Remove once works -namespace dolfinx { +namespace dolfinx +{ namespace { /// @todo Is it un-documented that the owning rank must come first in @@ -191,7 +192,7 @@ graph::AdjacencyList compute_destination_ranks( return g; } } // namespace -} +} // namespace dolfinx namespace dolfinx::refinement { @@ -203,10 +204,11 @@ create_identity_partitioner(const mesh::Mesh& parent_mesh, { // TODO: optimize for non ghosted mesh? - return - [&](MPI_Comm comm, int /*nparts*/, std::vector cell_types, - std::vector> cells) - -> graph::AdjacencyList + return [&parent_mesh, + parent_cell](MPI_Comm comm, int /*nparts*/, + std::vector cell_types, + std::vector> cells) + -> graph::AdjacencyList { auto parent_top = parent_mesh.topology(); auto parent_cell_im = parent_top->index_map(parent_top->dim()); @@ -215,7 +217,6 @@ create_identity_partitioner(const mesh::Mesh& parent_mesh, std::int32_t num_cells = cells.front().size() / num_cell_vertices; std::vector destinations(num_cells); - std::vector dest_offsets(num_cells + 1); int rank = dolfinx::MPI::rank(comm); for (std::int32_t i = 0; i < destinations.size(); i++) { @@ -223,7 +224,8 @@ create_identity_partitioner(const mesh::Mesh& parent_mesh, if (ghost_parent_cell) { destinations[i] - = parent_cell_im->owners()[parent_cell[i] - parent_cell_im->size_local()]; + = parent_cell_im + ->owners()[parent_cell[i] - parent_cell_im->size_local()]; } else { @@ -234,7 +236,7 @@ create_identity_partitioner(const mesh::Mesh& parent_mesh, auto dual_graph = mesh::build_dual_graph(comm, cell_types, cells); std::vector node_disp; node_disp = std::vector(MPI::size(comm) + 1, 0); - std::int32_t local_size = cells.front().size(); + std::int32_t local_size = dual_graph.num_nodes(); MPI_Allgather(&local_size, 1, dolfinx::MPI::mpi_t, node_disp.data() + 1, 1, dolfinx::MPI::mpi_t, comm); From 9b8314b34317b5d9ed26abc111127643e2b7c7ea Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Tue, 11 Mar 2025 22:28:42 +0100 Subject: [PATCH 34/44] Move compute_destination_ranks out of anonymous namespace --- cpp/dolfinx/graph/partitioners.cpp | 10 +- cpp/dolfinx/graph/partitioners.h | 6 + cpp/dolfinx/refinement/refine.h | 174 +---------------------------- 3 files changed, 11 insertions(+), 179 deletions(-) diff --git a/cpp/dolfinx/graph/partitioners.cpp b/cpp/dolfinx/graph/partitioners.cpp index fade5e160cf..e86b6ad6f2b 100644 --- a/cpp/dolfinx/graph/partitioners.cpp +++ b/cpp/dolfinx/graph/partitioners.cpp @@ -35,8 +35,6 @@ extern "C" using namespace dolfinx; -namespace -{ /// @todo Is it un-documented that the owning rank must come first in /// reach list of edges? /// @@ -49,7 +47,7 @@ namespace /// is the destination of the node with local index `i`. /// @return Destination ranks for each local node. template -graph::AdjacencyList compute_destination_ranks( +graph::AdjacencyList dolfinx::graph::compute_destination_ranks( MPI_Comm comm, const graph::AdjacencyList& graph, const std::vector& node_disp, const std::vector& part) { @@ -202,6 +200,7 @@ graph::AdjacencyList compute_destination_ranks( return g; } + //----------------------------------------------------------------------------- #ifdef HAS_PARMETIS template @@ -304,7 +303,6 @@ std::vector refine(MPI_Comm comm, const graph::AdjacencyList& adj_graph) //----------------------------------------------------------------------------- } #endif -} // namespace //----------------------------------------------------------------------------- #ifdef HAS_PTSCOTCH @@ -587,7 +585,7 @@ graph::partition_fn graph::parmetis::partitioner(double imbalance, { // FIXME: Is it implicit that the first entry is the owner? graph::AdjacencyList dest - = compute_destination_ranks(pcomm, graph, node_disp, part); + = graph::compute_destination_ranks(pcomm, graph, node_disp, part); if (split_comm) MPI_Comm_free(&pcomm); return dest; @@ -653,7 +651,7 @@ graph::partition_fn graph::kahip::partitioner(int mode, int seed, timer2.stop(); if (ghosting) - return compute_destination_ranks(comm, graph, node_disp, part); + return graph::compute_destination_ranks(comm, graph, node_disp, part); else { return regular_adjacency_list(std::vector(part.begin(), part.end()), diff --git a/cpp/dolfinx/graph/partitioners.h b/cpp/dolfinx/graph/partitioners.h index aba51a8430a..29aa62af97f 100644 --- a/cpp/dolfinx/graph/partitioners.h +++ b/cpp/dolfinx/graph/partitioners.h @@ -11,6 +11,12 @@ namespace dolfinx::graph { + +template +graph::AdjacencyList compute_destination_ranks( + MPI_Comm comm, const graph::AdjacencyList& graph, + const std::vector& node_disp, const std::vector& part); + namespace scotch { #ifdef HAS_PTSCOTCH diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h index 22e9ed49c5d..39cd82045a1 100644 --- a/cpp/dolfinx/refinement/refine.h +++ b/cpp/dolfinx/refinement/refine.h @@ -8,6 +8,7 @@ #include "dolfinx/common/MPI.h" #include "dolfinx/graph/AdjacencyList.h" +#include "dolfinx/graph/partitioners.h" #include "dolfinx/mesh/Mesh.h" #include "dolfinx/mesh/Topology.h" #include "dolfinx/mesh/cell_types.h" @@ -21,179 +22,6 @@ #include #include -// TODO: Remove once works -namespace dolfinx -{ -namespace -{ -/// @todo Is it un-documented that the owning rank must come first in -/// reach list of edges? -/// -/// @param[in] comm The communicator -/// @param[in] graph Graph, using global indices for graph edges -/// @param[in] node_disp The distribution of graph nodes across MPI -/// ranks. The global index `gidx` of local index `lidx` is `lidx + -/// node_disp[my_rank]`. -/// @param[in] part The destination rank for owned nodes, i.e. `dest[i]` -/// is the destination of the node with local index `i`. -/// @return Destination ranks for each local node. -template -graph::AdjacencyList compute_destination_ranks( - MPI_Comm comm, const graph::AdjacencyList& graph, - const std::vector& node_disp, const std::vector& part) -{ - common::Timer timer("Extend graph destination ranks for halo"); - - const int rank = dolfinx::MPI::rank(comm); - const std::int64_t range0 = node_disp[rank]; - const std::int64_t range1 = node_disp[rank + 1]; - assert(static_cast(range1 - range0) == graph.num_nodes()); - - // Wherever an owned 'node' goes, so must the nodes connected to it by - // an edge ('node1'). Task is to let the owner of node1 know the extra - // ranks that it needs to send node1 to. - std::vector> node_to_dest; - for (int node0 = 0; node0 < graph.num_nodes(); ++node0) - { - // Wherever 'node' goes to, so must the attached 'node1' - for (auto node1 : graph.links(node0)) - { - if (node1 < range0 or node1 >= range1) - { - auto it = std::ranges::upper_bound(node_disp, node1); - int remote_rank = std::distance(node_disp.begin(), it) - 1; - node_to_dest.push_back( - {remote_rank, node1, static_cast(part[node0])}); - } - else - node_to_dest.push_back( - {rank, node1, static_cast(part[node0])}); - } - } - - std::ranges::sort(node_to_dest); - auto [unique_end, range_end] = std::ranges::unique(node_to_dest); - node_to_dest.erase(unique_end, range_end); - - // Build send data and buffer - std::vector dest, send_sizes; - std::vector send_buffer; - { - auto it = node_to_dest.begin(); - while (it != node_to_dest.end()) - { - // Current destination rank - dest.push_back(it->front()); - - // Find iterator to next destination rank and pack send data - auto it1 - = std::find_if(it, node_to_dest.end(), [r0 = dest.back()](auto& idx) - { return idx[0] != r0; }); - send_sizes.push_back(2 * std::distance(it, it1)); - for (auto itx = it; itx != it1; ++itx) - { - send_buffer.push_back(itx->at(1)); - send_buffer.push_back(itx->at(2)); - } - - it = it1; - } - } - - // Prepare send displacements - std::vector send_disp(send_sizes.size() + 1, 0); - std::partial_sum(send_sizes.begin(), send_sizes.end(), - std::next(send_disp.begin())); - - // Discover src ranks. ParMETIS/KaHIP are not scalable (holding an - // array of size equal to the comm size), so no extra harm in using - // non-scalable neighbourhood detection (which might be faster for - // small rank counts). - const std::vector src - = dolfinx::MPI::compute_graph_edges_pcx(comm, dest); - - // Create neighbourhood communicator - MPI_Comm neigh_comm; - MPI_Dist_graph_create_adjacent(comm, src.size(), src.data(), MPI_UNWEIGHTED, - dest.size(), dest.data(), MPI_UNWEIGHTED, - MPI_INFO_NULL, false, &neigh_comm); - - // Determine receives sizes - std::vector recv_sizes(dest.size()); - send_sizes.reserve(1); - recv_sizes.reserve(1); - MPI_Neighbor_alltoall(send_sizes.data(), 1, MPI_INT, recv_sizes.data(), 1, - MPI_INT, neigh_comm); - - // Prepare receive displacements - std::vector recv_disp(recv_sizes.size() + 1, 0); - std::partial_sum(recv_sizes.begin(), recv_sizes.end(), - std::next(recv_disp.begin())); - - // Send/receive data - std::vector recv_buffer(recv_disp.back()); - MPI_Neighbor_alltoallv(send_buffer.data(), send_sizes.data(), - send_disp.data(), MPI_INT64_T, recv_buffer.data(), - recv_sizes.data(), recv_disp.data(), MPI_INT64_T, - neigh_comm); - MPI_Comm_free(&neigh_comm); - - // Prepare (local node index, destination rank) array. Add local data, - // then add the received data, and the make unique. - std::vector> local_node_to_dest; - for (auto d : part) - { - local_node_to_dest.push_back( - {static_cast(local_node_to_dest.size()), static_cast(d)}); - } - for (std::size_t i = 0; i < recv_buffer.size(); i += 2) - { - std::int64_t idx = recv_buffer[i]; - int d = recv_buffer[i + 1]; - assert(idx >= range0 and idx < range1); - std::int32_t idx_local = idx - range0; - local_node_to_dest.push_back({idx_local, d}); - } - - { - std::ranges::sort(local_node_to_dest); - auto [unique_end, range_end] = std::ranges::unique(local_node_to_dest); - local_node_to_dest.erase(unique_end, range_end); - } - // Compute offsets - std::vector offsets(graph.num_nodes() + 1, 0); - { - std::vector num_dests(graph.num_nodes(), 0); - for (auto x : local_node_to_dest) - ++num_dests[x[0]]; - std::partial_sum(num_dests.begin(), num_dests.end(), - std::next(offsets.begin())); - } - - // Fill data array - std::vector data(offsets.back()); - { - std::vector pos = offsets; - for (auto [x0, x1] : local_node_to_dest) - data[pos[x0]++] = x1; - } - - graph::AdjacencyList g(std::move(data), std::move(offsets)); - - // Make sure the owning rank comes first for each node - for (std::int32_t i = 0; i < g.num_nodes(); ++i) - { - auto d = g.links(i); - auto it = std::find(d.begin(), d.end(), part[i]); - assert(it != d.end()); - std::iter_swap(d.begin(), it); - } - - return g; -} -} // namespace -} // namespace dolfinx - namespace dolfinx::refinement { From c4766dcad705d0e12a97b8f56e44c87262ee5f00 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Wed, 12 Mar 2025 09:08:16 +0100 Subject: [PATCH 35/44] Add docstring --- cpp/dolfinx/refinement/refine.h | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h index 39cd82045a1..077d5bd5a1c 100644 --- a/cpp/dolfinx/refinement/refine.h +++ b/cpp/dolfinx/refinement/refine.h @@ -25,6 +25,17 @@ namespace dolfinx::refinement { +/** + * @brief Create a cell partitioner which maintains the partition of a coarse + * mesh. + * + * @tparam T floating point type of mesh geometry. + * @param parent_mesh mesh indicating the partition scheme to use, i.e. the + * coarse mesh. + * @param parent_cell list of cell indices mapping cells of the new refined mesh + * into the coarse mesh, usually output of `refinement::refine`. + * @return The created cell partition function. + */ template mesh::CellPartitionFunction create_identity_partitioner(const mesh::Mesh& parent_mesh, @@ -71,10 +82,6 @@ create_identity_partitioner(const mesh::Mesh& parent_mesh, std::partial_sum(node_disp.begin(), node_disp.end(), node_disp.begin()); return compute_destination_ranks(comm, dual_graph, node_disp, destinations); - - // std::iota(dest_offsets.begin(), dest_offsets.end(), 0); - // return graph::AdjacencyList(std::move(destinations), - // std::move(dest_offsets)); }; } From bb79adacd9a5f923e578546792672691de12b0e8 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Wed, 12 Mar 2025 09:23:40 +0100 Subject: [PATCH 36/44] Improve sequential code path --- cpp/dolfinx/refinement/refine.h | 22 ++++++++++++---------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h index 077d5bd5a1c..eae18da556f 100644 --- a/cpp/dolfinx/refinement/refine.h +++ b/cpp/dolfinx/refinement/refine.h @@ -18,6 +18,7 @@ #include "plaza.h" #include #include +#include #include #include #include @@ -49,22 +50,20 @@ create_identity_partitioner(const mesh::Mesh& parent_mesh, std::vector> cells) -> graph::AdjacencyList { - auto parent_top = parent_mesh.topology(); - auto parent_cell_im = parent_top->index_map(parent_top->dim()); + auto cell_im + = parent_mesh.topology()->index_map(parent_mesh.topology()->dim()); - int num_cell_vertices = mesh::num_cell_vertices(cell_types.front()); - std::int32_t num_cells = cells.front().size() / num_cell_vertices; + std::int32_t num_cells = cells.front().size() / mesh::num_cell_vertices(cell_types.front()); std::vector destinations(num_cells); int rank = dolfinx::MPI::rank(comm); for (std::int32_t i = 0; i < destinations.size(); i++) { - bool ghost_parent_cell = parent_cell[i] > parent_cell_im->size_local(); + bool ghost_parent_cell = parent_cell[i] > cell_im->size_local(); if (ghost_parent_cell) { destinations[i] - = parent_cell_im - ->owners()[parent_cell[i] - parent_cell_im->size_local()]; + = cell_im->owners()[parent_cell[i] - cell_im->size_local()]; } else { @@ -72,15 +71,18 @@ create_identity_partitioner(const mesh::Mesh& parent_mesh, } } + if (comm == MPI_COMM_NULL) + { + return graph::regular_adjacency_list(std::move(destinations), 1); + } + auto dual_graph = mesh::build_dual_graph(comm, cell_types, cells); - std::vector node_disp; - node_disp = std::vector(MPI::size(comm) + 1, 0); + std::vector node_disp(MPI::size(comm) + 1, 0); std::int32_t local_size = dual_graph.num_nodes(); MPI_Allgather(&local_size, 1, dolfinx::MPI::mpi_t, node_disp.data() + 1, 1, dolfinx::MPI::mpi_t, comm); std::partial_sum(node_disp.begin(), node_disp.end(), node_disp.begin()); - return compute_destination_ranks(comm, dual_graph, node_disp, destinations); }; } From d2ac3336561e80dffcd3f29be4e0800ff8b9dc00 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Wed, 12 Mar 2025 09:28:59 +0100 Subject: [PATCH 37/44] Add explicit instantiations --- cpp/dolfinx/graph/partitioners.cpp | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/cpp/dolfinx/graph/partitioners.cpp b/cpp/dolfinx/graph/partitioners.cpp index e86b6ad6f2b..a013068afaf 100644 --- a/cpp/dolfinx/graph/partitioners.cpp +++ b/cpp/dolfinx/graph/partitioners.cpp @@ -201,6 +201,15 @@ graph::AdjacencyList dolfinx::graph::compute_destination_ranks( return g; } +template graph::AdjacencyList dolfinx::graph::compute_destination_ranks( + MPI_Comm comm, const graph::AdjacencyList& graph, + const std::vector& node_disp, const std::vector& part); + +template graph::AdjacencyList dolfinx::graph::compute_destination_ranks( + MPI_Comm comm, const graph::AdjacencyList& graph, + const std::vector& node_disp, + const std::vector& part); + //----------------------------------------------------------------------------- #ifdef HAS_PARMETIS template From b59a1675e253b86ea08e8a40430c5221c9606eb8 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Wed, 12 Mar 2025 10:01:27 +0100 Subject: [PATCH 38/44] Doxygen fix explicit instantiation --- cpp/dolfinx/graph/partitioners.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/cpp/dolfinx/graph/partitioners.cpp b/cpp/dolfinx/graph/partitioners.cpp index a013068afaf..6d568fc5102 100644 --- a/cpp/dolfinx/graph/partitioners.cpp +++ b/cpp/dolfinx/graph/partitioners.cpp @@ -201,6 +201,7 @@ graph::AdjacencyList dolfinx::graph::compute_destination_ranks( return g; } +/// @cond template graph::AdjacencyList dolfinx::graph::compute_destination_ranks( MPI_Comm comm, const graph::AdjacencyList& graph, const std::vector& node_disp, const std::vector& part); @@ -209,6 +210,7 @@ template graph::AdjacencyList dolfinx::graph::compute_destination_ranks( MPI_Comm comm, const graph::AdjacencyList& graph, const std::vector& node_disp, const std::vector& part); +/// @endcond //----------------------------------------------------------------------------- #ifdef HAS_PARMETIS From 5419dee23c89e4db3cfc1a69ed8c3c5043c54659 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Wed, 12 Mar 2025 10:02:09 +0100 Subject: [PATCH 39/44] Move docstring to header --- cpp/dolfinx/graph/partitioners.h | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/cpp/dolfinx/graph/partitioners.h b/cpp/dolfinx/graph/partitioners.h index 29aa62af97f..26a78309790 100644 --- a/cpp/dolfinx/graph/partitioners.h +++ b/cpp/dolfinx/graph/partitioners.h @@ -12,6 +12,17 @@ namespace dolfinx::graph { +/// @todo Is it un-documented that the owning rank must come first in +/// reach list of edges? +/// +/// @param[in] comm The communicator +/// @param[in] graph Graph, using global indices for graph edges +/// @param[in] node_disp The distribution of graph nodes across MPI +/// ranks. The global index `gidx` of local index `lidx` is `lidx + +/// node_disp[my_rank]`. +/// @param[in] part The destination rank for owned nodes, i.e. `dest[i]` +/// is the destination of the node with local index `i`. +/// @return Destination ranks for each local node. template graph::AdjacencyList compute_destination_ranks( MPI_Comm comm, const graph::AdjacencyList& graph, From 6f1fcbcfb9694af5f07aaee6c16cb846f3285653 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Wed, 12 Mar 2025 10:07:48 +0100 Subject: [PATCH 40/44] Remove cpp docstring --- cpp/dolfinx/graph/partitioners.cpp | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/cpp/dolfinx/graph/partitioners.cpp b/cpp/dolfinx/graph/partitioners.cpp index 6d568fc5102..425ca4e2910 100644 --- a/cpp/dolfinx/graph/partitioners.cpp +++ b/cpp/dolfinx/graph/partitioners.cpp @@ -35,17 +35,6 @@ extern "C" using namespace dolfinx; -/// @todo Is it un-documented that the owning rank must come first in -/// reach list of edges? -/// -/// @param[in] comm The communicator -/// @param[in] graph Graph, using global indices for graph edges -/// @param[in] node_disp The distribution of graph nodes across MPI -/// ranks. The global index `gidx` of local index `lidx` is `lidx + -/// node_disp[my_rank]`. -/// @param[in] part The destination rank for owned nodes, i.e. `dest[i]` -/// is the destination of the node with local index `i`. -/// @return Destination ranks for each local node. template graph::AdjacencyList dolfinx::graph::compute_destination_ranks( MPI_Comm comm, const graph::AdjacencyList& graph, From 701e39906b90610adf22ab3032a214bc7677acd0 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Wed, 12 Mar 2025 15:40:36 +0100 Subject: [PATCH 41/44] Change defaults and add special case of config --- cpp/dolfinx/refinement/refine.h | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h index eae18da556f..1b7f44fa91a 100644 --- a/cpp/dolfinx/refinement/refine.h +++ b/cpp/dolfinx/refinement/refine.h @@ -21,6 +21,7 @@ #include #include #include +#include #include namespace dolfinx::refinement @@ -53,7 +54,8 @@ create_identity_partitioner(const mesh::Mesh& parent_mesh, auto cell_im = parent_mesh.topology()->index_map(parent_mesh.topology()->dim()); - std::int32_t num_cells = cells.front().size() / mesh::num_cell_vertices(cell_types.front()); + std::int32_t num_cells + = cells.front().size() / mesh::num_cell_vertices(cell_types.front()); std::vector destinations(num_cells); int rank = dolfinx::MPI::rank(comm); @@ -123,7 +125,7 @@ std::tuple, std::optional>, refine(const mesh::Mesh& mesh, std::optional> edges, mesh::CellPartitionFunction partitioner = nullptr, - Option option = Option::none) + Option option = Option::parent_cell) { auto topology = mesh.topology(); assert(topology); @@ -137,6 +139,11 @@ refine(const mesh::Mesh& mesh, if (!partitioner) { + if (!parent_cell) + { + throw std::runtime_error( + "Identity partitioner relies on parent cell computation"); + } assert(parent_cell); partitioner = create_identity_partitioner(mesh, parent_cell.value()); } From 1d19a2128f0012fe3c35172b083d4609d4d9381f Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Wed, 12 Mar 2025 15:51:37 +0100 Subject: [PATCH 42/44] Switch to optional to handle cases correctly --- cpp/dolfinx/refinement/refine.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h index 1b7f44fa91a..8fa4e61432b 100644 --- a/cpp/dolfinx/refinement/refine.h +++ b/cpp/dolfinx/refinement/refine.h @@ -124,7 +124,7 @@ std::tuple, std::optional>, std::optional>> refine(const mesh::Mesh& mesh, std::optional> edges, - mesh::CellPartitionFunction partitioner = nullptr, + std::optional partitioner = std::nullopt, Option option = Option::parent_cell) { auto topology = mesh.topology(); @@ -137,7 +137,7 @@ refine(const mesh::Mesh& mesh, ? interval::compute_refinement_data(mesh, edges, option) : plaza::compute_refinement_data(mesh, edges, option); - if (!partitioner) + if (!partitioner.has_value()) { if (!parent_cell) { @@ -150,7 +150,7 @@ refine(const mesh::Mesh& mesh, mesh::Mesh mesh1 = mesh::create_mesh( mesh.comm(), mesh.comm(), cell_adj.array(), mesh.geometry().cmap(), - mesh.comm(), new_vertex_coords, xshape, partitioner); + mesh.comm(), new_vertex_coords, xshape, *partitioner); // Report the number of refined cells const int D = topology->dim(); From 5629a0e5e82495f612971fe36cc7deb5806e6276 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Wed, 12 Mar 2025 15:54:48 +0100 Subject: [PATCH 43/44] Update docstring --- cpp/dolfinx/refinement/refine.h | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h index 8fa4e61432b..8675f16876b 100644 --- a/cpp/dolfinx/refinement/refine.h +++ b/cpp/dolfinx/refinement/refine.h @@ -103,20 +103,17 @@ create_identity_partitioner(const mesh::Mesh& parent_mesh, /// refined mesh will **not** include ghosts cells (cells connected by /// facet to an owned cell) even if the parent mesh is ghosted. /// -/// @warning Passing `nullptr` for `partitioner`, the refined mesh will -/// **not** have ghosts cells even if the parent mesh is ghosted. The -/// possibility to not re-partition the refined mesh and include ghost -/// cells in the refined mesh will be added in a future release. /// /// @param[in] mesh Input mesh to be refined. /// @param[in] edges Indices of the edges that should be split in the /// refinement. If not provided (`std::nullopt`), uniform refinement is /// performed. -/// @param[in] partitioner Partitioner to be used to distribute the -/// refined mesh. If not callable, refined cells will be on the same +/// @param[in] partitioner (Optional) partitioner to be used to distribute the +/// refined mesh. If not provided (`std::nullopt`) the partition of the coarse +/// mesh will be maintained. If not callable, refined cells will be on the same /// process as the parent cell. /// @param[in] option Control the computation of parent facets, parent -/// cells. If an option is not selected, an empty list is returned. +/// cells. /// @return New mesh, and optional parent cell indices and parent facet /// indices. template From e86a7ba49e2a83ccaa35c66f3fa1d843d14d9672 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Tue, 18 Mar 2025 10:47:27 +0100 Subject: [PATCH 44/44] Remove duplicated 'inclusion_mapping' definition --- cpp/dolfinx/multigrid/transfer_matrix.h | 68 +---------- cpp/test/multigrid/transfer_matrix.cpp | 144 +++++++++++++++++------- 2 files changed, 102 insertions(+), 110 deletions(-) diff --git a/cpp/dolfinx/multigrid/transfer_matrix.h b/cpp/dolfinx/multigrid/transfer_matrix.h index bb68fc49fa3..fc828875666 100644 --- a/cpp/dolfinx/multigrid/transfer_matrix.h +++ b/cpp/dolfinx/multigrid/transfer_matrix.h @@ -9,7 +9,6 @@ #include #include #include -#include #include #include @@ -19,75 +18,10 @@ #include "dolfinx/fem/FunctionSpace.h" #include "dolfinx/la/SparsityPattern.h" #include "dolfinx/la/utils.h" -#include "dolfinx/mesh/Mesh.h" namespace dolfinx::multigrid { -template -std::vector -inclusion_mapping(const dolfinx::mesh::Mesh& mesh_from, - const dolfinx::mesh::Mesh& mesh_to) -{ - - const common::IndexMap& im_from = *mesh_from.topology()->index_map(0); - const common::IndexMap& im_to = *mesh_to.topology()->index_map(0); - - std::vector map(im_from.size_global(), -1); - - std::span x_from = mesh_from.geometry().x(); - std::span x_to = mesh_to.geometry().x(); - - auto build_global_to_local = [&](const auto& im) - { - return [&](std::int32_t idx) - { - std::array tmp; - im.local_to_global(std::vector{idx}, tmp); - return tmp[0]; - }; - }; - - auto to_global_to = build_global_to_local(im_to); - auto to_global_from = build_global_to_local(im_from); - - for (std::int32_t i = 0; i < im_from.size_local(); i++) - { - std::ranges::subrange vertex_from(std::next(x_from.begin(), 3 * i), - std::next(x_from.begin(), 3 * (i + 1))); - for (std::int64_t j = 0; j < im_to.size_local() + im_to.num_ghosts(); j++) - { - std::ranges::subrange vertex_to(std::next(x_to.begin(), 3 * j), - std::next(x_to.begin(), 3 * (j + 1))); - - if (std::ranges::equal( - vertex_from, vertex_to, [](T a, T b) - { return std::abs(a - b) <= std::numeric_limits::epsilon(); })) - { - map[to_global_from(i)] = to_global_to(j); - break; - } - } - } - - if (dolfinx::MPI::size(mesh_to.comm()) == 1) - { - // no communication required - assert(std::ranges::all_of(map, [](auto e) { return e >= 0; })); - return map; - } - - // map holds at this point for every original local index the corresponding - // mapped global index. All other entries are still -1, but are available on - // other processes. - std::vector result(map.size(), -1); - MPI_Allreduce(map.data(), result.data(), map.size(), - dolfinx::MPI::mpi_t, MPI_MAX, mesh_from.comm()); - - assert(std::ranges::all_of(result, [](auto e) { return e >= 0; })); - return result; -} - template la::SparsityPattern create_sparsity_pattern(const dolfinx::fem::FunctionSpace& V_from, @@ -256,4 +190,4 @@ void assemble_transfer_matrix(la::MatSet auto mat_set, } } -} // namespace dolfinx::transfer \ No newline at end of file +} // namespace dolfinx::multigrid \ No newline at end of file diff --git a/cpp/test/multigrid/transfer_matrix.cpp b/cpp/test/multigrid/transfer_matrix.cpp index 60fa6dd3b04..d5066d98fb8 100644 --- a/cpp/test/multigrid/transfer_matrix.cpp +++ b/cpp/test/multigrid/transfer_matrix.cpp @@ -28,9 +28,10 @@ #include #include #include +#include +#include #include #include -#include using namespace dolfinx; using namespace Catch::Matchers; @@ -81,8 +82,8 @@ TEMPLATE_TEST_CASE("Transfer Matrix 1D", "[transfer_matrix]", double, float) la::MatrixCSR transfer_matrix(std::move(sp), la::BlockMode::compact); multigrid::assemble_transfer_matrix(transfer_matrix.mat_set_values(), - *V_coarse, *V_fine, inclusion_map, - weight); + *V_coarse, *V_fine, inclusion_map, + weight); // clang-format off std::vector expected{1.0, 0.5, 0.0, 0.0, 0.0, @@ -106,7 +107,8 @@ TEMPLATE_TEST_CASE("Transfer Matrix 1D", "[transfer_matrix]", double, float) // int rank = dolfinx::MPI::rank(MPI_COMM_WORLD); // // TODO: see https://github.com/FEniCS/dolfinx/issues/3358 -// // auto part = mesh::create_cell_partitioner(mesh::GhostMode::shared_facet); +// // auto part = +// mesh::create_cell_partitioner(mesh::GhostMode::shared_facet); // mesh::CellPartitionFunction part // = [&](MPI_Comm /* comm */, int /* nparts */, // const std::vector& /* cell_types */, @@ -193,7 +195,8 @@ TEMPLATE_TEST_CASE("Transfer Matrix 1D", "[transfer_matrix]", double, float) // RangeEquals(inclusion_map)); // la::SparsityPattern sp -// = multigrid::create_sparsity_pattern(*V_coarse, *V_fine, inclusion_map); +// = multigrid::create_sparsity_pattern(*V_coarse, *V_fine, +// inclusion_map); // la::MatrixCSR transfer_matrix(std::move(sp), la::BlockMode::compact); // multigrid::assemble_transfer_matrix(transfer_matrix.mat_set_values(), @@ -217,21 +220,34 @@ TEMPLATE_TEST_CASE("Transfer Matrix 1D", "[transfer_matrix]", double, float) // { // // clang-format off // expected[0] = { -// /* row_0 */ 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -// /* row_1 */ 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -// /* row_2 */ 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -// /* row_3 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -// /* row_4 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -// /* row_5 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -// /* row_6 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 +// /* row_0 */ 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* row_1 */ 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* row_2 */ 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* row_3 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* row_4 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, +// 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* row_5 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, +// 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* row_6 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 // }; // expected[1] = { -// /* row_0 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, -// /* row_1 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, -// /* row_2 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, -// /* row_3 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, -// /* row_4 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -// /* row_5 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 +// /* row_0 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* row_1 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, +// /* row_2 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, +// /* row_3 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, +// /* row_4 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* row_5 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 // }; // // clang-format on // } @@ -239,31 +255,73 @@ TEMPLATE_TEST_CASE("Transfer Matrix 1D", "[transfer_matrix]", double, float) // { // // clang-format off // expected[0] = { -// /* row_0 */ 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -// /* row_1 */ 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -// /* row_2 */ 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -// /* row_3 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -// /* row_4 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -// /* row_5 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -// /* row_6 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* row_0 */ 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, +// /* row_1 */ 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, +// /* row_2 */ 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, +// /* row_3 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, +// /* row_4 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, +// /* row_5 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, +// /* row_6 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, // }; // expected[1] = { -// /* row_0 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -// /* row_1 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -// /* row_2 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -// /* row_3 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -// /* row_4 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -// /* row_5 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, -// /* row_6 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* row_0 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.5, +// 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, +// /* row_1 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, 0.0, +// /* row_2 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, +// /* row_3 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, +// /* row_4 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, +// /* row_5 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.5, +// /* row_6 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, // }; // expected[2] = { -// /* row_0 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, -// /* row_1 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, -// /* row_2 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, -// /* row_3 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -// /* row_4 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -// /* row_5 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -// /* row_6 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 +// /* row_0 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.5, 1.0, 0.5, +// /* row_1 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.5, 1.0, 0.5, 0.0, 0.0, +// /* row_2 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, +// 0.0, 0.0, 0.0, 0.0, +// /* row_3 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.5, 0.5, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, +// /* row_4 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, +// /* row_5 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, +// /* row_6 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0 // }; // // clang-format on // } @@ -311,8 +369,8 @@ TEMPLATE_TEST_CASE("Transfer Matrix 2D", "[transfer_matrix]", double) = multigrid::create_sparsity_pattern(*V_coarse, *V_fine, inclusion_map); la::MatrixCSR transfer_matrix(std::move(sp), la::BlockMode::compact); multigrid::assemble_transfer_matrix(transfer_matrix.mat_set_values(), - *V_coarse, *V_fine, inclusion_map, - weight); + *V_coarse, *V_fine, inclusion_map, + weight); transfer_matrix.scatter_rev(); std::vector expected{0.5, 0.0, 0.5, 0.0, 1.0, 0.0, 0.5, 0.0, 0.0, 0.5, 1.0, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, @@ -358,8 +416,8 @@ TEMPLATE_TEST_CASE("Transfer Matrix 3D", "[transfer_matrix]", double) = multigrid::create_sparsity_pattern(*V_coarse, *V_fine, inclusion_map); la::MatrixCSR transfer_matrix(std::move(sp), la::BlockMode::compact); multigrid::assemble_transfer_matrix(transfer_matrix.mat_set_values(), - *V_coarse, *V_fine, inclusion_map, - weight); + *V_coarse, *V_fine, inclusion_map, + weight); std::vector expected{ 1.0, 0.5, 0.5, 0.5, 0.5, 0.5, 0.0, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,