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(
255 typename decltype(new_metric_prod_)::non_indices_domain_t::
256 discrete_element_type elem) {
258 new_metric_prod_[elem],
259 metric_prod_[elem],
261 });
262 }
263
264
265 return FillMetricProd<
266 ddc::detail::TypeSeq<TailIndex1...>,
267 ddc::detail::TypeSeq<TailIndex2...>>::
268 run(exec_space, metric_prod, metric, new_metric_prod_);
269 }
270};
271} // namespace detail
272
273template <
274 TensorIndex MetricIndex,
278 class ExecSpace>
280 typename MetricType::non_indices_domain_t,
281 MetricIndex,
282 Indices1,
283 Indices2,
284 typename MetricType::layout_type,
285 typename MetricType::memory_space>
287 ExecSpace const& exec_space,
289 typename MetricType::non_indices_domain_t,
290 MetricIndex,
291 Indices1,
292 Indices2,
293 typename MetricType::layout_type,
294 typename MetricType::memory_space> metric_prod,
295 MetricType metric)
296{
297 typename MetricType::non_indices_domain_t dom_(metric.domain());
298 ddc::Chunk metric_prod_alloc_(
299 dom_,
300 ddc::KokkosAllocator<double, typename ExecSpace::memory_space>());
301 tensor::Tensor metric_prod_(metric_prod_alloc_);
302 ddc::parallel_fill(exec_space, metric_prod_, 1.);
303
304 return detail::FillMetricProd<Indices1, Indices2>::
305 run(exec_space, metric_prod, metric, metric_prod_);
306}
307
308namespace detail {
309
310template <class Dom>
311using non_primes = ddc::type_seq_remove_t<ddc::to_type_seq_t<Dom>, primes<ddc::to_type_seq_t<Dom>>>;
312
313} // namespace detail
314
315// Apply metrics inplace (like g_mu_muprime*T^muprime^nu)
316template <
317 misc::Specialization<Tensor> MetricType,
318 misc::Specialization<Tensor> TensorType,
319 class ExecSpace>
321 TensorType,
323 detail::non_primes<typename MetricType::accessor_t::natural_domain_t>>
324inplace_apply_metric(ExecSpace const& exec_space, TensorType tensor, MetricType metric_prod)
325{
326 ddc::Chunk result_alloc(
329 detail::non_primes<typename MetricType::accessor_t::natural_domain_t>>,
330 detail::non_primes<typename MetricType::accessor_t::natural_domain_t>>(
331 tensor.domain()),
332 ddc::KokkosAllocator<double, typename ExecSpace::memory_space>());
334 TensorType,
336 detail::non_primes<typename MetricType::accessor_t::natural_domain_t>>
337 result(result_alloc);
338 ddc::parallel_for_each(
339 exec_space,
340 tensor.non_indices_domain(),
341 KOKKOS_LAMBDA(typename TensorType::non_indices_domain_t::discrete_element_type elem) {
343 result[elem],
344 metric_prod[elem],
346 swap_character_t<detail::non_primes<
347 typename MetricType::accessor_t::natural_domain_t>>,
348 primes<swap_character_t<detail::non_primes<
349 typename MetricType::accessor_t::natural_domain_t>>>>(
350 tensor)[elem]);
351 });
352 misc::detail::unsecure_parallel_deepcopy(exec_space, tensor, result);
353
356 detail::non_primes<typename MetricType::accessor_t::natural_domain_t>>(tensor);
357}
358
359template <
360 TensorIndex MetricIndex,
361 TensorNatIndex... Index1,
364 class ExecSpace>
366 TensorType,
367 swap_character_t<ddc::detail::TypeSeq<Index1...>>,
368 ddc::detail::TypeSeq<Index1...>>
369inplace_apply_metric(ExecSpace const& exec_space, TensorType tensor, MetricType metric)
370{
372 MetricIndex,
373 ddc::detail::TypeSeq<Index1...>,
374 primes<ddc::detail::TypeSeq<Index1...>>>>
375 metric_prod_accessor;
376 ddc::Chunk metric_prod_alloc(
377 ddc::cartesian_prod_t<
378 typename TensorType::non_indices_domain_t,
380 MetricIndex,
381 ddc::detail::TypeSeq<Index1...>,
382 primes<ddc::detail::TypeSeq<Index1...>>>>(
383 tensor.non_indices_domain(),
384 metric_prod_accessor.domain()),
385 ddc::KokkosAllocator<double, typename ExecSpace::memory_space>());
386 tensor::Tensor metric_prod(metric_prod_alloc);
387
389 MetricIndex,
390 ddc::detail::TypeSeq<Index1...>,
391 primes<ddc::detail::TypeSeq<Index1...>>>(exec_space, metric_prod, metric);
392
393 return inplace_apply_metric(exec_space, tensor, metric_prod);
394}
396template <misc::Specialization<Tensor> MetricType>
398 MetricType,
399 ddc::to_type_seq_t<typename MetricType::accessor_t::natural_domain_t>,
401
402// Compute invert metric (g_mu_nu for gmunu or gmunu for g_mu_nu)
403template <TensorIndex MetricIndex, misc::Specialization<Tensor> MetricType, class ExecSpace>
405 ExecSpace const& exec_space,
407 MetricType metric)
408{
409 if constexpr (
413 ddc::parallel_for_each(
414 exec_space,
415 inv_metric.domain(),
417 inv_metric.mem(elem)
418 = 1.
419 / metric.mem(
421 swap_character_t<ddc::to_type_seq_t<
422 typename MetricType::accessor_t::
424 ddc::to_type_seq_t<typename MetricType::accessor_t::
425 natural_domain_t>>(elem));
426 });
428 // Allocate a buffer mirroring the metric as a full matrix
429 ddc::Chunk buffer_alloc(
430 ddc::cartesian_prod_t<
431 typename MetricType::non_indices_domain_t,
432 ddc::DiscreteDomain<
435 metric.non_indices_domain(),
436 ddc::DiscreteDomain<
438 tensor::metric_index_2<MetricIndex>>(metric.natural_domain())),
439 ddc::KokkosAllocator<double, typename ExecSpace::memory_space>());
440 ddc::ChunkSpan buffer(buffer_alloc);
441 // Allocate a buffer for KokkosBatched::SerialInverseLU internal needs
442 ddc::Chunk buffer_alloc2(
443 ddc::cartesian_prod_t<
444 typename MetricType::non_indices_domain_t,
445 ddc::DiscreteDomain<tensor::metric_index_1<MetricIndex>>>(
446 metric.non_indices_domain(),
447 ddc::DiscreteDomain<tensor::metric_index_1<MetricIndex>>(
448 ddc::DiscreteElement<tensor::metric_index_1<MetricIndex>>(0),
449 ddc::DiscreteVector<tensor::metric_index_1<MetricIndex>>(
450 metric.natural_domain()
451 .template extent<
453 * metric.natural_domain()
454 .template extent<
456 ddc::KokkosAllocator<double, typename ExecSpace::memory_space>());
457 ddc::ChunkSpan buffer2(buffer_alloc2);
458
459 // process
460 ddc::parallel_for_each(
461 exec_space,
462 inv_metric.non_indices_domain(),
463 KOKKOS_LAMBDA(
464 typename invert_metric_t<
465 MetricType>::non_indices_domain_t::discrete_element_type elem) {
466 ddc::device_for_each(
467 ddc::DiscreteDomain<
469 tensor::metric_index_2<MetricIndex>>(metric.natural_domain()),
470 [&](ddc::DiscreteElement<
473 buffer(elem, index) = metric(metric.access_element(elem, index));
474 });
475
476 int err = KokkosBatched::SerialLU<KokkosBatched::Algo::SolveLU::Unblocked>::invoke(
477 buffer[elem]
478 .allocation_kokkos_view()); // Seems to require diagonal-dominant
479 err += KokkosBatched::SerialInverseLU<KokkosBatched::Algo::SolveLU::Unblocked>::
480 invoke(buffer[elem].allocation_kokkos_view(),
481 buffer2[elem].allocation_kokkos_view());
482 assert(err == 0 && "Kokkos-kernels failed at inverting metric tensor");
483 /*
484 Kokkos::View<double**, Kokkos::LayoutRight, Kokkos::HostSpace> const
485 sol("metric_inversion_sol", n, n);
486 for (std::size_t i=0; i<n; ++i) {
487 Kokkos::View<double*, Kokkos::LayoutRight, Kokkos::HostSpace> const
488 rhs("metric_inversion_rhs", n);
489 for (std::size_t j=0; j<n; ++j) {
490 rhs(j) = 0.;
491 }
492 rhs(i) = 1.;
493 Kokkos::View<double**, Kokkos::LayoutRight, Kokkos::HostSpace> const
494 tmp("metric_inversion_tmp", n, n+4);
495 const int err = KokkosBatched::SerialGesv<KokkosBatched::Gesv::NoPivoting>::invoke(metric_view, Kokkos::subview(sol, i, Kokkos::ALL), rhs, tmp);
496 assert(err==0 && "Kokkos-kernels failed at inverting metric tensor");
497 }
498 */
499
500 ddc::device_for_each(
501 tensor::swap_character_t<ddc::DiscreteDomain<MetricIndex>>(
502 inv_metric.domain()),
503 [&](tensor::swap_character_t<ddc::DiscreteElement<MetricIndex>>
504 mem_index) {
505 inv_metric.mem(elem, mem_index) = buffer(
506 elem,
507 tensor::relabelize_indices_in<
508 tensor::swap_character_t<ddc::detail::TypeSeq<
509 tensor::metric_index_1<MetricIndex>,
510 tensor::metric_index_2<MetricIndex>>>,
511 ddc::detail::TypeSeq<
512 tensor::metric_index_1<MetricIndex>,
513 tensor::metric_index_2<MetricIndex>>>(
514 inv_metric.accessor().canonical_natural_element(
515 mem_index)));
516 });
517 });
518 }
519
520 return inv_metric;
521}
522
523} // namespace tensor
524
525} // 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:395
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:402
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:284
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:322
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