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