SimiLie
Loading...
Searching...
No Matches
hodge_star.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/factorial.hpp>
9#include <similie/misc/specialization.hpp>
10#include <similie/misc/type_seq_conversion.hpp>
11#include <similie/tensor/antisymmetric_tensor.hpp>
12#include <similie/tensor/character.hpp>
13#include <similie/tensor/determinant.hpp>
14#include <similie/tensor/full_tensor.hpp>
15#include <similie/tensor/levi_civita_tensor.hpp>
16#include <similie/tensor/metric.hpp>
17#include <similie/tensor/prime.hpp>
18
19namespace sil {
20
21namespace exterior {
22
23template <
24 misc::Specialization<ddc::detail::TypeSeq> Indices1,
25 misc::Specialization<ddc::detail::TypeSeq> Indices2>
26using hodge_star_domain_t
27 = ddc::detail::convert_type_seq_to_discrete_domain_t<ddc::type_seq_merge_t<
28 ddc::detail::TypeSeq<misc::convert_type_seq_to_t<
29 tensor::TensorFullIndex,
30 Indices1>>, // TODO clarify Antisymmetric
31 std::conditional_t<
32 (ddc::type_seq_size_v<Indices2> == 0),
33 ddc::detail::TypeSeq<>,
34 ddc::detail::TypeSeq<misc::convert_type_seq_to_t<
35 tensor::TensorAntisymmetricIndex,
36 Indices2>>>>>;
37
38template <
39 misc::Specialization<ddc::detail::TypeSeq> Indices1,
40 misc::Specialization<ddc::detail::TypeSeq> Indices2,
41 misc::Specialization<tensor::Tensor> HodgeStarType,
42 class MetricDeterminantType,
43 misc::Specialization<tensor::Tensor> MetricProdType,
44 class ExecSpace>
45 requires(
46 misc::Specialization<MetricDeterminantType, ddc::ChunkSpan>
47 || misc::Specialization<MetricDeterminantType, tensor::Tensor>)
48HodgeStarType fill_hodge_star(
49 ExecSpace const& exec_space,
50 HodgeStarType hodge_star,
51 MetricDeterminantType metric_determinant,
52 MetricProdType metric_prod)
53{
54 sil::tensor::TensorAccessor<misc::convert_type_seq_to_t<
55 tensor::TensorLeviCivitaIndex,
56 ddc::type_seq_merge_t<tensor::primes<tensor::lower_t<Indices1>>, Indices2>>>
57 levi_civita_accessor;
58 ddc::DiscreteDomain<misc::convert_type_seq_to_t<
59 tensor::TensorLeviCivitaIndex,
60 ddc::type_seq_merge_t<tensor::primes<tensor::lower_t<Indices1>>, Indices2>>>
61 levi_civita_dom(levi_civita_accessor.domain());
62 ddc::Chunk levi_civita_alloc(
63 levi_civita_dom,
64 ddc::KokkosAllocator<
65 double,
66 typename ExecSpace::memory_space>()); // TODO consider avoid allocation
67 sil::tensor::Tensor levi_civita(levi_civita_alloc);
68
69 ddc::parallel_for_each(
70 exec_space,
71 hodge_star.non_indices_domain(),
72 KOKKOS_LAMBDA(
73 typename HodgeStarType::non_indices_domain_t::discrete_element_type elem) {
74 tensor_prod(hodge_star[elem], metric_prod[elem], levi_civita);
75 hodge_star[elem] *= Kokkos::sqrt(Kokkos::abs(metric_determinant(elem)))
76 / misc::factorial(ddc::type_seq_size_v<Indices1>);
77 });
78 return hodge_star;
79}
80
81template <
82 tensor::TensorIndex MetricIndex,
83 misc::Specialization<ddc::detail::TypeSeq> Indices1,
84 misc::Specialization<ddc::detail::TypeSeq> Indices2,
85 misc::Specialization<tensor::Tensor> HodgeStarType,
86 misc::Specialization<tensor::Tensor> MetricType,
87 class ExecSpace>
88HodgeStarType fill_hodge_star(
89 ExecSpace const& exec_space,
90 HodgeStarType hodge_star,
91 MetricType metric)
92{
93 static_assert(tensor::are_contravariant_v<
94 ddc::to_type_seq_t<typename MetricIndex::subindices_domain_t>>);
95 static_assert(tensor::are_contravariant_v<
96 ddc::to_type_seq_t<typename MetricType::accessor_t::natural_domain_t>>);
97 static_assert(tensor::are_contravariant_v<Indices1>);
98 static_assert(tensor::are_covariant_v<Indices2>);
99
100 // Allocate metric_det to receive metric field determinant values
101 ddc::Chunk metric_det_alloc(
102 metric.non_indices_domain(),
103 ddc::KokkosAllocator<double, typename ExecSpace::memory_space>());
104 ddc::ChunkSpan metric_det(metric_det_alloc);
105 // Allocate a buffer mirroring the metric as a full matrix, it will be overwritten by tensor::determinant() which involves a LU decomposition
106 ddc::Chunk buffer_alloc(
107 ddc::cartesian_prod_t<
108 typename MetricType::non_indices_domain_t,
109 ddc::DiscreteDomain<
110 tensor::metric_index_1<MetricIndex>,
111 tensor::metric_index_2<MetricIndex>>>(
112 metric.non_indices_domain(),
113 ddc::DiscreteDomain<
114 tensor::metric_index_1<MetricIndex>,
115 tensor::metric_index_2<MetricIndex>>(metric.natural_domain())),
116 ddc::KokkosAllocator<double, typename ExecSpace::memory_space>());
117 ddc::ChunkSpan buffer(buffer_alloc);
118 // Compute determinants
119 ddc::parallel_for_each(
120 exec_space,
121 metric.non_indices_domain(),
122 KOKKOS_LAMBDA(typename MetricType::non_indices_domain_t::discrete_element_type elem) {
123 ddc::annotated_for_each(
124 ddc::DiscreteDomain<
125 tensor::metric_index_1<MetricIndex>,
126 tensor::metric_index_2<MetricIndex>>(metric.natural_domain()),
127 [&](ddc::DiscreteElement<
128 tensor::metric_index_1<MetricIndex>,
129 tensor::metric_index_2<MetricIndex>> index) {
130 buffer(elem, index) = metric.get(metric.access_element(
131 elem,
132 index)); // TODO: triggers a "nvlink warning : Stack size for entry function cannot be statically determined"
133 });
134 metric_det(elem) = 1. / tensor::determinant(buffer[elem].allocation_kokkos_view());
135 });
136
137 // Allocate & compute the product of metrics
138 tensor::tensor_accessor_for_domain_t<
139 tensor::metric_prod_domain_t<MetricIndex, Indices1, tensor::primes<Indices1>>>
140 metric_prod_accessor;
141 ddc::Chunk metric_prod_alloc(
142 ddc::cartesian_prod_t<
143 typename MetricType::non_indices_domain_t,
144 tensor::metric_prod_domain_t<MetricIndex, Indices1, tensor::primes<Indices1>>>(
145 metric.non_indices_domain(),
146 metric_prod_accessor.domain()),
147 ddc::KokkosAllocator<double, typename ExecSpace::memory_space>());
148 tensor::Tensor metric_prod(metric_prod_alloc);
149
150 fill_metric_prod<
151 MetricIndex,
152 Indices1,
153 tensor::primes<Indices1>>(exec_space, metric_prod, metric);
154
155 // Compute Hodge star
156 return fill_hodge_star<Indices1, Indices2>(exec_space, hodge_star, metric_det, metric_prod);
157}
158
159} // namespace exterior
160
161} // namespace sil
sil::exterior::fill_hodge_star
HodgeStarType fill_hodge_star(ExecSpace const &exec_space, HodgeStarType hodge_star, MetricDeterminantType metric_determinant, MetricProdType metric_prod)
Definition hodge_star.hpp:47
sil::exterior::hodge_star_domain_t
ddc::detail::convert_type_seq_to_discrete_domain_t< ddc::type_seq_merge_t< ddc::detail::TypeSeq< misc::convert_type_seq_to_t< tensor::TensorFullIndex, Indices1 > >, std::conditional_t<(ddc::type_seq_size_v< Indices2 >==0), ddc::detail::TypeSeq<>, ddc::detail::TypeSeq< misc::convert_type_seq_to_t< tensor::TensorAntisymmetricIndex, Indices2 > > > > > hodge_star_domain_t
Definition hodge_star.hpp:26
sil::misc::factorial
constexpr std::size_t factorial(std::size_t k) noexcept
Definition factorial.hpp:13
sil::misc::convert_type_seq_to_t
typename detail::ConvertTypeSeqTo< T, Seq >::type convert_type_seq_to_t
Definition type_seq_conversion.hpp:45
sil::tensor::Tensor
Tensor(ddc::Chunk< ElementType, SupportType, Allocator >) -> Tensor< ElementType, SupportType, Kokkos::layout_right, typename Allocator::memory_space >
sil::tensor::metric_index_1
ddc::type_seq_element_t< 0, ddc::to_type_seq_t< typename MetricIndex::subindices_domain_t > > metric_index_1
Definition metric.hpp:34
sil::tensor::primes
detail::Primes< Indices, I >::type primes
Definition prime.hpp:50
sil::tensor::are_covariant_v
bool constexpr are_covariant_v
Definition character.hpp:230
sil::tensor::determinant
KOKKOS_FUNCTION ViewType::value_type determinant(const ViewType &matrix)
Definition determinant.hpp:18
sil::tensor::metric_index_2
ddc::type_seq_element_t< 1, ddc::to_type_seq_t< typename MetricIndex::subindices_domain_t > > metric_index_2
Definition metric.hpp:42
sil::tensor::metric_prod_domain_t
typename detail::MetricProdDomainType< MetricIndex, Indices1, Indices2 >::type metric_prod_domain_t
Definition metric.hpp:142
sil::tensor::tensor_accessor_for_domain_t
detail::TensorAccessorForDomain< Dom >::type tensor_accessor_for_domain_t
Definition tensor_impl.hpp:330
sil::tensor::are_contravariant_v
bool constexpr are_contravariant_v
Definition character.hpp:269
sil
The top-level namespace of SimiLie.
Definition csr.hpp:14