8#include <similie/misc/macros.hpp>
9#include <similie/misc/unsecure_parallel_deepcopy.hpp>
11#include <KokkosBatched_InverseLU_Decl.hpp>
12#include <KokkosBatched_LU_Decl.hpp>
14#include "character.hpp"
15#include "diagonal_tensor.hpp"
16#include "identity_tensor.hpp"
17#include "lorentzian_sign_tensor.hpp"
19#include "relabelization.hpp"
20#include "symmetric_tensor.hpp"
21#include "tensor_prod.hpp"
29template <
class... CDim>
34template <
class MetricIndex>
36 = ddc::type_seq_element_t<0, ddc::to_type_seq_t<typename MetricIndex::subindices_domain_t>>;
38template <
class... CDim>
43template <
class MetricIndex>
45 = ddc::type_seq_element_t<1, ddc::to_type_seq_t<typename MetricIndex::subindices_domain_t>>;
49template <
class TypeSeqCDim>
50struct ConvertTypeSeqToMetricIndex1;
52template <
class... CDim>
53struct ConvertTypeSeqToMetricIndex1<ddc::detail::TypeSeq<CDim...>>
58template <
class TypeSeqCDim>
59struct ConvertTypeSeqToMetricIndex2;
61template <
class... CDim>
62struct ConvertTypeSeqToMetricIndex2<ddc::detail::TypeSeq<CDim...>>
77 typename detail::ConvertTypeSeqToMetricIndex1<
78 typename Index1::type_seq_dimensions>::type,
79 typename detail::ConvertTypeSeqToMetricIndex2<
80 typename Index2::type_seq_dimensions>::type>,
83template <TensorNatIndex Index1, TensorNatIndex Index2,
class Dom>
89 typename detail::ConvertTypeSeqToMetricIndex1<
90 typename Index1::type_seq_dimensions>::type,
91 typename detail::ConvertTypeSeqToMetricIndex2<
92 typename Index2::type_seq_dimensions>::type>,
96template <misc::Specialization<Tensor> TensorType, TensorNatIndex Index1, TensorNatIndex Index2>
100 typename detail::ConvertTypeSeqToMetricIndex1<
101 typename Index1::type_seq_dimensions>::type,
102 typename detail::ConvertTypeSeqToMetricIndex2<
103 typename Index2::type_seq_dimensions>::type>,
106template <TensorNatIndex Index1, TensorNatIndex Index2, misc::Specialization<Tensor> TensorType>
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>,
121template <
class MetricIndex,
class Indices1,
class Indices2>
122struct MetricProdDomainType;
124template <
class MetricIndex,
class... Index1,
class... Index2>
125struct MetricProdDomainType<
127 ddc::detail::TypeSeq<Index1...>,
128 ddc::detail::TypeSeq<Index2...>>
130 static_assert(
sizeof...(Index1) ==
sizeof...(Index2));
132 ddc::DiscreteDomain<MetricIndex>,
134 ddc::type_seq_element_t<
135 ddc::type_seq_rank_v<Index1, ddc::detail::TypeSeq<Index1...>>,
136 ddc::detail::TypeSeq<Index2...>>>...>;
142 TensorIndex MetricIndex,
146 typename detail::MetricProdDomainType<MetricIndex, Indices1, Indices2>::type;
156 class LayoutStridedPolicy,
158struct MetricProdType;
165 class LayoutStridedPolicy,
167struct MetricProdType<
170 ddc::detail::TypeSeq<Index1...>,
171 ddc::detail::TypeSeq<Index2...>,
177 ddc::cartesian_prod_t<
181 ddc::detail::TypeSeq<Index1...>,
182 ddc::detail::TypeSeq<Index2...>>>,
191 TensorIndex MetricIndex,
194 class LayoutStridedPolicy,
206template <
class Indices1,
class Indices2>
207struct FillMetricProd;
210struct FillMetricProd<ddc::detail::TypeSeq<>, ddc::detail::TypeSeq<>>
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_)
219 misc::detail::unsecure_parallel_deepcopy(exec_space, metric_prod, metric_prod_);
224template <
class HeadIndex1,
class... TailIndex1,
class HeadIndex2,
class... TailIndex2>
225struct FillMetricProd<
226 ddc::detail::TypeSeq<HeadIndex1, TailIndex1...>,
227 ddc::detail::TypeSeq<HeadIndex2, TailIndex2...>>
229 template <
class MetricProdType,
class MetricType,
class MetricProdType_,
class ExecSpace>
230 static MetricProdType run(
231 ExecSpace
const& exec_space,
232 MetricProdType metric_prod,
234 MetricProdType_ metric_prod_)
236 ddc::cartesian_prod_t<
237 typename MetricProdType_::discrete_domain_type,
239 typename MetricType::indices_domain_t,
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>());
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",
256 new_metric_prod_.non_indices_domain(),
258 typename decltype(new_metric_prod_)::non_indices_domain_t::
259 discrete_element_type elem) {
261 new_metric_prod_[elem],
268 return FillMetricProd<
269 ddc::detail::TypeSeq<TailIndex1...>,
270 ddc::detail::TypeSeq<TailIndex2...>>::
271 run(exec_space, metric_prod, metric, new_metric_prod_);
277 TensorIndex MetricIndex,
283 typename MetricType::non_indices_domain_t,
287 typename MetricType::layout_type,
288 typename MetricType::memory_space>
290 ExecSpace
const& exec_space,
292 typename MetricType::non_indices_domain_t,
296 typename MetricType::layout_type,
297 typename MetricType::memory_space> metric_prod,
300 typename MetricType::non_indices_domain_t dom_(metric.domain());
301 ddc::Chunk metric_prod_alloc_(
303 ddc::KokkosAllocator<double, typename ExecSpace::memory_space>());
305 ddc::parallel_fill(exec_space, metric_prod_, 1.);
307 return detail::FillMetricProd<Indices1, Indices2>::
308 run(exec_space, metric_prod, metric, metric_prod_);
320 misc::Specialization<Tensor> MetricType,
321 misc::Specialization<Tensor> TensorType,
326 detail::non_primes<typename MetricType::accessor_t::natural_domain_t>>
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>>(
335 ddc::KokkosAllocator<double, typename ExecSpace::memory_space>());
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",
345 tensor.non_indices_domain(),
346 KOKKOS_LAMBDA(
typename TensorType::non_indices_domain_t::discrete_element_type elem) {
352 typename MetricType::accessor_t::natural_domain_t>>,
354 typename MetricType::accessor_t::natural_domain_t>>>>(
357 misc::detail::unsecure_parallel_deepcopy(exec_space, tensor, result);
361 detail::non_primes<typename MetricType::accessor_t::natural_domain_t>>(tensor);
365 TensorIndex MetricIndex,
366 TensorNatIndex... Index1,
373 ddc::detail::TypeSeq<Index1...>>
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,
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>());
395 ddc::detail::TypeSeq<Index1...>,
396 primes<ddc::detail::TypeSeq<Index1...>>>(exec_space, metric_prod, metric);
401template <misc::Specialization<Tensor> MetricType>
404 ddc::to_type_seq_t<typename MetricType::accessor_t::natural_domain_t>,
408template <TensorIndex MetricIndex, misc::Specialization<Tensor> MetricType,
class ExecSpace>
410 ExecSpace
const& exec_space,
418 SIMILIE_DEBUG_LOG(
"similie_invert_diagonal_metric");
419 ddc::parallel_for_each(
420 "similie_invert_diagonal_metric",
429 typename MetricType::accessor_t::
431 ddc::to_type_seq_t<
typename MetricType::accessor_t::
436 ddc::Chunk buffer_alloc(
437 ddc::cartesian_prod_t<
438 typename MetricType::non_indices_domain_t,
442 metric.non_indices_domain(),
446 ddc::KokkosAllocator<double, typename ExecSpace::memory_space>());
447 ddc::ChunkSpan buffer(buffer_alloc);
449 ddc::Chunk buffer_alloc2(
450 ddc::cartesian_prod_t<
451 typename MetricType::non_indices_domain_t,
453 metric.non_indices_domain(),
457 metric.natural_domain()
460 * metric.natural_domain()
463 ddc::KokkosAllocator<double, typename ExecSpace::memory_space>());
464 ddc::ChunkSpan buffer2(buffer_alloc2);
467 SIMILIE_DEBUG_LOG(
"similie_invert_metric");
468 ddc::parallel_for_each(
469 "similie_invert_metric",
471 inv_metric.non_indices_domain(),
474 MetricType>::non_indices_domain_t::discrete_element_type elem) {
475 ddc::device_for_each(
479 [&](ddc::DiscreteElement<
482 buffer(elem, index) = metric(metric.access_element(elem, index));
485 int err = KokkosBatched::SerialLU<KokkosBatched::Algo::SolveLU::Unblocked>::invoke(
487 .allocation_kokkos_view());
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");
509 ddc::device_for_each(
511 inv_metric.domain()),
514 inv_metric.mem(elem, mem_index) = buffer(
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(
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
typename detail::MetricProdType< NonMetricDom, MetricIndex, Indices1, Indices2, LayoutStridedPolicy, MemorySpace >::type metric_prod_t
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
detail::Primes< Indices, I >::type primes
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)
constexpr relabelize_metric_in_domain_t< Dom, Index1, Index2 > relabelize_metric_in_domain(Dom metric_dom)
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
constexpr relabelize_metric_t< TensorType, Index1, Index2 > relabelize_metric(TensorType tensor)
ddc::type_seq_element_t< 1, ddc::to_type_seq_t< typename MetricIndex::subindices_domain_t > > metric_index_2
typename detail::MetricProdDomainType< MetricIndex, Indices1, Indices2 >::type metric_prod_domain_t
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)
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)
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
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.