Skip to content

Commit

Permalink
Intrepid2: Implementation of team-level Basis::getValues
Browse files Browse the repository at this point in the history 
- Implemented team-level getValues for classic Lagrangian basis functions.
- Modified/added tests to compare the team-level getValues with host getValues
- Modified impelementation of JacobiPolynomial to reduce FAD temporaries

Signed-off-by: Mauro Perego <[email protected]>
  • Loading branch information
@mperego
mperego committed Oct 17, 2024
1 parent ffcdd92 commit f714191
Show file tree
Hide file tree
Showing 224 changed files with 16,969 additions and 1,320 deletions.
55 changes: 55 additions & 0 deletions packages/intrepid2/src/Discretization/Basis/Intrepid2_Basis.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -379,6 +379,61 @@ using HostBasisPtr = BasisPtr<typename Kokkos::HostSpace::device_type, OutputTyp
}
}


/** \brief Return the size of the scratch space, in bytes, needed for using the team-level implementation of getValues.
Warning, <var>inputPoints</var> is only used to deduce the type of the points where to evaluate basis functions.
The rank of </var>inputPoints</var> and its size are not relevant, however,
when using DFAD types, </var>inputPoints</var> cannot be empty,
otherwise the size of the scracth space needed won't be deduced correctly.
\param space [in] - inputPoints
\param perTeamSpaceSize [out] - size of the scratch space needed per team
\param perThreadeSize [out] - size of the scratch space beeded per thread
*/
virtual
void getScratchSpaceSize( ordinal_type& perTeamSpaceSize,
ordinal_type& perThreadSpaceSize,
const PointViewType inputPoints,
const EOperator operatorType = OPERATOR_VALUE) const {
INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE( true, std::logic_error,
">>> ERROR (Basis::getValuesScratchSpace): this method is not supported or should be overridden accordingly by derived classes.");
}


/** \brief Team-level evaluation of basis functions on a <strong>reference cell</strong>.
Returns values of <var>operatorType</var> acting on basis functions for a set of
points in the <strong>reference cell</strong> for which the basis is defined.
The interface allow also to select basis functions associated to a particular entity.
As an example, if <var>subcellDim==1</var> (edges) and <var>subcellOrdinal==0</var>, <var>outputValues</var> will contain all the basis functions associated with the first edge.
<var>outputValues</var> will contain all the cell basis functions when the default value (-1) is used for <var>subcellDim</var> and <var>subcellOrdinal</var>
\param outputValues [out] - variable rank array with the basis values
\param inputPoints [in] - rank-2 array (P,D) with the evaluation points
\param operatorType [in] - the operator acting on the basis functions
\param teamMember [in] - team member of the Kokkos::TemaPolicy
\param scratchStorage [in] - scratch space to use by each team
\param subcellDim [in] - the dimension of the subcells, the default values of -1 returns basis functions associated to subcells of all dimensions
\param subcellOrdinal [in] - the ordinal of the subcell, the default values of -1 returns basis functions associated to subcells of all ordinals
\remark This function is supposed to be called within a TeamPolicy kernel.
The size of the required scratch space is determined by the getScratchSpaceSize function.
*/
KOKKOS_INLINE_FUNCTION
virtual
void getValues( OutputViewType /* outputValues */,
const PointViewType /* inputPoints */,
const EOperator /* operatorType */,
const typename Kokkos::TeamPolicy<ExecutionSpace>::member_type& teamMember,
const typename ExecutionSpace::scratch_memory_space &scratchStorage,
const ordinal_type subcellDim=-1,
const ordinal_type subcellOrdinal=-1) const {
INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE( true, std::logic_error,
">>> ERROR (Basis::getValues): this method is not supported or should be overridden accordingly by derived classes.");
}

/** \brief Evaluation of a FEM basis on a <strong>reference cell</strong>.
Returns values of <var>operatorType</var> acting on FEM basis functions for a set of
Expand Down
17 changes: 17 additions & 0 deletions packages/intrepid2/src/Discretization/Basis/Intrepid2_HCURL_HEX_I1_FEM.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -185,6 +185,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
51 changes: 51 additions & 0 deletions packages/intrepid2/src/Discretization/Basis/Intrepid2_HCURL_HEX_I1_FEMDef.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -330,6 +330,57 @@ namespace Intrepid2 {

}

template<typename DT, typename OT, typename PT>
void
Basis_HCURL_HEX_I1_FEM<DT,OT,PT>::getScratchSpaceSize(
ordinal_type& perTeamSpaceSize,
ordinal_type& perThreadSpaceSize,
const PointViewType inputPoints,
const EOperator operatorType) const {
perTeamSpaceSize = 0;
perThreadSpaceSize = 0;
}

template<typename DT, typename OT, typename PT>
KOKKOS_INLINE_FUNCTION
void
Basis_HCURL_HEX_I1_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_HCURL_HEX_I1_FEM::getValues), The capability of selecting subsets of basis functions has not been implemented yet.");

(void) scratchStorage; //avoid unused variable warning

const int numPoints = inputPoints.extent(0);

switch(operatorType) {
case OPERATOR_VALUE:
Kokkos::parallel_for (Kokkos::TeamThreadRange (team_member, numPoints), [=] (ordinal_type& pt) {
auto output = Kokkos::subview( outputValues, Kokkos::ALL(), pt, Kokkos::ALL() );
const auto input = Kokkos::subview( inputPoints, pt, Kokkos::ALL() );
Impl::Basis_HCURL_HEX_I1_FEM::Serial<OPERATOR_VALUE>::getValues( output, input);
});
break;
case OPERATOR_CURL:
Kokkos::parallel_for (Kokkos::TeamThreadRange (team_member, numPoints), [=] (ordinal_type& pt) {
auto output = Kokkos::subview( outputValues, Kokkos::ALL(), pt, Kokkos::ALL() );
const auto input = Kokkos::subview( inputPoints, pt, Kokkos::ALL() );
Impl::Basis_HCURL_HEX_I1_FEM::Serial<OPERATOR_CURL>::getValues( output, input);
});
break;
default: {
INTREPID2_TEST_FOR_ABORT( true, ">>> ERROR: (Intrepid2::Basis_HCURL_HEX_I1_FEM::getValues), Operator Type not supported.");
}
}
}

}// namespace Intrepid2

#endif
32 changes: 25 additions & 7 deletions packages/intrepid2/src/Discretization/Basis/Intrepid2_HCURL_HEX_In_FEM.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -148,20 +148,21 @@ namespace Intrepid2 {
class Basis_HCURL_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_HCURL_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 @@ -184,6 +185,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
112 changes: 88 additions & 24 deletions packages/intrepid2/src/Discretization/Basis/Intrepid2_HCURL_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_HCURL_HEX_In_FEM::Serial<opType>::
Basis_HCURL_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,22 +44,22 @@ 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() + 3*cardLine*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_A(Kokkos::view_wrap(ptr1, vcprop), cardLine, npts);
viewType outputLine_B(Kokkos::view_wrap(ptr2, vcprop), cardLine, npts);
viewType outputBubble(Kokkos::view_wrap(ptr3, vcprop), cardBubble, npts);
ViewType workLine(Kokkos::view_wrap(ptr0, vcprop), cardLine, npts);
ViewType outputLine_A(Kokkos::view_wrap(ptr1, vcprop), cardLine, npts);
ViewType outputLine_B(Kokkos::view_wrap(ptr2, vcprop), cardLine, npts);
ViewType outputBubble(Kokkos::view_wrap(ptr3, vcprop), cardBubble, npts);

// tensor product
ordinal_type idx = 0;
Expand Down Expand Up @@ -142,12 +142,12 @@ namespace Intrepid2 {
auto ptr4 = work.data() + 4*cardLine*npts*dim_s;
auto ptr5 = work.data() + 5*cardLine*npts*dim_s;

viewType workLine(Kokkos::view_wrap(ptr0, vcprop), cardLine, npts);
viewType outputLine_A(Kokkos::view_wrap(ptr1, vcprop), cardLine, npts);
viewType outputLine_B(Kokkos::view_wrap(ptr2, vcprop), cardLine, npts);
viewType outputLine_DA(Kokkos::view_wrap(ptr3, vcprop), cardLine, npts, 1);
viewType outputLine_DB(Kokkos::view_wrap(ptr4, vcprop), cardLine, npts, 1);
viewType outputBubble(Kokkos::view_wrap(ptr5, vcprop), cardBubble, npts);
ViewType workLine(Kokkos::view_wrap(ptr0, vcprop), cardLine, npts);
ViewType outputLine_A(Kokkos::view_wrap(ptr1, vcprop), cardLine, npts);
ViewType outputLine_B(Kokkos::view_wrap(ptr2, vcprop), cardLine, npts);
ViewType outputLine_DA(Kokkos::view_wrap(ptr3, vcprop), cardLine, npts, 1);
ViewType outputLine_DB(Kokkos::view_wrap(ptr4, vcprop), cardLine, npts, 1);
ViewType outputBubble(Kokkos::view_wrap(ptr5, vcprop), cardBubble, npts);

// tensor product
ordinal_type idx = 0;
Expand Down Expand Up @@ -588,6 +588,70 @@ 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_HCURL_HEX_In_FEM<DT,OT,PT>::getScratchSpaceSize(
ordinal_type& perTeamSpaceSize,
ordinal_type& perThreadSpaceSize,
const PointViewType inputPoints,
const EOperator operatorType) const {
perTeamSpaceSize = 0;
ordinal_type scalarWorkViewExtent = (operatorType == OPERATOR_VALUE) ?
3*this->vinvLine_.extent(0)+this->vinvBubble_.extent(0):
5*this->vinvLine_.extent(0)+this->vinvBubble_.extent(0);
perThreadSpaceSize = scalarWorkViewExtent*get_dimension_scalar(inputPoints)*sizeof(typename BasisBase::scalarType);
}

template<typename DT, typename OT, typename PT>
KOKKOS_INLINE_FUNCTION
void
Basis_HCURL_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_HCURL_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 scalarSizePerPoint = (operatorType == OPERATOR_VALUE) ?
3*this->vinvLine_.extent(0)+this->vinvBubble_.extent(0):
5*this->vinvLine_.extent(0)+this->vinvBubble_.extent(0);
ordinal_type sizePerPoint = scalarSizePerPoint*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_HCURL_HEX_In_FEM::Serial<OPERATOR_VALUE>::getValues( output, input, work, this->vinvLine_, this->vinvBubble_ );
});
break;
case OPERATOR_CURL:
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_HCURL_HEX_In_FEM::Serial<OPERATOR_CURL>::getValues( output, input, work, this->vinvLine_, this->vinvBubble_ );
});
break;
default: {
INTREPID2_TEST_FOR_ABORT( true,
">>> ERROR (Basis_HCURL_HEX_In_FEM): getValues not implemented for this operator");
}
}
}

} // namespace Intrepid2

#endif
Loading

0 comments on commit f714191

Please sign in to comment.