SimiLie
Loading...
Searching...
No Matches
metric.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/unsecure_parallel_deepcopy.hpp>
9
10#include <KokkosBatched_InverseLU_Decl.hpp>
11#include <KokkosBatched_LU_Decl.hpp>
12
13#include "character.hpp"
14#include "diagonal_tensor.hpp"
15#include "identity_tensor.hpp"
16#include "lorentzian_sign_tensor.hpp"
17#include "prime.hpp"
18#include "relabelization.hpp"
19#include "symmetric_tensor.hpp"
20#include "tensor_prod.hpp"
21
22namespace sil {
23
24namespace tensor {
25
26// Abstract indices used to characterize a generic metric
27
28template <class... CDim>
30{
31};
32
33template <class MetricIndex>
35 = ddc::type_seq_element_t<0, ddc::to_type_seq_t<typename MetricIndex::subindices_domain_t>>;
36
37template <class... CDim>
38struct MetricIndex2 : TensorNaturalIndex<CDim...>
39{
40};
41
42template <class MetricIndex>
44 = ddc::type_seq_element_t<1, ddc::to_type_seq_t<typename MetricIndex::subindices_domain_t>>;
45
46namespace detail {
47
48template <class TypeSeqCDim>
49struct ConvertTypeSeqToMetricIndex1;
50
51template <class... CDim>
52struct ConvertTypeSeqToMetricIndex1<ddc::detail::TypeSeq<CDim...>>
53{
54 using type = MetricIndex1<CDim...>;
55};
56
57template <class TypeSeqCDim>
58struct ConvertTypeSeqToMetricIndex2;
59
60template <class... CDim>
61struct ConvertTypeSeqToMetricIndex2<ddc::detail::TypeSeq<CDim...>>
62{
63 using type = MetricIndex2<CDim...>;
64};
65
66} // namespace detail
67
68// Relabelize metric
69template <
72 TensorNatIndex Index2>
74 Dom,
75 ddc::detail::TypeSeq<
76 typename detail::ConvertTypeSeqToMetricIndex1<
77 typename Index1::type_seq_dimensions>::type,
78 typename detail::ConvertTypeSeqToMetricIndex2<
79 typename Index2::type_seq_dimensions>::type>,
80 ddc::detail::TypeSeq<uncharacterize_t<Index1>, uncharacterize_t<Index2>>>;
82template <TensorNatIndex Index1, TensorNatIndex Index2, class Dom>
84 Dom metric_dom)
85{
87 ddc::detail::TypeSeq<
88 typename detail::ConvertTypeSeqToMetricIndex1<
89 typename Index1::type_seq_dimensions>::type,
90 typename detail::ConvertTypeSeqToMetricIndex2<
91 typename Index2::type_seq_dimensions>::type>,
92 ddc::detail::TypeSeq<uncharacterize_t<Index1>, uncharacterize_t<Index2>>>(metric_dom);
93}
95template <misc::Specialization<Tensor> TensorType, TensorNatIndex Index1, TensorNatIndex Index2>
97 TensorType,
98 ddc::detail::TypeSeq<
99 typename detail::ConvertTypeSeqToMetricIndex1<
100 typename Index1::type_seq_dimensions>::type,
101 typename detail::ConvertTypeSeqToMetricIndex2<
102 typename Index2::type_seq_dimensions>::type>,
103 ddc::detail::TypeSeq<uncharacterize_t<Index1>, uncharacterize_t<Index2>>>;
105template <TensorNatIndex Index1, TensorNatIndex Index2, misc::Specialization<Tensor> TensorType>
107{
109 ddc::detail::TypeSeq<
110 typename detail::ConvertTypeSeqToMetricIndex1<
111 typename Index1::type_seq_dimensions>::type,
112 typename detail::ConvertTypeSeqToMetricIndex2<
113 typename Index2::type_seq_dimensions>::type>,
114 ddc::detail::TypeSeq<uncharacterize_t<Index1>, uncharacterize_t<Index2>>>(tensor);
115}
116
117// Compute domain for a tensor product of metrics (ie. g_mu_muprime*g_nu_nuprime*...)
118namespace detail {
119
120template <class MetricIndex, class Indices1, class Indices2>
121struct MetricProdDomainType;
122
123template <class MetricIndex, class... Index1, class... Index2>
124struct MetricProdDomainType<
125 MetricIndex,
126 ddc::detail::TypeSeq<Index1...>,
127 ddc::detail::TypeSeq<Index2...>>
128{
129 static_assert(sizeof...(Index1) == sizeof...(Index2));
130 using type = ddc::cartesian_prod_t<relabelize_metric_in_domain_t<
131 ddc::DiscreteDomain<MetricIndex>,
132 Index1,
133 ddc::type_seq_element_t<
134 ddc::type_seq_rank_v<Index1, ddc::detail::TypeSeq<Index1...>>,
135 ddc::detail::TypeSeq<Index2...>>>...>;
136};
137
138} // namespace detail
139
140template <
141 TensorIndex MetricIndex,
145 typename detail::MetricProdDomainType<MetricIndex, Indices1, Indices2>::type;
146
147namespace detail {
148
149// Compute tensor product of metrics (ie. g_mu_muprime*g_nu_nuprime*...)
150template <
151 class NonMetricDom,
152 class MetricIndex,
153 class Indices1,
154 class Indices2,
155 class LayoutStridedPolicy,
156 class MemorySpace>
157struct MetricProdType;
158
159template <
160 class NonMetricDom,
161 class MetricIndex,
162 class... Index1,
163 class... Index2,
164 class LayoutStridedPolicy,
165 class MemorySpace>
166struct MetricProdType<
167 NonMetricDom,
168 MetricIndex,
169 ddc::detail::TypeSeq<Index1...>,
170 ddc::detail::TypeSeq<Index2...>,
171 LayoutStridedPolicy,
172 MemorySpace>
173{
174 using type = tensor::Tensor<
175 double,
176 ddc::cartesian_prod_t<
177 NonMetricDom,
179 MetricIndex,
180 ddc::detail::TypeSeq<Index1...>,
181 ddc::detail::TypeSeq<Index2...>>>,
182 LayoutStridedPolicy,
183 MemorySpace>;
184};
185
186} // namespace detail
187
188template <
190 TensorIndex MetricIndex,
193 class LayoutStridedPolicy,
194 class MemorySpace>
195using metric_prod_t = typename detail::MetricProdType<
196 NonMetricDom,
197 MetricIndex,
198 Indices1,
199 Indices2,
200 LayoutStridedPolicy,
201 MemorySpace>::type;
202
203namespace detail {
204
205template <class Indices1, class Indices2>
206struct FillMetricProd;
207
208template <>
209struct FillMetricProd<ddc::detail::TypeSeq<>, ddc::detail::TypeSeq<>>
210{
211 template <class MetricProdType, class MetricType, class MetricProdType_, class ExecSpace>
212 static MetricProdType run(
213 ExecSpace const& exec_space,
214 MetricProdType metric_prod,
215 [[maybe_unused]] MetricType metric,
216 MetricProdType_ metric_prod_)
217 {
218 misc::detail::unsecure_parallel_deepcopy(exec_space, metric_prod, metric_prod_);
219 return metric_prod;
220 }
221};
222
223template <class HeadIndex1, class... TailIndex1, class HeadIndex2, class... TailIndex2>
224struct FillMetricProd<
225 ddc::detail::TypeSeq<HeadIndex1, TailIndex1...>,
226 ddc::detail::TypeSeq<HeadIndex2, TailIndex2...>>
227{
228 template <class MetricProdType, class MetricType, class MetricProdType_, class ExecSpace>
229 static MetricProdType run(
230 ExecSpace const& exec_space,
231 MetricProdType metric_prod,
232 MetricType metric,
233 MetricProdType_ metric_prod_)
234 {
235 ddc::cartesian_prod_t<
236 typename MetricProdType_::discrete_domain_type,
238 typename MetricType::indices_domain_t,
239 HeadIndex1,
240 HeadIndex2>>
241 new_metric_prod_dom_(
242 metric_prod_.domain(),
244 metric.indices_domain()));
245 ddc::Chunk new_metric_prod_alloc_(
246 new_metric_prod_dom_,
247 ddc::KokkosAllocator<double, typename ExecSpace::memory_space>());
248 tensor::Tensor new_metric_prod_(new_metric_prod_alloc_);
249
250 if (new_metric_prod_dom_.size() != 0) {
251 ddc::parallel_for_each(
252 exec_space,
253 new_metric_prod_.non_indices_domain(),
254 KOKKOS_LAMBDA(typename decltype(new_metric_prod_)::non_indices_domain_t::
255 discrete_element_type elem) {
257 new_metric_prod_[elem],
258 metric_prod_[elem],
260 });
261 }
262
263
264 return FillMetricProd<
265 ddc::detail::TypeSeq<TailIndex1...>,
266 ddc::detail::TypeSeq<TailIndex2...>>::
267 run(exec_space, metric_prod, metric, new_metric_prod_);
268 }
269};
270} // namespace detail
271
272template <
273 TensorIndex MetricIndex,
277 class ExecSpace>
279 typename MetricType::non_indices_domain_t,
280 MetricIndex,
281 Indices1,
282 Indices2,
283 typename MetricType::layout_type,
284 typename MetricType::memory_space>
286 ExecSpace const& exec_space,
288 typename MetricType::non_indices_domain_t,
289 MetricIndex,
290 Indices1,
291 Indices2,
292 typename MetricType::layout_type,
293 typename MetricType::memory_space> metric_prod,
294 MetricType metric)
295{
296 typename MetricType::non_indices_domain_t dom_(metric.domain());
297 ddc::Chunk metric_prod_alloc_(
298 dom_,
299 ddc::KokkosAllocator<double, typename ExecSpace::memory_space>());
300 tensor::Tensor metric_prod_(metric_prod_alloc_);
301 ddc::parallel_fill(exec_space, metric_prod_, 1.);
302
303 return detail::FillMetricProd<Indices1, Indices2>::
304 run(exec_space, metric_prod, metric, metric_prod_);
305}
306
307namespace detail {
308
309template <class Dom>
310using non_primes = ddc::type_seq_remove_t<ddc::to_type_seq_t<Dom>, primes<ddc::to_type_seq_t<Dom>>>;
311
312} // namespace detail
313
314// Apply metrics inplace (like g_mu_muprime*T^muprime^nu)
315template <
316 misc::Specialization<Tensor> MetricType,
317 misc::Specialization<Tensor> TensorType,
318 class ExecSpace>
320 TensorType,
322 detail::non_primes<typename MetricType::accessor_t::natural_domain_t>>
323inplace_apply_metric(ExecSpace const& exec_space, TensorType tensor, MetricType metric_prod)
324{
325 ddc::Chunk result_alloc(
328 detail::non_primes<typename MetricType::accessor_t::natural_domain_t>>,
329 detail::non_primes<typename MetricType::accessor_t::natural_domain_t>>(
330 tensor.domain()),
331 ddc::KokkosAllocator<double, typename ExecSpace::memory_space>());
333 TensorType,
335 detail::non_primes<typename MetricType::accessor_t::natural_domain_t>>
336 result(result_alloc);
337 ddc::parallel_for_each(
338 exec_space,
339 tensor.non_indices_domain(),
340 KOKKOS_LAMBDA(typename TensorType::non_indices_domain_t::discrete_element_type elem) {
342 result[elem],
343 metric_prod[elem],
345 swap_character_t<detail::non_primes<
346 typename MetricType::accessor_t::natural_domain_t>>,
347 primes<swap_character_t<detail::non_primes<
348 typename MetricType::accessor_t::natural_domain_t>>>>(
349 tensor)[elem]);
350 });
351 misc::detail::unsecure_parallel_deepcopy(exec_space, tensor, result);
352
355 detail::non_primes<typename MetricType::accessor_t::natural_domain_t>>(tensor);
356}
357
358template <
359 TensorIndex MetricIndex,
360 TensorNatIndex... Index1,
363 class ExecSpace>
365 TensorType,
366 swap_character_t<ddc::detail::TypeSeq<Index1...>>,
367 ddc::detail::TypeSeq<Index1...>>
368inplace_apply_metric(ExecSpace const& exec_space, TensorType tensor, MetricType metric)
369{
371 MetricIndex,
372 ddc::detail::TypeSeq<Index1...>,
373 primes<ddc::detail::TypeSeq<Index1...>>>>
374 metric_prod_accessor;
375 ddc::Chunk metric_prod_alloc(
376 ddc::cartesian_prod_t<
377 typename TensorType::non_indices_domain_t,
379 MetricIndex,
380 ddc::detail::TypeSeq<Index1...>,
381 primes<ddc::detail::TypeSeq<Index1...>>>>(
382 tensor.non_indices_domain(),
383 metric_prod_accessor.domain()),
384 ddc::KokkosAllocator<double, typename ExecSpace::memory_space>());
385 tensor::Tensor metric_prod(metric_prod_alloc);
386
388 MetricIndex,
389 ddc::detail::TypeSeq<Index1...>,
390 primes<ddc::detail::TypeSeq<Index1...>>>(exec_space, metric_prod, metric);
391
392 return inplace_apply_metric(exec_space, tensor, metric_prod);
393}
395template <misc::Specialization<Tensor> MetricType>
397 MetricType,
398 ddc::to_type_seq_t<typename MetricType::accessor_t::natural_domain_t>,
400
401// Compute invert metric (g_mu_nu for gmunu or gmunu for g_mu_nu)
402template <TensorIndex MetricIndex, misc::Specialization<Tensor> MetricType, class ExecSpace>
404 ExecSpace const& exec_space,
406 MetricType metric)
407{
408 if constexpr (
412 ddc::parallel_for_each(
413 exec_space,
414 inv_metric.domain(),
416 inv_metric.mem(elem)
417 = 1.
418 / metric.mem(relabelize_indices_in<
419 swap_character_t<ddc::to_type_seq_t<
420 typename MetricType::accessor_t::
422 ddc::to_type_seq_t<typename MetricType::accessor_t::
423 natural_domain_t>>(elem));
424 });
426 // Allocate a buffer mirroring the metric as a full matrix
427 ddc::Chunk buffer_alloc(
428 ddc::cartesian_prod_t<
429 typename MetricType::non_indices_domain_t,
430 ddc::DiscreteDomain<
433 metric.non_indices_domain(),
434 ddc::DiscreteDomain<
436 tensor::metric_index_2<MetricIndex>>(metric.natural_domain())),
437 ddc::KokkosAllocator<double, typename ExecSpace::memory_space>());
438 ddc::ChunkSpan buffer(buffer_alloc);
439 // Allocate a buffer for KokkosBatched::SerialInverseLU internal needs
440 ddc::Chunk buffer_alloc2(
441 ddc::cartesian_prod_t<
442 typename MetricType::non_indices_domain_t,
443 ddc::DiscreteDomain<tensor::metric_index_1<MetricIndex>>>(
444 metric.non_indices_domain(),
445 ddc::DiscreteDomain<tensor::metric_index_1<MetricIndex>>(
446 ddc::DiscreteElement<tensor::metric_index_1<MetricIndex>>(0),
447 ddc::DiscreteVector<tensor::metric_index_1<MetricIndex>>(
448 metric.natural_domain()
449 .template extent<
451 * metric.natural_domain()
452 .template extent<
454 ddc::KokkosAllocator<double, typename ExecSpace::memory_space>());
455 ddc::ChunkSpan buffer2(buffer_alloc2);
456
457 // process
458 ddc::parallel_for_each(
459 exec_space,
460 inv_metric.non_indices_domain(),
461 KOKKOS_LAMBDA(typename invert_metric_t<
462 MetricType>::non_indices_domain_t::discrete_element_type elem) {
463 ddc::annotated_for_each(
464 ddc::DiscreteDomain<
466 tensor::metric_index_2<MetricIndex>>(metric.natural_domain()),
467 [&](ddc::DiscreteElement<
470 buffer(elem, index) = metric(metric.access_element(elem, index));
471 });
472
473 int err = KokkosBatched::SerialLU<KokkosBatched::Algo::SolveLU::Unblocked>::invoke(
474 buffer[elem]
475 .allocation_kokkos_view()); // Seems to require diagonal-dominant
476 err += KokkosBatched::SerialInverseLU<KokkosBatched::Algo::SolveLU::Unblocked>::
477 invoke(buffer[elem].allocation_kokkos_view(),
478 buffer2[elem].allocation_kokkos_view());
479 assert(err == 0 && "Kokkos-kernels failed at inverting metric tensor");
480 /*
481 Kokkos::View<double**, Kokkos::LayoutRight, Kokkos::HostSpace> const
482 sol("metric_inversion_sol", n, n);
483 for (std::size_t i=0; i<n; ++i) {
484 Kokkos::View<double*, Kokkos::LayoutRight, Kokkos::HostSpace> const
485 rhs("metric_inversion_rhs", n);
486 for (std::size_t j=0; j<n; ++j) {
487 rhs(j) = 0.;
488 }
489 rhs(i) = 1.;
490 Kokkos::View<double**, Kokkos::LayoutRight, Kokkos::HostSpace> const
491 tmp("metric_inversion_tmp", n, n+4);
492 const int err = KokkosBatched::SerialGesv<KokkosBatched::Gesv::NoPivoting>::invoke(metric_view, Kokkos::subview(sol, i, Kokkos::ALL), rhs, tmp);
493 assert(err==0 && "Kokkos-kernels failed at inverting metric tensor");
494 }
495 */
496
497 ddc::annotated_for_each(
498 tensor::swap_character_t<ddc::DiscreteDomain<MetricIndex>>(
499 inv_metric.domain()),
500 [&](tensor::swap_character_t<ddc::DiscreteElement<MetricIndex>>
501 mem_index) {
502 inv_metric.mem(elem, mem_index) = buffer(
503 elem,
504 tensor::relabelize_indices_in<
505 tensor::swap_character_t<ddc::detail::TypeSeq<
506 tensor::metric_index_1<MetricIndex>,
507 tensor::metric_index_2<MetricIndex>>>,
508 ddc::detail::TypeSeq<
509 tensor::metric_index_1<MetricIndex>,
510 tensor::metric_index_2<MetricIndex>>>(
511 inv_metric.accessor().canonical_natural_element(
512 mem_index)));
513 });
514 });
515 }
516
517 return inv_metric;
518}
519
520} // namespace tensor
521
522} // namespace sil
Tensor(ddc::Chunk< ElementType, SupportType, Allocator >) -> Tensor< ElementType, SupportType, Kokkos::layout_right, typename Allocator::memory_space >
relabelize_indices_of_t< MetricType, ddc::to_type_seq_t< typename MetricType::accessor_t::natural_domain_t >, swap_character_t< ddc::to_type_seq_t< typename MetricType::accessor_t::natural_domain_t > > > invert_metric_t
Definition metric.hpp:394
typename detail::MetricProdType< NonMetricDom, MetricIndex, Indices1, Indices2, LayoutStridedPolicy, MemorySpace >::type metric_prod_t
Definition metric.hpp:193
detail::RelabelizeIndicesOfType< TensorType, OldIndices, NewIndices >::type relabelize_indices_of_t
ddc::type_seq_element_t< 0, ddc::to_type_seq_t< typename MetricIndex::subindices_domain_t > > metric_index_1
Definition metric.hpp:34
detail::Primes< Indices, I >::type primes
Definition prime.hpp:50
constexpr relabelize_indices_in_t< T, OldIndices, NewIndices > relabelize_indices_in(T t)
invert_metric_t< MetricType > fill_inverse_metric(ExecSpace const &exec_space, invert_metric_t< MetricType > inv_metric, MetricType metric)
Definition metric.hpp:401
constexpr relabelize_metric_in_domain_t< Dom, Index1, Index2 > relabelize_metric_in_domain(Dom metric_dom)
Definition metric.hpp:81
relabelize_indices_of_t< TensorType, ddc::detail::TypeSeq< typename detail::ConvertTypeSeqToMetricIndex1< typename Index1::type_seq_dimensions >::type, typename detail::ConvertTypeSeqToMetricIndex2< typename Index2::type_seq_dimensions >::type >, ddc::detail::TypeSeq< uncharacterize_t< Index1 >, uncharacterize_t< Index2 > > > relabelize_metric_t
Definition metric.hpp:94
constexpr relabelize_metric_t< TensorType, Index1, Index2 > relabelize_metric(TensorType tensor)
Definition metric.hpp:104
ddc::type_seq_element_t< 1, ddc::to_type_seq_t< typename MetricIndex::subindices_domain_t > > metric_index_2
Definition metric.hpp:42
typename detail::MetricProdDomainType< MetricIndex, Indices1, Indices2 >::type metric_prod_domain_t
Definition metric.hpp:142
metric_prod_t< typename MetricType::non_indices_domain_t, MetricIndex, Indices1, Indices2, typename MetricType::layout_type, typename MetricType::memory_space > fill_metric_prod(ExecSpace const &exec_space, metric_prod_t< typename MetricType::non_indices_domain_t, MetricIndex, Indices1, Indices2, typename MetricType::layout_type, typename MetricType::memory_space > metric_prod, MetricType metric)
Definition metric.hpp:283
Tensor< ElementType, ddc::DiscreteDomain< ProdDDim... >, LayoutStridedPolicy, MemorySpace > tensor_prod(Tensor< ElementType, ddc::DiscreteDomain< ProdDDim... >, LayoutStridedPolicy, MemorySpace > prod_tensor, Tensor< ElementType, ddc::DiscreteDomain< DDim1... >, LayoutStridedPolicy, MemorySpace > tensor1, Tensor< ElementType, ddc::DiscreteDomain< DDim2... >, LayoutStridedPolicy, MemorySpace > tensor2)
detail::SwapCharacter< T >::type swap_character_t
relabelize_indices_of_t< TensorType, swap_character_t< detail::non_primes< typename MetricType::accessor_t::natural_domain_t > >, detail::non_primes< typename MetricType::accessor_t::natural_domain_t > > inplace_apply_metric(ExecSpace const &exec_space, TensorType tensor, MetricType metric_prod)
Definition metric.hpp:321
relabelize_indices_in_t< Dom, ddc::detail::TypeSeq< typename detail::ConvertTypeSeqToMetricIndex1< typename Index1::type_seq_dimensions >::type, typename detail::ConvertTypeSeqToMetricIndex2< typename Index2::type_seq_dimensions >::type >, ddc::detail::TypeSeq< uncharacterize_t< Index1 >, uncharacterize_t< Index2 > > > relabelize_metric_in_domain_t
Definition metric.hpp:71
detail::Uncharacterize< Index >::type uncharacterize_t
constexpr relabelize_indices_of_t< Tensor, OldIndices, NewIndices > relabelize_indices_of(Tensor tensor)
typename detail::NaturalDomainType< Index >::type natural_domain_t
detail::TensorAccessorForDomain< Dom >::type tensor_accessor_for_domain_t
typename detail::RelabelizeIndicesInType< T, OldIndices, NewIndices >::type relabelize_indices_in_t
The top-level namespace of SimiLie.
Definition csr.hpp:14