1#ifndef _COMPADRE_FUNCTORS_HPP_
2#define _COMPADRE_FUNCTORS_HPP_
12#include "KokkosBatched_Gemm_Decl.hpp"
114 KOKKOS_INLINE_FUNCTION
117 const int local_index = teamMember.league_rank();
142 KOKKOS_INLINE_FUNCTION
148 const int local_index = teamMember.league_rank();
179 KOKKOS_INLINE_FUNCTION
187 const int local_index = teamMember.league_rank();
203 dimensions-1, dimensions);
208 for (
int j = 0; j < delta.extent_int(0); ++j) {
211 for (
int j = 0; j < thread_workspace.extent_int(0); ++j) {
212 thread_workspace(j) = 0;
224 dimensions, dimensions);
234 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
235 Kokkos::single(Kokkos::PerThread(teamMember), [&] () {
241 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
242 Kokkos::single(Kokkos::PerThread(teamMember), [&] () {
243 for (
int j=0; j<dimensions; ++j) {
244 for (
int k=0; k<dimensions-1; ++k) {
252 &&
"StaggeredEdgeIntegralSample prestencil weight only written for manifolds.");
253 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember,
255 Kokkos::single(Kokkos::PerThread(teamMember), [&] () {
256 for (int quadrature = 0; quadrature<_data._qm.getNumberOfQuadraturePoints(); ++quadrature) {
257 XYZ tangent_quadrature_coord_2d;
258 for (int j=0; j<dimensions-1; ++j) {
259 tangent_quadrature_coord_2d[j] = _data._pc.getTargetCoordinate(target_index, j, &T);
260 tangent_quadrature_coord_2d[j] -= _data._pc.getNeighborCoordinate(target_index, m, j, &T);
262 double tangent_vector[3];
263 tangent_vector[0] = tangent_quadrature_coord_2d[0]*T(0,0) + tangent_quadrature_coord_2d[1]*T(1,0);
264 tangent_vector[1] = tangent_quadrature_coord_2d[0]*T(0,1) + tangent_quadrature_coord_2d[1]*T(1,1);
265 tangent_vector[2] = tangent_quadrature_coord_2d[0]*T(0,2) + tangent_quadrature_coord_2d[1]*T(1,2);
267 for (int j=0; j<dimensions; ++j) {
268 _data._prestencil_weights(0,target_index,m,0,j) += (1-_data._qm.getSite(quadrature,0))*tangent_vector[j]*_data._qm.getWeight(quadrature);
269 _data._prestencil_weights(1,target_index,m,0,j) += _data._qm.getSite(quadrature,0)*tangent_vector[j]*_data._qm.getWeight(quadrature);
276 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember,
279 Kokkos::single(Kokkos::PerThread(teamMember), [&] () {
280 calcGradientPij(_data, teamMember, delta.data(), thread_workspace.data(), target_index,
281 m, 0 , 0 , dimensions-1, _data._curvature_poly_order,
282 false , &T, ReconstructionSpace::ScalarTaylorPolynomial,
286 double grad_xi1 = 0, grad_xi2 = 0;
287 Kokkos::parallel_reduce(Kokkos::ThreadVectorRange(teamMember,
290 for (int l=0; l<_data.manifold_NP; ++l) {
291 alpha_ij += delta(l)*Q(l,i);
294 double normal_coordinate = rel_coord[dimensions-1];
297 t_grad_xi1 += alpha_ij * normal_coordinate;
301 Kokkos::single(Kokkos::PerThread(teamMember), [&] () {
307 Kokkos::parallel_reduce(Kokkos::ThreadVectorRange(teamMember,
310 for (int l=0; l<_data.manifold_NP; ++l) {
311 alpha_ij += delta(l)*Q(l,i);
314 double normal_coordinate = rel_coord[dimensions-1];
317 if (dimensions>2) t_grad_xi2 += alpha_ij * normal_coordinate;
322 teamMember.team_barrier();
324 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember,
327 Kokkos::single(Kokkos::PerThread(teamMember), [&] () {
328 for (int j=0; j<dimensions; ++j) {
329 tangent(0,j) = t1(m)*T(dimensions-1,j) + T(0,j);
330 tangent(1,j) = t2(m)*T(dimensions-1,j) + T(1,j);
335 for (int j=0; j<dimensions; ++j) {
336 norm += tangent(0,j)*tangent(0,j);
340 norm = std::sqrt(norm);
341 for (int j=0; j<dimensions; ++j) {
342 tangent(0,j) /= norm;
346 if (dimensions-1 == 2) {
347 double dot_product = tangent(0,0)*tangent(1,0)
348 + tangent(0,1)*tangent(1,1)
349 + tangent(0,2)*tangent(1,2);
350 for (int j=0; j<dimensions; ++j) {
351 tangent(1,j) -= dot_product*tangent(0,j);
355 for (int j=0; j<dimensions; ++j) {
356 norm += tangent(1,j)*tangent(1,j);
358 norm = std::sqrt(norm);
359 for (int j=0; j<dimensions; ++j) {
360 tangent(1,j) /= norm;
365 for (int j=0; j<dimensions; ++j) {
366 for (int k=0; k<dimensions-1; ++k) {
367 _data._prestencil_weights(0,target_index,m,k,j) = tangent(k,j);
373 teamMember.team_barrier();
384 KOKKOS_INLINE_FUNCTION
392 const int local_index = teamMember.league_rank();
426 double * rhs_data = RHS.data();
427 Kokkos::parallel_for(Kokkos::TeamVectorRange(teamMember, this_num_rows), [&] (
const int i) {
428 rhs_data[i] = std::sqrt(w(i));
436 KokkosBatched::TeamVectorGemm<member_type,KokkosBatched::Trans::Transpose,KokkosBatched::Trans::NoTranspose,KokkosBatched::Algo::Gemm::Unblocked>
443 teamMember.team_barrier();
446 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, _data.
max_num_rows), [&] (
const int i) {
447 for (int j=0; j < _data.this_num_cols; j++) {
448 PsqrtW(i,j) = PsqrtW(i,j)*std::sqrt(w(i));
451 teamMember.team_barrier();
462 double cutoff_p = _data.
_epsilons(target_index);
468 teamMember.team_barrier();
485 KOKKOS_INLINE_FUNCTION
493 const int local_index = teamMember.league_rank();
508 dimensions, dimensions);
511 dimensions, dimensions);
527 KokkosBatched::TeamVectorGemm<member_type,KokkosBatched::Trans::Transpose,KokkosBatched::Trans::NoTranspose,KokkosBatched::Algo::Gemm::Unblocked>
534 teamMember.team_barrier();
537 GMLS_LinearAlgebra::largestTwoEigenvectorsThreeByThreeSymmetric(teamMember, T, PTP, dimensions,
538 const_cast<pool_type&
>(_random_number_pool));
540 teamMember.team_barrier();
551 KOKKOS_INLINE_FUNCTION
559 const int local_index = teamMember.league_rank();
596 double * rhs_data = RHS.data();
597 Kokkos::parallel_for(Kokkos::TeamVectorRange(teamMember, this_num_neighbors), [&] (
const int i) {
598 rhs_data[i] = std::sqrt(w(i));
607 KokkosBatched::TeamVectorGemm<
member_type,KokkosBatched::Trans::Transpose,
608 KokkosBatched::Trans::NoTranspose,KokkosBatched::Algo::Gemm::Unblocked>
615 teamMember.team_barrier();
618 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, this_num_neighbors), [&] (
const int i) {
620 CurvaturePsqrtW(i, j) = CurvaturePsqrtW(i, j)*std::sqrt(w(i));
624 teamMember.team_barrier();
635 KOKKOS_INLINE_FUNCTION
643 const int local_index = teamMember.league_rank();
655 dimensions, dimensions);
679 teamMember.team_barrier();
681 double grad_xi1 = 0, grad_xi2 = 0;
683 for (
int k=0; k<dimensions-1; ++k) {
686 Kokkos::parallel_reduce(Kokkos::TeamThreadRange(teamMember,
687 _data.
manifold_NP), [&] (
const int l,
double &talpha_ij) {
688 Kokkos::single(Kokkos::PerThread(teamMember), [&] () {
689 talpha_ij += P_target_row(offset,l)*Q(l,i);
692 teamMember.team_barrier();
693 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
695 manifold_gradient(i*(dimensions-1) + k) = alpha_ij;
698 teamMember.team_barrier();
701 double normal_coordinate = rel_coord[dimensions-1];
704 grad_xi1 += manifold_gradient(i*(dimensions-1)) * normal_coordinate;
705 if (dimensions>2) grad_xi2 += manifold_gradient(i*(dimensions-1)+1) * normal_coordinate;
706 teamMember.team_barrier();
710 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
712 double grad_xi[2] = {grad_xi1, grad_xi2};
716 for (
int i=0; i<dimensions-1; ++i) {
717 for (
int j=0; j<dimensions; ++j) {
721 for (
int j=0; j<dimensions; ++j) {
722 T(i,j) = grad_xi[i]*T(dimensions-1,j);
729 for (
int j=0; j<dimensions; ++j) {
730 norm += T(0,j)*T(0,j);
734 norm = std::sqrt(norm);
735 for (
int j=0; j<dimensions; ++j) {
740 if (dimensions-1 == 2) {
741 double dot_product = T(0,0)*T(1,0) + T(0,1)*T(1,1) + T(0,2)*T(1,2);
742 for (
int j=0; j<dimensions; ++j) {
743 T(1,j) -= dot_product*T(0,j);
747 for (
int j=0; j<dimensions; ++j) {
748 norm += T(1,j)*T(1,j);
750 norm = std::sqrt(norm);
751 for (
int j=0; j<dimensions; ++j) {
757 double norm_t_normal = 0;
759 T(dimensions-1,0) = T(0,1)*T(1,2) - T(1,1)*T(0,2);
760 norm_t_normal += T(dimensions-1,0)*T(dimensions-1,0);
761 T(dimensions-1,1) = -(T(0,0)*T(1,2) - T(1,0)*T(0,2));
762 norm_t_normal += T(dimensions-1,1)*T(dimensions-1,1);
763 T(dimensions-1,2) = T(0,0)*T(1,1) - T(1,0)*T(0,1);
764 norm_t_normal += T(dimensions-1,2)*T(dimensions-1,2);
766 T(dimensions-1,0) = T(1,1) - T(0,1);
767 norm_t_normal += T(dimensions-1,0)*T(dimensions-1,0);
768 T(dimensions-1,1) = T(0,0) - T(1,0);
769 norm_t_normal += T(dimensions-1,1)*T(dimensions-1,1);
771 norm_t_normal = std::sqrt(norm_t_normal);
772 for (
int i=0; i<dimensions-1; ++i) {
773 T(dimensions-1,i) /= norm_t_normal;
776 teamMember.team_barrier();
788 KOKKOS_INLINE_FUNCTION
808 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
810 &&
"FixTangentDirectionOrder called on manifold with a dimension of 0.");
811 double dot_product = 0;
812 for (
int i=0; i<dimensions; ++i) {
813 dot_product += T(dimensions-1,i) * N(i);
816 if (dot_product < 0) {
818 for (
int i=0; i<dimensions; ++i) {
825 for (
int i=0; i<dimensions; ++i) {
827 T(dimensions-1,i) *= -1;
832 teamMember.team_barrier();
843 KOKKOS_INLINE_FUNCTION
851 const int local_index = teamMember.league_rank();
864 dimensions, dimensions);
886 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
888 manifold_coeffs(j) = 0;
891 teamMember.team_barrier();
894 double normal_coordinate = rel_coord[dimensions-1];
896 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
897 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
900 manifold_coeffs(j) += Q(j,i) * normal_coordinate;
904 teamMember.team_barrier();
916 KOKKOS_INLINE_FUNCTION
924 const int local_index = teamMember.league_rank();
953 createWeightsAndP(_data, teamMember, delta, thread_workspace, PsqrtW, w, dimensions-1,
959 double * Q_data = Q.data();
960 Kokkos::parallel_for(Kokkos::TeamVectorRange(teamMember,this_num_rows), [&] (
const int i) {
961 Q_data[i] = std::sqrt(w(i));
971 KokkosBatched::TeamVectorGemm<member_type,KokkosBatched::Trans::Transpose,KokkosBatched::Trans::NoTranspose,KokkosBatched::Algo::Gemm::Unblocked>
978 teamMember.team_barrier();
981 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, _data.
max_num_rows), [&] (
const int i) {
982 Kokkos::parallel_for(Kokkos::ThreadVectorRange(teamMember, _data.this_num_cols), [&] (const int j) {
983 PsqrtW(i, j) = PsqrtW(i, j)*std::sqrt(w(i));
987 teamMember.team_barrier();
998 KOKKOS_INLINE_FUNCTION
1006 const int local_index = teamMember.league_rank();
#define compadre_kernel_assert_debug(condition)
Kokkos::View< double *, Kokkos::MemoryTraits< Kokkos::Unmanaged > > scratch_vector_type
Kokkos::Random_XorShift64_Pool pool_type
team_policy::member_type member_type
#define TO_GLOBAL(variable)
Kokkos::View< double **, layout_right, Kokkos::MemoryTraits< Kokkos::Unmanaged > > scratch_matrix_right_type
KOKKOS_INLINE_FUNCTION int getThreadScratchLevel(const int level) const
KOKKOS_INLINE_FUNCTION int getTeamScratchLevel(const int level) const
KOKKOS_INLINE_FUNCTION void computeCurvatureFunctionals(const TargetData &data, const member_type &teamMember, scratch_vector_type delta, scratch_vector_type thread_workspace, scratch_matrix_right_type P_target_row, const scratch_matrix_right_type *V, const local_index_type local_neighbor_index=-1)
Evaluates a polynomial basis for the curvature with a gradient target functional applied.
ProblemType
Problem type, that optionally can handle manifolds.
@ MANIFOLD
Solve GMLS problem on a manifold (will use QR or SVD to solve the resultant GMLS problem dependent on...
constexpr SamplingFunctional VaryingManifoldVectorPointSample
For integrating polynomial dotted with normal over an edge.
KOKKOS_INLINE_FUNCTION void computeTargetFunctionals(const TargetData &data, const member_type &teamMember, scratch_vector_type delta, scratch_vector_type thread_workspace, scratch_matrix_right_type P_target_row)
Evaluates a polynomial basis with a target functional applied to each member of the basis.
ConstraintType
Constraint type.
constexpr SamplingFunctional StaggeredEdgeIntegralSample
Samples consist of the result of integrals of a vector dotted with the tangent along edges between ne...
KOKKOS_INLINE_FUNCTION void createWeightsAndP(const BasisData &data, const member_type &teamMember, scratch_vector_type delta, scratch_vector_type thread_workspace, scratch_matrix_right_type P, scratch_vector_type w, const int dimension, int polynomial_order, bool weight_p=false, scratch_matrix_right_type *V=NULL, const ReconstructionSpace reconstruction_space=ReconstructionSpace::ScalarTaylorPolynomial, const SamplingFunctional polynomial_sampling_functional=PointSample)
Fills the _P matrix with either P or P*sqrt(w)
constexpr SamplingFunctional PointSample
Available sampling functionals.
constexpr SamplingFunctional StaggeredEdgeAnalyticGradientIntegralSample
Analytical integral of a gradient source vector is just a difference of the scalar source at neighbor...
KOKKOS_INLINE_FUNCTION void createWeightsAndPForCurvature(const BasisData &data, const member_type &teamMember, scratch_vector_type delta, scratch_vector_type thread_workspace, scratch_matrix_right_type P, scratch_vector_type w, const int dimension, bool only_specific_order, scratch_matrix_right_type *V=NULL)
Fills the _P matrix with P*sqrt(w) for use in solving for curvature.
DenseSolverType
Dense solver type.
WeightingFunctionType
Available weighting kernel function types.
constexpr SamplingFunctional ManifoldVectorPointSample
Point evaluations of the entire vector source function (but on a manifold, so it includes a transform...
KOKKOS_INLINE_FUNCTION void applyTargetsToCoefficients(const SolutionData &data, const member_type &teamMember, scratch_matrix_right_type Q, scratch_matrix_right_type P_target_row)
For applying the evaluations from a target functional to the polynomial coefficients.
KOKKOS_INLINE_FUNCTION void evaluateConstraints(scratch_matrix_right_type M, scratch_matrix_right_type PsqrtW, const ConstraintType constraint_type, const ReconstructionSpace reconstruction_space, const int NP, const double cutoff_p, const int dimension, const int num_neighbors=0, scratch_matrix_right_type *T=NULL)
KOKKOS_INLINE_FUNCTION void calcGradientPij(const BasisData &data, const member_type &teamMember, double *delta, double *thread_workspace, const int target_index, int neighbor_index, const double alpha, const int partial_direction, const int dimension, const int poly_order, bool specific_order_only, const scratch_matrix_right_type *V, const ReconstructionSpace reconstruction_space, const SamplingFunctional polynomial_sampling_functional, const int evaluation_site_local_index=0)
Evaluates the gradient of a polynomial basis under the Dirac Delta (pointwise) sampling function.
KOKKOS_INLINE_FUNCTION void computeTargetFunctionalsOnManifold(const TargetData &data, const member_type &teamMember, scratch_vector_type delta, scratch_vector_type thread_workspace, scratch_matrix_right_type P_target_row, scratch_matrix_right_type V, scratch_vector_type curvature_coefficients)
Evaluates a polynomial basis with a target functional applied, using information from the manifold cu...
ReconstructionSpace
Space in which to reconstruct polynomial.
@ ScalarTaylorPolynomial
Scalar polynomial basis centered at the target site and scaled by sum of basis powers e....
Functor to evaluate curvature targets and apply to coefficients of curvature reconstruction.
KOKKOS_INLINE_FUNCTION void operator()(const member_type &teamMember) const
ApplyCurvatureTargets(GMLSBasisData data)
Functor to apply target evaluation to polynomial coefficients to store in _alphas.
KOKKOS_INLINE_FUNCTION void operator()(const member_type &teamMember) const
ApplyTargets(GMLSSolutionData data)
Functor to assemble the P*sqrt(weights) matrix and construct sqrt(weights)*Identity for curvature.
AssembleCurvaturePsqrtW(GMLSBasisData data)
KOKKOS_INLINE_FUNCTION void operator()(const member_type &teamMember) const
Functor to assemble the P*sqrt(weights) matrix and construct sqrt(weights)*Identity.
AssembleManifoldPsqrtW(GMLSBasisData data)
KOKKOS_INLINE_FUNCTION void operator()(const member_type &teamMember) const
Functor to assemble the P*sqrt(weights) matrix and construct sqrt(weights)*Identity.
KOKKOS_INLINE_FUNCTION void operator()(const member_type &teamMember) const
AssembleStandardPsqrtW(GMLSBasisData data)
Functor to create a coarse tangent approximation from a given neighborhood of points.
ComputeCoarseTangentPlane(GMLSBasisData data)
KOKKOS_INLINE_FUNCTION void operator()(const member_type &teamMember) const
pool_type _random_number_pool
Functor to calculate prestencil weights to apply to data to transform into a format expected by a GML...
KOKKOS_INLINE_FUNCTION void operator()(const member_type &teamMember) const
ComputePrestencilWeights(GMLSBasisData data)
Functor to evaluate targets on a manifold.
EvaluateManifoldTargets(GMLSBasisData data)
KOKKOS_INLINE_FUNCTION void operator()(const member_type &teamMember) const
Functor to evaluate targets operations on the basis.
KOKKOS_INLINE_FUNCTION void operator()(const member_type &teamMember) const
EvaluateStandardTargets(GMLSBasisData data)
Functor to determine if tangent directions need reordered, and to reorder them if needed We require t...
KOKKOS_INLINE_FUNCTION void operator()(const member_type &teamMember) const
FixTangentDirectionOrdering(GMLSBasisData data)
double * P_target_row_data
int _reconstruction_space_rank
int manifold_gradient_dim
GMLS::point_connections_type _pc
SamplingFunctional _polynomial_sampling_functional
int _initial_index_for_batch
ReconstructionSpace _reconstruction_space
DenseSolverType _dense_solver_type
SolutionSet< device_memory_space > _d_ss
Kokkos::View< double **, layout_right > _source_extra_data
double * manifold_curvature_coefficients_data
ProblemType _problem_type
int _data_sampling_multiplier
WeightingFunctionType _curvature_weighting_type
int _curvature_weighting_n
int _order_of_quadrature_points
GMLS::point_connections_type _additional_pc
Kokkos::View< double *****, layout_right > _prestencil_weights
Kokkos::View< TargetOperation * > _curvature_support_operations
Kokkos::View< double * > _epsilons
Kokkos::View< TargetOperation * > _operations
WeightingFunctionType _weighting_type
Kokkos::View< double **, layout_right > _target_extra_data
SamplingFunctional _data_sampling_functional
int _curvature_poly_order
int _curvature_weighting_p
ConstraintType _constraint_type
int _dimension_of_quadrature_points
double * P_target_row_data
Kokkos::View< int * > additional_number_of_neighbors_list
Kokkos::View< int * > number_of_neighbors_list
SolutionSet< device_memory_space > _d_ss
int _initial_index_for_batch
Functor to evaluate curvature targets and construct accurate tangent direction approximation for mani...
KOKKOS_INLINE_FUNCTION void operator()(const member_type &teamMember) const
GetAccurateTangentDirections(GMLSBasisData data)
KOKKOS_INLINE_FUNCTION int getNumberOfNeighborsDevice(int target_index) const
Get number of neighbors for a given target (device)
KOKKOS_INLINE_FUNCTION XYZ getRelativeCoord(const int target_index, const int neighbor_list_num, const int dimension, const scratch_matrix_right_type *V=NULL) const
Returns the relative coordinate as a vector between the target site and the neighbor site.
All vairables and functionality related to the layout and storage of GMLS solutions (alpha values)
KOKKOS_INLINE_FUNCTION int getTargetOffsetIndex(const int lro_num, const int input_component, const int output_component, const int evaluation_site_local_index=0) const
Handles offset from operation input/output + extra evaluation sites.