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