Skip to content

Commit

Permalink
Intrepid2: adding some HDIV team-level getValues functions and tests
Browse files Browse the repository at this point in the history
  • Loading branch information
@mperego
mperego committed Sep 26, 2024
1 parent ed10272 commit 5e164b6
Show file tree
Hide file tree
Showing 23 changed files with 1,812 additions and 243 deletions.
4 changes: 2 additions & 2 deletions packages/intrepid2/src/Discretization/Basis/Intrepid2_HCURL_TET_In_FEMDef.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -571,8 +571,8 @@ Basis_HCURL_TET_In_FEM<DT,OT,PT>::getScratchSpaceSize(
const PointViewType inputPoints,
const EOperator operatorType) const {
perTeamSpaceSize = 0;
ordinal_type scalarViewExtent = (operatorType == OPERATOR_VALUE) ? this->basisCardinality_ : 7*this->basisCardinality_;
perThreadSpaceSize = scalarViewExtent*get_dimension_scalar(inputPoints)*sizeof(typename BasisBase::scalarType);
ordinal_type scalarWorkViewExtent = (operatorType == OPERATOR_VALUE) ? this->basisCardinality_ : 7*this->basisCardinality_;
perThreadSpaceSize = scalarWorkViewExtent*get_dimension_scalar(inputPoints)*sizeof(typename BasisBase::scalarType);
}

template<typename DT, typename OT, typename PT>
Expand Down
4 changes: 2 additions & 2 deletions packages/intrepid2/src/Discretization/Basis/Intrepid2_HCURL_TRI_In_FEMDef.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -460,8 +460,8 @@ namespace Intrepid2 {
const PointViewType inputPoints,
const EOperator operatorType) const {
perTeamSpaceSize = 0;
ordinal_type scalarViewExtent = (operatorType == OPERATOR_VALUE) ? this->basisCardinality_ : 5*this->basisCardinality_;
perThreadSpaceSize = scalarViewExtent*get_dimension_scalar(inputPoints)*sizeof(typename BasisBase::scalarType);
ordinal_type scalarWorkViewExtent = (operatorType == OPERATOR_VALUE) ? this->basisCardinality_ : 5*this->basisCardinality_;
perThreadSpaceSize = scalarWorkViewExtent*get_dimension_scalar(inputPoints)*sizeof(typename BasisBase::scalarType);
}

template<typename DT, typename OT, typename PT>
Expand Down
34 changes: 25 additions & 9 deletions packages/intrepid2/src/Discretization/Basis/Intrepid2_HDIV_HEX_In_FEM.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -138,20 +138,21 @@ namespace Intrepid2 {
class Basis_HDIV_HEX_In_FEM
: public Basis<DeviceType,outputValueType,pointValueType> {
public:
using OrdinalTypeArray1DHost = typename Basis<DeviceType,outputValueType,pointValueType>::OrdinalTypeArray1DHost;
using OrdinalTypeArray2DHost = typename Basis<DeviceType,outputValueType,pointValueType>::OrdinalTypeArray2DHost;
using OrdinalTypeArray3DHost = typename Basis<DeviceType,outputValueType,pointValueType>::OrdinalTypeArray3DHost;
using BasisBase = Basis<DeviceType, outputValueType, pointValueType>;
using OrdinalTypeArray1DHost = typename BasisBase::OrdinalTypeArray1DHost;
using OrdinalTypeArray2DHost = typename BasisBase::OrdinalTypeArray2DHost;
using OrdinalTypeArray3DHost = typename BasisBase::OrdinalTypeArray3DHost;

/** \brief Constructor.
*/
Basis_HDIV_HEX_In_FEM(const ordinal_type order,
const EPointType pointType = POINTTYPE_EQUISPACED);

using OutputViewType = typename Basis<DeviceType,outputValueType,pointValueType>::OutputViewType;
using PointViewType = typename Basis<DeviceType,outputValueType,pointValueType>::PointViewType;
using ScalarViewType = typename Basis<DeviceType,outputValueType,pointValueType>::ScalarViewType;
using OutputViewType = typename BasisBase::OutputViewType;
using PointViewType = typename BasisBase::PointViewType;
using ScalarViewType = typename BasisBase::ScalarViewType;

using Basis<DeviceType,outputValueType,pointValueType>::getValues;
using BasisBase::getValues;

virtual
void
Expand All @@ -174,6 +175,23 @@ namespace Intrepid2 {
operatorType );
}

virtual void
getScratchSpaceSize( ordinal_type& perTeamSpaceSize,
ordinal_type& perThreadSpaceSize,
const PointViewType inputPointsconst,
const EOperator operatorType = OPERATOR_VALUE) const override;

KOKKOS_INLINE_FUNCTION
virtual void
getValues(
OutputViewType outputValues,
const PointViewType inputPoints,
const EOperator operatorType,
const typename Kokkos::TeamPolicy<typename DeviceType::execution_space>::member_type& team_member,
const typename DeviceType::execution_space::scratch_memory_space & scratchStorage,
const ordinal_type subcellDim = -1,
const ordinal_type subcellOrdinal = -1) const override;

virtual
void
getDofCoords( ScalarViewType dofCoords ) const override {
Expand Down Expand Up @@ -254,8 +272,6 @@ namespace Intrepid2 {

}// namespace Intrepid2



#include "Intrepid2_HDIV_HEX_In_FEMDef.hpp"

#endif
102 changes: 80 additions & 22 deletions packages/intrepid2/src/Discretization/Basis/Intrepid2_HDIV_HEX_In_FEMDef.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,19 +21,19 @@ namespace Intrepid2 {
// -------------------------------------------------------------------------------------
namespace Impl {

template<EOperator opType>
template<EOperator OpType>
template<typename OutputViewType,
typename inputViewType,
typename workViewType,
typename vinvViewType>
typename InputViewType,
typename WorkViewType,
typename VinvViewType>
KOKKOS_INLINE_FUNCTION
void
Basis_HDIV_HEX_In_FEM::Serial<opType>::
Basis_HDIV_HEX_In_FEM::Serial<OpType>::
getValues( OutputViewType output,
const inputViewType input,
workViewType work,
const vinvViewType vinvLine,
const vinvViewType vinvBubble) {
const InputViewType input,
WorkViewType work,
const VinvViewType vinvLine,
const VinvViewType vinvBubble) {
const ordinal_type cardLine = vinvLine.extent(0);
const ordinal_type cardBubble = vinvBubble.extent(0);

Expand All @@ -44,21 +44,21 @@ namespace Intrepid2 {
const auto input_y = Kokkos::subview(input, Kokkos::ALL(), range_type(1,2));
const auto input_z = Kokkos::subview(input, Kokkos::ALL(), range_type(2,3));

const ordinal_type dim_s = get_dimension_scalar(work);
const ordinal_type dim_s = get_dimension_scalar(input);
auto ptr0 = work.data();
auto ptr1 = work.data()+cardLine*npts*dim_s;
auto ptr2 = work.data()+2*cardLine*npts*dim_s;
auto ptr3 = work.data()+(2*cardLine+cardBubble)*npts*dim_s;

typedef typename Kokkos::DynRankView<typename workViewType::value_type, typename workViewType::memory_space> viewType;
auto vcprop = Kokkos::common_view_alloc_prop(work);
typedef typename Kokkos::DynRankView<typename InputViewType::value_type, typename WorkViewType::memory_space> ViewType;
auto vcprop = Kokkos::common_view_alloc_prop(input);

switch (opType) {
switch (OpType) {
case OPERATOR_VALUE: {
viewType workLine(Kokkos::view_wrap(ptr0, vcprop), cardLine, npts);
viewType outputLine(Kokkos::view_wrap(ptr1, vcprop), cardLine, npts);
viewType outputBubble_A(Kokkos::view_wrap(ptr2, vcprop), cardBubble, npts);
viewType outputBubble_B(Kokkos::view_wrap(ptr3, vcprop), cardBubble, npts);
ViewType workLine(Kokkos::view_wrap(ptr0, vcprop), cardLine, npts);
ViewType outputLine(Kokkos::view_wrap(ptr1, vcprop), cardLine, npts);
ViewType outputBubble_A(Kokkos::view_wrap(ptr2, vcprop), cardBubble, npts);
ViewType outputBubble_B(Kokkos::view_wrap(ptr3, vcprop), cardBubble, npts);

// tensor product
ordinal_type idx = 0;
Expand Down Expand Up @@ -138,13 +138,13 @@ namespace Intrepid2 {
break;
}
case OPERATOR_DIV: {
viewType workLine(Kokkos::view_wrap(ptr0, vcprop), cardLine, npts);
ViewType workLine(Kokkos::view_wrap(ptr0, vcprop), cardLine, npts);
// A line value
viewType outputBubble_A(Kokkos::view_wrap(ptr2, vcprop), cardBubble, npts);
ViewType outputBubble_A(Kokkos::view_wrap(ptr2, vcprop), cardBubble, npts);
// B line value
viewType outputBubble_B(Kokkos::view_wrap(ptr3, vcprop), cardBubble, npts);
ViewType outputBubble_B(Kokkos::view_wrap(ptr3, vcprop), cardBubble, npts);
// Line grad
viewType outputLine(Kokkos::view_wrap(ptr1, vcprop), cardLine, npts, 1);
ViewType outputLine(Kokkos::view_wrap(ptr1, vcprop), cardLine, npts, 1);

// tensor product
ordinal_type idx = 0;
Expand Down Expand Up @@ -508,6 +508,64 @@ namespace Intrepid2 {
this->dofCoeffs_ = Kokkos::create_mirror_view(typename DT::memory_space(), dofCoeffsHost);
Kokkos::deep_copy(this->dofCoeffs_, dofCoeffsHost);
}
}

template<typename DT, typename OT, typename PT>
void
Basis_HDIV_HEX_In_FEM<DT,OT,PT>::getScratchSpaceSize(
ordinal_type& perTeamSpaceSize,
ordinal_type& perThreadSpaceSize,
const PointViewType inputPoints,
const EOperator operatorType) const {
perTeamSpaceSize = 0;
perThreadSpaceSize = (2*this->vinvLine_.extent(0)+2*this->vinvBubble_.extent(0))*get_dimension_scalar(inputPoints)*sizeof(typename BasisBase::scalarType);
}

template<typename DT, typename OT, typename PT>
KOKKOS_INLINE_FUNCTION
void
Basis_HDIV_HEX_In_FEM<DT,OT,PT>::getValues(
OutputViewType outputValues,
const PointViewType inputPoints,
const EOperator operatorType,
const typename Kokkos::TeamPolicy<typename DT::execution_space>::member_type& team_member,
const typename DT::execution_space::scratch_memory_space & scratchStorage,
const ordinal_type subcellDim,
const ordinal_type subcellOrdinal) const {

INTREPID2_TEST_FOR_ABORT( !((subcellDim == -1) && (subcellOrdinal == -1)),
">>> ERROR: (Intrepid2::Basis_HDIV_HEX_In_FEM::getValues), The capability of selecting subsets of basis functions has not been implemented yet.");

const int numPoints = inputPoints.extent(0);
using ScalarType = typename ScalarTraits<typename PointViewType::value_type>::scalar_type;
using WorkViewType = Kokkos::DynRankView< ScalarType,typename DT::execution_space::scratch_memory_space,Kokkos::MemoryTraits<Kokkos::Unmanaged> >;
ordinal_type sizePerPoint = (2*this->vinvLine_.extent(0)+2*this->vinvBubble_.extent(0))*get_dimension_scalar(inputPoints);
WorkViewType workView(scratchStorage, sizePerPoint*team_member.team_size());
using range_type = Kokkos::pair<ordinal_type,ordinal_type>;

switch(operatorType) {
case OPERATOR_VALUE:
Kokkos::parallel_for (Kokkos::TeamThreadRange (team_member, numPoints), [=] (ordinal_type& pt) {
auto output = Kokkos::subview( outputValues, Kokkos::ALL(), range_type (pt,pt+1), Kokkos::ALL() );
const auto input = Kokkos::subview( inputPoints, range_type(pt, pt+1), Kokkos::ALL() );
WorkViewType work(workView.data() + sizePerPoint*team_member.team_rank(), sizePerPoint);
Impl::Basis_HDIV_HEX_In_FEM::Serial<OPERATOR_VALUE>::getValues( output, input, work, this->vinvLine_, this->vinvBubble_ );
});
break;
case OPERATOR_DIV:
Kokkos::parallel_for (Kokkos::TeamThreadRange (team_member, numPoints), [=] (ordinal_type& pt) {
auto output = Kokkos::subview( outputValues, Kokkos::ALL(), range_type(pt,pt+1), Kokkos::ALL() );
const auto input = Kokkos::subview( inputPoints, range_type(pt,pt+1), Kokkos::ALL() );
WorkViewType work(workView.data() + sizePerPoint*team_member.team_rank(), sizePerPoint);
Impl::Basis_HDIV_HEX_In_FEM::Serial<OPERATOR_DIV>::getValues( output, input, work, this->vinvLine_, this->vinvBubble_ );
});
break;
default: {
INTREPID2_TEST_FOR_ABORT( true,
">>> ERROR (Basis_HDIV_HEX_In_FEM): getValues not implemented for this operator");
}
}
}

} // namespace Intrepid2

#endif
33 changes: 26 additions & 7 deletions packages/intrepid2/src/Discretization/Basis/Intrepid2_HDIV_QUAD_In_FEM.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -135,20 +135,21 @@ namespace Intrepid2 {
class Basis_HDIV_QUAD_In_FEM
: public Basis<DeviceType,outputValueType,pointValueType> {
public:
using OrdinalTypeArray1DHost = typename Basis<DeviceType,outputValueType,pointValueType>::OrdinalTypeArray1DHost;
using OrdinalTypeArray2DHost = typename Basis<DeviceType,outputValueType,pointValueType>::OrdinalTypeArray2DHost;
using OrdinalTypeArray3DHost = typename Basis<DeviceType,outputValueType,pointValueType>::OrdinalTypeArray3DHost;
using BasisBase = Basis<DeviceType, outputValueType, pointValueType>;
using OrdinalTypeArray1DHost = typename BasisBase::OrdinalTypeArray1DHost;
using OrdinalTypeArray2DHost = typename BasisBase::OrdinalTypeArray2DHost;
using OrdinalTypeArray3DHost = typename BasisBase::OrdinalTypeArray3DHost;

/** \brief Constructor.
*/
Basis_HDIV_QUAD_In_FEM(const ordinal_type order,
const EPointType pointType = POINTTYPE_EQUISPACED);

using OutputViewType = typename Basis<DeviceType,outputValueType,pointValueType>::OutputViewType;
using PointViewType = typename Basis<DeviceType,outputValueType,pointValueType>::PointViewType;
using ScalarViewType = typename Basis<DeviceType,outputValueType,pointValueType>::ScalarViewType;
using OutputViewType = typename BasisBase::OutputViewType;
using PointViewType = typename BasisBase::PointViewType;
using ScalarViewType = typename BasisBase::ScalarViewType;

using Basis<DeviceType,outputValueType,pointValueType>::getValues;
using BasisBase::getValues;

virtual
void
Expand All @@ -170,6 +171,24 @@ namespace Intrepid2 {
this->vinvBubble_,
operatorType );
}

virtual void
getScratchSpaceSize( ordinal_type& perTeamSpaceSize,
ordinal_type& perThreadSpaceSize,
const PointViewType inputPointsconst,
const EOperator operatorType = OPERATOR_VALUE) const override;

KOKKOS_INLINE_FUNCTION
virtual void
getValues(
OutputViewType outputValues,
const PointViewType inputPoints,
const EOperator operatorType,
const typename Kokkos::TeamPolicy<typename DeviceType::execution_space>::member_type& team_member,
const typename DeviceType::execution_space::scratch_memory_space & scratchStorage,
const ordinal_type subcellDim = -1,
const ordinal_type subcellOrdinal = -1) const override;

virtual
void
getDofCoords( ScalarViewType dofCoords ) const override {
Expand Down
Loading

0 comments on commit 5e164b6

Please sign in to comment.