From 9d549d00f0f03c1a60cf94302175299199dd8a6b Mon Sep 17 00:00:00 2001 From: blegouix Date: Sun, 1 Mar 2026 11:12:51 +0100 Subject: [PATCH] wip --- AGENTS.md | 6 +- src/similie/exterior/coboundary.hpp | 307 ++++++++++++++++++---------- 2 files changed, 209 insertions(+), 104 deletions(-) diff --git a/AGENTS.md b/AGENTS.md index 0c9eca0..cd988b3 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -9,7 +9,11 @@ SPDX-License-Identifier: MIT - Clone with `git clone --recurse-submodules git@github.com:blegouix/similie.git` or any similar commande which suits the network configuration. Be sure to pull the submodules recursively. # Build environment -- **Requirement:** Build and use the `docker/similie/Dockerfile` image before compiling or running tests. Do not build on the host directly. If `docker` is not installed, install it first (refer to the online documentation). +- **Requirement:** Build and use the `docker/similie/Dockerfile` image only if mandatory host tools are missing for the selected build mode. + - For CPU builds, required host tools include `cmake` and `g++-13`. + - For CUDA builds, required host tools include `cmake`, `nvcc`, and a working CUDA toolchain compatible with Kokkos. + - If all mandatory tools are available on host, host builds/tests are allowed. + - If any mandatory tool is missing, use the Docker image before compiling or running tests. If `docker` is not installed, install it first (refer to the online documentation). - Configure `CMake` for CPU with `cmake -DCMAKE_CXX_COMPILER=g++-13 -DCMAKE_BUILD_TYPE=Debug -DKokkos_ENABLE_OPENMP=ON -B build -S .` or for CUDA with `cmake -DCMAKE_CXX_COMPILER=$PWD/vendor/ddc/vendor/kokkos/bin/nvcc_wrapper -DCMAKE_BUILD_TYPE=Debug -DKokkos_ENABLE_CUDA=ON -DKokkos_ARCH_=ON -B build -S .`. Of course replace `` with the arch namecode which suits the GPU. - Compile as usual once configured (`make -jN_PROC` from `build`). diff --git a/src/similie/exterior/coboundary.hpp b/src/similie/exterior/coboundary.hpp index b817ec1..6d50f01 100644 --- a/src/similie/exterior/coboundary.hpp +++ b/src/similie/exterior/coboundary.hpp @@ -3,9 +3,12 @@ #pragma once +#include + #include #include +#include #include #include #include @@ -184,10 +187,206 @@ struct NonSpectatorDimension> ddc::DiscreteDomain<>>...>; }; -struct CoboundaryDummyIndex +template +struct CoboundaryRank; + +template +struct CoboundaryRank> { + static constexpr std::size_t value = sizeof...(DDim); }; +template +inline constexpr std::size_t coboundary_rank_v = CoboundaryRank::value; + +template +KOKKOS_FUNCTION void fill_simplex_boundary_half_( + std::array, SimplexType::dimension()>& out, + typename boundary_t::discrete_element_type elem, + typename SimplexType::discrete_vector_type vect, + bool negative = false) +{ + auto array = ddc::detail::array(vect); + int id_dist = -1; + for (std::size_t i = 0; i < SimplexType::dimension(); ++i) { + auto array_ = array; + auto id = array_.begin() + id_dist + 1; + while (id != array_.end() && *id == 0) { + ++id; + } + assert(id != array_.end()); + id_dist = static_cast(id - array_.begin()); + *id = 0; + typename SimplexType::discrete_vector_type vect_; + ddc::detail::array(vect_) = array_; + out[i] = boundary_t(elem, vect_, (negative + i) % 2); + } +} + +template +KOKKOS_FUNCTION void fill_simplex_boundary_( + std::array, 2 * SimplexType::dimension()>& out, + SimplexType simplex) +{ + std::array, SimplexType::dimension()> lower_half; + fill_simplex_boundary_half_( + lower_half, + simplex.discrete_element(), + simplex.discrete_vector(), + SimplexType::dimension() % 2); + for (std::size_t i = 0; i < SimplexType::dimension(); ++i) { + out[i] = lower_half[i]; + } + + std::array, SimplexType::dimension()> upper_half; + fill_simplex_boundary_half_( + upper_half, + simplex.discrete_element() + simplex.discrete_vector(), + -simplex.discrete_vector()); + for (std::size_t i = 0; i < SimplexType::dimension(); ++i) { + out[i + SimplexType::dimension()] = upper_half[i]; + } + + if constexpr (SimplexType::dimension() % 2 == 0) { + if (!simplex.negative()) { + for (std::size_t i = 0; i < 2 * SimplexType::dimension(); ++i) { + out[i].negative() = !out[i].negative(); + } + } + } else { + if (simplex.negative()) { + for (std::size_t i = 0; i < 2 * SimplexType::dimension(); ++i) { + out[i].negative() = !out[i].negative(); + } + } + } +} + +template +KOKKOS_FUNCTION void fill_tangent_basis_(std::array& basis) +{ + if constexpr (NBasis == 0) { + return; + } else if constexpr (K == 0) { + VectType vect; + auto array = ddc::detail::array(vect); + for (std::size_t j = 0; j < NDims; ++j) { + array[j] = 0; + } + basis[0] = vect; + return; + } else { + std::array indices {}; + for (std::size_t i = 0; i < K; ++i) { + indices[i] = static_cast(i); + } + + std::size_t out_idx = 0; + while (true) { + VectType vect; + auto array = ddc::detail::array(vect); + for (std::size_t j = 0; j < NDims; ++j) { + array[j] = 0; + } + for (std::size_t j = 0; j < K; ++j) { + array[static_cast(indices[j])] = 1; + } + basis[out_idx++] = vect; + + int pos = static_cast(K) - 1; + while (pos >= 0 && indices[static_cast(pos)] == static_cast(NDims - K + static_cast(pos))) { + --pos; + } + if (pos < 0) { + break; + } + ++indices[static_cast(pos)]; + for (std::size_t j = static_cast(pos) + 1; j < K; ++j) { + indices[j] = indices[j - 1] + 1; + } + } + } +} + +template < + tensor::TensorNatIndex TagToAddToCochain, + tensor::TensorIndex CochainTag, + misc::Specialization TensorType> +KOKKOS_FUNCTION void coboundary( + typename ddc::remove_dims_of_t< + typename coboundary_tensor_t< + TagToAddToCochain, + CochainTag, + TensorType>::discrete_domain_type, + coboundary_index_t>::discrete_element_type elem, + coboundary_tensor_t coboundary_tensor, + TensorType tensor) +{ + using full_domain_t = ddc::remove_dims_of_t< + typename coboundary_tensor_t< + TagToAddToCochain, + CochainTag, + TensorType>::discrete_domain_type, + coboundary_index_t>; + using non_spectator_domain_t = typename detail::NonSpectatorDimension< + TagToAddToCochain, + typename TensorType::non_indices_domain_t>::type; + static constexpr std::size_t n_dims = coboundary_rank_v; + static constexpr std::size_t k = CochainTag::rank(); + static constexpr std::size_t chain_size = misc::binomial_coefficient(n_dims, k + 1); + static constexpr std::size_t lower_chain_size = misc::binomial_coefficient(n_dims, k); + static constexpr std::size_t boundary_size = 2 * (k + 1); + + using simplex_top_t = simplex_for_domain_t; + using simplex_boundary_t = boundary_t; + using non_spectator_simplex_t = simplex_for_domain_t; + using non_spectator_discrete_vector_type = typename non_spectator_simplex_t::discrete_vector_type; + + std::array chain; + std::array lower_chain; + fill_tangent_basis_(chain); + fill_tangent_basis_( + lower_chain); + + std::array simplex_boundary; + std::array boundary_values; + + for (std::size_t i = 0; i < chain_size; ++i) { + simplex_top_t simplex( + std::integral_constant {}, + elem, + chain[i]); + + fill_simplex_boundary_(simplex_boundary, simplex); + + for (std::size_t j = 0; j < boundary_size; ++j) { + auto const boundary_simplex = simplex_boundary[j]; + auto lower_it = misc::detail::find( + lower_chain.begin(), + lower_chain.end(), + boundary_simplex.discrete_vector()); + assert(lower_it != lower_chain.end()); + + boundary_values[j] = tensor.mem( + misc::domain_contains(tensor.domain(), boundary_simplex.discrete_element()) + ? boundary_simplex.discrete_element() + : elem, + ddc::DiscreteElement(static_cast( + lower_it - lower_chain.begin()))); + } + + double integral = 0.; + for (std::size_t j = 0; j < boundary_size; ++j) { + integral += (simplex_boundary[j].negative() ? -1. : 1.) * boundary_values[j]; + } + + coboundary_tensor.mem( + elem, + ddc::DiscreteElement>(i)) + = integral; + } +} + } // namespace detail template < @@ -204,65 +403,6 @@ coboundary_tensor_t coboundary( = ddc::remove_dims_of>( coboundary_tensor.domain()); - // buffer to store the K-chain containing the boundary of each K+1-simplex of the mesh - ddc::Chunk simplex_boundary_alloc( - ddc::cartesian_prod_t< - ddc::remove_dims_of_t< - typename coboundary_tensor_t< - TagToAddToCochain, - CochainTag, - TensorType>::discrete_domain_type, - coboundary_index_t>, - ddc::DiscreteDomain>( - batch_dom, - ddc::DiscreteDomain( - ddc::DiscreteElement(0), - ddc::DiscreteVector( - 2 * (CochainTag::rank() + 1)))), - ddc::KokkosAllocator< - simplex_for_domain_t< - CochainTag::rank(), - ddc::remove_dims_of_t< - typename coboundary_tensor_t< - TagToAddToCochain, - CochainTag, - TensorType>::discrete_domain_type, - coboundary_index_t>>, - typename ExecSpace::memory_space>()); - ddc::ChunkSpan simplex_boundary(simplex_boundary_alloc); - - // buffer to store the values of the K-cochain on the boundary of each K+1-cosimplex of the mesh - ddc::Chunk boundary_values_alloc( - ddc::cartesian_prod_t< - ddc::remove_dims_of_t< - typename coboundary_tensor_t< - TagToAddToCochain, - CochainTag, - TensorType>::discrete_domain_type, - coboundary_index_t>, - ddc::DiscreteDomain>( - batch_dom, - ddc::DiscreteDomain( - ddc::DiscreteElement(0), - ddc::DiscreteVector( - 2 * (CochainTag::rank() + 1)))), - ddc::KokkosAllocator()); - ddc::ChunkSpan boundary_values(boundary_values_alloc); - - // compute the tangent K+1-basis for each node of the mesh. This is a local K+1-chain. - auto chain = tangent_basis< - CochainTag::rank() + 1, - typename detail::NonSpectatorDimension< - TagToAddToCochain, - typename TensorType::non_indices_domain_t>::type>(exec_space); - - // compute the tangent K-basis for each node of the mesh. This is a local K-chain. - auto lower_chain = tangent_basis< - CochainTag::rank(), - typename detail::NonSpectatorDimension< - TagToAddToCochain, - typename TensorType::non_indices_domain_t>::type>(exec_space); - // iterate over every node, we will work inside the tangent space associated to each of them SIMILIE_DEBUG_LOG("similie_compute_coboundary"); ddc::parallel_for_each( @@ -270,49 +410,10 @@ coboundary_tensor_t coboundary( exec_space, batch_dom, KOKKOS_LAMBDA(typename decltype(batch_dom)::discrete_element_type elem) { - // declare a K+1-cochain storing the K+1-cosimplices of the output cochain for the current tangent space and iterate over them - auto cochain = Cochain(chain, coboundary_tensor[elem]); - for (auto i = cochain.begin(); i < cochain.end(); ++i) { - // extract the K+1-simplex from the current K+1-cosimplex (this is not absolutly trivial because the cochain is based on a LocalChain) - Simplex simplex = Simplex( - std::integral_constant {}, - elem, - (*i).discrete_vector()); - // compute its boundary as a K-chain in the simplex_boundary buffer - Chain boundary_chain - = boundary(simplex_boundary[elem].allocation_kokkos_view(), simplex); - // iterate over every K-simplex forming the boundary - for (auto j = boundary_chain.begin(); j < boundary_chain.end(); ++j) { - // extract from the input K-cochain the values associated to every K-simplex of the boundary and fill the boundary_values buffer - boundary_values[elem].allocation_kokkos_view()( - Kokkos::Experimental::distance(boundary_chain.begin(), j)) - = tensor.mem( - misc::domain_contains( - tensor.domain(), - (*j).discrete_element()) - ? (*j).discrete_element() - : elem, // TODO this is an assumption on boundary condition (free boundary), needs to be generalized - ddc::DiscreteElement( - Kokkos::Experimental::distance( - lower_chain.begin(), - misc::detail:: - find(lower_chain.begin(), - lower_chain.end(), - (*j).discrete_vector())))); - } - // build the cochain of the boundary - Cochain cochain_boundary( - boundary_chain, - boundary_values[elem].allocation_kokkos_view()); - // integrate over the cochain forming the boundary to compute the coboundary - // (*i).value() = cochain_boundary.integrate(); // Cannot be used because CochainIterator::operator* does not return a reference - coboundary_tensor - .mem(elem, - ddc::DiscreteElement< - coboundary_index_t>( - Kokkos::Experimental::distance(cochain.begin(), i))) - = cochain_boundary.integrate(); - } + detail::coboundary( + elem, + coboundary_tensor, + tensor); }); return coboundary_tensor;