SimiLie
Loading...
Searching...
No Matches
coboundary.hpp
1// SPDX-FileCopyrightText: 2024 Baptiste Legouix
2// SPDX-License-Identifier: AGPL-3.0-or-later
3
4#pragma once
5
6#include <array>
7#include <utility>
8
9#include <ddc/ddc.hpp>
10
11#include <similie/misc/are_all_same.hpp>
12#include <similie/misc/filled_struct.hpp>
13#include <similie/misc/macros.hpp>
14#include <similie/misc/portable_stl.hpp>
15#include <similie/misc/select_from_type_seq.hpp>
16#include <similie/misc/specialization.hpp>
17#include <similie/misc/type_seq_conversion.hpp>
18#include <similie/tensor/antisymmetric_tensor.hpp>
19#include <similie/tensor/dummy_index.hpp>
20#include <similie/tensor/owning_tensor.hpp>
21#include <similie/tensor/tensor_impl.hpp>
22
23#include <Kokkos_StdAlgorithms.hpp>
24
25#include "boundary.hpp"
26#include "cochain.hpp"
27#include "cosimplex.hpp"
28#include "evaluators.hpp"
29
30
31namespace sil {
32
33namespace exterior {
34
35namespace detail {
36
37template <class T>
38struct CoboundaryType;
39
40template <
41 std::size_t K,
42 class... Tag,
43 class ElementType,
44 class LayoutStridedPolicy1,
45 class LayoutStridedPolicy2,
46 class ExecSpace>
47struct CoboundaryType<
48 Cochain<Chain<Simplex<K, Tag...>, LayoutStridedPolicy1, ExecSpace>,
49 ElementType,
50 LayoutStridedPolicy2>>
51{
52 using type = Cosimplex<Simplex<K + 1, Tag...>, ElementType>;
53};
54
55} // namespace detail
56
57template <misc::Specialization<Cochain> CochainType>
58using coboundary_t = typename detail::CoboundaryType<CochainType>::type;
59
60namespace detail {
61
62template <class TagToAddToCochain, class CochainTag>
63struct CoboundaryIndex;
64
65template <tensor::TensorNatIndex TagToAddToCochain, tensor::TensorNatIndex CochainTag>
66 requires(CochainTag::rank() == 0)
67struct CoboundaryIndex<TagToAddToCochain, CochainTag>
68{
69 using type = TagToAddToCochain;
70};
71
72template <tensor::TensorNatIndex TagToAddToCochain, tensor::TensorNatIndex CochainTag>
73 requires(CochainTag::rank() == 1)
74struct CoboundaryIndex<TagToAddToCochain, CochainTag>
75{
76 using type = tensor::TensorAntisymmetricIndex<TagToAddToCochain, CochainTag>;
77};
78
79template <tensor::TensorNatIndex TagToAddToCochain, tensor::TensorNatIndex... Tag>
80struct CoboundaryIndex<TagToAddToCochain, tensor::TensorAntisymmetricIndex<Tag...>>
81{
82 using type = tensor::TensorAntisymmetricIndex<TagToAddToCochain, Tag...>;
83};
84
85} // namespace detail
86
87template <class TagToAddToCochain, class CochainTag>
88using coboundary_index_t = typename detail::CoboundaryIndex<TagToAddToCochain, CochainTag>::type;
89
90namespace detail {
91
92template <
93 tensor::TensorNatIndex TagToAddToCochain,
94 tensor::TensorIndex CochainTag,
95 misc::Specialization<tensor::Tensor> TensorType>
96struct CoboundaryTensorType;
97
98template <
99 tensor::TensorNatIndex TagToAddToCochain,
100 tensor::TensorIndex CochainIndex,
101 class ElementType,
102 class... DDim,
103 class SupportType,
104 class MemorySpace>
105struct CoboundaryTensorType<
106 TagToAddToCochain,
107 CochainIndex,
108 tensor::Tensor<ElementType, ddc::DiscreteDomain<DDim...>, SupportType, MemorySpace>>
109{
110 static_assert(ddc::type_seq_contains_v<
111 ddc::detail::TypeSeq<CochainIndex>,
112 ddc::detail::TypeSeq<DDim...>>);
113 using type = tensor::Tensor<
114 ElementType,
115 ddc::replace_dim_of_t<
116 ddc::DiscreteDomain<DDim...>,
117 CochainIndex,
118 coboundary_index_t<TagToAddToCochain, CochainIndex>>,
119 SupportType,
120 MemorySpace>;
121};
122
123} // namespace detail
124
125template <
126 tensor::TensorNatIndex TagToAddToCochain,
127 tensor::TensorIndex CochainTag,
128 misc::Specialization<tensor::Tensor> TensorType>
129using coboundary_tensor_t =
130 typename detail::CoboundaryTensorType<TagToAddToCochain, CochainTag, TensorType>::type;
131
132namespace detail {
133
134template <misc::Specialization<Chain> ChainType>
135struct ComputeSimplex;
136
137template <std::size_t K, class... Tag, class LayoutStridedPolicy, class ExecSpace>
138struct ComputeSimplex<Chain<Simplex<K, Tag...>, LayoutStridedPolicy, ExecSpace>>
139{
140 KOKKOS_FUNCTION static Simplex<K + 1, Tag...> run(
141 Chain<Simplex<K, Tag...>, LayoutStridedPolicy, ExecSpace> const& chain)
142 {
143 ddc::DiscreteVector<Tag...> vect {
144 0 * ddc::type_seq_rank_v<Tag, ddc::detail::TypeSeq<Tag...>>...};
145 for (auto i = chain.begin(); i < chain.end(); ++i) {
146 vect = ddc::DiscreteVector<Tag...> {
147 (static_cast<bool>(vect.template get<Tag>())
148 || static_cast<bool>((*i).discrete_vector().template get<Tag>()))...};
149 }
150 // This assumes that the chain has been produced using boundary().
151 return Simplex(
152 std::integral_constant<std::size_t, K + 1> {},
153 (*chain.begin())
154 .discrete_element(), // This is an assumption on the structure of the chain, which is satisfied if it has been produced using boundary()
155 vect);
156 }
157};
158
159} // namespace detail
160
161template <class... Args>
162struct Coboundary;
163
164template <class... Args>
165struct TransposedCoboundary;
166
167template <misc::Specialization<Cochain> CochainType>
168struct Coboundary<CochainType>
169{
170 KOKKOS_FUNCTION static coboundary_t<CochainType>
171 run(CochainType
172 cochain) // Warning: only cochain.chain() produced using boundary() are supported
173 {
174 assert(cochain.size() == 2 * (cochain.dimension() + 1)
175 && "only cochain over the boundary of a single simplex is supported");
176
177 /* Commented because would require an additional buffer
178 assert(boundary(cochain.chain()) == boundary_t<typename CochainType::chain_type> {}
179 && "only cochain over the boundary of a single simplex is supported");
180 */
181
182 return coboundary_t<CochainType>(
183 detail::ComputeSimplex<typename CochainType::chain_type>::run(cochain.chain()),
184 cochain.integrate());
185 }
186};
187
188namespace detail {
189
190template <tensor::TensorNatIndex Index, class Dom>
191struct NonSpectatorDimension;
192
193template <tensor::TensorNatIndex Index, class... DDim>
194struct NonSpectatorDimension<Index, ddc::DiscreteDomain<DDim...>>
195{
196 using type = ddc::cartesian_prod_t<std::conditional_t<
197 ddc::type_seq_contains_v<
198 ddc::detail::TypeSeq<typename DDim::continuous_dimension_type>,
199 typename Index::type_seq_dimensions>,
200 ddc::DiscreteDomain<DDim>,
201 ddc::DiscreteDomain<>>...>;
202};
203
204struct CoboundaryDummyIndex
205{
206};
207
208template <class MemorySpace = Kokkos::HostSpace, class TensorIndex, class... DDims>
209KOKKOS_FUNCTION auto make_stencil(ddc::DiscreteElement<DDims...> front)
210 -> sil::tensor::OwningTensor<
211 sil::tensor::Tensor<
212 double,
213 ddc::DiscreteDomain<DDims..., TensorIndex>,
214 Kokkos::layout_right,
215 MemorySpace>,
216 std::array<
217 double,
218 (sizeof...(DDims) == 0 ? 1UL : (1UL << sizeof...(DDims)))
219 * TensorIndex::mem_size()>>
220{
221 std::array<
222 double,
223 (sizeof...(DDims) == 0 ? 1UL : (1UL << sizeof...(DDims))) * TensorIndex::mem_size()>
224 storage {};
225 ddc::DiscreteDomain<DDims..., TensorIndex> const
226 domain(ddc::DiscreteDomain<DDims...>(
227 front,
228 typename ddc::DiscreteDomain<DDims...>::discrete_vector_type(
229 ddc::DiscreteVector<DDims>(2)...)),
230 ddc::DiscreteDomain<TensorIndex>(
231 ddc::DiscreteElement<TensorIndex>(0),
232 ddc::DiscreteVector<TensorIndex>(TensorIndex::mem_size())));
233 ddc::ChunkSpan<
234 double,
235 ddc::DiscreteDomain<DDims..., TensorIndex>,
236 Kokkos::layout_right,
237 MemorySpace>
238 span(storage.data(), domain);
239 sil::tensor::Tensor<
240 double,
241 ddc::DiscreteDomain<DDims..., TensorIndex>,
242 Kokkos::layout_right,
243 MemorySpace>
244 tensor(span);
245 return sil::tensor::OwningTensor<
246 sil::tensor::Tensor<
247 double,
248 ddc::DiscreteDomain<DDims..., TensorIndex>,
249 Kokkos::layout_right,
250 MemorySpace>,
251 std::array<
252 double,
253 (sizeof...(DDims) == 0 ? 1UL : (1UL << sizeof...(DDims)))
254 * TensorIndex::mem_size()>>(tensor, std::move(storage));
255}
256
257template <class Elem>
258KOKKOS_FUNCTION Elem decrement_all(Elem elem)
259{
260 Elem shifted = elem;
261 constexpr std::size_t RANK = ddc::type_seq_size_v<ddc::to_type_seq_t<Elem>>;
262 for (std::size_t dim_id = 0; dim_id < RANK; ++dim_id) {
263 ddc::detail::array(shifted)[dim_id] -= 1;
264 }
265 return shifted;
266}
267
268template <class LowerChainType, class VectorType>
269KOKKOS_FUNCTION auto find_discrete_vector(
270 LowerChainType const& lower_chain,
271 VectorType const& vector)
272{
273 for (auto i = lower_chain.begin(); i < lower_chain.end(); ++i) {
274 if (i->discrete_vector() == vector) {
275 return i;
276 }
277 }
278 return lower_chain.end();
279}
280
281} // namespace detail
282
283template <tensor::TensorNatIndex TagToAddToCochain, tensor::TensorIndex CochainTag>
284struct Coboundary<TagToAddToCochain, CochainTag>
285{
286 template <class Evaluator, class ChainType, class LowerChainType, class Elem, class NaturalElem>
287 KOKKOS_FUNCTION static auto value(
288 Evaluator evaluator,
289 ChainType chain,
290 LowerChainType lower_chain,
291 Elem elem,
292 NaturalElem natural_elem)
293 {
294 std::size_t const stored_component_id
295 = coboundary_index_t<TagToAddToCochain, CochainTag>::access_id_to_mem_id(
296 natural_elem
297 .template uid<coboundary_index_t<TagToAddToCochain, CochainTag>>());
298 typename ChainType::simplex_type
299 simplex(std::integral_constant<std::size_t, CochainTag::rank() + 1> {},
300 typename ChainType::simplex_type::discrete_element_type(elem),
301 chain[stored_component_id].discrete_vector());
302 auto boundary_chain = boundary<typename ChainType::memory_space>(simplex);
303
304 typename LowerChainType::simplex_type::discrete_element_type front(
305 boundary_chain.begin()->discrete_element());
306 for (auto j = boundary_chain.begin(); j < boundary_chain.end(); ++j) {
307 for (std::size_t dim_id = 0;
308 dim_id < ddc::type_seq_size_v<ddc::to_type_seq_t<
309 typename LowerChainType::simplex_type::discrete_element_type>>;
310 ++dim_id) {
311 ddc::detail::array(front)[dim_id] = std::
312 min(ddc::detail::array(front)[dim_id],
313 ddc::detail::array(j->discrete_element())[dim_id]);
314 }
315 }
316
317 auto stencil = detail::make_stencil<typename ChainType::memory_space, CochainTag>(front);
318 ddc::device_for_each(stencil.domain(), [&](auto stencil_elem) {
319 auto basis_stencil
320 = detail::make_stencil<typename ChainType::memory_space, CochainTag>(front);
321 basis_stencil.mem(stencil_elem) = 1.0;
322
323 [[maybe_unused]] sil::tensor::TensorAccessor<
324 coboundary_index_t<TagToAddToCochain, CochainTag>> accessor;
325 std::array<double, coboundary_index_t<TagToAddToCochain, CochainTag>::access_size()>
326 output_storage {};
327 ddc::ChunkSpan<
328 double,
329 ddc::DiscreteDomain<coboundary_index_t<TagToAddToCochain, CochainTag>>,
330 Kokkos::layout_right,
331 typename ChainType::memory_space>
332 output_span(output_storage.data(), accessor.domain());
333 sil::tensor::Tensor output_tensor(output_span);
334
335 auto basis_evaluator = [&](auto sampled_elem, auto cochain_elem) {
336 auto basis_sampler = [&](auto basis_elem, auto basis_cochain_elem) {
337 if (!basis_stencil.non_indices_domain().contains(basis_elem)) {
338 return 0.0;
339 }
340 return basis_stencil(
341 basis_stencil.access_element(basis_elem, basis_cochain_elem));
342 };
343 return evaluator.value(basis_sampler, sampled_elem, cochain_elem);
344 };
345
346 run(output_tensor, basis_evaluator, chain, lower_chain, elem);
347 stencil.mem(stencil_elem) = output_tensor(natural_elem);
348 });
349
350 return stencil;
351 }
352
353 template <
354 class CoboundaryTensorType,
355 class Evaluator,
356 class ChainType,
357 class LowerChainType,
358 class Elem>
359 KOKKOS_FUNCTION static void run(
360 CoboundaryTensorType coboundary_tensor,
361 Evaluator evaluator,
362 ChainType chain,
363 LowerChainType lower_chain,
364 Elem elem)
365 {
366 constexpr std::size_t BOUNDARY_SIZE = 2 * (CochainTag::rank() + 1);
367 std::array<double, BOUNDARY_SIZE> boundary_values_alloc {};
368 ddc::DiscreteDomain<detail::CoboundaryDummyIndex> const boundary_domain(
369 ddc::DiscreteElement<detail::CoboundaryDummyIndex>(0),
370 ddc::DiscreteVector<detail::CoboundaryDummyIndex>(BOUNDARY_SIZE));
371 ddc::ChunkSpan<
372 double,
373 ddc::DiscreteDomain<detail::CoboundaryDummyIndex>,
374 Kokkos::layout_right,
375 typename CoboundaryTensorType::memory_space>
376 boundary_values(boundary_values_alloc.data(), boundary_domain);
377
378 typename ChainType::simplex_type::discrete_element_type const elem_on_chain
379 = misc::select_from_type_seq<ddc::to_type_seq_t<
380 typename ChainType::simplex_type::discrete_element_type>>(elem);
381 for (auto i = chain.begin(); i < chain.end(); ++i) {
382 std::size_t const chain_id = Kokkos::Experimental::distance(chain.begin(), i);
383 typename ChainType::simplex_type
384 simplex(std::integral_constant<std::size_t, CochainTag::rank() + 1> {},
385 elem_on_chain,
386 i->discrete_vector());
387 auto boundary_chain = boundary<typename CoboundaryTensorType::memory_space>(simplex);
388 for (auto j = boundary_chain.begin(); j < boundary_chain.end(); ++j) {
389 std::size_t const boundary_id
390 = Kokkos::Experimental::distance(boundary_chain.begin(), j);
391 if constexpr (
392 ddc::type_seq_size_v<ddc::type_seq_remove_t<
393 ddc::to_type_seq_t<Elem>,
394 ddc::to_type_seq_t<typename LowerChainType::simplex_type::
395 discrete_element_type>>>
396 == 0) {
397 boundary_values(ddc::DiscreteElement<detail::CoboundaryDummyIndex>(boundary_id))
398 = evaluator(
399 Elem(j->discrete_element()),
400 ddc::DiscreteElement<CochainTag>(Kokkos::Experimental::distance(
401 lower_chain.begin(),
402 detail::find_discrete_vector(
403 lower_chain,
404 j->discrete_vector()))));
405 } else {
406 boundary_values(ddc::DiscreteElement<detail::CoboundaryDummyIndex>(boundary_id))
407 = evaluator(
408 Elem(j->discrete_element(),
409 misc::select_from_type_seq<ddc::type_seq_remove_t<
410 ddc::to_type_seq_t<Elem>,
411 ddc::to_type_seq_t<
412 typename LowerChainType::simplex_type::
413 discrete_element_type>>>(elem)),
414 ddc::DiscreteElement<CochainTag>(Kokkos::Experimental::distance(
415 lower_chain.begin(),
416 detail::find_discrete_vector(
417 lower_chain,
418 j->discrete_vector()))));
419 }
420 }
421 Cochain cochain_boundary(boundary_chain, boundary_values.allocation_kokkos_view());
422 coboundary_tensor.mem(
423 ddc::DiscreteElement<coboundary_index_t<TagToAddToCochain, CochainTag>>(
424 chain_id))
425 = cochain_boundary.integrate();
426 }
427 }
428};
429
430template <tensor::TensorNatIndex TagToAddToCochain, tensor::TensorIndex CochainTag>
431struct TransposedCoboundary<TagToAddToCochain, CochainTag>
432{
433 template <class Evaluator, class ChainType, class LowerChainType, class Elem, class NaturalElem>
434 KOKKOS_FUNCTION static auto value(
435 Evaluator evaluator,
436 ChainType chain,
437 LowerChainType lower_chain,
438 Elem elem,
439 NaturalElem natural_elem)
440 {
441 constexpr std::size_t BOUNDARY_SIZE = 2 * (CochainTag::rank() + 1);
442 std::size_t const stored_component_id
443 = coboundary_index_t<TagToAddToCochain, CochainTag>::access_id_to_mem_id(
444 natural_elem
445 .template uid<coboundary_index_t<TagToAddToCochain, CochainTag>>());
446 typename ChainType::simplex_type
447 simplex(std::integral_constant<std::size_t, CochainTag::rank() + 1> {},
448 typename ChainType::simplex_type::discrete_element_type(elem),
449 chain[stored_component_id].discrete_vector());
450 auto boundary_chain = boundary<typename ChainType::memory_space>(simplex);
451
452 typename LowerChainType::simplex_type::discrete_element_type front(
453 boundary_chain.begin()->discrete_element());
454 for (auto j = boundary_chain.begin(); j < boundary_chain.end(); ++j) {
455 auto sampled_face_elem = j->discrete_element();
456 std::size_t const boundary_id
457 = Kokkos::Experimental::distance(boundary_chain.begin(), j);
458 if (boundary_id >= CochainTag::rank() + 1) {
459 for (std::size_t dim_id = 0;
460 dim_id < ddc::type_seq_size_v<ddc::to_type_seq_t<
461 typename LowerChainType::simplex_type::discrete_element_type>>;
462 ++dim_id) {
463 if (ddc::detail::array(chain[stored_component_id].discrete_vector())[dim_id]
464 != 0
465 && ddc::detail::array(j->discrete_vector())[dim_id] == 0) {
466 ddc::detail::array(sampled_face_elem)[dim_id] -= 2;
467 break;
468 }
469 }
470 }
471 for (std::size_t dim_id = 0;
472 dim_id < ddc::type_seq_size_v<ddc::to_type_seq_t<
473 typename LowerChainType::simplex_type::discrete_element_type>>;
474 ++dim_id) {
475 ddc::detail::array(front)[dim_id] = std::
476 min(ddc::detail::array(front)[dim_id],
477 ddc::detail::array(sampled_face_elem)[dim_id]);
478 }
479 }
480
481 auto stencil = detail::make_stencil<typename ChainType::memory_space, CochainTag>(front);
482 ddc::device_for_each(stencil.domain(), [&](auto stencil_elem) {
483 auto basis_stencil
484 = detail::make_stencil<typename ChainType::memory_space, CochainTag>(front);
485 basis_stencil.mem(stencil_elem) = 1.0;
486
487 [[maybe_unused]] sil::tensor::TensorAccessor<
488 coboundary_index_t<TagToAddToCochain, CochainTag>> accessor;
489 std::array<double, coboundary_index_t<TagToAddToCochain, CochainTag>::access_size()>
490 output_storage {};
491 ddc::ChunkSpan<
492 double,
493 ddc::DiscreteDomain<coboundary_index_t<TagToAddToCochain, CochainTag>>,
494 Kokkos::layout_right,
495 typename ChainType::memory_space>
496 output_span(output_storage.data(), accessor.domain());
497 sil::tensor::Tensor output_tensor(output_span);
498
499 auto basis_evaluator = [&](auto sampled_elem, auto cochain_elem) {
500 auto basis_sampler = [&](auto basis_elem, auto basis_cochain_elem) {
501 if (!basis_stencil.non_indices_domain().contains(basis_elem)) {
502 return 0.0;
503 }
504 return basis_stencil(
505 basis_stencil.access_element(basis_elem, basis_cochain_elem));
506 };
507 return evaluator.value(basis_sampler, sampled_elem, cochain_elem);
508 };
509
510 run(output_tensor, basis_evaluator, chain, lower_chain, elem);
511 stencil.mem(stencil_elem) = output_tensor(natural_elem);
512 });
513
514 return stencil;
515 }
516
517 template <
518 class CoboundaryTensorType,
519 class Evaluator,
520 class ChainType,
521 class LowerChainType,
522 class Elem>
523 KOKKOS_FUNCTION static void run(
524 CoboundaryTensorType coboundary_tensor,
525 Evaluator evaluator,
526 ChainType chain,
527 LowerChainType lower_chain,
528 Elem elem)
529 {
530 constexpr std::size_t BOUNDARY_SIZE = 2 * (CochainTag::rank() + 1);
531 std::array<double, BOUNDARY_SIZE> boundary_values_alloc {};
532 ddc::DiscreteDomain<detail::CoboundaryDummyIndex> const boundary_domain(
533 ddc::DiscreteElement<detail::CoboundaryDummyIndex>(0),
534 ddc::DiscreteVector<detail::CoboundaryDummyIndex>(BOUNDARY_SIZE));
535 ddc::ChunkSpan<
536 double,
537 ddc::DiscreteDomain<detail::CoboundaryDummyIndex>,
538 Kokkos::layout_right,
539 typename CoboundaryTensorType::memory_space>
540 boundary_values(boundary_values_alloc.data(), boundary_domain);
541
542 typename ChainType::simplex_type::discrete_element_type const elem_on_chain
543 = misc::select_from_type_seq<ddc::to_type_seq_t<
544 typename ChainType::simplex_type::discrete_element_type>>(elem);
545 for (auto i = chain.begin(); i < chain.end(); ++i) {
546 std::size_t const chain_id = Kokkos::Experimental::distance(chain.begin(), i);
547 auto simplex_vector = i->discrete_vector();
548 typename ChainType::simplex_type
549 simplex(std::integral_constant<std::size_t, CochainTag::rank() + 1> {},
550 elem_on_chain,
551 simplex_vector);
552 auto boundary_chain = boundary<typename CoboundaryTensorType::memory_space>(simplex);
553 for (auto j = boundary_chain.begin(); j < boundary_chain.end(); ++j) {
554 std::size_t const boundary_id
555 = Kokkos::Experimental::distance(boundary_chain.begin(), j);
556 auto sampled_face_elem = j->discrete_element();
557 if (boundary_id >= CochainTag::rank() + 1) {
558 for (std::size_t dim_id = 0;
559 dim_id < ddc::type_seq_size_v<ddc::to_type_seq_t<
560 typename ChainType::simplex_type::discrete_element_type>>;
561 ++dim_id) {
562 if (ddc::detail::array(simplex_vector)[dim_id] != 0
563 && ddc::detail::array(j->discrete_vector())[dim_id] == 0) {
564 ddc::detail::array(sampled_face_elem)[dim_id] -= 2;
565 break;
566 }
567 }
568 }
569 if constexpr (
570 ddc::type_seq_size_v<ddc::type_seq_remove_t<
571 ddc::to_type_seq_t<Elem>,
572 ddc::to_type_seq_t<typename LowerChainType::simplex_type::
573 discrete_element_type>>>
574 == 0) {
575 boundary_values(ddc::DiscreteElement<detail::CoboundaryDummyIndex>(boundary_id))
576 = evaluator(
577 Elem(sampled_face_elem),
578 ddc::DiscreteElement<CochainTag>(Kokkos::Experimental::distance(
579 lower_chain.begin(),
580 detail::find_discrete_vector(
581 lower_chain,
582 j->discrete_vector()))));
583 } else {
584 boundary_values(ddc::DiscreteElement<detail::CoboundaryDummyIndex>(boundary_id))
585 = evaluator(
586 Elem(sampled_face_elem,
587 misc::select_from_type_seq<ddc::type_seq_remove_t<
588 ddc::to_type_seq_t<Elem>,
589 ddc::to_type_seq_t<
590 typename LowerChainType::simplex_type::
591 discrete_element_type>>>(elem)),
592 ddc::DiscreteElement<CochainTag>(Kokkos::Experimental::distance(
593 lower_chain.begin(),
594 detail::find_discrete_vector(
595 lower_chain,
596 j->discrete_vector()))));
597 }
598 }
599 Cochain cochain_boundary(boundary_chain, boundary_values.allocation_kokkos_view());
600 coboundary_tensor.mem(
601 ddc::DiscreteElement<coboundary_index_t<TagToAddToCochain, CochainTag>>(
602 chain_id))
603 = -cochain_boundary.integrate();
604 }
605 }
606};
607
608template <misc::Specialization<Cochain> CochainType>
609KOKKOS_FUNCTION coboundary_t<CochainType> coboundary(CochainType cochain)
610{
611 return Coboundary<CochainType>::run(cochain);
612}
613
614template <
615 tensor::TensorNatIndex TagToAddToCochain,
616 tensor::TensorIndex CochainTag,
617 misc::Specialization<tensor::Tensor> TensorType,
618 class ExecSpace>
619coboundary_tensor_t<TagToAddToCochain, CochainTag, TensorType> coboundary(
620 ExecSpace const& exec_space,
621 coboundary_tensor_t<TagToAddToCochain, CochainTag, TensorType> coboundary_tensor,
622 TensorType tensor)
623{
624 ddc::DiscreteDomain batch_dom
625 = ddc::remove_dims_of<coboundary_index_t<TagToAddToCochain, CochainTag>>(
626 coboundary_tensor.domain());
627
628 auto chain = tangent_basis<
629 CochainTag::rank() + 1,
630 typename detail::NonSpectatorDimension<
631 TagToAddToCochain,
632 typename TensorType::non_indices_domain_t>::type>(exec_space);
633 auto lower_chain = tangent_basis<
634 CochainTag::rank(),
635 typename detail::NonSpectatorDimension<
636 TagToAddToCochain,
637 typename TensorType::non_indices_domain_t>::type>(exec_space);
638
639 SIMILIE_DEBUG_LOG("similie_compute_coboundary");
640 ddc::parallel_for_each(
641 "similie_compute_coboundary",
642 exec_space,
643 batch_dom,
644 KOKKOS_LAMBDA(typename decltype(batch_dom)::discrete_element_type elem) {
645 Coboundary<TagToAddToCochain, CochainTag>::
646 run(coboundary_tensor[elem],
647 detail::ClampedTensorEvaluator<TensorType> {tensor},
648 chain,
649 lower_chain,
650 elem);
651 });
652
653 return coboundary_tensor;
654}
655
656template <
657 tensor::TensorNatIndex TagToAddToCochain,
658 tensor::TensorIndex CochainTag,
659 misc::Specialization<tensor::Tensor> TensorType,
660 class ExecSpace>
661coboundary_tensor_t<TagToAddToCochain, CochainTag, TensorType> transposed_coboundary(
662 ExecSpace const& exec_space,
663 coboundary_tensor_t<TagToAddToCochain, CochainTag, TensorType> coboundary_tensor,
664 TensorType tensor)
665{
666 ddc::DiscreteDomain batch_dom
667 = ddc::remove_dims_of<coboundary_index_t<TagToAddToCochain, CochainTag>>(
668 coboundary_tensor.domain());
669
670 auto chain = tangent_basis<
671 CochainTag::rank() + 1,
672 typename detail::NonSpectatorDimension<
673 TagToAddToCochain,
674 typename TensorType::non_indices_domain_t>::type>(exec_space);
675 auto lower_chain = tangent_basis<
676 CochainTag::rank(),
677 typename detail::NonSpectatorDimension<
678 TagToAddToCochain,
679 typename TensorType::non_indices_domain_t>::type>(exec_space);
680
681 SIMILIE_DEBUG_LOG("similie_compute_transposed_coboundary");
682 ddc::parallel_for_each(
683 "similie_compute_transposed_coboundary",
684 exec_space,
685 batch_dom,
686 KOKKOS_LAMBDA(typename decltype(batch_dom)::discrete_element_type elem) {
687 TransposedCoboundary<TagToAddToCochain, CochainTag>::
688 run(coboundary_tensor[elem],
689 detail::ZeroOutsideTensorEvaluator<TensorType> {tensor},
690 chain,
691 lower_chain,
692 elem);
693 });
694
695 return coboundary_tensor;
696}
697
698template <
699 tensor::TensorNatIndex TagToAddToCochain,
700 tensor::TensorIndex CochainTag,
701 misc::Specialization<tensor::Tensor> TensorType,
702 class ExecSpace>
703coboundary_tensor_t<TagToAddToCochain, CochainTag, TensorType> deriv(
704 ExecSpace const& exec_space,
705 coboundary_tensor_t<TagToAddToCochain, CochainTag, TensorType> coboundary_tensor,
706 TensorType tensor)
707{
708 return coboundary<TagToAddToCochain, CochainTag>(exec_space, coboundary_tensor, tensor);
709}
710
711} // namespace exterior
712
713} // namespace sil
sil::exterior::Chain
Chain class.
Definition chain.hpp:26
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::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::tangent_basis
KOKKOS_FUNCTION constexpr auto tangent_basis(Elem elem)
Definition local_chain.hpp:443
sil::exterior::coboundary
KOKKOS_FUNCTION coboundary_t< CochainType > coboundary(CochainType cochain)
Definition coboundary.hpp:609
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:129
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:703
sil::exterior::transposed_coboundary
coboundary_tensor_t< TagToAddToCochain, CochainTag, TensorType > transposed_coboundary(ExecSpace const &exec_space, coboundary_tensor_t< TagToAddToCochain, CochainTag, TensorType > coboundary_tensor, TensorType tensor)
Definition coboundary.hpp:661
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:169
sil::exterior::coboundary_index_t
typename detail::CoboundaryIndex< TagToAddToCochain, CochainTag >::type coboundary_index_t
Definition coboundary.hpp:88
sil::exterior::coboundary_t
typename detail::CoboundaryType< CochainType >::type coboundary_t
Definition coboundary.hpp:58
sil::misc::select_from_type_seq
KOKKOS_FUNCTION auto select_from_type_seq(T t)
Definition select_from_type_seq.hpp:30
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
sil::exterior::Coboundary< CochainType >::run
static KOKKOS_FUNCTION coboundary_t< CochainType > run(CochainType cochain)
Definition coboundary.hpp:171
sil::exterior::Coboundary< TagToAddToCochain, CochainTag >::run
static KOKKOS_FUNCTION void run(CoboundaryTensorType coboundary_tensor, Evaluator evaluator, ChainType chain, LowerChainType lower_chain, Elem elem)
Definition coboundary.hpp:359
sil::exterior::Coboundary< TagToAddToCochain, CochainTag >::value
static KOKKOS_FUNCTION auto value(Evaluator evaluator, ChainType chain, LowerChainType lower_chain, Elem elem, NaturalElem natural_elem)
Definition coboundary.hpp:287
sil::exterior::TransposedCoboundary< TagToAddToCochain, CochainTag >::value
static KOKKOS_FUNCTION auto value(Evaluator evaluator, ChainType chain, LowerChainType lower_chain, Elem elem, NaturalElem natural_elem)
Definition coboundary.hpp:434
sil::exterior::TransposedCoboundary< TagToAddToCochain, CochainTag >::run
static KOKKOS_FUNCTION void run(CoboundaryTensorType coboundary_tensor, Evaluator evaluator, ChainType chain, LowerChainType lower_chain, Elem elem)
Definition coboundary.hpp:523