8#include <similie/misc/unsecure_parallel_deepcopy.hpp>
10#include <KokkosBatched_InverseLU_Decl.hpp>
11#include <KokkosBatched_LU_Decl.hpp>
13#include "character.hpp"
14#include "diagonal_tensor.hpp"
15#include "identity_tensor.hpp"
16#include "lorentzian_sign_tensor.hpp"
18#include "relabelization.hpp"
19#include "symmetric_tensor.hpp"
20#include "tensor_prod.hpp"
28template <
class... CDim>
33template <
class MetricIndex>
35 = ddc::type_seq_element_t<0, ddc::to_type_seq_t<typename MetricIndex::subindices_domain_t>>;
37template <
class... CDim>
42template <
class MetricIndex>
44 = ddc::type_seq_element_t<1, ddc::to_type_seq_t<typename MetricIndex::subindices_domain_t>>;
48template <
class TypeSeqCDim>
49struct ConvertTypeSeqToMetricIndex1;
51template <
class... CDim>
52struct ConvertTypeSeqToMetricIndex1<ddc::detail::TypeSeq<CDim...>>
57template <
class TypeSeqCDim>
58struct ConvertTypeSeqToMetricIndex2;
60template <
class... CDim>
61struct ConvertTypeSeqToMetricIndex2<ddc::detail::TypeSeq<CDim...>>
76 typename detail::ConvertTypeSeqToMetricIndex1<
77 typename Index1::type_seq_dimensions>::type,
78 typename detail::ConvertTypeSeqToMetricIndex2<
79 typename Index2::type_seq_dimensions>::type>,
82template <TensorNatIndex Index1, TensorNatIndex Index2,
class Dom>
88 typename detail::ConvertTypeSeqToMetricIndex1<
89 typename Index1::type_seq_dimensions>::type,
90 typename detail::ConvertTypeSeqToMetricIndex2<
91 typename Index2::type_seq_dimensions>::type>,
95template <misc::Specialization<Tensor> TensorType, TensorNatIndex Index1, TensorNatIndex Index2>
99 typename detail::ConvertTypeSeqToMetricIndex1<
100 typename Index1::type_seq_dimensions>::type,
101 typename detail::ConvertTypeSeqToMetricIndex2<
102 typename Index2::type_seq_dimensions>::type>,
105template <TensorNatIndex Index1, TensorNatIndex Index2, misc::Specialization<Tensor> TensorType>
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>,
120template <
class MetricIndex,
class Indices1,
class Indices2>
121struct MetricProdDomainType;
123template <
class MetricIndex,
class... Index1,
class... Index2>
124struct MetricProdDomainType<
126 ddc::detail::TypeSeq<Index1...>,
127 ddc::detail::TypeSeq<Index2...>>
129 static_assert(
sizeof...(Index1) ==
sizeof...(Index2));
131 ddc::DiscreteDomain<MetricIndex>,
133 ddc::type_seq_element_t<
134 ddc::type_seq_rank_v<Index1, ddc::detail::TypeSeq<Index1...>>,
135 ddc::detail::TypeSeq<Index2...>>>...>;
141 TensorIndex MetricIndex,
145 typename detail::MetricProdDomainType<MetricIndex, Indices1, Indices2>::type;
155 class LayoutStridedPolicy,
157struct MetricProdType;
164 class LayoutStridedPolicy,
166struct MetricProdType<
169 ddc::detail::TypeSeq<Index1...>,
170 ddc::detail::TypeSeq<Index2...>,
176 ddc::cartesian_prod_t<
180 ddc::detail::TypeSeq<Index1...>,
181 ddc::detail::TypeSeq<Index2...>>>,
190 TensorIndex MetricIndex,
193 class LayoutStridedPolicy,
205template <
class Indices1,
class Indices2>
206struct FillMetricProd;
209struct FillMetricProd<ddc::detail::TypeSeq<>, ddc::detail::TypeSeq<>>
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_)
218 misc::detail::unsecure_parallel_deepcopy(exec_space, metric_prod, metric_prod_);
223template <
class HeadIndex1,
class... TailIndex1,
class HeadIndex2,
class... TailIndex2>
224struct FillMetricProd<
225 ddc::detail::TypeSeq<HeadIndex1, TailIndex1...>,
226 ddc::detail::TypeSeq<HeadIndex2, TailIndex2...>>
228 template <
class MetricProdType,
class MetricType,
class MetricProdType_,
class ExecSpace>
229 static MetricProdType run(
230 ExecSpace
const& exec_space,
231 MetricProdType metric_prod,
233 MetricProdType_ metric_prod_)
235 ddc::cartesian_prod_t<
236 typename MetricProdType_::discrete_domain_type,
238 typename MetricType::indices_domain_t,
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>());
250 if (new_metric_prod_dom_.size() != 0) {
251 ddc::parallel_for_each(
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],
264 return FillMetricProd<
265 ddc::detail::TypeSeq<TailIndex1...>,
266 ddc::detail::TypeSeq<TailIndex2...>>::
267 run(exec_space, metric_prod, metric, new_metric_prod_);
273 TensorIndex MetricIndex,
279 typename MetricType::non_indices_domain_t,
283 typename MetricType::layout_type,
284 typename MetricType::memory_space>
286 ExecSpace
const& exec_space,
288 typename MetricType::non_indices_domain_t,
292 typename MetricType::layout_type,
293 typename MetricType::memory_space> metric_prod,
296 typename MetricType::non_indices_domain_t dom_(metric.domain());
297 ddc::Chunk metric_prod_alloc_(
299 ddc::KokkosAllocator<double, typename ExecSpace::memory_space>());
301 ddc::parallel_fill(exec_space, metric_prod_, 1.);
303 return detail::FillMetricProd<Indices1, Indices2>::
304 run(exec_space, metric_prod, metric, metric_prod_);
316 misc::Specialization<Tensor> MetricType,
317 misc::Specialization<Tensor> TensorType,
322 detail::non_primes<typename MetricType::accessor_t::natural_domain_t>>
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>>(
331 ddc::KokkosAllocator<double, typename ExecSpace::memory_space>());
335 detail::non_primes<typename MetricType::accessor_t::natural_domain_t>>
336 result(result_alloc);
337 ddc::parallel_for_each(
339 tensor.non_indices_domain(),
340 KOKKOS_LAMBDA(
typename TensorType::non_indices_domain_t::discrete_element_type elem) {
346 typename MetricType::accessor_t::natural_domain_t>>,
348 typename MetricType::accessor_t::natural_domain_t>>>>(
351 misc::detail::unsecure_parallel_deepcopy(exec_space, tensor, result);
355 detail::non_primes<typename MetricType::accessor_t::natural_domain_t>>(tensor);
359 TensorIndex MetricIndex,
360 TensorNatIndex... Index1,
367 ddc::detail::TypeSeq<Index1...>>
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,
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>());
389 ddc::detail::TypeSeq<Index1...>,
390 primes<ddc::detail::TypeSeq<Index1...>>>(exec_space, metric_prod, metric);
395template <misc::Specialization<Tensor> MetricType>
398 ddc::to_type_seq_t<typename MetricType::accessor_t::natural_domain_t>,
402template <TensorIndex MetricIndex, misc::Specialization<Tensor> MetricType,
class ExecSpace>
404 ExecSpace
const& exec_space,
412 ddc::parallel_for_each(
420 typename MetricType::accessor_t::
422 ddc::to_type_seq_t<
typename MetricType::accessor_t::
427 ddc::Chunk buffer_alloc(
428 ddc::cartesian_prod_t<
429 typename MetricType::non_indices_domain_t,
433 metric.non_indices_domain(),
437 ddc::KokkosAllocator<double, typename ExecSpace::memory_space>());
438 ddc::ChunkSpan buffer(buffer_alloc);
440 ddc::Chunk buffer_alloc2(
441 ddc::cartesian_prod_t<
442 typename MetricType::non_indices_domain_t,
444 metric.non_indices_domain(),
448 metric.natural_domain()
451 * metric.natural_domain()
454 ddc::KokkosAllocator<double, typename ExecSpace::memory_space>());
455 ddc::ChunkSpan buffer2(buffer_alloc2);
458 ddc::parallel_for_each(
460 inv_metric.non_indices_domain(),
462 MetricType>::non_indices_domain_t::discrete_element_type elem) {
463 ddc::annotated_for_each(
467 [&](ddc::DiscreteElement<
470 buffer(elem, index) = metric(metric.access_element(elem, index));
473 int err = KokkosBatched::SerialLU<KokkosBatched::Algo::SolveLU::Unblocked>::invoke(
475 .allocation_kokkos_view());
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");
497 ddc::annotated_for_each(
499 inv_metric.domain()),
502 inv_metric.mem(elem, mem_index) = buffer(
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(
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.