SimiLie
Loading...
Searching...
No Matches
coboundary.hpp
1// SPDX-FileCopyrightText: 2024 Baptiste Legouix
2// SPDX-License-Identifier: MIT
3
4#pragma once
5
6#include <ddc/ddc.hpp>
7
8#include <similie/misc/are_all_same.hpp>
9#include <similie/misc/domain_contains.hpp>
10#include <similie/misc/filled_struct.hpp>
11#include <similie/misc/macros.hpp>
12#include <similie/misc/portable_stl.hpp>
13#include <similie/misc/specialization.hpp>
14#include <similie/misc/type_seq_conversion.hpp>
15#include <similie/tensor/antisymmetric_tensor.hpp>
16#include <similie/tensor/dummy_index.hpp>
17#include <similie/tensor/tensor_impl.hpp>
18
19#include <Kokkos_StdAlgorithms.hpp>
20
21#include "cochain.hpp"
22#include "cosimplex.hpp"
23
24
25namespace sil {
26
27namespace exterior {
28
29namespace detail {
30
31template <class T>
32struct CoboundaryType;
33
34template <
35 std::size_t K,
36 class... Tag,
37 class ElementType,
38 class LayoutStridedPolicy1,
39 class LayoutStridedPolicy2,
40 class ExecSpace>
41struct CoboundaryType<
42 Cochain<Chain<Simplex<K, Tag...>, LayoutStridedPolicy1, ExecSpace>,
43 ElementType,
44 LayoutStridedPolicy2>>
45{
46 using type = Cosimplex<Simplex<K + 1, Tag...>, ElementType>;
47};
48
49} // namespace detail
50
51template <misc::Specialization<Cochain> CochainType>
52using coboundary_t = typename detail::CoboundaryType<CochainType>::type;
53
54namespace detail {
55
56template <class TagToAddToCochain, class CochainTag>
57struct CoboundaryIndex;
58
59template <tensor::TensorNatIndex TagToAddToCochain, tensor::TensorNatIndex CochainTag>
60 requires(CochainTag::rank() == 0)
61struct CoboundaryIndex<TagToAddToCochain, CochainTag>
62{
63 using type = TagToAddToCochain;
64};
65
66template <tensor::TensorNatIndex TagToAddToCochain, tensor::TensorNatIndex CochainTag>
67 requires(CochainTag::rank() == 1)
68struct CoboundaryIndex<TagToAddToCochain, CochainTag>
69{
70 using type = tensor::TensorAntisymmetricIndex<TagToAddToCochain, CochainTag>;
71};
72
73template <tensor::TensorNatIndex TagToAddToCochain, tensor::TensorNatIndex... Tag>
74struct CoboundaryIndex<TagToAddToCochain, tensor::TensorAntisymmetricIndex<Tag...>>
75{
76 using type = tensor::TensorAntisymmetricIndex<TagToAddToCochain, Tag...>;
77};
78
79} // namespace detail
80
81template <class TagToAddToCochain, class CochainTag>
82using coboundary_index_t = typename detail::CoboundaryIndex<TagToAddToCochain, CochainTag>::type;
83
84namespace detail {
85
86template <
87 tensor::TensorNatIndex TagToAddToCochain,
88 tensor::TensorIndex CochainTag,
89 misc::Specialization<tensor::Tensor> TensorType>
90struct CoboundaryTensorType;
91
92template <
93 tensor::TensorNatIndex TagToAddToCochain,
94 tensor::TensorIndex CochainIndex,
95 class ElementType,
96 class... DDim,
97 class SupportType,
98 class MemorySpace>
99struct CoboundaryTensorType<
100 TagToAddToCochain,
101 CochainIndex,
102 tensor::Tensor<ElementType, ddc::DiscreteDomain<DDim...>, SupportType, MemorySpace>>
103{
104 static_assert(ddc::type_seq_contains_v<
105 ddc::detail::TypeSeq<CochainIndex>,
106 ddc::detail::TypeSeq<DDim...>>);
107 using type = tensor::Tensor<
108 ElementType,
109 ddc::replace_dim_of_t<
110 ddc::DiscreteDomain<DDim...>,
111 CochainIndex,
112 coboundary_index_t<TagToAddToCochain, CochainIndex>>,
113 SupportType,
114 MemorySpace>;
115};
116
117} // namespace detail
118
119template <
120 tensor::TensorNatIndex TagToAddToCochain,
121 tensor::TensorIndex CochainTag,
122 misc::Specialization<tensor::Tensor> TensorType>
123using coboundary_tensor_t =
124 typename detail::CoboundaryTensorType<TagToAddToCochain, CochainTag, TensorType>::type;
125
126namespace detail {
127
128template <misc::Specialization<Chain> ChainType>
129struct ComputeSimplex;
130
131template <std::size_t K, class... Tag, class LayoutStridedPolicy, class ExecSpace>
132struct ComputeSimplex<Chain<Simplex<K, Tag...>, LayoutStridedPolicy, ExecSpace>>
133{
134 KOKKOS_FUNCTION static Simplex<K + 1, Tag...> run(
135 Chain<Simplex<K, Tag...>, LayoutStridedPolicy, ExecSpace> const& chain)
136 {
137 ddc::DiscreteVector<Tag...> vect {
138 0 * ddc::type_seq_rank_v<Tag, ddc::detail::TypeSeq<Tag...>>...};
139 for (auto i = chain.begin(); i < chain.end(); ++i) {
140 vect = ddc::DiscreteVector<Tag...> {
141 (static_cast<bool>(vect.template get<Tag>())
142 || static_cast<bool>((*i).discrete_vector().template get<Tag>()))...};
143 }
144 return Simplex(
145 std::integral_constant<std::size_t, K + 1> {},
146 chain[0].discrete_element(), // This is an assumption on the structure of the chain, which is satisfied if it has been produced using boundary()
147 vect);
148 }
149};
150
151} // namespace detail
152
153template <misc::Specialization<Cochain> CochainType>
154KOKKOS_FUNCTION coboundary_t<CochainType> coboundary(
155 CochainType
156 cochain) // Warning: only cochain.chain() produced using boundary() are supported
157{
158 assert(cochain.size() == 2 * (cochain.dimension() + 1)
159 && "only cochain over the boundary of a single simplex is supported");
160
161 /* Commented because would require an additional buffer
162 assert(boundary(cochain.chain()) == boundary_t<typename CochainType::chain_type> {}
163 && "only cochain over the boundary of a single simplex is supported");
164 */
165
166 return coboundary_t<CochainType>(
167 detail::ComputeSimplex<typename CochainType::chain_type>::run(cochain.chain()),
168 cochain.integrate());
169}
170
171namespace detail {
172
173template <tensor::TensorNatIndex Index, class Dom>
174struct NonSpectatorDimension;
175
176template <tensor::TensorNatIndex Index, class... DDim>
177struct NonSpectatorDimension<Index, ddc::DiscreteDomain<DDim...>>
178{
179 using type = ddc::cartesian_prod_t<std::conditional_t<
180 ddc::type_seq_contains_v<
181 ddc::detail::TypeSeq<typename DDim::continuous_dimension_type>,
182 typename Index::type_seq_dimensions>,
183 ddc::DiscreteDomain<DDim>,
184 ddc::DiscreteDomain<>>...>;
185};
186
187struct CoboundaryDummyIndex
188{
189};
190
191} // namespace detail
192
193template <
194 tensor::TensorNatIndex TagToAddToCochain,
195 tensor::TensorIndex CochainTag,
196 misc::Specialization<tensor::Tensor> TensorType,
197 class ExecSpace>
198coboundary_tensor_t<TagToAddToCochain, CochainTag, TensorType> coboundary(
199 ExecSpace const& exec_space,
200 coboundary_tensor_t<TagToAddToCochain, CochainTag, TensorType> coboundary_tensor,
201 TensorType tensor)
202{
203 ddc::DiscreteDomain batch_dom
204 = ddc::remove_dims_of<coboundary_index_t<TagToAddToCochain, CochainTag>>(
205 coboundary_tensor.domain());
206
207 // buffer to store the K-chain containing the boundary of each K+1-simplex of the mesh
208 ddc::Chunk simplex_boundary_alloc(
209 ddc::cartesian_prod_t<
210 ddc::remove_dims_of_t<
211 typename coboundary_tensor_t<
212 TagToAddToCochain,
213 CochainTag,
214 TensorType>::discrete_domain_type,
215 coboundary_index_t<TagToAddToCochain, CochainTag>>,
216 ddc::DiscreteDomain<detail::CoboundaryDummyIndex>>(
217 batch_dom,
218 ddc::DiscreteDomain<detail::CoboundaryDummyIndex>(
219 ddc::DiscreteElement<detail::CoboundaryDummyIndex>(0),
220 ddc::DiscreteVector<detail::CoboundaryDummyIndex>(
221 2 * (CochainTag::rank() + 1)))),
222 ddc::KokkosAllocator<
223 simplex_for_domain_t<
224 CochainTag::rank(),
225 ddc::remove_dims_of_t<
226 typename coboundary_tensor_t<
227 TagToAddToCochain,
228 CochainTag,
229 TensorType>::discrete_domain_type,
230 coboundary_index_t<TagToAddToCochain, CochainTag>>>,
231 typename ExecSpace::memory_space>());
232 ddc::ChunkSpan simplex_boundary(simplex_boundary_alloc);
233
234 // buffer to store the values of the K-cochain on the boundary of each K+1-cosimplex of the mesh
235 ddc::Chunk boundary_values_alloc(
236 ddc::cartesian_prod_t<
237 ddc::remove_dims_of_t<
238 typename coboundary_tensor_t<
239 TagToAddToCochain,
240 CochainTag,
241 TensorType>::discrete_domain_type,
242 coboundary_index_t<TagToAddToCochain, CochainTag>>,
243 ddc::DiscreteDomain<detail::CoboundaryDummyIndex>>(
244 batch_dom,
245 ddc::DiscreteDomain<detail::CoboundaryDummyIndex>(
246 ddc::DiscreteElement<detail::CoboundaryDummyIndex>(0),
247 ddc::DiscreteVector<detail::CoboundaryDummyIndex>(
248 2 * (CochainTag::rank() + 1)))),
249 ddc::KokkosAllocator<double, typename ExecSpace::memory_space>());
250 ddc::ChunkSpan boundary_values(boundary_values_alloc);
251
252 // compute the tangent K+1-basis for each node of the mesh. This is a local K+1-chain.
253 auto chain = tangent_basis<
254 CochainTag::rank() + 1,
255 typename detail::NonSpectatorDimension<
256 TagToAddToCochain,
257 typename TensorType::non_indices_domain_t>::type>(exec_space);
258
259 // compute the tangent K-basis for each node of the mesh. This is a local K-chain.
260 auto lower_chain = tangent_basis<
261 CochainTag::rank(),
262 typename detail::NonSpectatorDimension<
263 TagToAddToCochain,
264 typename TensorType::non_indices_domain_t>::type>(exec_space);
265
266 // iterate over every node, we will work inside the tangent space associated to each of them
267 SIMILIE_DEBUG_LOG("similie_compute_coboundary");
268 ddc::parallel_for_each(
269 "similie_compute_coboundary",
270 exec_space,
271 batch_dom,
272 KOKKOS_LAMBDA(typename decltype(batch_dom)::discrete_element_type elem) {
273 // declare a K+1-cochain storing the K+1-cosimplices of the output cochain for the current tangent space and iterate over them
274 auto cochain = Cochain(chain, coboundary_tensor[elem]);
275 for (auto i = cochain.begin(); i < cochain.end(); ++i) {
276 // extract the K+1-simplex from the current K+1-cosimplex (this is not absolutly trivial because the cochain is based on a LocalChain)
277 Simplex simplex = Simplex(
278 std::integral_constant<std::size_t, CochainTag::rank() + 1> {},
279 elem,
280 (*i).discrete_vector());
281 // compute its boundary as a K-chain in the simplex_boundary buffer
282 Chain boundary_chain
283 = boundary(simplex_boundary[elem].allocation_kokkos_view(), simplex);
284 // iterate over every K-simplex forming the boundary
285 for (auto j = boundary_chain.begin(); j < boundary_chain.end(); ++j) {
286 // extract from the input K-cochain the values associated to every K-simplex of the boundary and fill the boundary_values buffer
287 boundary_values[elem].allocation_kokkos_view()(
288 Kokkos::Experimental::distance(boundary_chain.begin(), j))
289 = tensor.mem(
290 misc::domain_contains(
291 tensor.domain(),
292 (*j).discrete_element())
293 ? (*j).discrete_element()
294 : elem, // TODO this is an assumption on boundary condition (free boundary), needs to be generalized
295 ddc::DiscreteElement<CochainTag>(
296 Kokkos::Experimental::distance(
297 lower_chain.begin(),
298 misc::detail::
299 find(lower_chain.begin(),
300 lower_chain.end(),
301 (*j).discrete_vector()))));
302 }
303 // build the cochain of the boundary
304 Cochain cochain_boundary(
305 boundary_chain,
306 boundary_values[elem].allocation_kokkos_view());
307 // integrate over the cochain forming the boundary to compute the coboundary
308 // (*i).value() = cochain_boundary.integrate(); // Cannot be used because CochainIterator::operator* does not return a reference
309 coboundary_tensor
310 .mem(elem,
311 ddc::DiscreteElement<
312 coboundary_index_t<TagToAddToCochain, CochainTag>>(
313 Kokkos::Experimental::distance(cochain.begin(), i)))
314 = cochain_boundary.integrate();
315 }
316 });
317
318 return coboundary_tensor;
319}
320
321template <
322 tensor::TensorNatIndex TagToAddToCochain,
323 tensor::TensorIndex CochainTag,
324 misc::Specialization<tensor::Tensor> TensorType,
325 class ExecSpace>
326coboundary_tensor_t<TagToAddToCochain, CochainTag, TensorType> deriv(
327 ExecSpace const& exec_space,
328 coboundary_tensor_t<TagToAddToCochain, CochainTag, TensorType> coboundary_tensor,
329 TensorType tensor)
330{
331 return coboundary<
332 TagToAddToCochain,
333 CochainTag,
334 TensorType>(exec_space, coboundary_tensor, tensor);
335}
336
337} // namespace exterior
338
339} // namespace sil
sil::exterior::Chain
Chain class.
Definition chain.hpp:26
sil::exterior::Chain::end
KOKKOS_FUNCTION auto end()
Definition chain.hpp:165
sil::exterior::Chain::begin
KOKKOS_FUNCTION auto begin()
Definition chain.hpp:155
sil::exterior::Cochain
Cochain class.
Definition cochain.hpp:30
sil::exterior::Cochain::integrate
KOKKOS_FUNCTION element_type integrate() noexcept
Definition cochain.hpp:187
sil::exterior::Simplex
Simplex class.
Definition simplex.hpp:75
sil::exterior::tangent_basis
constexpr auto tangent_basis(ExecSpace const &exec_space)
Definition local_chain.hpp:410
sil::exterior::Cochain
Cochain(ChainType, TensorType) -> Cochain< ChainType, typename TensorType::value_type, ddc::detail::mdspan_to_kokkos_layout_t< typename TensorType::layout_type > >
sil::exterior::Simplex
Simplex(ddc::DiscreteElement< Tag... >, ddc::DiscreteVector< T... >) -> Simplex< sizeof...(T), Tag... >
sil::exterior::coboundary
KOKKOS_FUNCTION coboundary_t< CochainType > coboundary(CochainType cochain)
Definition coboundary.hpp:154
sil::exterior::Chain
Chain(Head, Tail...) -> Chain< typename Head::value_type, typename Head::array_layout, typename Head::memory_space >
sil::exterior::coboundary_tensor_t
typename detail::CoboundaryTensorType< TagToAddToCochain, CochainTag, TensorType >::type coboundary_tensor_t
Definition coboundary.hpp:123
sil::exterior::deriv
coboundary_tensor_t< TagToAddToCochain, CochainTag, TensorType > deriv(ExecSpace const &exec_space, coboundary_tensor_t< TagToAddToCochain, CochainTag, TensorType > coboundary_tensor, TensorType tensor)
Definition coboundary.hpp:326
sil::exterior::boundary
KOKKOS_FUNCTION Chain< boundary_t< SimplexType >, typename AllocationType::array_layout, typename AllocationType::memory_space > boundary(AllocationType allocation, SimplexType simplex)
Definition boundary.hpp:79
sil::exterior::simplex_for_domain_t
detail::SimplexForDomain< K, Dom >::type simplex_for_domain_t
Definition simplex.hpp:232
sil::exterior::coboundary_index_t
typename detail::CoboundaryIndex< TagToAddToCochain, CochainTag >::type coboundary_index_t
Definition coboundary.hpp:82
sil::exterior::coboundary_t
typename detail::CoboundaryType< CochainType >::type coboundary_t
Definition coboundary.hpp:52
sil::misc::domain_contains
KOKKOS_FUNCTION bool domain_contains(ddc::DiscreteDomain< DDim... > dom, ddc::DiscreteElement< ODDim... > elem)
Definition domain_contains.hpp:35
sil::tensor::Tensor
Tensor(ddc::Chunk< ElementType, SupportType, Allocator >) -> Tensor< ElementType, SupportType, Kokkos::layout_right, typename Allocator::memory_space >
sil
The top-level namespace of SimiLie.
Definition csr.hpp:15