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 <array>
7
8#include <ddc/ddc.hpp>
9
10#include <similie/misc/are_all_same.hpp>
11#include <similie/misc/clamp_to_domain.hpp>
12#include <similie/misc/domain_contains.hpp>
13#include <similie/misc/filled_struct.hpp>
14#include <similie/misc/macros.hpp>
15#include <similie/misc/portable_stl.hpp>
16#include <similie/misc/select_from_type_seq.hpp>
17#include <similie/misc/specialization.hpp>
18#include <similie/misc/type_seq_conversion.hpp>
19#include <similie/tensor/antisymmetric_tensor.hpp>
20#include <similie/tensor/dummy_index.hpp>
21#include <similie/tensor/tensor_impl.hpp>
22
23#include <Kokkos_StdAlgorithms.hpp>
24
25#include "cochain.hpp"
26#include "cosimplex.hpp"
27
28
29namespace sil {
30
31namespace exterior {
32
33namespace detail {
34
35template <class T>
36struct CoboundaryType;
37
38template <
39 std::size_t K,
40 class... Tag,
41 class ElementType,
42 class LayoutStridedPolicy1,
43 class LayoutStridedPolicy2,
44 class ExecSpace>
45struct CoboundaryType<
46 Cochain<Chain<Simplex<K, Tag...>, LayoutStridedPolicy1, ExecSpace>,
47 ElementType,
48 LayoutStridedPolicy2>>
49{
50 using type = Cosimplex<Simplex<K + 1, Tag...>, ElementType>;
51};
52
53} // namespace detail
54
55template <misc::Specialization<Cochain> CochainType>
56using coboundary_t = typename detail::CoboundaryType<CochainType>::type;
57
58namespace detail {
59
60template <class TagToAddToCochain, class CochainTag>
61struct CoboundaryIndex;
62
63template <tensor::TensorNatIndex TagToAddToCochain, tensor::TensorNatIndex CochainTag>
64 requires(CochainTag::rank() == 0)
65struct CoboundaryIndex<TagToAddToCochain, CochainTag>
66{
67 using type = TagToAddToCochain;
68};
69
70template <tensor::TensorNatIndex TagToAddToCochain, tensor::TensorNatIndex CochainTag>
71 requires(CochainTag::rank() == 1)
72struct CoboundaryIndex<TagToAddToCochain, CochainTag>
73{
75};
76
77template <tensor::TensorNatIndex TagToAddToCochain, tensor::TensorNatIndex... Tag>
78struct CoboundaryIndex<TagToAddToCochain, tensor::TensorAntisymmetricIndex<Tag...>>
79{
80 using type = tensor::TensorAntisymmetricIndex<TagToAddToCochain, Tag...>;
81};
82
83} // namespace detail
84
85template <class TagToAddToCochain, class CochainTag>
86using coboundary_index_t = typename detail::CoboundaryIndex<TagToAddToCochain, CochainTag>::type;
87
88namespace detail {
89
90template <
91 tensor::TensorNatIndex TagToAddToCochain,
92 tensor::TensorIndex CochainTag,
94struct CoboundaryTensorType;
95
96template <
97 tensor::TensorNatIndex TagToAddToCochain,
98 tensor::TensorIndex CochainIndex,
99 class ElementType,
100 class... DDim,
101 class SupportType,
102 class MemorySpace>
103struct CoboundaryTensorType<
104 TagToAddToCochain,
105 CochainIndex,
106 tensor::Tensor<ElementType, ddc::DiscreteDomain<DDim...>, SupportType, MemorySpace>>
107{
108 static_assert(ddc::type_seq_contains_v<
109 ddc::detail::TypeSeq<CochainIndex>,
110 ddc::detail::TypeSeq<DDim...>>);
111 using type = tensor::Tensor<
112 ElementType,
113 ddc::replace_dim_of_t<
114 ddc::DiscreteDomain<DDim...>,
115 CochainIndex,
117 SupportType,
118 MemorySpace>;
119};
120
121} // namespace detail
122
123template <
124 tensor::TensorNatIndex TagToAddToCochain,
125 tensor::TensorIndex CochainTag,
128 typename detail::CoboundaryTensorType<TagToAddToCochain, CochainTag, TensorType>::type;
129
130namespace detail {
131
132template <misc::Specialization<Chain> ChainType>
133struct ComputeSimplex;
134
135template <std::size_t K, class... Tag, class LayoutStridedPolicy, class ExecSpace>
136struct ComputeSimplex<Chain<Simplex<K, Tag...>, LayoutStridedPolicy, ExecSpace>>
137{
138 KOKKOS_FUNCTION static Simplex<K + 1, Tag...> run(
139 Chain<Simplex<K, Tag...>, LayoutStridedPolicy, ExecSpace> const& chain)
140 {
141 ddc::DiscreteVector<Tag...> vect {
142 0 * ddc::type_seq_rank_v<Tag, ddc::detail::TypeSeq<Tag...>>...};
143 for (auto i = chain.begin(); i < chain.end(); ++i) {
144 vect = ddc::DiscreteVector<Tag...> {
145 (static_cast<bool>(vect.template get<Tag>())
146 || static_cast<bool>((*i).discrete_vector().template get<Tag>()))...};
147 }
148 return Simplex(
149 std::integral_constant<std::size_t, K + 1> {},
150 chain[0].discrete_element(), // This is an assumption on the structure of the chain, which is satisfied if it has been produced using boundary()
151 vect);
152 }
153};
154
155} // namespace detail
156
157template <class... Args>
159
160template <class... Args>
162
163template <misc::Specialization<Cochain> CochainType>
164struct Coboundary<CochainType>
165{
166 KOKKOS_FUNCTION static coboundary_t<CochainType>
167 run(CochainType
168 cochain) // Warning: only cochain.chain() produced using boundary() are supported
169 {
170 assert(cochain.size() == 2 * (cochain.dimension() + 1)
171 && "only cochain over the boundary of a single simplex is supported");
172
173 /* Commented because would require an additional buffer
174 assert(boundary(cochain.chain()) == boundary_t<typename CochainType::chain_type> {}
175 && "only cochain over the boundary of a single simplex is supported");
176 */
177
179 detail::ComputeSimplex<typename CochainType::chain_type>::run(cochain.chain()),
180 cochain.integrate());
181 }
182};
183
184namespace detail {
185
186template <tensor::TensorNatIndex Index, class Dom>
187struct NonSpectatorDimension;
188
189template <tensor::TensorNatIndex Index, class... DDim>
190struct NonSpectatorDimension<Index, ddc::DiscreteDomain<DDim...>>
191{
192 using type = ddc::cartesian_prod_t<std::conditional_t<
193 ddc::type_seq_contains_v<
194 ddc::detail::TypeSeq<typename DDim::continuous_dimension_type>,
195 typename Index::type_seq_dimensions>,
196 ddc::DiscreteDomain<DDim>,
197 ddc::DiscreteDomain<>>...>;
198};
199
200struct CoboundaryDummyIndex
201{
202};
203
204} // namespace detail
205
206template <tensor::TensorNatIndex TagToAddToCochain, tensor::TensorIndex CochainTag>
207struct Coboundary<TagToAddToCochain, CochainTag>
208{
209 template <
210 class CoboundaryTensorType,
211 class Evaluator,
212 class ChainType,
213 class LowerChainType,
214 class Elem>
215 KOKKOS_FUNCTION static void run(
216 CoboundaryTensorType coboundary_tensor,
217 Evaluator evaluator,
218 ChainType chain,
219 LowerChainType lower_chain,
220 Elem elem)
221 {
222 constexpr std::size_t BOUNDARY_SIZE = 2 * (CochainTag::rank() + 1);
223 std::array<typename LowerChainType::simplex_type, BOUNDARY_SIZE> simplex_boundary_alloc {};
224 std::array<double, BOUNDARY_SIZE> boundary_values_alloc {};
225 ddc::DiscreteDomain<detail::CoboundaryDummyIndex> const boundary_domain(
226 ddc::DiscreteElement<detail::CoboundaryDummyIndex>(0),
227 ddc::DiscreteVector<detail::CoboundaryDummyIndex>(BOUNDARY_SIZE));
228 ddc::ChunkSpan<
229 typename LowerChainType::simplex_type,
230 ddc::DiscreteDomain<detail::CoboundaryDummyIndex>,
231 Kokkos::layout_right,
232 typename CoboundaryTensorType::memory_space>
233 simplex_boundary(simplex_boundary_alloc.data(), boundary_domain);
234 ddc::ChunkSpan<
235 double,
236 ddc::DiscreteDomain<detail::CoboundaryDummyIndex>,
237 Kokkos::layout_right,
238 typename CoboundaryTensorType::memory_space>
239 boundary_values(boundary_values_alloc.data(), boundary_domain);
240
241 auto cochain = Cochain(chain, coboundary_tensor);
242 for (auto i = cochain.begin(); i < cochain.end(); ++i) {
243 typename ChainType::simplex_type
244 simplex(std::integral_constant<std::size_t, CochainTag::rank() + 1> {},
245 typename ChainType::simplex_type::discrete_element_type(elem),
246 (*i).discrete_vector());
247 auto boundary_chain = boundary(simplex_boundary.allocation_kokkos_view(), simplex);
248 for (auto j = boundary_chain.begin(); j < boundary_chain.end(); ++j) {
249 std::size_t const boundary_id
250 = Kokkos::Experimental::distance(boundary_chain.begin(), j);
251 if constexpr (
252 ddc::type_seq_size_v<ddc::type_seq_remove_t<
253 ddc::to_type_seq_t<Elem>,
254 ddc::to_type_seq_t<typename LowerChainType::simplex_type::
255 discrete_element_type>>>
256 == 0) {
257 boundary_values(ddc::DiscreteElement<detail::CoboundaryDummyIndex>(boundary_id))
258 = evaluator(
259 Elem((*j).discrete_element()),
260 ddc::DiscreteElement<CochainTag>(Kokkos::Experimental::distance(
261 lower_chain.begin(),
262 misc::detail::
263 find(lower_chain.begin(),
264 lower_chain.end(),
265 (*j).discrete_vector()))));
266 } else {
267 boundary_values(ddc::DiscreteElement<detail::CoboundaryDummyIndex>(boundary_id))
268 = evaluator(
269 Elem((*j).discrete_element(),
270 misc::select_from_type_seq<ddc::type_seq_remove_t<
271 ddc::to_type_seq_t<Elem>,
272 ddc::to_type_seq_t<
273 typename LowerChainType::simplex_type::
274 discrete_element_type>>>(elem)),
275 ddc::DiscreteElement<CochainTag>(Kokkos::Experimental::distance(
276 lower_chain.begin(),
277 misc::detail::
278 find(lower_chain.begin(),
279 lower_chain.end(),
280 (*j).discrete_vector()))));
281 }
282 }
283 Cochain cochain_boundary(boundary_chain, boundary_values.allocation_kokkos_view());
284 coboundary_tensor.mem(
286 Kokkos::Experimental::distance(cochain.begin(), i)))
287 = cochain_boundary.integrate();
288 }
289 }
290};
291
292template <tensor::TensorNatIndex TagToAddToCochain, tensor::TensorIndex CochainTag>
293struct TransposedCoboundary<TagToAddToCochain, CochainTag>
294{
295 template <
296 class CoboundaryTensorType,
297 class Evaluator,
298 class ChainType,
299 class LowerChainType,
300 class Elem>
301 KOKKOS_FUNCTION static void run(
302 CoboundaryTensorType coboundary_tensor,
303 Evaluator evaluator,
304 ChainType chain,
305 LowerChainType lower_chain,
306 Elem elem)
307 {
308 constexpr std::size_t BOUNDARY_SIZE = 2 * (CochainTag::rank() + 1);
309 std::array<typename LowerChainType::simplex_type, BOUNDARY_SIZE> simplex_boundary_alloc {};
310 std::array<double, BOUNDARY_SIZE> boundary_values_alloc {};
311 ddc::DiscreteDomain<detail::CoboundaryDummyIndex> const boundary_domain(
312 ddc::DiscreteElement<detail::CoboundaryDummyIndex>(0),
313 ddc::DiscreteVector<detail::CoboundaryDummyIndex>(BOUNDARY_SIZE));
314 ddc::ChunkSpan<
315 typename LowerChainType::simplex_type,
316 ddc::DiscreteDomain<detail::CoboundaryDummyIndex>,
317 Kokkos::layout_right,
318 typename CoboundaryTensorType::memory_space>
319 simplex_boundary(simplex_boundary_alloc.data(), boundary_domain);
320 ddc::ChunkSpan<
321 double,
322 ddc::DiscreteDomain<detail::CoboundaryDummyIndex>,
323 Kokkos::layout_right,
324 typename CoboundaryTensorType::memory_space>
325 boundary_values(boundary_values_alloc.data(), boundary_domain);
326
327 auto cochain = Cochain(chain, coboundary_tensor);
328 for (auto i = cochain.begin(); i < cochain.end(); ++i) {
329 typename ChainType::simplex_type
330 simplex(std::integral_constant<std::size_t, CochainTag::rank() + 1> {},
331 typename ChainType::simplex_type::discrete_element_type(elem),
332 (*i).discrete_vector());
333 auto boundary_chain = boundary(simplex_boundary.allocation_kokkos_view(), simplex);
334 for (auto j = boundary_chain.begin(); j < boundary_chain.end(); ++j) {
335 std::size_t const boundary_id
336 = Kokkos::Experimental::distance(boundary_chain.begin(), j);
337 auto sampled_face_vector = (*j).discrete_vector();
338 auto sampled_face_elem = (*j).discrete_element();
339 if (boundary_id >= CochainTag::rank() + 1) {
340 auto simplex_vector = (*i).discrete_vector();
341 for (std::size_t dim_id = 0;
342 dim_id < ddc::type_seq_size_v<ddc::to_type_seq_t<
343 typename ChainType::simplex_type::discrete_element_type>>;
344 ++dim_id) {
345 if (ddc::detail::array(simplex_vector)[dim_id] != 0
346 && ddc::detail::array(sampled_face_vector)[dim_id] == 0) {
347 ddc::detail::array(sampled_face_elem)[dim_id] -= 2;
348 break;
349 }
350 }
351 }
352 if constexpr (
353 ddc::type_seq_size_v<ddc::type_seq_remove_t<
354 ddc::to_type_seq_t<Elem>,
355 ddc::to_type_seq_t<typename LowerChainType::simplex_type::
356 discrete_element_type>>>
357 == 0) {
358 boundary_values(ddc::DiscreteElement<detail::CoboundaryDummyIndex>(boundary_id))
359 = evaluator(
360 Elem(sampled_face_elem),
361 ddc::DiscreteElement<CochainTag>(Kokkos::Experimental::distance(
362 lower_chain.begin(),
363 misc::detail::
364 find(lower_chain.begin(),
365 lower_chain.end(),
366 sampled_face_vector))));
367 } else {
368 boundary_values(ddc::DiscreteElement<detail::CoboundaryDummyIndex>(boundary_id))
369 = evaluator(
370 Elem(sampled_face_elem,
371 misc::select_from_type_seq<ddc::type_seq_remove_t<
372 ddc::to_type_seq_t<Elem>,
373 ddc::to_type_seq_t<
374 typename LowerChainType::simplex_type::
375 discrete_element_type>>>(elem)),
376 ddc::DiscreteElement<CochainTag>(Kokkos::Experimental::distance(
377 lower_chain.begin(),
378 misc::detail::
379 find(lower_chain.begin(),
380 lower_chain.end(),
381 sampled_face_vector))));
382 }
383 }
384 Cochain cochain_boundary(boundary_chain, boundary_values.allocation_kokkos_view());
385 coboundary_tensor.mem(
387 Kokkos::Experimental::distance(cochain.begin(), i)))
388 = -cochain_boundary.integrate();
389 }
390 }
391};
392
393template <misc::Specialization<Cochain> CochainType>
394KOKKOS_FUNCTION coboundary_t<CochainType> coboundary(CochainType cochain)
395{
396 return Coboundary<CochainType>::run(cochain);
397}
398
399template <
400 tensor::TensorNatIndex TagToAddToCochain,
401 tensor::TensorIndex CochainTag,
403 class ExecSpace>
405 ExecSpace const& exec_space,
407 TensorType tensor)
408{
409 ddc::DiscreteDomain batch_dom
410 = ddc::remove_dims_of<coboundary_index_t<TagToAddToCochain, CochainTag>>(
411 coboundary_tensor.domain());
412
413 auto chain = tangent_basis<
414 CochainTag::rank() + 1,
415 typename detail::NonSpectatorDimension<
416 TagToAddToCochain,
417 typename TensorType::non_indices_domain_t>::type>(exec_space);
418 auto lower_chain = tangent_basis<
419 CochainTag::rank(),
420 typename detail::NonSpectatorDimension<
421 TagToAddToCochain,
422 typename TensorType::non_indices_domain_t>::type>(exec_space);
423
424 SIMILIE_DEBUG_LOG("similie_compute_coboundary");
425 ddc::parallel_for_each(
426 "similie_compute_coboundary",
427 exec_space,
428 batch_dom,
429 KOKKOS_LAMBDA(typename decltype(batch_dom)::discrete_element_type elem) {
431 coboundary_tensor[elem],
432 // TODO this is an assumption on boundary condition (free boundary), needs to be generalized.
433 [&](auto sampled_elem, auto cochain_elem) {
434 auto const clamped_elem = misc::
435 clamp_to_domain(tensor.non_indices_domain(), sampled_elem);
436 return tensor.mem(clamped_elem, cochain_elem);
437 },
438 chain,
439 lower_chain,
440 elem);
441 });
442
443 return coboundary_tensor;
444}
445
446template <
447 tensor::TensorNatIndex TagToAddToCochain,
448 tensor::TensorIndex CochainTag,
450 class ExecSpace>
452 ExecSpace const& exec_space,
454 TensorType tensor)
455{
456 ddc::DiscreteDomain batch_dom
457 = ddc::remove_dims_of<coboundary_index_t<TagToAddToCochain, CochainTag>>(
458 coboundary_tensor.domain());
459
460 auto chain = tangent_basis<
461 CochainTag::rank() + 1,
462 typename detail::NonSpectatorDimension<
463 TagToAddToCochain,
464 typename TensorType::non_indices_domain_t>::type>(exec_space);
465 auto lower_chain = tangent_basis<
466 CochainTag::rank(),
467 typename detail::NonSpectatorDimension<
468 TagToAddToCochain,
469 typename TensorType::non_indices_domain_t>::type>(exec_space);
470
471 SIMILIE_DEBUG_LOG("similie_compute_transposed_coboundary");
472 ddc::parallel_for_each(
473 "similie_compute_transposed_coboundary",
474 exec_space,
475 batch_dom,
476 KOKKOS_LAMBDA(typename decltype(batch_dom)::discrete_element_type elem) {
478 coboundary_tensor[elem],
479 [&](auto sampled_elem, auto cochain_elem) {
480 if (!misc::domain_contains(tensor.non_indices_domain(), sampled_elem)) {
481 return 0.0;
482 }
483 return tensor.mem(sampled_elem, cochain_elem);
484 },
485 chain,
486 lower_chain,
487 elem);
488 });
489
490 return coboundary_tensor;
491}
492
493template <
494 tensor::TensorNatIndex TagToAddToCochain,
495 tensor::TensorIndex CochainTag,
497 class ExecSpace>
499 ExecSpace const& exec_space,
501 TensorType tensor)
502{
503 return coboundary<TagToAddToCochain, CochainTag>(exec_space, coboundary_tensor, tensor);
504}
505
506} // namespace exterior
507
508} // namespace sil
Chain class.
Definition chain.hpp:26
Cochain class.
Definition cochain.hpp:30
KOKKOS_FUNCTION element_type integrate() noexcept
Definition cochain.hpp:187
Simplex class.
Definition simplex.hpp:75
constexpr auto tangent_basis(ExecSpace const &exec_space)
Cochain(ChainType, TensorType) -> Cochain< ChainType, typename TensorType::value_type, ddc::detail::mdspan_to_kokkos_layout_t< typename TensorType::layout_type > >
Simplex(ddc::DiscreteElement< Tag... >, ddc::DiscreteVector< T... >) -> Simplex< sizeof...(T), Tag... >
KOKKOS_FUNCTION coboundary_t< CochainType > coboundary(CochainType cochain)
Chain(Head, Tail...) -> Chain< typename Head::value_type, typename Head::array_layout, typename Head::memory_space >
typename detail::CoboundaryTensorType< TagToAddToCochain, CochainTag, TensorType >::type coboundary_tensor_t
coboundary_tensor_t< TagToAddToCochain, CochainTag, TensorType > deriv(ExecSpace const &exec_space, coboundary_tensor_t< TagToAddToCochain, CochainTag, TensorType > coboundary_tensor, TensorType tensor)
coboundary_tensor_t< TagToAddToCochain, CochainTag, TensorType > transposed_coboundary(ExecSpace const &exec_space, coboundary_tensor_t< TagToAddToCochain, CochainTag, TensorType > coboundary_tensor, TensorType tensor)
KOKKOS_FUNCTION Chain< boundary_t< SimplexType >, typename AllocationType::array_layout, typename AllocationType::memory_space > boundary(AllocationType allocation, SimplexType simplex)
Definition boundary.hpp:142
typename detail::CoboundaryIndex< TagToAddToCochain, CochainTag >::type coboundary_index_t
typename detail::CoboundaryType< CochainType >::type coboundary_t
KOKKOS_FUNCTION Elem clamp_to_domain(BatchDomain const &batch_domain, Elem const &elem)
KOKKOS_FUNCTION bool domain_contains(ddc::DiscreteDomain< DDim... > dom, ddc::DiscreteElement< ODDim... > elem)
KOKKOS_FUNCTION auto select_from_type_seq(T t)
Tensor(ddc::Chunk< ElementType, SupportType, Allocator >) -> Tensor< ElementType, SupportType, Kokkos::layout_right, typename Allocator::memory_space >
The top-level namespace of SimiLie.
Definition csr.hpp:15
static KOKKOS_FUNCTION coboundary_t< CochainType > run(CochainType cochain)
static KOKKOS_FUNCTION void run(CoboundaryTensorType coboundary_tensor, Evaluator evaluator, ChainType chain, LowerChainType lower_chain, Elem elem)
static KOKKOS_FUNCTION void run(CoboundaryTensorType coboundary_tensor, Evaluator evaluator, ChainType chain, LowerChainType lower_chain, Elem elem)