47 #ifndef __INTREPID2_PROJECTIONTOOLS_HPP__
48 #define __INTREPID2_PROJECTIONTOOLS_HPP__
50 #include "Intrepid2_ConfigDefs.hpp"
54 #include "Shards_CellTopology.hpp"
55 #include "Shards_BasicTopologies.hpp"
78 #include "Teuchos_LAPACK.hpp"
83 #ifdef HAVE_INTREPID2_KOKKOSKERNELS
84 #include "KokkosBatched_QR_Serial_Internal.hpp"
85 #include "KokkosBatched_ApplyQ_Serial_Internal.hpp"
86 #include "KokkosBatched_Trsv_Serial_Internal.hpp"
87 #include "KokkosBatched_Util.hpp"
144 template<
typename DeviceType>
147 using ExecSpaceType =
typename DeviceType::execution_space;
148 using MemSpaceType =
typename DeviceType::memory_space;
149 using EvalPointsType =
typename ProjectionStruct<DeviceType, double>::EvalPointsType;
171 template<
typename basisCoeffsValueType,
class ...basisCoeffsProperties,
172 typename funValsValueType,
class ...funValsProperties,
174 typename ortValueType,
class ...ortProperties>
176 getL2BasisCoeffs(Kokkos::DynRankView<basisCoeffsValueType,basisCoeffsProperties...> basisCoeffs,
177 const Kokkos::DynRankView<funValsValueType,funValsProperties...> targetAtEvalPoints,
178 const Kokkos::DynRankView<ortValueType, ortProperties...> cellOrientations,
179 const BasisType* cellBasis,
206 template<
typename basisCoeffsValueType,
class ...basisCoeffsProperties,
207 typename funValsValueType,
class ...funValsProperties,
209 typename ortValueType,
class ...ortProperties>
211 getL2DGBasisCoeffs(Kokkos::DynRankView<basisCoeffsValueType,basisCoeffsProperties...> basisCoeffs,
212 const Kokkos::DynRankView<funValsValueType,funValsProperties...> targetAtEvalPoints,
213 const Kokkos::DynRankView<ortValueType, ortProperties...> cellOrientations,
214 const BasisType* cellBasis,
240 template<
typename basisViewType,
typename targetViewType,
typename BasisType>
243 const targetViewType targetAtTargetEPoints,
244 const BasisType* cellBasis,
268 template<
class BasisCoeffsViewType,
class TargetValueViewType,
class TargetGradViewType,
269 class BasisType,
class OrientationViewType>
272 const TargetValueViewType targetAtEvalPoints,
273 const TargetGradViewType targetGradAtGradEvalPoints,
274 const OrientationViewType cellOrientations,
275 const BasisType* cellBasis,
302 template<
typename basisCoeffsValueType,
class ...basisCoeffsProperties,
303 typename funValsValueType,
class ...funValsProperties,
305 typename ortValueType,
class ...ortProperties>
307 getHCurlBasisCoeffs(Kokkos::DynRankView<basisCoeffsValueType,basisCoeffsProperties...> basisCoeffs,
308 const Kokkos::DynRankView<funValsValueType,funValsProperties...> targetAtEvalPoints,
309 const Kokkos::DynRankView<funValsValueType,funValsProperties...> targetCurlAtCurlEvalPoints,
310 const Kokkos::DynRankView<ortValueType, ortProperties...> cellOrientations,
311 const BasisType* cellBasis,
335 template<
typename basisCoeffsValueType,
class ...basisCoeffsProperties,
336 typename funValsValueType,
class ...funValsProperties,
338 typename ortValueType,
class ...ortProperties>
340 getHDivBasisCoeffs(Kokkos::DynRankView<basisCoeffsValueType,basisCoeffsProperties...> basisCoeffs,
341 const Kokkos::DynRankView<funValsValueType,funValsProperties...> targetAtEvalPoints,
342 const Kokkos::DynRankView<funValsValueType,funValsProperties...> targetDivAtDivEvalPoints,
343 const Kokkos::DynRankView<ortValueType, ortProperties...> cellOrientations,
344 const BasisType* cellBasis,
365 template<
typename basisCoeffsValueType,
class ...basisCoeffsProperties,
366 typename funValsValueType,
class ...funValsProperties,
368 typename ortValueType,
class ...ortProperties>
370 getHVolBasisCoeffs(Kokkos::DynRankView<basisCoeffsValueType,basisCoeffsProperties...> basisCoeffs,
371 const Kokkos::DynRankView<funValsValueType,funValsProperties...> targetAtEvalPoints,
372 [[maybe_unused]]
const Kokkos::DynRankView<ortValueType, ortProperties...> cellOrientations,
373 const BasisType* cellBasis,
393 template<
typename dstViewType,
394 typename dstBasisType,
395 typename srcViewType,
396 typename srcBasisType,
397 typename OrientationViewType>
400 const dstBasisType* dstBasis,
401 const srcViewType srcCoeffs,
402 const srcBasisType* srcBasis,
403 const OrientationViewType cellOrientations){
406 INTREPID2_TEST_FOR_EXCEPTION(dstBasis->getFunctionSpace() != srcBasis->getFunctionSpace(), std::runtime_error,
407 "The source and destination bases are not compatible. They need to belong to the same function space");
408 INTREPID2_TEST_FOR_EXCEPTION(dstBasis->getBaseCellTopology().getKey() != srcBasis->getBaseCellTopology().getKey(), std::runtime_error,
409 "The source and destination bases are not compatible. They do not have the same basic cell topology");
415 ordinal_type numCells = cellOrientations.extent(0);
416 ordinal_type dim = srcBasis->getBaseCellTopology().getDimension();
417 ordinal_type srcBasisCardinality = srcBasis->getCardinality();
418 ordinal_type fieldDimension = (srcBasis->getFunctionSpace() == Intrepid2::FUNCTION_SPACE_HCURL || srcBasis->getFunctionSpace() == Intrepid2::FUNCTION_SPACE_HDIV) ? dim : 1;
420 auto evaluationPoints = projStruct.getAllEvalPoints();
421 ordinal_type numPoints = evaluationPoints.extent(0);
423 using outViewType = Kokkos::DynRankView<typename srcBasisType::OutputValueType, DeviceType>;
424 outViewType srcAtEvalPoints, refBasisAtEvalPoints, basisAtEvalPoints;
425 if(fieldDimension == dim) {
426 srcAtEvalPoints = outViewType(
"srcAtEvalPoints", numCells, numPoints, dim);
427 refBasisAtEvalPoints = outViewType(
"refBasisAtEvalPoints", srcBasisCardinality, numPoints, dim);
428 basisAtEvalPoints = outViewType(
"basisAtEvalPoints", numCells, srcBasisCardinality, numPoints, dim);
430 srcAtEvalPoints = outViewType(
"srcAtEvalPoints", numCells, numPoints);
431 refBasisAtEvalPoints = outViewType(
"refBasisAtEvalPoints", srcBasisCardinality, numPoints);
432 basisAtEvalPoints = outViewType(
"basisAtEvalPoints", numCells, srcBasisCardinality, numPoints);
435 srcBasis->getValues(refBasisAtEvalPoints,evaluationPoints);
439 refBasisAtEvalPoints,
443 Kokkos::parallel_for(Kokkos::RangePolicy<typename DeviceType::execution_space>(0,numCells),
444 KOKKOS_LAMBDA (
const int &ic) {
445 for(
int j=0; j<numPoints; ++j) {
446 for(
int k=0; k<srcBasisCardinality; ++k) {
447 for(
int d=0; d<fieldDimension; ++d)
448 srcAtEvalPoints.access(ic,j,d) += srcCoeffs(ic,k)*basisAtEvalPoints.access(ic,k,j,d);
452 ExecSpaceType().fence();
475 std::string systemName_;
476 bool matrixIndependentOfCell_;
485 ElemSystem (std::string systemName,
bool matrixIndependentOfCell) :
486 systemName_(systemName), matrixIndependentOfCell_(matrixIndependentOfCell){};
515 template<
typename ViewType1,
typename ViewType2,
typename ViewType3,
typename ViewType4>
516 void solve(ViewType1 basisCoeffs, ViewType2 elemMat, ViewType2 elemRhs, ViewType2 tau,
517 ViewType3 w,
const ViewType4 elemDof, ordinal_type n, ordinal_type m=0) {
518 #ifdef HAVE_INTREPID2_KOKKOSKERNELS
519 solveDevice(basisCoeffs, elemMat, elemRhs, tau,
522 solveHost(basisCoeffs, elemMat, elemRhs, tau,
530 #ifdef HAVE_INTREPID2_KOKKOSKERNELS
531 template<
typename ViewType1,
typename ViewType2,
typename ViewType3,
typename ViewType4>
532 void solveDevice(ViewType1 basisCoeffs, ViewType2 elemMat, ViewType2 elemRhs, ViewType2 taul,
533 ViewType3 work,
const ViewType4 elemDof, ordinal_type n, ordinal_type m) {
534 using HostDeviceType = Kokkos::Device<Kokkos::DefaultHostExecutionSpace,Kokkos::HostSpace>;
536 ordinal_type numCells = basisCoeffs.extent(0);
538 if(matrixIndependentOfCell_) {
539 auto A0 = Kokkos::subview(elemMat, 0, Kokkos::ALL(), Kokkos::ALL());
540 auto tau0 = Kokkos::subview(taul, 0, Kokkos::ALL());
542 Kokkos::DynRankView<typename ViewType2::value_type, HostDeviceType> A0_host(
"A0_host", A0.extent(0),A0.extent(1));
543 auto A0_device = Kokkos::create_mirror_view(
typename DeviceType::memory_space(), A0_host);
544 Kokkos::deep_copy(A0_device, A0);
545 Kokkos::deep_copy(A0_host, A0_device);
547 for(ordinal_type i=n; i<n+m; ++i)
548 for(ordinal_type j=0; j<n; ++j)
549 A0_host(i,j) = A0_host(j,i);
551 Kokkos::DynRankView<typename ViewType2::value_type, HostDeviceType> tau0_host(
"A0_host", tau0.extent(0));
552 auto tau0_device = Kokkos::create_mirror_view(
typename DeviceType::memory_space(), tau0_host);
553 auto w0_host = Kokkos::create_mirror_view(Kokkos::subview(work, 0, Kokkos::ALL()));
556 KokkosBatched::SerialQR_Internal::invoke(A0_host.extent(0), A0_host.extent(1),
557 A0_host.data(), A0_host.stride_0(), A0_host.stride_1(),
558 tau0_host.data(), tau0_host.stride_0(), w0_host.data());
560 Kokkos::deep_copy(A0_device, A0_host);
561 Kokkos::deep_copy(A0, A0_device);
562 Kokkos::deep_copy(tau0_device, tau0_host);
563 Kokkos::deep_copy(tau0, tau0_device);
565 Kokkos::parallel_for (systemName_,
566 Kokkos::RangePolicy<ExecSpaceType, int> (0, numCells),
567 KOKKOS_LAMBDA (
const size_t ic) {
568 auto w = Kokkos::subview(work, ic, Kokkos::ALL());
570 auto b = Kokkos::subview(elemRhs, ic, Kokkos::ALL());
573 KokkosBatched::SerialApplyQ_RightForwardInternal::invoke(
574 1, A0.extent(0), A0.extent(1),
575 A0.data(), A0.stride_0(), A0.stride_1(),
576 tau0.data(), tau0.stride_0(),
577 b.data(), 1, b.stride_0(),
581 KokkosBatched::SerialTrsvInternalUpper<KokkosBatched::Algo::Trsv::Unblocked>::invoke(
false,
584 A0.data(), A0.stride_0(), A0.stride_1(),
585 b.data(), b.stride_0());
588 for(ordinal_type i=0; i<n; ++i){
589 basisCoeffs(ic,elemDof(i)) = b(i);
595 Kokkos::parallel_for (systemName_,
596 Kokkos::RangePolicy<ExecSpaceType, int> (0, numCells),
597 KOKKOS_LAMBDA (
const size_t ic) {
599 auto A = Kokkos::subview(elemMat, ic, Kokkos::ALL(), Kokkos::ALL());
600 auto tau = Kokkos::subview(taul, ic, Kokkos::ALL());
601 auto w = Kokkos::subview(work, ic, Kokkos::ALL());
603 for(ordinal_type i=n; i<n+m; ++i)
604 for(ordinal_type j=0; j<n; ++j)
608 KokkosBatched::SerialQR_Internal::invoke(A.extent(0), A.extent(1),
609 A.data(), A.stride_0(), A.stride_1(), tau.data(), tau.stride_0(), w.data());
611 auto b = Kokkos::subview(elemRhs, ic, Kokkos::ALL());
614 KokkosBatched::SerialApplyQ_RightForwardInternal::invoke(
615 1, A.extent(0), A.extent(1),
616 A.data(), A.stride_0(), A.stride_1(),
617 tau.data(), tau.stride_0(),
618 b.data(), 1, b.stride_0(),
622 KokkosBatched::SerialTrsvInternalUpper<KokkosBatched::Algo::Trsv::Unblocked>::invoke(
false,
625 A.data(), A.stride_0(), A.stride_1(),
626 b.data(), b.stride_0());
629 for(ordinal_type i=0; i<n; ++i){
630 basisCoeffs(ic,elemDof(i)) = b(i);
640 template<
typename ViewType1,
typename ViewType2,
typename ViewType3,
typename ViewType4>
641 void solveHost(ViewType1 basisCoeffs, ViewType2 elemMat, ViewType2 elemRhs, ViewType2 ,
642 ViewType3,
const ViewType4 elemDof, ordinal_type n, ordinal_type m) {
643 using value_type =
typename ViewType2::value_type;
644 using device_type = DeviceType;
645 using host_exec_space = Kokkos::DefaultHostExecutionSpace;
646 using host_memory_space = Kokkos::HostSpace;
647 using host_device_type = Kokkos::Device<host_exec_space,host_memory_space>;
648 using vector_host_type = Kokkos::View<value_type*,host_device_type>;
649 using scratch_host_type = Kokkos::View<value_type*,host_exec_space::scratch_memory_space>;
650 using matrix_host_type = Kokkos::View<value_type**,Kokkos::LayoutLeft,host_device_type>;
651 using do_not_init_tag = Kokkos::ViewAllocateWithoutInitializing;
652 using host_team_policy_type = Kokkos::TeamPolicy<host_exec_space>;
653 using host_range_policy_type = Kokkos::RangePolicy<host_exec_space>;
659 const ordinal_type numCells = basisCoeffs.extent(0);
660 const ordinal_type numRows = m+n, numCols = n;
663 Teuchos::LAPACK<ordinal_type,value_type> lapack;
666 Kokkos::View<ordinal_type*,host_device_type> elemDof_host(do_not_init_tag(
"elemDof_host"), elemDof.extent(0));
668 auto elemDof_device = Kokkos::create_mirror_view(
typename device_type::memory_space(), elemDof_host);
669 Kokkos::deep_copy(elemDof_device, elemDof); Kokkos::fence();
670 Kokkos::deep_copy(elemDof_host, elemDof_device);
674 auto elemRhs_host = Kokkos::create_mirror_view_and_copy(host_memory_space(), elemRhs);
675 auto elemMat_host = Kokkos::create_mirror_view_and_copy(host_memory_space(), elemMat);
678 auto basisCoeffs_host = Kokkos::create_mirror_view_and_copy(host_memory_space(), basisCoeffs);
680 if (matrixIndependentOfCell_) {
682 matrix_host_type A(do_not_init_tag(
"A"), numRows, numRows);
684 for (ordinal_type j=0;j<numRows;++j)
685 for (ordinal_type i=0;i<numRows;++i)
686 A(i, j) = elemMat_host(0, i, j);
688 for (ordinal_type j=0;j<numCols;++j)
689 for (ordinal_type i=numCols;i<numRows;++i)
693 ordinal_type lwork(-1);
695 ordinal_type info(0);
698 numRows, numRows, numCells,
699 nullptr, std::max(1,numRows),
700 nullptr, std::max(1,numRows),
706 matrix_host_type C(do_not_init_tag(
"C"), numRows, numCells);
708 host_range_policy_type policy(0, numCells);
711 (
"ProjectionTools::solveHost::matrixIndependentOfCell::pack",
712 policy, [=](
const ordinal_type & ic) {
713 for (ordinal_type i=0;i<numRows;++i)
714 C(i,ic) = elemRhs_host(ic, i);
719 vector_host_type work(do_not_init_tag(
"work"), lwork);
720 ordinal_type info(0);
722 numRows, numRows, numCells,
723 A.data(), std::max(1,numRows),
724 C.data(), std::max(1,numRows),
727 INTREPID2_TEST_FOR_EXCEPTION
728 (info != 0, std::runtime_error,
"GELS return non-zero info code");
732 (
"ProjectionTools::solveHost::matrixIndependentOfCell::unpack",
733 policy, [=](
const ordinal_type & ic) {
734 for (ordinal_type i=0;i<numCols;++i)
735 basisCoeffs_host(ic,elemDof_host(i)) = C(i,ic);
739 const ordinal_type level(0);
740 host_team_policy_type policy(numCells, 1, 1);
743 ordinal_type lwork(-1);
745 ordinal_type info(0);
749 nullptr, std::max(1,numRows),
750 nullptr, std::max(1,numRows),
756 const ordinal_type per_team_extent = numRows*numRows + numRows + lwork;
757 const ordinal_type per_team_scratch = scratch_host_type::shmem_size(per_team_extent);
758 policy.set_scratch_size(level, Kokkos::PerTeam(per_team_scratch));
762 (
"ProjectionTools::solveHost::matrixDependentOfCell",
763 policy, [=](
const typename host_team_policy_type::member_type& member) {
764 const ordinal_type ic = member.league_rank();
766 scratch_host_type scratch(member.team_scratch(level), per_team_extent);
767 value_type * sptr = scratch.data();
770 matrix_host_type A(sptr, numRows, numRows); sptr += A.span();
771 for (ordinal_type j=0;j<numRows;++j)
772 for (ordinal_type i=0;i<numRows;++i)
773 A(i, j) = elemMat_host(ic, i, j);
775 for (ordinal_type j=0;j<numCols;++j)
776 for (ordinal_type i=numCols;i<numRows;++i)
779 vector_host_type c(sptr, numRows); sptr += c.span();
780 for (ordinal_type i=0;i<numRows;++i)
781 c(i) = elemRhs_host(ic, i);
783 vector_host_type work(sptr, lwork); sptr += work.span();
784 ordinal_type info(0);
787 A.data(), std::max(1,numRows),
788 c.data(), std::max(1,numRows),
791 INTREPID2_TEST_FOR_EXCEPTION
792 (info != 0, std::runtime_error,
"GELS return non-zero info code");
795 for (ordinal_type i=0;i<numCols;++i)
796 basisCoeffs_host(ic,elemDof_host(i)) = c(i);
799 Kokkos::deep_copy(basisCoeffs, basisCoeffs_host);
Header file for the Intrepid2::Basis_HDIV_TET_I1_FEM class.
Header file for the Intrepid2::Basis_HDIV_HEX_I1_FEM class.
Header file for the Intrepid2::Basis_HDIV_WEDGE_I1_FEM class.
Header function for Intrepid2::Util class and other utility functions.
Header file for the Intrepid2::Basis_HDIV_QUAD_I1_FEM class.
Header file for the Intrepid2::Basis_HDIV_TRI_I1_FEM class.
Contains definitions of custom data types in Intrepid2.
Header file for the Intrepid2::ProjectionStruct.
void createL2ProjectionStruct(const BasisPtrType cellBasis, const ordinal_type targetCubDegree)
Initialize the ProjectionStruct for L2 projections.
Header file for the Intrepid2::Basis_HCURL_QUAD_I1_FEM class.
Header file for the Intrepid2::Basis_HCURL_WEDGE_I1_FEM class.
Header file for the Intrepid2::Basis_HCURL_TET_I1_FEM class.
Header file for the Intrepid2::Basis_HCURL_HEX_I1_FEM class.
Header file for the Intrepid2::Basis_HCURL_TRI_I1_FEM class.
An helper class to compute the evaluation points and weights needed for performing projections...
Stateless class that acts as a factory for a family of nodal bases (hypercube topologies only at this...
Header file for the abstract base class Intrepid2::Basis.