Compadre  1.2.0
GMLS_Manifold_Multiple_Evaluation_Sites.cpp
Go to the documentation of this file.
1 #include <iostream>
2 #include <string>
3 #include <vector>
4 #include <map>
5 #include <stdlib.h>
6 #include <cstdio>
7 #include <random>
8 
9 #include <Compadre_Config.h>
10 #include <Compadre_GMLS.hpp>
11 #include <Compadre_Evaluator.hpp>
13 
14 #include "GMLS_Manifold.hpp"
15 
16 #ifdef COMPADRE_USE_MPI
17 #include <mpi.h>
18 #endif
19 
20 #include <Kokkos_Timer.hpp>
21 #include <Kokkos_Core.hpp>
22 
23 using namespace Compadre;
24 //! [Ambient to Local Back To Ambient Helper Function]
25 void AmbientLocalAmbient(XYZ& u, double* T_data, double* P_data) {
26  // reconstructions with vector output on a manifold are defined
27  // in their local chart (2D). Next, they are mapped to 3D using
28  // the tangent vector calculated at the target site.
29  //
30  // We are interested in a mapping using the tangent vector at
31  // additional evaluation site instead, so we map back to the
32  // local chart using the tangent vector defined at the target
33  // site (T), then mapping from the local chart to ambient space
34  // using the tangent vector calculated at the extra site (P).
35 
36 
37  host_scratch_matrix_right_type T(T_data, 3, 3);
38  host_scratch_matrix_right_type P(P_data, 3, 3);
39 
40  // first get back to local chart
41  double local_vec[3] = {0,0};
42  for (int j=0; j<3; ++j) {
43  local_vec[0] += T(0,j) * u[j];
44  local_vec[1] += T(1,j) * u[j];
45  }
46  // second go to ambient space using tangent for first neighbor
47  for (int j=0; j<3; ++j) u[j] = 0;
48  for (int j=0; j<3; ++j) {
49  u[j] += P(0, j) * local_vec[0];
50  u[j] += P(1, j) * local_vec[1];
51  }
52 
53 }
54 
55 //! [Ambient to Local Back To Ambient Helper Function]
56 
57 // called from command line
58 int main (int argc, char* args[]) {
59 
60 // initializes MPI (if available) with command line arguments given
61 #ifdef COMPADRE_USE_MPI
62 MPI_Init(&argc, &args);
63 #endif
64 
65 // initializes Kokkos with command line arguments given
66 Kokkos::initialize(argc, args);
67 
68 // code block to reduce scope for all Kokkos View allocations
69 // otherwise, Views may be deallocating when we call Kokkos::finalize() later
70 {
71  // check if 8 arguments are given from the command line, the first being the program name
72  // constraint_type used in solving each GMLS problem:
73  // 0 - No constraints used in solving each GMLS problem
74  // 1 - Neumann Gradient Scalar used in solving each GMLS problem
75  int constraint_type = 0; // No constraints by default
76  if (argc >= 8) {
77  int arg8toi = atoi(args[7]);
78  if (arg8toi > 0) {
79  constraint_type = arg8toi;
80  }
81  }
82 
83  // check if 7 arguments are given from the command line, the first being the program name
84  // problem_type used in solving each GMLS problem:
85  // 0 - Standard GMLS problem
86  // 1 - Manifold GMLS problem
87  int problem_type = 1; // Manifold for this example
88  if (argc >= 7) {
89  int arg7toi = atoi(args[6]);
90  if (arg7toi > 0) {
91  problem_type = arg7toi;
92  }
93  }
94 
95  // check if 6 arguments are given from the command line, the first being the program name
96  // solver_type used for factorization in solving each GMLS problem:
97  // 0 - SVD used for factorization in solving each GMLS problem
98  // 1 - QR used for factorization in solving each GMLS problem
99  // 2 - LU used for factorization in solving each GMLS problem
100  int solver_type = 1; // QR by default
101  if (argc >= 6) {
102  int arg6toi = atoi(args[5]);
103  if (arg6toi >= 0) {
104  solver_type = arg6toi;
105  }
106  }
107 
108  // check if 5 arguments are given from the command line, the first being the program name
109  // N_pts_on_sphere used to determine spatial resolution
110  int N_pts_on_sphere = 1000; // 1000 points by default
111  if (argc >= 5) {
112  int arg5toi = atoi(args[4]);
113  if (arg5toi > 0) {
114  N_pts_on_sphere = arg5toi;
115  }
116  }
117 
118  // check if 4 arguments are given from the command line
119  // dimension for the coordinates and the solution
120  int dimension = 3; // dimension 3 by default
121  if (argc >= 4) {
122  int arg4toi = atoi(args[3]);
123  if (arg4toi > 0) {
124  dimension = arg4toi;
125  }
126  }
127 
128  // check if 3 arguments are given from the command line
129  // set the number of target sites where we will reconstruct the target functionals at
130  int number_target_coords = 200; // 200 target sites by default
131  if (argc >= 3) {
132  int arg3toi = atoi(args[2]);
133  if (arg3toi > 0) {
134  number_target_coords = arg3toi;
135  }
136  }
137 
138  // check if 2 arguments are given from the command line
139  // set the number of target sites where we will reconstruct the target functionals at
140  int order = 3; // 3rd degree polynomial basis by default
141  if (argc >= 2) {
142  int arg2toi = atoi(args[1]);
143  if (arg2toi > 0) {
144  order = arg2toi;
145  }
146  }
147 
148  // minimum neighbors for unisolvency is the same as the size of the polynomial basis
149  // dimension has one subtracted because it is a D-1 manifold represented in D dimensions
150  const int min_neighbors = Compadre::GMLS::getNP(order, dimension-1);
151 
152  //! [Parse Command Line Arguments]
153  Kokkos::Timer timer;
154  Kokkos::Profiling::pushRegion("Setup Point Data");
155  //! [Setting Up The Point Cloud]
156 
157 
158  // coordinates of source sites, bigger than needed then resized later
159  Kokkos::View<double**, Kokkos::DefaultExecutionSpace> source_coords_device("source coordinates",
160  1.25*N_pts_on_sphere, 3);
161  Kokkos::View<double**>::HostMirror source_coords = Kokkos::create_mirror_view(source_coords_device);
162 
163  double r = 1.0;
164 
165  // set number of source coordinates from what was calculated
166  int number_source_coords;
167 
168  {
169  // fill source coordinates from surface of a sphere with quasiuniform points
170  // algorithm described at https://www.cmu.edu/biolphys/deserno/pdf/sphere_equi.pdf
171  int N_count = 0;
172  double a = 4*PI*r*r/N_pts_on_sphere;
173  double d = std::sqrt(a);
174  int M_theta = std::round(PI/d);
175  double d_theta = PI/M_theta;
176  double d_phi = a/d_theta;
177  for (int i=0; i<M_theta; ++i) {
178  double theta = PI*(i + 0.5)/M_theta;
179  int M_phi = std::round(2*PI*std::sin(theta)/d_phi);
180  for (int j=0; j<M_phi; ++j) {
181  double phi = 2*PI*j/M_phi;
182  source_coords(N_count, 0) = theta;
183  source_coords(N_count, 1) = phi;
184  N_count++;
185  }
186  }
187 
188  for (int i=0; i<N_count; ++i) {
189  double theta = source_coords(i,0);
190  double phi = source_coords(i,1);
191  source_coords(i,0) = r*std::sin(theta)*std::cos(phi);
192  source_coords(i,1) = r*std::sin(theta)*std::sin(phi);
193  source_coords(i,2) = r*cos(theta);
194  //printf("%f %f %f\n", source_coords(i,0), source_coords(i,1), source_coords(i,2));
195  }
196  number_source_coords = N_count;
197  }
198 
199  // resize source_coords to the size actually needed
200  Kokkos::resize(source_coords, number_source_coords, 3);
201  Kokkos::resize(source_coords_device, number_source_coords, 3);
202 
203  // coordinates of target sites
204  Kokkos::View<double**, Kokkos::DefaultExecutionSpace> target_coords_device ("target coordinates",
205  number_target_coords, 3);
206  Kokkos::View<double**>::HostMirror target_coords = Kokkos::create_mirror_view(target_coords_device);
207 
208  {
209  bool enough_pts_found = false;
210  // fill target coordinates from surface of a sphere with quasiuniform points
211  // stop when enough points are found
212  int N_count = 0;
213  double a = 4.0*PI*r*r/((double)(5.0*number_target_coords)); // 5.0 is in case number_target_coords is set to something very small (like 1)
214  double d = std::sqrt(a);
215  int M_theta = std::round(PI/d);
216  double d_theta = PI/((double)M_theta);
217  double d_phi = a/d_theta;
218  for (int i=0; i<M_theta; ++i) {
219  double theta = PI*(i + 0.5)/M_theta;
220  int M_phi = std::round(2*PI*std::sin(theta)/d_phi);
221  for (int j=0; j<M_phi; ++j) {
222  double phi = 2*PI*j/M_phi;
223  target_coords(N_count, 0) = theta;
224  target_coords(N_count, 1) = phi;
225  N_count++;
226  if (N_count == number_target_coords) {
227  enough_pts_found = true;
228  break;
229  }
230  }
231  if (enough_pts_found) break;
232  }
233 
234  for (int i=0; i<N_count; ++i) {
235  double theta = target_coords(i,0);
236  double phi = target_coords(i,1);
237  target_coords(i,0) = r*std::sin(theta)*std::cos(phi);
238  target_coords(i,1) = r*std::sin(theta)*std::sin(phi);
239  target_coords(i,2) = r*cos(theta);
240  //printf("%f %f %f\n", target_coords(i,0), target_coords(i,1), target_coords(i,2));
241  }
242  }
243 
244 
245  //! [Setting Up The Point Cloud]
246 
247  Kokkos::Profiling::popRegion();
248  Kokkos::fence();
249  Kokkos::Profiling::pushRegion("Creating Data");
250 
251  //! [Creating The Data]
252 
253 
254  // source coordinates need copied to device before using to construct sampling data
255  Kokkos::deep_copy(source_coords_device, source_coords);
256 
257  // target coordinates copied next, because it is a convenient time to send them to device
258  Kokkos::deep_copy(target_coords_device, target_coords);
259 
260  // ensure that source coordinates are sent to device before evaluating sampling data based on them
261  Kokkos::fence();
262 
263 
264  // need Kokkos View storing true solution (for samples)
265  Kokkos::View<double*, Kokkos::DefaultExecutionSpace> sampling_data_device("samples of true solution",
266  source_coords_device.extent(0));
267  Kokkos::View<double*, Kokkos::DefaultExecutionSpace> ones_data_device("samples of ones",
268  source_coords_device.extent(0));
269  Kokkos::deep_copy(ones_data_device, 1.0);
270 
271  // need Kokkos View storing true vector solution (for samples)
272  Kokkos::View<double**, Kokkos::DefaultExecutionSpace> sampling_vector_data_device("samples of vector true solution",
273  source_coords_device.extent(0), 3);
274 
275  Kokkos::parallel_for("Sampling Manufactured Solutions", Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>
276  (0,source_coords.extent(0)), KOKKOS_LAMBDA(const int i) {
277 
278  // coordinates of source site i
279  double xval = source_coords_device(i,0);
280  double yval = (dimension>1) ? source_coords_device(i,1) : 0;
281  double zval = (dimension>2) ? source_coords_device(i,2) : 0;
282 
283  // data for targets with scalar input
284  sampling_data_device(i) = sphere_harmonic54(xval, yval, zval);
285 
286  for (int j=0; j<3; ++j) {
287  double gradient[3] = {0,0,0};
288  gradient_sphereHarmonic54_ambient(gradient, xval, yval, zval);
289  sampling_vector_data_device(i,j) = gradient[j];
290  }
291  });
292 
293 
294  //! [Creating The Data]
295 
296  Kokkos::Profiling::popRegion();
297  Kokkos::Profiling::pushRegion("Neighbor Search");
298 
299  //! [Performing Neighbor Search]
300 
301 
302  // Point cloud construction for neighbor search
303  // CreatePointCloudSearch constructs an object of type PointCloudSearch, but deduces the templates for you
304  auto point_cloud_search(CreatePointCloudSearch(source_coords, dimension));
305 
306  // each row is a neighbor list for a target site, with the first column of each row containing
307  // the number of neighbors for that rows corresponding target site
308  double epsilon_multiplier = 1.7;
309  int estimated_upper_bound_number_neighbors =
310  point_cloud_search.getEstimatedNumberNeighborsUpperBound(min_neighbors, dimension, epsilon_multiplier);
311 
312  Kokkos::View<int**, Kokkos::DefaultExecutionSpace> neighbor_lists_device("neighbor lists",
313  number_target_coords, estimated_upper_bound_number_neighbors); // first column is # of neighbors
314  Kokkos::View<int**>::HostMirror neighbor_lists = Kokkos::create_mirror_view(neighbor_lists_device);
315 
316  // each target site has a window size
317  Kokkos::View<double*, Kokkos::DefaultExecutionSpace> epsilon_device("h supports", number_target_coords);
318  Kokkos::View<double*>::HostMirror epsilon = Kokkos::create_mirror_view(epsilon_device);
319 
320  // query the point cloud to generate the neighbor lists using a kdtree to produce the n nearest neighbor
321  // to each target site, adding (epsilon_multiplier-1)*100% to whatever the distance away the further neighbor used is from
322  // each target to the view for epsilon
323  point_cloud_search.generate2DNeighborListsFromKNNSearch(false /*not dry run*/, target_coords, neighbor_lists,
324  epsilon, min_neighbors, epsilon_multiplier);
325 
326  //! [Performing Neighbor Search]
327 
328  Kokkos::Profiling::popRegion();
329  Kokkos::fence(); // let call to build neighbor lists complete before copying back to device
330  timer.reset();
331 
332  //! [Setting Up The GMLS Object]
333 
334 
335  // Copy data back to device (they were filled on the host)
336  // We could have filled Kokkos Views with memory space on the host
337  // and used these instead, and then the copying of data to the device
338  // would be performed in the GMLS class
339  Kokkos::deep_copy(neighbor_lists_device, neighbor_lists);
340  Kokkos::deep_copy(epsilon_device, epsilon);
341 
342  // solver name for passing into the GMLS class
343  std::string solver_name;
344  if (solver_type == 0) { // SVD
345  solver_name = "SVD";
346  } else if (solver_type == 1) { // QR
347  solver_name = "QR";
348  } else if (solver_type == 2) { // LU
349  solver_name = "LU";
350  }
351 
352  // problem name for passing into the GMLS class
353  std::string problem_name;
354  if (problem_type == 0) { // Standard
355  problem_name = "STANDARD";
356  } else if (problem_type == 1) { // Manifold
357  problem_name = "MANIFOLD";
358  }
359 
360  // boundary name for passing into the GMLS class
361  std::string constraint_name;
362  if (constraint_type == 0) { // No constraints
363  constraint_name = "NO_CONSTRAINT";
364  } else if (constraint_type == 1) { // Neumann Gradient Scalar
365  constraint_name = "NEUMANN_GRAD_SCALAR";
366  }
367 
368  // initialize an instance of the GMLS class for problems with a scalar basis and
369  // traditional point sampling as the sampling functional
370  GMLS my_GMLS_scalar(order, dimension,
371  solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
372  order /*manifold order*/);
373 
374  // pass in neighbor lists, source coordinates, target coordinates, and window sizes
375  //
376  // neighbor lists have the format:
377  // dimensions: (# number of target sites) X (# maximum number of neighbors for any given target + 1)
378  // the first column contains the number of neighbors for that rows corresponding target index
379  //
380  // source coordinates have the format:
381  // dimensions: (# number of source sites) X (dimension)
382  // entries in the neighbor lists (integers) correspond to rows of this 2D array
383  //
384  // target coordinates have the format:
385  // dimensions: (# number of target sites) X (dimension)
386  // # of target sites is same as # of rows of neighbor lists
387  //
388  my_GMLS_scalar.setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
389 
390  // set up additional sites to evaluate target operators
391  // these sites will be the neighbors of the target site
392  my_GMLS_scalar.setAdditionalEvaluationSitesData(neighbor_lists_device, source_coords_device);
393 
394  // set a reference outward normal direction, causing a surface orientation after
395  // the GMLS instance computes an approximate tangent bundle
396  // on a sphere, the ambient coordinates are the outward normal direction
397  my_GMLS_scalar.setReferenceOutwardNormalDirection(target_coords_device, true /* use to orient surface */);
398 
399  // create a vector of target operations
400  std::vector<TargetOperation> lro_scalar(3);
401  lro_scalar[0] = ScalarPointEvaluation;
402  lro_scalar[1] = GaussianCurvaturePointEvaluation;
403  lro_scalar[2] = CurlOfVectorPointEvaluation;
404 
405  // and then pass them to the GMLS class
406  my_GMLS_scalar.addTargets(lro_scalar);
407 
408  // sets the weighting kernel function from WeightingFunctionType for curvature
409  my_GMLS_scalar.setCurvatureWeightingType(WeightingFunctionType::Power);
410 
411  // power to use in the weighting kernel function for curvature coefficients
412  my_GMLS_scalar.setCurvatureWeightingPower(2);
413 
414  // sets the weighting kernel function from WeightingFunctionType
415  my_GMLS_scalar.setWeightingType(WeightingFunctionType::Power);
416 
417  // power to use in that weighting kernel function
418  my_GMLS_scalar.setWeightingPower(2);
419 
420  // generate the alphas that to be combined with data for each target operation requested in lro
421  my_GMLS_scalar.generateAlphas();
422 
423  Kokkos::Profiling::pushRegion("Full Polynomial Basis GMLS Solution");
424  // initialize another instance of the GMLS class for problems with a vector basis on a manifold and point
425  // evaluation of that vector as the sampling functional
426  // VectorTaylorPolynomial indicates that the basis will be a polynomial with as many components as the
427  // dimension of the manifold. This differs from another possibility, which follows this class.
429  order, dimension,
430  solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
431  order /*manifold order*/);
432 
433  my_GMLS_vector.setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
434  my_GMLS_vector.setAdditionalEvaluationSitesData(neighbor_lists_device, source_coords_device);
435  std::vector<TargetOperation> lro_vector(2);
436  lro_vector[0] = VectorPointEvaluation;
437  //lro_vector[1] = DivergenceOfVectorPointEvaluation;
438  my_GMLS_vector.addTargets(lro_vector);
439  my_GMLS_vector.setCurvatureWeightingType(WeightingFunctionType::Power);
440  my_GMLS_vector.setCurvatureWeightingPower(2);
441  my_GMLS_vector.setWeightingType(WeightingFunctionType::Power);
442  my_GMLS_vector.setWeightingPower(2);
443  my_GMLS_vector.generateAlphas();
444  Kokkos::Profiling::popRegion();
445 
446  Kokkos::Profiling::pushRegion("Scalar Polynomial Basis Repeated to Form a Vector GMLS Solution");
447  // initialize another instance of the GMLS class for problems with a vector basis on a manifold and point
448  // evaluation of that vector as the sampling functional
449  // VectorOfScalarClonesTaylorPolynomial indicates a scalar polynomial will be solved for, since
450  // each component of the reconstructed vector are independent. This results in solving a smaller system
451  // for each GMLS problem, and is the suggested way to do vector reconstructions when sampling functionals
452  // acting on the basis would not create non-zero offdiagonal blocks.
453  //
454  // As an example, consider a 2D manifold in 3D ambient space. The GMLS problem is posed in the local chart,
455  // meaning that the problem being solved looks like
456  //
457  // [P_0 0 | where P_0 has dimension #number of neighbors for a target X #dimension of a scalar basis
458  // | 0 P_1]
459  //
460  // P_1 is similarly shaped, and for sampling functional that is a point evaluation, P_0 and P_1 are
461  // identical and their degrees of freedom in this system are disjoint, allowing us to solve for the
462  // degrees of freedom for either block independently. Additionally, the will produce the exact
463  // same polynomial coefficients for the point sampling functional, therefore it makes sense to use
464  // VectorOfScalarClonesTaylorPolynomial.
465  //
466  // In the print-out for this program, we include the timings and errors on this and VectorTaylorPolynomial
467  // in order to demonstrate that they produce exactly the same answer, but that one is much more efficient.
468  //
470  order, dimension,
471  solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
472  order /*manifold order*/);
473 
474  my_GMLS_vector_of_scalar_clones.setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
475  my_GMLS_vector_of_scalar_clones.setAdditionalEvaluationSitesData(neighbor_lists_device, source_coords_device);
476  std::vector<TargetOperation> lro_vector_of_scalar_clones(2);
477  lro_vector_of_scalar_clones[0] = VectorPointEvaluation;
478  //lro_vector_of_scalar_clones[1] = DivergenceOfVectorPointEvaluation;
479  my_GMLS_vector_of_scalar_clones.addTargets(lro_vector_of_scalar_clones);
480  my_GMLS_vector_of_scalar_clones.setCurvatureWeightingType(WeightingFunctionType::Power);
481  my_GMLS_vector_of_scalar_clones.setCurvatureWeightingPower(2);
482  my_GMLS_vector_of_scalar_clones.setWeightingType(WeightingFunctionType::Power);
483  my_GMLS_vector_of_scalar_clones.setWeightingPower(2);
484  my_GMLS_vector_of_scalar_clones.generateAlphas();
485  Kokkos::Profiling::popRegion();
486 
487 
488  //! [Setting Up The GMLS Object]
489 
490  double instantiation_time = timer.seconds();
491  std::cout << "Took " << instantiation_time << "s to complete alphas generation." << std::endl;
492  Kokkos::fence(); // let generateAlphas finish up before using alphas
493  Kokkos::Profiling::pushRegion("Apply Alphas to Data");
494 
495  //! [Apply GMLS Alphas To Data]
496 
497 
498  // it is important to note that if you expect to use the data as a 1D view, then you should use double*
499  // however, if you know that the target operation will result in a 2D view (vector or matrix output),
500  // then you should template with double** as this is something that can not be infered from the input data
501  // or the target operator at compile time. Additionally, a template argument is required indicating either
502  // Kokkos::HostSpace or Kokkos::DefaultExecutionSpace::memory_space()
503 
504  // The Evaluator class takes care of handling input data views as well as the output data views.
505  // It uses information from the GMLS class to determine how many components are in the input
506  // as well as output for any choice of target functionals and then performs the contactions
507  // on the data using the alpha coefficients generated by the GMLS class, all on the device.
508  Evaluator scalar_gmls_evaluator(&my_GMLS_scalar);
509  Evaluator vector_gmls_evaluator(&my_GMLS_vector);
510  Evaluator vector_gmls_evaluator_of_scalar_clones(&my_GMLS_vector_of_scalar_clones);
511 
512  auto output_value = scalar_gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double*, Kokkos::HostSpace>
513  (sampling_data_device, ScalarPointEvaluation, PointSample,
514  true /*scalar_as_vector_if_needed*/, 1 /*evaluation site index*/);
515 
516  auto output_gaussian_curvature = scalar_gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double*, Kokkos::HostSpace>
517  (ones_data_device, GaussianCurvaturePointEvaluation, PointSample,
518  true /*scalar_as_vector_if_needed*/, 1 /*evaluation site index*/);
519 
520  auto output_curl = scalar_gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double**, Kokkos::HostSpace>
521  (sampling_data_device, CurlOfVectorPointEvaluation, PointSample,
522  true /*scalar_as_vector_if_needed*/, 1 /*evaluation site index*/);
523 
524  //auto output_laplacian = scalar_gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double*, Kokkos::HostSpace>
525  // (sampling_data_device, LaplacianOfScalarPointEvaluation, PointSample,
526  // true /*scalar_as_vector_if_needed*/, 1 /*evaluation site index*/);
527 
528  //auto output_gradient = scalar_gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double**, Kokkos::HostSpace>
529  // (sampling_data_device, GradientOfScalarPointEvaluation, PointSample,
530  // true /*scalar_as_vector_if_needed*/, 1 /*evaluation site index*/);
531 
532  auto output_vector = vector_gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double**, Kokkos::HostSpace>
533  (sampling_vector_data_device, VectorPointEvaluation, VaryingManifoldVectorPointSample,
534  true /*scalar_as_vector_if_needed*/, 1 /*evaluation site index*/);
535 
536  //auto output_divergence = vector_gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double*, Kokkos::HostSpace>
537  // (sampling_vector_data_device, DivergenceOfVectorPointEvaluation, ManifoldVectorPointSample,
538  // true /*scalar_as_vector_if_needed*/, 1 /*evaluation site index*/);
539 
540  auto output_vector_of_scalar_clones =
541  vector_gmls_evaluator_of_scalar_clones.applyAlphasToDataAllComponentsAllTargetSites<double**,
542  Kokkos::HostSpace>(sampling_vector_data_device, VectorPointEvaluation, VaryingManifoldVectorPointSample,
543  true /*scalar_as_vector_if_needed*/, 1 /*evaluation site index*/);
544 
545  //auto output_divergence_of_scalar_clones =
546  // vector_gmls_evaluator_of_scalar_clones.applyAlphasToDataAllComponentsAllTargetSites<double*,
547  // Kokkos::HostSpace>(sampling_vector_data_device, DivergenceOfVectorPointEvaluation, ManifoldVectorPointSample,
548  // true /*scalar_as_vector_if_needed*/, 1 /*evaluation site index*/);
549 
550 
551  // Kokkos::fence(); // let application of alphas to data finish before using results
552  //
553  //// move gradient data to device so that it can be transformed into velocity
554  //auto output_gradient_device_mirror = Kokkos::create_mirror(Kokkos::DefaultExecutionSpace::memory_space(), output_gradient);
555  //Kokkos::deep_copy(output_gradient_device_mirror, output_gradient);
556  //Kokkos::parallel_for("Create Velocity From Surface Gradient", Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>
557  // (0,target_coords.extent(0)), KOKKOS_LAMBDA(const int i) {
558  //
559  // // coordinates of target site i
560  // double xval = target_coords_device(i,0);
561  // double yval = (dimension>1) ? target_coords_device(i,1) : 0;
562  // double zval = (dimension>2) ? target_coords_device(i,2) : 0;
563 
564  // double gradx = output_gradient_device_mirror(i,0);
565  // double grady = output_gradient_device_mirror(i,1);
566  // double gradz = output_gradient_device_mirror(i,2);
567  //
568  // // overwrites gradient with velocity
569  // output_gradient_device_mirror(i,0) = (grady*zval - yval*gradz);
570  // output_gradient_device_mirror(i,1) = (-(gradx*zval - xval*gradz));
571  // output_gradient_device_mirror(i,2) = (gradx*yval - xval*grady);
572  //
573  //});
574  //Kokkos::deep_copy(output_gradient, output_gradient_device_mirror);
575 
576 
577  //! [Apply GMLS Alphas To Data]
578 
579  Kokkos::fence(); // let application of alphas to data finish before using results
580  Kokkos::Profiling::popRegion();
581  // times the Comparison in Kokkos
582  Kokkos::Profiling::pushRegion("Comparison");
583 
584  //! [Check That Solutions Are Correct]
585 
586  double tangent_bundle_error = 0;
587  double tangent_bundle_norm = 0;
588  double values_error = 0;
589  double values_norm = 0;
590  double gc_error = 0;
591  double gc_norm = 0;
592  double curl_ambient_error = 0;
593  double curl_ambient_norm = 0;
594  //double laplacian_error = 0;
595  //double laplacian_norm = 0;
596  //double gradient_ambient_error = 0;
597  //double gradient_ambient_norm = 0;
598  double vector_ambient_error = 0;
599  double vector_ambient_norm = 0;
600  //double divergence_ambient_error = 0;
601  //double divergence_ambient_norm = 0;
602  double vector_of_scalar_clones_ambient_error = 0;
603  double vector_of_scalar_clones_ambient_norm = 0;
604  //double divergence_of_scalar_clones_ambient_error = 0;
605  //double divergence_of_scalar_clones_ambient_norm = 0;
606 
607  // tangent vectors for each source coordinate are stored here
608  auto d_prestencil_weights = my_GMLS_vector_of_scalar_clones.getPrestencilWeights();
609  auto prestencil_weights = Kokkos::create_mirror_view(d_prestencil_weights);
610  Kokkos::deep_copy(prestencil_weights, d_prestencil_weights);
611 
612  // tangent vector at target sites are stored here
613  auto d_tangent_directions = my_GMLS_vector_of_scalar_clones.getTangentDirections();
614  auto tangent_directions = Kokkos::create_mirror_view(d_tangent_directions);
615  Kokkos::deep_copy(tangent_directions, d_tangent_directions);
616 
617  // loop through the target sites
618  for (int i=0; i<number_target_coords; i++) {
619 
621  (tangent_directions.data() + TO_GLOBAL(i)*TO_GLOBAL(3)*TO_GLOBAL(3), 3, 3);
622  XYZ u;
623 
624 
625  // load value from output
626  double GMLS_value = output_value(i);
627  //printf("GMLS val: %f, %d\n", GMLS_value, i);
628 
629  double GMLS_gc = output_gaussian_curvature(i);
630  //printf("GMLS gc: %f, %d\n", GMLS_gc, i);
631 
632  // // load laplacian from output
633  // double GMLS_Laplacian = output_laplacian(i);
634 
635  // target site i's coordinate
636  //double xval = target_coords(i,0);
637  //double yval = (dimension>1) ? target_coords(i,1) : 0;
638  //double zval = (dimension>2) ? target_coords(i,2) : 0;
639  double xval = source_coords(neighbor_lists(i,1),0);
640  double yval = (dimension>1) ? source_coords(neighbor_lists(i,1),1) : 0;
641  double zval = (dimension>2) ? source_coords(neighbor_lists(i,1),2) : 0;
642  double coord[3] = {xval, yval, zval};
643 
644  // get tangent vector and see if orthgonal to coordinate (it should be on a sphere)
645  for (int j=0; j<dimension-1; ++j) {
646  double tangent_inner_prod = 0;
647  for (int k=0; k<dimension; ++k) {
648  tangent_inner_prod += coord[k] * prestencil_weights(0, i, 0 /* local neighbor index */, j, k);
649  }
650  tangent_bundle_error += tangent_inner_prod * tangent_inner_prod;
651  }
652  double normal_inner_prod = 0;
653  for (int k=0; k<dimension; ++k) {
654  normal_inner_prod += coord[k] * my_GMLS_scalar.getTangentBundle(i, dimension-1, k);
655  }
656  // inner product could be plus or minus 1 (depends on tangent direction ordering)
657  double normal_error_to_sum = (normal_inner_prod > 0) ? normal_inner_prod - 1 : normal_inner_prod + 1;
658  tangent_bundle_error += normal_error_to_sum * normal_error_to_sum;
659  tangent_bundle_norm += 1;
660 
661  // evaluation of various exact solutions
662  double actual_value = sphere_harmonic54(xval, yval, zval);
663  double actual_gc = 1.0; // Gaussian curvature is constant
664  // double actual_Laplacian = laplace_beltrami_sphere_harmonic54(xval, yval, zval);
665  double actual_Gradient_ambient[3] = {0,0,0}; // initialized for 3, but only filled up to dimension
666  gradient_sphereHarmonic54_ambient(actual_Gradient_ambient, xval, yval, zval);
667  // //velocity_sphereHarmonic54_ambient(actual_Gradient_ambient, xval, yval, zval);
668  double actual_Curl_ambient[3] = {0,0,0};
669  curl_sphere_harmonic54(actual_Curl_ambient, xval, yval, zval);
670 
671  values_error += (GMLS_value - actual_value)*(GMLS_value - actual_value);
672  values_norm += actual_value*actual_value;
673 
674  gc_error += (GMLS_gc - actual_gc)*(GMLS_gc - actual_gc);
675  gc_norm += 1;
676 
677  for (int j=0; j<dimension; ++j) u[j] = output_curl(i,j);
678  AmbientLocalAmbient(u, T.data(), &prestencil_weights(0, i, 0 /* local neighbor index */, 0, 0));
679  for (int j=0; j<dimension; ++j) output_curl(i,j) = u[j];
680 
681  for (int j=0; j<dimension; ++j) {
682  curl_ambient_error += (output_curl(i,j) - actual_Curl_ambient[j])*(output_curl(i,j) - actual_Curl_ambient[j]);
683  curl_ambient_norm += actual_Curl_ambient[j]*actual_Curl_ambient[j];
684  }
685 
686  // laplacian_error += (GMLS_Laplacian - actual_Laplacian)*(GMLS_Laplacian - actual_Laplacian);
687  // laplacian_norm += actual_Laplacian*actual_Laplacian;
688 
689  // for (int j=0; j<dimension; ++j) {
690  // gradient_ambient_error += (output_gradient(i,j) - actual_Gradient_ambient[j])*(output_gradient(i,j) - actual_Gradient_ambient[j]);
691  // gradient_ambient_norm += actual_Gradient_ambient[j]*actual_Gradient_ambient[j];
692  // }
693 
694  for (int j=0; j<dimension; ++j) u[j] = output_vector(i,j);
695  AmbientLocalAmbient(u, T.data(), &prestencil_weights(0, i, 0 /* local neighbor index */, 0, 0));
696  for (int j=0; j<dimension; ++j) output_vector(i,j) = u[j];
697 
698  for (int j=0; j<dimension; ++j) {
699  vector_ambient_error += (output_vector(i,j) - actual_Gradient_ambient[j])*(output_vector(i,j) - actual_Gradient_ambient[j]);
700  vector_ambient_norm += actual_Gradient_ambient[j]*actual_Gradient_ambient[j];
701  }
702 
703  // divergence_ambient_error += (output_divergence(i) - actual_Laplacian)*(output_divergence(i) - actual_Laplacian);
704  // divergence_ambient_norm += actual_Laplacian*actual_Laplacian;
705 
706 
707  for (int j=0; j<dimension; ++j) u[j] = output_vector_of_scalar_clones(i,j);
708  AmbientLocalAmbient(u, T.data(), &prestencil_weights(0, i, 0 /* local neighbor index */, 0, 0));
709  for (int j=0; j<dimension; ++j) output_vector_of_scalar_clones(i,j) = u[j];
710  //// first get back to local chart
711  //double local_vec[3] = {0,0};
712  //for (int j=0; j<dimension; ++j) {
713  // local_vec[0] += T(0,j) * output_vector_of_scalar_clones(i,j);
714  // local_vec[1] += T(1,j) * output_vector_of_scalar_clones(i,j);
715  //}
716  //// second go to ambient space using tangent for first neighbor
717  //for (int j=0; j<dimension; ++j) output_vector_of_scalar_clones(i,j) = 0;
718  //for (int j=0; j<dimension; ++j) {
719  // output_vector_of_scalar_clones(i,j) += prestencil_weights(0, i, 0 /* local neighbor index */, 0, j) * local_vec[0];
720  // output_vector_of_scalar_clones(i,j) += prestencil_weights(0, i, 0 /* local neighbor index */, 1, j) * local_vec[1];
721  //}
722 
723 
724  for (int j=0; j<dimension; ++j) {
725  vector_of_scalar_clones_ambient_error += (output_vector_of_scalar_clones(i,j) - actual_Gradient_ambient[j])*(output_vector_of_scalar_clones(i,j) - actual_Gradient_ambient[j]);
726  vector_of_scalar_clones_ambient_norm += actual_Gradient_ambient[j]*actual_Gradient_ambient[j];
727  }
728 
729  // divergence_of_scalar_clones_ambient_error += (output_divergence_of_scalar_clones(i) - actual_Laplacian)*(output_divergence_of_scalar_clones(i) - actual_Laplacian);
730  // divergence_of_scalar_clones_ambient_norm += actual_Laplacian*actual_Laplacian;
731 
732  }
733 
734  tangent_bundle_error /= number_target_coords;
735  tangent_bundle_error = std::sqrt(tangent_bundle_error);
736  tangent_bundle_norm /= number_target_coords;
737  tangent_bundle_norm = std::sqrt(tangent_bundle_norm);
738 
739  values_error /= number_target_coords;
740  values_error = std::sqrt(values_error);
741  values_norm /= number_target_coords;
742  values_norm = std::sqrt(values_norm);
743 
744  gc_error /= number_target_coords;
745  gc_error = std::sqrt(gc_error);
746  gc_norm /= number_target_coords;
747  gc_norm = std::sqrt(gc_norm);
748 
749  curl_ambient_error /= number_target_coords;
750  curl_ambient_error = std::sqrt(curl_ambient_error);
751  curl_ambient_norm /= number_target_coords;
752  curl_ambient_norm = std::sqrt(curl_ambient_norm);
753 
754  //laplacian_error /= number_target_coords;
755  //laplacian_error = std::sqrt(laplacian_error);
756  //laplacian_norm /= number_target_coords;
757  //laplacian_norm = std::sqrt(laplacian_norm);
758 
759  //gradient_ambient_error /= number_target_coords;
760  //gradient_ambient_error = std::sqrt(gradient_ambient_error);
761  //gradient_ambient_norm /= number_target_coords;
762  //gradient_ambient_norm = std::sqrt(gradient_ambient_norm);
763 
764  vector_ambient_error /= number_target_coords;
765  vector_ambient_error = std::sqrt(vector_ambient_error);
766  vector_ambient_norm /= number_target_coords;
767  vector_ambient_norm = std::sqrt(vector_ambient_norm);
768 
769  //divergence_ambient_error /= number_target_coords;
770  //divergence_ambient_error = std::sqrt(divergence_ambient_error);
771  //divergence_ambient_norm /= number_target_coords;
772  //divergence_ambient_norm = std::sqrt(divergence_ambient_norm);
773 
774  vector_of_scalar_clones_ambient_error /= number_target_coords;
775  vector_of_scalar_clones_ambient_error = std::sqrt(vector_of_scalar_clones_ambient_error);
776  vector_of_scalar_clones_ambient_norm /= number_target_coords;
777  vector_of_scalar_clones_ambient_norm = std::sqrt(vector_of_scalar_clones_ambient_norm);
778 
779  //divergence_of_scalar_clones_ambient_error /= number_target_coords;
780  //divergence_of_scalar_clones_ambient_error = std::sqrt(divergence_of_scalar_clones_ambient_error);
781  //divergence_of_scalar_clones_ambient_norm /= number_target_coords;
782  //divergence_of_scalar_clones_ambient_norm = std::sqrt(divergence_of_scalar_clones_ambient_norm);
783 
784  printf("Tangent Bundle Error: %g\n", tangent_bundle_error / tangent_bundle_norm);
785  printf("Point Value Error: %g\n", values_error / values_norm);
786  printf("Gaussian Curvature Error: %g\n", gc_error / gc_norm);
787  printf("Surface Curl (Ambient) Error: %g\n", curl_ambient_error / curl_ambient_norm);
788  //printf("Laplace-Beltrami Error: %g\n", laplacian_error / laplacian_norm);
789  //printf("Surface Gradient (Ambient) Error: %g\n", gradient_ambient_error / gradient_ambient_norm);
790  printf("Surface Vector (VectorBasis) Error: %g\n", vector_ambient_error / vector_ambient_norm);
791  //printf("Surface Divergence (VectorBasis) Error: %g\n", divergence_ambient_error / divergence_ambient_norm);
792  printf("Surface Vector (ScalarClones) Error: %g\n",
793  vector_of_scalar_clones_ambient_error / vector_of_scalar_clones_ambient_norm);
794  //printf("Surface Divergence (ScalarClones) Error: %g\n",
795  // divergence_of_scalar_clones_ambient_error / divergence_of_scalar_clones_ambient_norm);
796  //! [Check That Solutions Are Correct]
797  // popRegion hidden from tutorial
798  // stop timing comparison loop
799  Kokkos::Profiling::popRegion();
800  //! [Finalize Program]
801 
802 
803 } // end of code block to reduce scope, causing Kokkos View de-allocations
804 // otherwise, Views may be deallocating when we call Kokkos::finalize() later
805 
806 // finalize Kokkos and MPI (if available)
807 Kokkos::finalize();
808 #ifdef COMPADRE_USE_MPI
809 MPI_Finalize();
810 #endif
811 
812 return 0;
813 
814 } // main
815 
816 
817 //! [Finalize Program]
KOKKOS_INLINE_FUNCTION double sphere_harmonic54(double x, double y, double z)
Point evaluation of the curl of a vector (results in a vector)
Lightweight Evaluator Helper This class is a lightweight wrapper for extracting and applying all rele...
Point evaluation of a vector (reconstructs entire vector at once, requiring a ReconstructionSpace hav...
PointCloudSearch< view_type > CreatePointCloudSearch(view_type src_view, const local_index_type dimensions=-1, const local_index_type max_leaf=-1)
CreatePointCloudSearch allows for the construction of an object of type PointCloudSearch with templat...
constexpr SamplingFunctional VaryingManifoldVectorPointSample
For integrating polynomial dotted with normal over an edge.
#define PI
KOKKOS_INLINE_FUNCTION void curl_sphere_harmonic54(double *curl, double x, double y, double z)
Point evaluation of Gaussian curvature.
void AmbientLocalAmbient(XYZ &u, double *T_data, double *P_data)
[Ambient to Local Back To Ambient Helper Function]
#define TO_GLOBAL(variable)
Point evaluation of a scalar.
Kokkos::View< double **, layout_right, host_execution_space, Kokkos::MemoryTraits< Kokkos::Unmanaged > > host_scratch_matrix_right_type
static KOKKOS_INLINE_FUNCTION int getNP(const int m, const int dimension=3, const ReconstructionSpace r_space=ReconstructionSpace::ScalarTaylorPolynomial)
Returns size of the basis for a given polynomial order and dimension General to dimension 1...
Scalar basis reused as many times as there are components in the vector resulting in a much cheaper p...
Kokkos::View< output_data_type, output_array_layout, output_memory_space > applyAlphasToDataAllComponentsAllTargetSites(view_type_input_data sampling_data, TargetOperation lro, const SamplingFunctional sro_in=PointSample, bool scalar_as_vector_if_needed=true, const int evaluation_site_local_index=0) const
Transformation of data under GMLS.
int main(int argc, char *args[])
[Ambient to Local Back To Ambient Helper Function]
KOKKOS_INLINE_FUNCTION void gradient_sphereHarmonic54_ambient(double *gradient, double x, double y, double z)
Generalized Moving Least Squares (GMLS)
void setProblemData(view_type_1 neighbor_lists, view_type_2 source_coordinates, view_type_3 target_coordinates, view_type_4 epsilons)
Sets basic problem data (neighbor lists, source coordinates, and target coordinates) ...
constexpr SamplingFunctional PointSample
Available sampling functionals.
Vector polynomial basis having # of components _dimensions, or (_dimensions-1) in the case of manifol...