Device-callable coboundary

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_<YOUR_ARCH_NAMECODE>=ON -B build -S .`. Of course replace `<YOUR_ARCH_NAMECODE>` with the arch namecode which suits the GPU.
 - Compile as usual once configured (`make -jN_PROC` from `build`).
 

src/similie/exterior/coboundary.hpp

@@ -3,9 +3,12 @@
 
 #pragma once
 
+#include <array>
+
 #include <ddc/ddc.hpp>
 
 #include <similie/misc/are_all_same.hpp>
+#include <similie/misc/binomial_coefficient.hpp>
 #include <similie/misc/domain_contains.hpp>
 #include <similie/misc/filled_struct.hpp>
 #include <similie/misc/macros.hpp>
@@ -184,10 +187,206 @@ struct NonSpectatorDimension<Index, ddc::DiscreteDomain<DDim...>>
             ddc::DiscreteDomain<>>...>;
 };
 
-struct CoboundaryDummyIndex
+template <class DDim>
+struct CoboundaryRank;
+
+template <class... DDim>
+struct CoboundaryRank<ddc::DiscreteDomain<DDim...>>
 {
+    static constexpr std::size_t value = sizeof...(DDim);
 };
 
+template <class DDim>
+inline constexpr std::size_t coboundary_rank_v = CoboundaryRank<DDim>::value;
+
+template <class SimplexType>
+KOKKOS_FUNCTION void fill_simplex_boundary_half_(
+        std::array<boundary_t<SimplexType>, SimplexType::dimension()>& out,
+        typename boundary_t<SimplexType>::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<int>(id - array_.begin());
+        *id = 0;
+        typename SimplexType::discrete_vector_type vect_;
+        ddc::detail::array(vect_) = array_;
+        out[i] = boundary_t<SimplexType>(elem, vect_, (negative + i) % 2);
+    }
+}
+
+template <class SimplexType>
+KOKKOS_FUNCTION void fill_simplex_boundary_(
+        std::array<boundary_t<SimplexType>, 2 * SimplexType::dimension()>& out,
+        SimplexType simplex)
+{
+    std::array<boundary_t<SimplexType>, SimplexType::dimension()> lower_half;
+    fill_simplex_boundary_half_<SimplexType>(
+            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<boundary_t<SimplexType>, SimplexType::dimension()> upper_half;
+    fill_simplex_boundary_half_<SimplexType>(
+            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 <class VectType, std::size_t NDims, std::size_t K, std::size_t NBasis>
+KOKKOS_FUNCTION void fill_tangent_basis_(std::array<VectType, NBasis>& 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<int, K> indices {};
+        for (std::size_t i = 0; i < K; ++i) {
+            indices[i] = static_cast<int>(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<std::size_t>(indices[j])] = 1;
+            }
+            basis[out_idx++] = vect;
+
+            int pos = static_cast<int>(K) - 1;
+            while (pos >= 0 && indices[static_cast<std::size_t>(pos)] == static_cast<int>(NDims - K + static_cast<std::size_t>(pos))) {
+                --pos;
+            }
+            if (pos < 0) {
+                break;
+            }
+            ++indices[static_cast<std::size_t>(pos)];
+            for (std::size_t j = static_cast<std::size_t>(pos) + 1; j < K; ++j) {
+                indices[j] = indices[j - 1] + 1;
+            }
+        }
+    }
+}
+
+template <
+        tensor::TensorNatIndex TagToAddToCochain,
+        tensor::TensorIndex CochainTag,
+        misc::Specialization<tensor::Tensor> TensorType>
+KOKKOS_FUNCTION void coboundary(
+        typename ddc::remove_dims_of_t<
+                typename coboundary_tensor_t<
+                        TagToAddToCochain,
+                        CochainTag,
+                        TensorType>::discrete_domain_type,
+                coboundary_index_t<TagToAddToCochain, CochainTag>>::discrete_element_type elem,
+        coboundary_tensor_t<TagToAddToCochain, CochainTag, TensorType> 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<TagToAddToCochain, CochainTag>>;
+    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<non_spectator_domain_t>;
+    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<k + 1, full_domain_t>;
+    using simplex_boundary_t = boundary_t<simplex_top_t>;
+    using non_spectator_simplex_t = simplex_for_domain_t<k + 1, non_spectator_domain_t>;
+    using non_spectator_discrete_vector_type = typename non_spectator_simplex_t::discrete_vector_type;
+
+    std::array<non_spectator_discrete_vector_type, chain_size> chain;
+    std::array<non_spectator_discrete_vector_type, lower_chain_size> lower_chain;
+    fill_tangent_basis_<non_spectator_discrete_vector_type, n_dims, k + 1, chain_size>(chain);
+    fill_tangent_basis_<non_spectator_discrete_vector_type, n_dims, k, lower_chain_size>(
+            lower_chain);
+
+    std::array<simplex_boundary_t, boundary_size> simplex_boundary;
+    std::array<double, boundary_size> boundary_values;
+
+    for (std::size_t i = 0; i < chain_size; ++i) {
+        simplex_top_t simplex(
+                std::integral_constant<std::size_t, k + 1> {},
+                elem,
+                chain[i]);
+
+        fill_simplex_boundary_<simplex_top_t>(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<CochainTag>(static_cast<std::size_t>(
+                            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<coboundary_index_t<TagToAddToCochain, CochainTag>>(i))
+                = integral;
+    }
+}
+
 } // namespace detail
 
 template <
@@ -204,115 +403,17 @@ coboundary_tensor_t<TagToAddToCochain, CochainTag, TensorType> coboundary(
             = ddc::remove_dims_of<coboundary_index_t<TagToAddToCochain, CochainTag>>(
                     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<TagToAddToCochain, CochainTag>>,
-                    ddc::DiscreteDomain<detail::CoboundaryDummyIndex>>(
-                    batch_dom,
-                    ddc::DiscreteDomain<detail::CoboundaryDummyIndex>(
-                            ddc::DiscreteElement<detail::CoboundaryDummyIndex>(0),
-                            ddc::DiscreteVector<detail::CoboundaryDummyIndex>(
-                                    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<TagToAddToCochain, CochainTag>>>,
-                    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<TagToAddToCochain, CochainTag>>,
-                    ddc::DiscreteDomain<detail::CoboundaryDummyIndex>>(
-                    batch_dom,
-                    ddc::DiscreteDomain<detail::CoboundaryDummyIndex>(
-                            ddc::DiscreteElement<detail::CoboundaryDummyIndex>(0),
-                            ddc::DiscreteVector<detail::CoboundaryDummyIndex>(
-                                    2 * (CochainTag::rank() + 1)))),
-            ddc::KokkosAllocator<double, typename ExecSpace::memory_space>());
-    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(
             "similie_compute_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<std::size_t, CochainTag::rank() + 1> {},
-                            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<CochainTag>(
-                                                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<TagToAddToCochain, CochainTag>>(
-                                         Kokkos::Experimental::distance(cochain.begin(), i)))
-                            = cochain_boundary.integrate();
-                }
+                detail::coboundary<TagToAddToCochain, CochainTag, TensorType>(
+                        elem,
+                        coboundary_tensor,
+                        tensor);
             });
 
     return coboundary_tensor;