Compadre 1.5.5
Loading...
Searching...
No Matches
Compadre_Functors.hpp
Go to the documentation of this file.
1#ifndef _COMPADRE_FUNCTORS_HPP_
2#define _COMPADRE_FUNCTORS_HPP_
3
6#include "Compadre_Basis.hpp"
10#include "Compadre_Functors.hpp"
12#include "KokkosBatched_Gemm_Decl.hpp"
13#include "Compadre_GMLS.hpp"
14
15namespace Compadre {
16
18
19 Kokkos::View<double**, layout_right> _source_extra_data;
20 Kokkos::View<double**, layout_right> _target_extra_data;
21 Kokkos::View<double*> _epsilons;
22 Kokkos::View<double*****, layout_right> _prestencil_weights;
23 Kokkos::View<TargetOperation*> _curvature_support_operations;
24 Kokkos::View<TargetOperation*> _operations;
25
28 int _NP;
43
57
58 // convenience variables (not from GMLS class)
60 double * RHS_data;
62 double * P_data;
66 double * Coeffs_data;
67 double * w_data;
74 double * T_data;
76 double * ref_N_data;
80
82};
83
101
102/** @name Functors
103 * Member functions that perform operations on the entire batch
104 */
105///@{
106
107//! Functor to apply target evaluation to polynomial coefficients to store in _alphas
109
111
113
114 KOKKOS_INLINE_FUNCTION
115 void operator()(const member_type& teamMember) const {
116
117 const int local_index = teamMember.league_rank();
118
119 /*
120 * Data
121 */
122
123 // Coefficients for polynomial basis have overwritten _data._RHS
130
131 applyTargetsToCoefficients(_data, teamMember, Coeffs, P_target_row);
132 }
133};
134
135//! Functor to evaluate targets operations on the basis
137
139
141
142 KOKKOS_INLINE_FUNCTION
143 void operator()(const member_type& teamMember) const {
144 /*
145 * Dimensions
146 */
147
148 const int local_index = teamMember.league_rank();
149
150 /*
151 * Data
152 */
153
157
158 scratch_vector_type delta(teamMember.thread_scratch(_data._pm.getThreadScratchLevel(0)),
160 scratch_vector_type thread_workspace(teamMember.thread_scratch(_data._pm.getThreadScratchLevel(0)),
162
163 /*
164 * Evaluate Standard Targets
165 */
166
167 computeTargetFunctionals(_data, teamMember, delta, thread_workspace, P_target_row);
168
169 }
170};
171
172//! Functor to calculate prestencil weights to apply to data to transform into a format expected by a GMLS stencil
174
176
178
179 KOKKOS_INLINE_FUNCTION
180 void operator()(const member_type& teamMember) const {
181
182 /*
183 * Dimensions
184 */
185
186 const int target_index = _data._initial_index_for_batch + teamMember.league_rank();
187 const int local_index = teamMember.league_rank();
188 const int dimensions = _data._dimensions;
189
190 /*
191 * Data
192 */
193
194 scratch_vector_type delta(teamMember.thread_scratch(_data._pm.getThreadScratchLevel(0)),
196 scratch_vector_type thread_workspace(teamMember.thread_scratch(_data._pm.getThreadScratchLevel(0)),
198
200 scratch_vector_type t1, t2;
202 tangent = scratch_matrix_right_type(teamMember.thread_scratch(_data._pm.getThreadScratchLevel(1)),
203 dimensions-1, dimensions);
204 t1 = scratch_vector_type(teamMember.team_scratch(_data._pm.getTeamScratchLevel(1)),
206 t2 = scratch_vector_type(teamMember.team_scratch(_data._pm.getTeamScratchLevel(1)),
208 for (int j = 0; j < delta.extent_int(0); ++j) {
209 delta(j) = 0;
210 }
211 for (int j = 0; j < thread_workspace.extent_int(0); ++j) {
212 thread_workspace(j) = 0;
213 }
214 }
215
216
217 // holds polynomial coefficients for curvature reconstruction
221
223 + TO_GLOBAL(target_index)*TO_GLOBAL(dimensions)*TO_GLOBAL(dimensions),
224 dimensions, dimensions);
225
226 scratch_vector_type manifold_gradient(teamMember.team_scratch(_data._pm.getTeamScratchLevel(1)),
228
229 /*
230 * Prestencil Weight Calculations
231 */
232
234 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
235 Kokkos::single(Kokkos::PerThread(teamMember), [&] () {
236 _data._prestencil_weights(0,0,0,0,0) = -1;
237 _data._prestencil_weights(1,0,0,0,0) = 1;
238 });
239 });
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) {
245 _data._prestencil_weights(0,target_index,0,k,j) = T(k,j);
246 }
247 }
248 });
249 });
252 && "StaggeredEdgeIntegralSample prestencil weight only written for manifolds.");
253 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember,
254 _data._pc._nla.getNumberOfNeighborsDevice(target_index)), [&] (const int m) {
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);
261 }
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);
266
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);
270 }
271 }
272 });
273 });
275
276 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember,
277 _data._pc._nla.getNumberOfNeighborsDevice(target_index)), [&] (const int m) {
278
279 Kokkos::single(Kokkos::PerThread(teamMember), [&] () {
280 calcGradientPij(_data, teamMember, delta.data(), thread_workspace.data(), target_index,
281 m, 0 /*alpha*/, 0 /*partial_direction*/, dimensions-1, _data._curvature_poly_order,
282 false /*specific order only*/, &T, ReconstructionSpace::ScalarTaylorPolynomial,
283 PointSample);
284 });
285 // reconstructs gradient at local neighbor index m
286 double grad_xi1 = 0, grad_xi2 = 0;
287 Kokkos::parallel_reduce(Kokkos::ThreadVectorRange(teamMember,
288 _data._pc._nla.getNumberOfNeighborsDevice(target_index)), [&] (const int i, double &t_grad_xi1) {
289 double alpha_ij = 0;
290 for (int l=0; l<_data.manifold_NP; ++l) {
291 alpha_ij += delta(l)*Q(l,i);
292 }
293 XYZ rel_coord = _data._pc.getRelativeCoord(target_index, i, dimensions, &T);
294 double normal_coordinate = rel_coord[dimensions-1];
295
296 // apply coefficients to sample data
297 t_grad_xi1 += alpha_ij * normal_coordinate;
298 }, grad_xi1);
299 t1(m) = grad_xi1;
300
301 Kokkos::single(Kokkos::PerThread(teamMember), [&] () {
302 calcGradientPij(_data, teamMember, delta.data(), thread_workspace.data(), target_index,
303 m, 0 /*alpha*/, 1 /*partial_direction*/, dimensions-1, _data._curvature_poly_order,
304 false /*specific order only*/, &T, ReconstructionSpace::ScalarTaylorPolynomial,
306 });
307 Kokkos::parallel_reduce(Kokkos::ThreadVectorRange(teamMember,
308 _data._pc._nla.getNumberOfNeighborsDevice(target_index)), [&] (const int i, double &t_grad_xi2) {
309 double alpha_ij = 0;
310 for (int l=0; l<_data.manifold_NP; ++l) {
311 alpha_ij += delta(l)*Q(l,i);
312 }
313 XYZ rel_coord = _data._pc.getRelativeCoord(target_index, i, dimensions, &T);
314 double normal_coordinate = rel_coord[dimensions-1];
315
316 // apply coefficients to sample data
317 if (dimensions>2) t_grad_xi2 += alpha_ij * normal_coordinate;
318 }, grad_xi2);
319 t2(m) = grad_xi2;
320 });
321
322 teamMember.team_barrier();
323
324 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember,
325 _data._pc._nla.getNumberOfNeighborsDevice(target_index)), [&] (const int m) {
326 // constructs tangent vector at neighbor site
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);
331 }
332
333 // calculate norm
334 double norm = 0;
335 for (int j=0; j<dimensions; ++j) {
336 norm += tangent(0,j)*tangent(0,j);
337 }
338
339 // normalize first vector
340 norm = std::sqrt(norm);
341 for (int j=0; j<dimensions; ++j) {
342 tangent(0,j) /= norm;
343 }
344
345 // orthonormalize next vector
346 if (dimensions-1 == 2) { // 2d manifold
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);
352 }
353 // normalize second vector
354 norm = 0;
355 for (int j=0; j<dimensions; ++j) {
356 norm += tangent(1,j)*tangent(1,j);
357 }
358 norm = std::sqrt(norm);
359 for (int j=0; j<dimensions; ++j) {
360 tangent(1,j) /= norm;
361 }
362 }
363
364 // stores matrix of tangent and normal directions as a prestencil weight
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);
368 }
369 }
370 });
371 });
372 }
373 teamMember.team_barrier();
374 }
375};
376
377//! Functor to assemble the P*sqrt(weights) matrix and construct sqrt(weights)*Identity
379
381
383
384 KOKKOS_INLINE_FUNCTION
385 void operator()(const member_type& teamMember) const {
386
387 /*
388 * Dimensions
389 */
390
391 const int target_index = _data._initial_index_for_batch + teamMember.league_rank();
392 const int local_index = teamMember.league_rank();
393 const int this_num_rows = _data._sampling_multiplier*_data._pc._nla.getNumberOfNeighborsDevice(target_index);
394
395 /*
396 * Data
397 */
398
400 + TO_GLOBAL(local_index)*TO_GLOBAL(_data.P_dim_0*_data.P_dim_1),
401 _data.P_dim_0, _data.P_dim_1);
403 + TO_GLOBAL(local_index)*TO_GLOBAL(_data.RHS_dim_0*_data.RHS_dim_1),
404 _data.RHS_dim_0, _data.RHS_dim_1);
406 + TO_GLOBAL(local_index)*TO_GLOBAL(_data.max_num_rows),
407 _data.max_num_rows);
408
409 // delta, used for each thread
410 scratch_vector_type delta(teamMember.thread_scratch(_data._pm.getThreadScratchLevel(0)),
411 _data.this_num_cols);
412 scratch_vector_type thread_workspace(teamMember.thread_scratch(_data._pm.getThreadScratchLevel(0)),
414
415 /*
416 * Assemble P*sqrt(W) and sqrt(w)*Identity
417 */
418
419 // creates the matrix sqrt(W)*P
420 createWeightsAndP(_data, teamMember, delta, thread_workspace, PsqrtW, w, _data._dimensions,
421 _data._poly_order, true /*weight_p*/, NULL /*&V*/, _data._reconstruction_space,
423
424 if ((_data._constraint_type == ConstraintType::NO_CONSTRAINT) && (_data._dense_solver_type != DenseSolverType::LU)) {
425 // fill in RHS with Identity * sqrt(weights)
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));
429 });
430 } else {
431 // create global memory for matrix M = PsqrtW^T*PsqrtW
432 // don't need to cast into scratch_matrix_left_type since the matrix is symmetric
434 + TO_GLOBAL(local_index)*TO_GLOBAL(_data.RHS_dim_0*_data.RHS_dim_1),
435 _data.RHS_dim_0, _data.RHS_dim_1);
436 KokkosBatched::TeamVectorGemm<member_type,KokkosBatched::Trans::Transpose,KokkosBatched::Trans::NoTranspose,KokkosBatched::Algo::Gemm::Unblocked>
437 ::invoke(teamMember,
438 1.0,
439 PsqrtW,
440 PsqrtW,
441 0.0,
442 M);
443 teamMember.team_barrier();
444
445 // Multiply PsqrtW with sqrt(W) to get PW
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));
449 }
450 });
451 teamMember.team_barrier();
452
453 // conditionally fill in rows determined by constraint type
454 if (_data._constraint_type == ConstraintType::NEUMANN_GRAD_SCALAR) {
455 // normal vector is contained in last row of T
457 + TO_GLOBAL(target_index)*TO_GLOBAL(_data._dimensions*_data._dimensions),
458 _data._dimensions, _data._dimensions);
459
460 // Get the number of neighbors for target index
461 int num_neigh_target = _data._pc._nla.getNumberOfNeighborsDevice(target_index);
462 double cutoff_p = _data._epsilons(target_index);
463
465 _data._NP, cutoff_p, _data._dimensions, num_neigh_target, &T);
466 }
467 }
468 teamMember.team_barrier();
469 }
470};
471
472//! Functor to create a coarse tangent approximation from a given neighborhood of points
474
476
477 // random number generator pool
479
481 // seed random number generator pool
482 _random_number_pool = pool_type(1);
483 }
484
485 KOKKOS_INLINE_FUNCTION
486 void operator()(const member_type& teamMember) const {
487
488 /*
489 * Dimensions
490 */
491
492 const int target_index = _data._initial_index_for_batch + teamMember.league_rank();
493 const int local_index = teamMember.league_rank();
494 const int dimensions = _data._dimensions;
495
496 /*
497 * Data
498 */
499
501 + TO_GLOBAL(local_index)*TO_GLOBAL(_data.P_dim_0*_data.P_dim_1),
502 _data.P_dim_0, _data.P_dim_1);
504 + TO_GLOBAL(local_index)*TO_GLOBAL(_data.max_num_rows),
505 _data.max_num_rows);
507 + TO_GLOBAL(target_index)*TO_GLOBAL(dimensions*dimensions),
508 dimensions, dimensions);
509
510 scratch_matrix_right_type PTP(teamMember.team_scratch(_data._pm.getTeamScratchLevel(1)),
511 dimensions, dimensions);
512
513 // delta, used for each thread
514 scratch_vector_type delta(teamMember.thread_scratch(_data._pm.getThreadScratchLevel(0)),
515 _data.this_num_cols);
516 scratch_vector_type thread_workspace(teamMember.thread_scratch(_data._pm.getThreadScratchLevel(0)),
518
519 /*
520 * Determine Coarse Approximation of Manifold Tangent Plane
521 */
522
523 // getting x y and z from which to derive a manifold
524 createWeightsAndPForCurvature(_data, teamMember, delta, thread_workspace, PsqrtW, w, dimensions, true /* only specific order */);
525
526 // create PsqrtW^T*PsqrtW
527 KokkosBatched::TeamVectorGemm<member_type,KokkosBatched::Trans::Transpose,KokkosBatched::Trans::NoTranspose,KokkosBatched::Algo::Gemm::Unblocked>
528 ::invoke(teamMember,
529 1.0,
530 PsqrtW,
531 PsqrtW,
532 0.0,
533 PTP);
534 teamMember.team_barrier();
535
536 // create coarse approximation of tangent plane in first two rows of T, with normal direction in third column
537 GMLS_LinearAlgebra::largestTwoEigenvectorsThreeByThreeSymmetric(teamMember, T, PTP, dimensions,
538 const_cast<pool_type&>(_random_number_pool));
539
540 teamMember.team_barrier();
541 }
542};
543
544//! Functor to assemble the P*sqrt(weights) matrix and construct sqrt(weights)*Identity for curvature
546
548
550
551 KOKKOS_INLINE_FUNCTION
552 void operator()(const member_type& teamMember) const {
553
554 /*
555 * Dimensions
556 */
557
558 const int target_index = _data._initial_index_for_batch + teamMember.league_rank();
559 const int local_index = teamMember.league_rank();
560 const int this_num_neighbors = _data._pc._nla.getNumberOfNeighborsDevice(target_index);
561
562 /*
563 * Data
564 */
565
566 scratch_matrix_right_type CurvaturePsqrtW(_data.P_data
567 + TO_GLOBAL(local_index)*TO_GLOBAL(_data.P_dim_0*_data.P_dim_1),
568 _data.P_dim_0, _data.P_dim_1);
570 + TO_GLOBAL(local_index)*TO_GLOBAL(_data.RHS_dim_0*_data.RHS_dim_1),
571 _data.RHS_dim_0, _data.RHS_dim_1);
573 + TO_GLOBAL(local_index)*TO_GLOBAL(_data.max_num_rows), _data.max_num_rows);
575 + TO_GLOBAL(target_index)*TO_GLOBAL(_data._dimensions*_data._dimensions),
576 _data._dimensions, _data._dimensions);
577
578 // delta, used for each thread
579 scratch_vector_type delta(teamMember.thread_scratch(_data._pm.getThreadScratchLevel(0)),
580 _data.this_num_cols);
581 scratch_vector_type thread_workspace(teamMember.thread_scratch(_data._pm.getThreadScratchLevel(0)),
583
584 //
585 // RECONSTRUCT ON THE TANGENT PLANE USING LOCAL COORDINATES
586 //
587
588 // creates the matrix sqrt(W)*P
589 createWeightsAndPForCurvature(_data, teamMember, delta, thread_workspace, CurvaturePsqrtW, w,
590 _data._dimensions-1, false /* only specific order */, &T);
591
592 // CurvaturePsqrtW is sized according to max_num_rows x this_num_cols of which in this case
593 // we are only using this_num_neighbors x manifold_NP
594 if (_data._dense_solver_type != DenseSolverType::LU) {
595 // fill in RHS with Identity * sqrt(weights)
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));
599 });
600 } else {
601 // create global memory for matrix M = PsqrtW^T*PsqrtW
602 // don't need to cast into scratch_matrix_left_type since the matrix is symmetric
604 + TO_GLOBAL(local_index)*TO_GLOBAL(_data.RHS_dim_0*_data.RHS_dim_1),
605 _data.RHS_dim_0, _data.RHS_dim_1);
606 // Assemble matrix M
607 KokkosBatched::TeamVectorGemm<member_type,KokkosBatched::Trans::Transpose,
608 KokkosBatched::Trans::NoTranspose,KokkosBatched::Algo::Gemm::Unblocked>
609 ::invoke(teamMember,
610 1.0,
611 CurvaturePsqrtW,
612 CurvaturePsqrtW,
613 0.0,
614 M);
615 teamMember.team_barrier();
616
617 // Multiply PsqrtW with sqrt(W) to get PW
618 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, this_num_neighbors), [&] (const int i) {
619 for (int j=0; j < _data.manifold_NP; j++) {
620 CurvaturePsqrtW(i, j) = CurvaturePsqrtW(i, j)*std::sqrt(w(i));
621 }
622 });
623 }
624 teamMember.team_barrier();
625 }
626};
627
628//! Functor to evaluate curvature targets and construct accurate tangent direction approximation for manifolds
630
632
634
635 KOKKOS_INLINE_FUNCTION
636 void operator()(const member_type& teamMember) const {
637
638 /*
639 * Dimensions
640 */
641
642 const int target_index = _data._initial_index_for_batch + teamMember.league_rank();
643 const int local_index = teamMember.league_rank();
644 auto dimensions = _data._dimensions;
645
646 /*
647 * Data
648 */
649
651 + TO_GLOBAL(local_index)*TO_GLOBAL(_data.Coeffs_dim_0*_data.Coeffs_dim_1),
652 _data.Coeffs_dim_0, _data.Coeffs_dim_1);
654 + TO_GLOBAL(target_index)*TO_GLOBAL(dimensions*dimensions),
655 dimensions, dimensions);
657 + TO_GLOBAL(local_index)*TO_GLOBAL(_data.P_target_row_dim_0*_data.P_target_row_dim_1),
659
660 scratch_vector_type manifold_gradient(teamMember.team_scratch(_data._pm.getTeamScratchLevel(1)),
662
663 // delta, used for each thread
664 scratch_vector_type delta(teamMember.thread_scratch(_data._pm.getThreadScratchLevel(0)),
665 _data.this_num_cols);
666 scratch_vector_type thread_workspace(teamMember.thread_scratch(_data._pm.getThreadScratchLevel(0)),
668
669 /*
670 * Manifold
671 */
672
673
674 //
675 // GET TARGET COEFFICIENTS RELATED TO GRADIENT TERMS
676 //
677 // reconstruct grad_xi1 and grad_xi2, not used for manifold_coeffs
678 computeCurvatureFunctionals(_data, teamMember, delta, thread_workspace, P_target_row, &T);
679 teamMember.team_barrier();
680
681 double grad_xi1 = 0, grad_xi2 = 0;
682 for (int i=0; i<_data._pc._nla.getNumberOfNeighborsDevice(target_index); ++i) {
683 for (int k=0; k<dimensions-1; ++k) {
684 double alpha_ij = 0;
685 int offset = _data._d_ss.getTargetOffsetIndex(0, 0, k, 0);
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);
690 });
691 }, alpha_ij);
692 teamMember.team_barrier();
693 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
694 // stored staggered, grad_xi1, grad_xi2, grad_xi1, grad_xi2, ....
695 manifold_gradient(i*(dimensions-1) + k) = alpha_ij;
696 });
697 }
698 teamMember.team_barrier();
699
700 XYZ rel_coord = _data._pc.getRelativeCoord(target_index, i, dimensions, &T);
701 double normal_coordinate = rel_coord[dimensions-1];
702
703 // apply coefficients to sample data
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();
707 }
708
709 // Constructs high order orthonormal tangent space T and inverse of metric tensor
710 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
711
712 double grad_xi[2] = {grad_xi1, grad_xi2};
713 double T_row[3];
714
715 // Construct T (high order approximation of orthonormal tangent vectors)
716 for (int i=0; i<dimensions-1; ++i) {
717 for (int j=0; j<dimensions; ++j) {
718 T_row[j] = T(i,j);
719 }
720 // build
721 for (int j=0; j<dimensions; ++j) {
722 T(i,j) = grad_xi[i]*T(dimensions-1,j);
723 T(i,j) += T_row[j];
724 }
725 }
726
727 // calculate norm
728 double norm = 0;
729 for (int j=0; j<dimensions; ++j) {
730 norm += T(0,j)*T(0,j);
731 }
732
733 // normalize first vector
734 norm = std::sqrt(norm);
735 for (int j=0; j<dimensions; ++j) {
736 T(0,j) /= norm;
737 }
738
739 // orthonormalize next vector
740 if (dimensions-1 == 2) { // 2d manifold
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);
744 }
745 // normalize second vector
746 norm = 0;
747 for (int j=0; j<dimensions; ++j) {
748 norm += T(1,j)*T(1,j);
749 }
750 norm = std::sqrt(norm);
751 for (int j=0; j<dimensions; ++j) {
752 T(1,j) /= norm;
753 }
754 }
755
756 // get normal vector to first two rows of T
757 double norm_t_normal = 0;
758 if (dimensions>2) {
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);
765 } else {
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);
770 }
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;
774 }
775 });
776 teamMember.team_barrier();
777 }
778};
779
780//! Functor to determine if tangent directions need reordered, and to reorder them if needed
781//! We require that the normal is consistent with a right hand rule on the tangent vectors
783
785
787
788 KOKKOS_INLINE_FUNCTION
789 void operator()(const member_type& teamMember) const {
790
791 /*
792 * Dimensions
793 */
794
795 const int target_index = _data._initial_index_for_batch + teamMember.league_rank();
796 auto dimensions = _data._dimensions;
797
798 /*
799 * Data
800 */
801
802 scratch_matrix_right_type T(_data.T_data + target_index*dimensions*dimensions, dimensions, dimensions);
803 scratch_vector_type N(_data.ref_N_data + target_index*dimensions, dimensions);
804
805 // take the dot product of the calculated normal in the tangent bundle with a given reference outward normal
806 // direction provided by the user. if the dot product is negative, flip the tangent vector ordering
807 // and flip the sign on the normal direction.
808 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
809 compadre_kernel_assert_debug(dimensions > 1
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);
814
815 }
816 if (dot_product < 0) {
817 if (dimensions==3) {
818 for (int i=0; i<dimensions; ++i) {
819 // swap tangent directions
820 double tmp = T(0,i);
821 T(0,i) = T(1,i);
822 T(1,i) = tmp;
823 }
824 }
825 for (int i=0; i<dimensions; ++i) {
826 // flip the sign of the normal vector
827 T(dimensions-1,i) *= -1;
828
829 }
830 }
831 });
832 teamMember.team_barrier();
833 }
834};
835
836//! Functor to evaluate curvature targets and apply to coefficients of curvature reconstruction
838
840
842
843 KOKKOS_INLINE_FUNCTION
844 void operator()(const member_type& teamMember) const {
845
846 /*
847 * Dimensions
848 */
849
850 const int target_index = _data._initial_index_for_batch + teamMember.league_rank();
851 const int local_index = teamMember.league_rank();
852 auto dimensions = _data._dimensions;
853
854 /*
855 * Data
856 */
857
859 + TO_GLOBAL(local_index)*TO_GLOBAL(_data.Coeffs_dim_0*_data.Coeffs_dim_1),
860 _data.Coeffs_dim_0, _data.Coeffs_dim_1);
861
863 + TO_GLOBAL(target_index)*TO_GLOBAL(dimensions*dimensions),
864 dimensions, dimensions);
865
867 + target_index*TO_GLOBAL(_data.manifold_NP), _data.manifold_NP);
868
870 + TO_GLOBAL(local_index)*TO_GLOBAL(_data.P_target_row_dim_0*_data.P_target_row_dim_1),
872
873
874 // delta, used for each thread
875 scratch_vector_type delta(teamMember.thread_scratch(_data._pm.getThreadScratchLevel(0)),
876 _data.this_num_cols);
877 scratch_vector_type thread_workspace(teamMember.thread_scratch(_data._pm.getThreadScratchLevel(0)),
879
880 /*
881 * Manifold
882 */
883
884 computeCurvatureFunctionals(_data, teamMember, delta, thread_workspace, P_target_row, &T);
885
886 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
887 for (int j=0; j<_data.manifold_NP; ++j) { // set to zero
888 manifold_coeffs(j) = 0;
889 }
890 });
891 teamMember.team_barrier();
892 for (int i=0; i<_data._pc._nla.getNumberOfNeighborsDevice(target_index); ++i) {
893 XYZ rel_coord = _data._pc.getRelativeCoord(target_index, i, dimensions, &T);
894 double normal_coordinate = rel_coord[dimensions-1];
895
896 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
897 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
898 // coefficients without a target premultiplied
899 for (int j=0; j<_data.manifold_NP; ++j) {
900 manifold_coeffs(j) += Q(j,i) * normal_coordinate;
901 }
902 });
903 });
904 teamMember.team_barrier();
905 }
906 }
907};
908
909//! Functor to assemble the P*sqrt(weights) matrix and construct sqrt(weights)*Identity
911
913
915
916 KOKKOS_INLINE_FUNCTION
917 void operator()(const member_type& teamMember) const {
918
919 /*
920 * Dimensions
921 */
922
923 const int target_index = _data._initial_index_for_batch + teamMember.league_rank();
924 const int local_index = teamMember.league_rank();
925 auto dimensions = _data._dimensions;
926 const int this_num_rows = _data._sampling_multiplier*_data._pc._nla.getNumberOfNeighborsDevice(target_index);
927
928 /*
929 * Data
930 */
931
933 + TO_GLOBAL(local_index)*TO_GLOBAL(_data.P_dim_0)*TO_GLOBAL(_data.P_dim_1),
934 _data.P_dim_0, _data.P_dim_1);
936 + TO_GLOBAL(local_index)*TO_GLOBAL(_data.RHS_dim_0)*TO_GLOBAL(_data.RHS_dim_1),
937 _data.RHS_dim_0, _data.RHS_dim_1);
939 + TO_GLOBAL(local_index)*TO_GLOBAL(_data.max_num_rows), _data.max_num_rows);
941 + TO_GLOBAL(target_index)*TO_GLOBAL(dimensions)*TO_GLOBAL(dimensions), dimensions, dimensions);
942
943 // delta, used for each thread
944 scratch_vector_type delta(teamMember.thread_scratch(_data._pm.getThreadScratchLevel(0)),
945 _data.this_num_cols);
946 scratch_vector_type thread_workspace(teamMember.thread_scratch(_data._pm.getThreadScratchLevel(0)),
948
949 /*
950 * Manifold
951 */
952
953 createWeightsAndP(_data, teamMember, delta, thread_workspace, PsqrtW, w, dimensions-1,
954 _data._poly_order, true /* weight with W*/, &T, _data._reconstruction_space,
956
957 if (_data._dense_solver_type != DenseSolverType::LU) {
958 // fill in RHS with Identity * sqrt(weights)
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));
962 });
963 } else {
964 // create global memory for matrix M = PsqrtW^T*PsqrtW
965 // don't need to cast into scratch_matrix_left_type since the matrix is symmetric
967 + TO_GLOBAL(local_index)*TO_GLOBAL(_data.RHS_dim_0*_data.RHS_dim_1),
968 _data.RHS_dim_0, _data.RHS_dim_1);
969
970 // Assemble matrix M
971 KokkosBatched::TeamVectorGemm<member_type,KokkosBatched::Trans::Transpose,KokkosBatched::Trans::NoTranspose,KokkosBatched::Algo::Gemm::Unblocked>
972 ::invoke(teamMember,
973 1.0,
974 PsqrtW,
975 PsqrtW,
976 0.0,
977 M);
978 teamMember.team_barrier();
979
980 // Multiply PsqrtW with sqrt(W) to get PW
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));
984 });
985 });
986 }
987 teamMember.team_barrier();
988 }
989};
990
991//! Functor to evaluate targets on a manifold
993
995
997
998 KOKKOS_INLINE_FUNCTION
999 void operator()(const member_type& teamMember) const {
1000
1001 /*
1002 * Dimensions
1003 */
1004
1005 const int target_index = _data._initial_index_for_batch + teamMember.league_rank();
1006 const int local_index = teamMember.league_rank();
1007 auto dimensions = _data._dimensions;
1008
1009 /*
1010 * Data
1011 */
1012
1014 + TO_GLOBAL(target_index)*TO_GLOBAL(dimensions*dimensions), dimensions, dimensions);
1016 + target_index*TO_GLOBAL(_data.manifold_NP), _data.manifold_NP);
1018 + TO_GLOBAL(local_index)*TO_GLOBAL(_data.P_target_row_dim_0*_data.P_target_row_dim_1),
1020
1021 // delta, used for each thread
1022 scratch_vector_type delta(teamMember.thread_scratch(_data._pm.getThreadScratchLevel(0)),
1023 _data.this_num_cols);
1024 scratch_vector_type thread_workspace(teamMember.thread_scratch(_data._pm.getThreadScratchLevel(0)),
1025 _data.thread_workspace_dim);
1026
1027 /*
1028 * Apply Standard Target Evaluations to Polynomial Coefficients
1029 */
1030
1031 computeTargetFunctionalsOnManifold(_data, teamMember, delta, thread_workspace, P_target_row, T, manifold_coeffs);
1032
1033 }
1034};
1035
1036///@}
1037
1038} // Compadre
1039
1040#endif
#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
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.
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
Functor to assemble the P*sqrt(weights) matrix and construct sqrt(weights)*Identity.
KOKKOS_INLINE_FUNCTION void operator()(const member_type &teamMember) const
Functor to create a coarse tangent approximation from a given neighborhood of points.
KOKKOS_INLINE_FUNCTION void operator()(const member_type &teamMember) const
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
Functor to evaluate targets on a manifold.
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
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
GMLS::point_connections_type _pc
SamplingFunctional _polynomial_sampling_functional
ReconstructionSpace _reconstruction_space
DenseSolverType _dense_solver_type
SolutionSet< device_memory_space > _d_ss
Kokkos::View< double **, layout_right > _source_extra_data
WeightingFunctionType _curvature_weighting_type
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
Kokkos::View< int * > additional_number_of_neighbors_list
Kokkos::View< int * > number_of_neighbors_list
SolutionSet< device_memory_space > _d_ss
Functor to evaluate curvature targets and construct accurate tangent direction approximation for mani...
KOKKOS_INLINE_FUNCTION void operator()(const member_type &teamMember) const
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.