Compadre  1.3.3
GMLS_DivergenceFree.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_Tutorial.hpp"
15 #include "CommandLineProcessor.hpp"
16 
17 #ifdef COMPADRE_USE_MPI
18 #include <mpi.h>
19 #endif
20 
21 #include <Kokkos_Timer.hpp>
22 #include <Kokkos_Core.hpp>
23 
24 using namespace Compadre;
25 
26 //! [Parse Command Line Arguments]
27 
28 // called from command line
29 int main (int argc, char* args[]) {
30 
31 // initializes MPI (if available) with command line arguments given
32 #ifdef COMPADRE_USE_MPI
33 MPI_Init(&argc, &args);
34 #endif
35 
36 // initializes Kokkos with command line arguments given
37 Kokkos::initialize(argc, args);
38 
39 // becomes false if the computed solution not within the failure_threshold of the actual solution
40 bool all_passed = true;
41 
42 // code block to reduce scope for all Kokkos View allocations
43 // otherwise, Views may be deallocating when we call Kokkos::finalize() later
44 {
45 
46  CommandLineProcessor clp(argc, args);
47  auto order = clp.order;
48  auto dimension = clp.dimension;
49  auto number_target_coords = clp.number_target_coords;
50  auto constraint_name = clp.constraint_name;
51  auto solver_name = clp.solver_name;
52  auto problem_name = clp.problem_name;
53 
54  // the functions we will be seeking to reconstruct are in the span of the basis
55  // of the reconstruction space we choose for GMLS, so the error should be very small
56  const double failure_tolerance = 1e-9;
57 
58  // minimum neighbors for unisolvency is the same as the size of the polynomial basis
59  const int min_neighbors = Compadre::GMLS::getNP(order, dimension, DivergenceFreeVectorTaylorPolynomial);
60 
61  //! [Parse Command Line Arguments]
62  Kokkos::Timer timer;
63  Kokkos::Profiling::pushRegion("Setup Point Data");
64  //! [Setting Up The Point Cloud]
65 
66  // approximate spacing of source sites
67  double h_spacing = 0.05;
68  int n_neg1_to_1 = 2*(1/h_spacing) + 1; // always odd
69 
70  // number of source coordinate sites that will fill a box of [-1,1]x[-1,1]x[-1,1] with a spacing approximately h
71  const int number_source_coords = std::pow(n_neg1_to_1, dimension);
72 
73  // coordinates of source sites
74  Kokkos::View<double**, Kokkos::DefaultExecutionSpace> source_coords_device("source coordinates",
75  number_source_coords, 3);
76  Kokkos::View<double**>::HostMirror source_coords = Kokkos::create_mirror_view(source_coords_device);
77 
78  // coordinates of target sites
79  Kokkos::View<double**, Kokkos::DefaultExecutionSpace> target_coords_device ("target coordinates", number_target_coords, 3);
80  Kokkos::View<double**>::HostMirror target_coords = Kokkos::create_mirror_view(target_coords_device);
81 
82 
83  // fill source coordinates with a uniform grid
84  int source_index = 0;
85  double this_coord[3] = {0,0,0};
86  for (int i=-n_neg1_to_1/2; i<n_neg1_to_1/2+1; ++i) {
87  this_coord[0] = i*h_spacing;
88  for (int j=-n_neg1_to_1/2; j<n_neg1_to_1/2+1; ++j) {
89  this_coord[1] = j*h_spacing;
90  for (int k=-n_neg1_to_1/2; k<n_neg1_to_1/2+1; ++k) {
91  this_coord[2] = k*h_spacing;
92  if (dimension==3) {
93  source_coords(source_index,0) = this_coord[0];
94  source_coords(source_index,1) = this_coord[1];
95  source_coords(source_index,2) = this_coord[2];
96  source_index++;
97  }
98  }
99  if (dimension==2) {
100  source_coords(source_index,0) = this_coord[0];
101  source_coords(source_index,1) = this_coord[1];
102  source_coords(source_index,2) = 0;
103  source_index++;
104  }
105  }
106  if (dimension==1) {
107  source_coords(source_index,0) = this_coord[0];
108  source_coords(source_index,1) = 0;
109  source_coords(source_index,2) = 0;
110  source_index++;
111  }
112  }
113 
114  // Generate target points - these are random permutation from available source points
115  // Note that this is assuming that the number of targets in this test will not exceed
116  // the number of source coords, which is 41^3 = 68921
117  // seed random number generator
118  std::mt19937 rng(50);
119  // generate random integers in [0..number_source_coords-1] (used to pick target sites)
120  std::uniform_int_distribution<int> gen_num_neighbours(0, source_index);
121  // fill target sites with random selections from source sites
122  for (int i=0; i<number_target_coords; ++i) {
123  const int source_site_to_copy = gen_num_neighbours(rng);
124  for (int j=0; j<3; j++) {
125  target_coords(i, j) = source_coords(source_site_to_copy, j);
126  }
127  }
128 
129  //! [Setting Up The Point Cloud]
130 
131  Kokkos::Profiling::popRegion();
132  Kokkos::Profiling::pushRegion("Creating Data");
133 
134  //! [Creating The Data]
135 
136 
137  // source coordinates need copied to device before using to construct sampling data
138  Kokkos::deep_copy(source_coords_device, source_coords);
139 
140  // target coordinates copied next, because it is a convenient time to send them to device
141  Kokkos::deep_copy(target_coords_device, target_coords);
142 
143  // need Kokkos View storing true solution
144  Kokkos::View<double**, Kokkos::DefaultExecutionSpace> vector_sampling_span_basis_data_device("samples of true vector solution",
145  source_coords_device.extent(0), dimension);
146  Kokkos::View<double**, Kokkos::DefaultExecutionSpace> vector_sampling_single_polynomial_data_device("samples of true vector solution",
147  source_coords_device.extent(0), dimension);
148 
149  Kokkos::parallel_for("Sampling Manufactured Solutions", Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>
150  (0,source_coords.extent(0)), KOKKOS_LAMBDA(const int i) {
151  // coordinates of source site i
152  double xval = source_coords_device(i,0);
153  double yval = (dimension>1) ? source_coords_device(i,1) : 0;
154  double zval = (dimension>2) ? source_coords_device(i,2) : 0;
155 
156  // data for targets with vector input
157  for (int j=0; j<dimension; ++j) {
158  vector_sampling_span_basis_data_device(i, j) = divfreeTestSolution_span_basis(xval, yval, zval, j, dimension, order);
159  vector_sampling_single_polynomial_data_device(i, j) = divfreeTestSolution_single_polynomial(xval, yval, zval, j, dimension);
160  }
161  });
162 
163 
164  //! [Creating The Data]
165 
166  Kokkos::Profiling::popRegion();
167  Kokkos::Profiling::pushRegion("Neighbor Search");
168 
169  //! [Performing Neighbor Search]
170 
171 
172  // Point cloud construction for neighbor search
173  // CreatePointCloudSearch constructs an object of type PointCloudSearch, but deduces the templates for you
174  auto point_cloud_search(CreatePointCloudSearch(source_coords, dimension));
175 
176  // each row is a neighbor list for a target site, with the first column of each row containing
177  // the number of neighbors for that rows corresponding target site
178  // for the default values in this test, the multiplier is suggested to be 2.2
179  double epsilon_multiplier = 2.2;
180  int estimated_upper_bound_number_neighbors =
181  point_cloud_search.getEstimatedNumberNeighborsUpperBound(min_neighbors, dimension, epsilon_multiplier);
182 
183  Kokkos::View<int**, Kokkos::DefaultExecutionSpace> neighbor_lists_device("neighbor lists",
184  number_target_coords, estimated_upper_bound_number_neighbors); // first column is # of neighbors
185  Kokkos::View<int**>::HostMirror neighbor_lists = Kokkos::create_mirror_view(neighbor_lists_device);
186 
187  // each target site has a window size
188  Kokkos::View<double*, Kokkos::DefaultExecutionSpace> epsilon_device("h supports", number_target_coords);
189  Kokkos::View<double*>::HostMirror epsilon = Kokkos::create_mirror_view(epsilon_device);
190 
191  // query the point cloud to generate the neighbor lists using a kdtree to produce the n nearest neighbor
192  // to each target site, adding (epsilon_multiplier-1)*100% to whatever the distance away the further neighbor used is from
193  // each target to the view for epsilon
194  point_cloud_search.generate2DNeighborListsFromKNNSearch(false /*not dry run*/, target_coords, neighbor_lists,
195  epsilon, min_neighbors, epsilon_multiplier);
196 
197  //! [Performing Neighbor Search]
198 
199  Kokkos::Profiling::popRegion();
200  Kokkos::fence(); // let call to build neighbor lists complete before copying back to device
201  timer.reset();
202 
203  //! [Setting Up The GMLS Object]
204 
205 
206  // Copy data back to device (they were filled on the host)
207  // We could have filled Kokkos Views with memory space on the host
208  // and used these instead, and then the copying of data to the device
209  // would be performed in the GMLS class
210  Kokkos::deep_copy(neighbor_lists_device, neighbor_lists);
211  Kokkos::deep_copy(epsilon_device, epsilon);
212 
213  // initialize an instance of the GMLS class
214  // NULL in manifold order indicates non-manifold case
215  // Vector basis to perform GMLS on divergence free basis
216  GMLS vector_divfree_basis_gmls(DivergenceFreeVectorTaylorPolynomial,
218  order, dimension,
219  solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
220  0 /*manifold order*/);
221 
222  // pass in neighbor lists, source coordinates, target coordinates, and window sizes
223  //
224  // neighbor lists have the format:
225  // dimensions: (# number of target sites) X (# maximum number of neighbors for any given target + 1)
226  // the first column contains the number of neighbors for that rows corresponding target index
227  //
228  // source coordinates have the format:
229  // dimensions: (# number of source sites) X (dimension)
230  // entries in the neighbor lists (integers) correspond to rows of this 2D array
231  //
232  // target coordinates have the format:
233  // dimensions: (# number of target sites) X (dimension)
234  // # of target sites is same as # of rows of neighbor lists
235  //
236  vector_divfree_basis_gmls.setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
237 
238  // create a vector of target operations
239  std::vector<TargetOperation> lro(4);
240  lro[0] = VectorPointEvaluation;
244 
245  // and then pass them to the GMLS class
246  vector_divfree_basis_gmls.addTargets(lro);
247 
248  // sets the weighting kernel function from WeightingFunctionType
249  vector_divfree_basis_gmls.setWeightingType(WeightingFunctionType::Power);
250 
251  // power to use in that weighting kernel function
252  vector_divfree_basis_gmls.setWeightingPower(2);
253 
254  // generate the alphas that to be combined with data for each target operation requested in lro
255  vector_divfree_basis_gmls.generateAlphas(15 /* # batches */);
256 
257  //! [Setting Up The GMLS Object]
258 
259  double instantiation_time = timer.seconds();
260  std::cout << "Took " << instantiation_time << "s to complete alphas generation." << std::endl;
261  Kokkos::fence(); // let generateAlphas finish up before using alphas
262  Kokkos::Profiling::pushRegion("Apply Alphas to Data");
263 
264  //! [Apply GMLS Alphas To Data]
265 
266  // it is important to note that if you expect to use the data as a 1D view, then you should use double*
267  // however, if you know that the target operation will result in a 2D view (vector or matrix output),
268  // then you should template with double** as this is something that can not be infered from the input data
269  // or the target operator at compile time. Additionally, a template argument is required indicating either
270  // Kokkos::HostSpace or Kokkos::DefaultExecutionSpace::memory_space()
271 
272  // The Evaluator class takes care of handling input data views as well as the output data views.
273  // It uses information from the GMLS class to determine how many components are in the input
274  // as well as output for any choice of target functionals and then performs the contactions
275  // on the data using the alpha coefficients generated by the GMLS class, all on the device.
276  Evaluator gmls_evaluator_vector(&vector_divfree_basis_gmls);
277 
278  auto output_vector_evaluation = gmls_evaluator_vector.applyAlphasToDataAllComponentsAllTargetSites<double**, Kokkos::HostSpace>
279  (vector_sampling_span_basis_data_device, VectorPointEvaluation, VectorPointSample);
280  auto output_curl_vector_evaluation = gmls_evaluator_vector.applyAlphasToDataAllComponentsAllTargetSites<double**, Kokkos::HostSpace>
281  (vector_sampling_span_basis_data_device, CurlOfVectorPointEvaluation, VectorPointSample);
282  auto output_curlcurl_vector_evaluation = gmls_evaluator_vector.applyAlphasToDataAllComponentsAllTargetSites<double**, Kokkos::HostSpace>
283  (vector_sampling_single_polynomial_data_device, CurlCurlOfVectorPointEvaluation, VectorPointSample);
284  auto output_gradient_vector_evaluation = gmls_evaluator_vector.applyAlphasToDataAllComponentsAllTargetSites<double**, Kokkos::HostSpace>
285  (vector_sampling_span_basis_data_device, GradientOfVectorPointEvaluation, VectorPointSample);
286 
287  //! [Apply GMLS Alphas To Data]
288 
289  Kokkos::fence(); // let application of alphas to data finish before using results
290  Kokkos::Profiling::popRegion();
291  // times the Comparison in Kokkos
292  Kokkos::Profiling::pushRegion("Comparison");
293 
294  //! [Check That Solutions Are Correct]
295 
296  // loop through the target sites
297  for (int i=0; i<number_target_coords; i++) {
298  // load vector components from output
299  double GMLS_DivFree_VectorX = output_vector_evaluation(i, 0);
300  double GMLS_DivFree_VectorY = (dimension>1) ? output_vector_evaluation(i, 1) : 0;
301  double GMLS_DivFree_VectorZ = (dimension>2) ? output_vector_evaluation(i, 2) : 0;
302 
303  // load curl of vector components from output
304  double GMLS_Curl_DivFree_VectorX = output_curl_vector_evaluation(i, 0);
305  double GMLS_Curl_DivFree_VectorY = (dimension>2) ? output_curl_vector_evaluation(i, 1) : 0;
306  double GMLS_Curl_DivFree_VectorZ = (dimension>2) ? output_curl_vector_evaluation(i, 2) : 0;
307 
308  // load curl curl of vector components from output
309  double GMLS_CurlCurl_DivFree_VectorX = output_curlcurl_vector_evaluation(i, 0);
310  double GMLS_CurlCurl_DivFree_VectorY = (dimension>1) ? output_curlcurl_vector_evaluation(i, 1) : 0;
311  double GMLS_CurlCurl_DivFree_VectorZ = (dimension>2) ? output_curlcurl_vector_evaluation(i, 2) : 0;
312 
313  // load gradient of vector components from output
314  double GMLS_Gradient_DivFree_VectorXX = output_gradient_vector_evaluation(i, 0);
315  double GMLS_Gradient_DivFree_VectorXY = output_gradient_vector_evaluation(i, 1);
316  double GMLS_Gradient_DivFree_VectorXZ = (dimension==3) ? output_gradient_vector_evaluation(i, 2) : 0.0;
317  double GMLS_Gradient_DivFree_VectorYX = (dimension==3) ? output_gradient_vector_evaluation(i, 3) : output_gradient_vector_evaluation(i, 2);
318  double GMLS_Gradient_DivFree_VectorYY = (dimension==3) ? output_gradient_vector_evaluation(i, 4) : output_gradient_vector_evaluation(i, 3);
319  double GMLS_Gradient_DivFree_VectorYZ = (dimension==3) ? output_gradient_vector_evaluation(i, 5) : 0.0;
320  double GMLS_Gradient_DivFree_VectorZX = (dimension==3) ? output_gradient_vector_evaluation(i, 6) : 0.0;
321  double GMLS_Gradient_DivFree_VectorZY = (dimension==3) ? output_gradient_vector_evaluation(i, 7) : 0.0;
322  double GMLS_Gradient_DivFree_VectorZZ = (dimension==3) ? output_gradient_vector_evaluation(i, 8) : 0.0;
323 
324  // target site i's coordinate
325  double xval = target_coords(i,0);
326  double yval = (dimension>1) ? target_coords(i,1) : 0;
327  double zval = (dimension>2) ? target_coords(i,2) : 0;
328 
329  // evaluation of vector exact solutions
330  double actual_vector[3] = {0, 0, 0};
331  if (dimension>=1) {
332  actual_vector[0] = divfreeTestSolution_span_basis(xval, yval, zval, 0, dimension, order);
333  if (dimension>=2) {
334  actual_vector[1] = divfreeTestSolution_span_basis(xval, yval, zval, 1, dimension, order);
335  if (dimension==3) {
336  actual_vector[2] = divfreeTestSolution_span_basis(xval, yval, zval, 2, dimension, order);
337  }
338  }
339  }
340 
341  // evaluation of curl of vector exact solutions
342  double actual_curl_vector[3] = {0, 0, 0};
343  if (dimension>=1) {
344  actual_curl_vector[0] = curldivfreeTestSolution_span_basis(xval, yval, zval, 0, dimension, order);
345  if (dimension==3) {
346  actual_curl_vector[1] = curldivfreeTestSolution_span_basis(xval, yval, zval, 1, dimension, order);
347  actual_curl_vector[2] = curldivfreeTestSolution_span_basis(xval, yval, zval, 2, dimension, order);
348  }
349  }
350 
351  // evaluation of curl of curl of vector exact solutions
352  double actual_curlcurl_vector[3] = {0, 0, 0};
353  if (dimension>=1) {
354  actual_curlcurl_vector[0] = curlcurldivfreeTestSolution_single_polynomial(xval, yval, zval, 0, dimension);
355  if (dimension>=2) {
356  actual_curlcurl_vector[1] = curlcurldivfreeTestSolution_single_polynomial(xval, yval, zval, 1, dimension);
357  if (dimension==3) {
358  actual_curlcurl_vector[2] = curlcurldivfreeTestSolution_single_polynomial(xval, yval, zval, 2, dimension);
359  }
360  }
361  }
362 
363  double actual_gradient_vector[9] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
364  if (dimension==3) {
365  for (int axes = 0; axes < 9; ++axes)
366  actual_gradient_vector[axes] = gradientdivfreeTestSolution_span_basis(xval, yval, zval, axes/dimension, axes%dimension, dimension, order);
367  }
368  if (dimension==2) {
369  for (int axes = 0; axes < 4; ++axes)
370  actual_gradient_vector[axes] = gradientdivfreeTestSolution_span_basis(xval, yval, zval, axes/dimension, axes%dimension, dimension, order);
371  }
372 
373  // check vector evaluation
374  if(std::abs(actual_vector[0] - GMLS_DivFree_VectorX) > failure_tolerance) {
375  all_passed = false;
376  std::cout << i << " Failed VectorX by: " << std::abs(actual_vector[0] - GMLS_DivFree_VectorX) << std::endl;
377  if (dimension>1) {
378  if(std::abs(actual_vector[1] - GMLS_DivFree_VectorY) > failure_tolerance) {
379  all_passed = false;
380  std::cout << i << " Failed VectorY by: " << std::abs(actual_vector[1] - GMLS_DivFree_VectorY) << std::endl;
381  }
382  }
383  if (dimension>2) {
384  if(std::abs(actual_vector[2] - GMLS_DivFree_VectorZ) > failure_tolerance) {
385  all_passed = false;
386  std::cout << i << " Failed VectorZ by: " << std::abs(actual_vector[2] - GMLS_DivFree_VectorZ) << std::endl;
387  }
388  }
389  }
390 
391  // check curl of vector evaluation
392  if (dimension==2) {
393  if(std::abs(actual_curl_vector[0] - GMLS_Curl_DivFree_VectorX) > failure_tolerance) {
394  all_passed = false;
395  std::cout << i << " Failed curl VectorX by: " << std::abs(actual_curl_vector[2] - GMLS_Curl_DivFree_VectorX) << std::endl;
396  }
397  } else if (dimension>2) {
398  if(std::abs(actual_curl_vector[0] - GMLS_Curl_DivFree_VectorX) > failure_tolerance) {
399  all_passed = false;
400  std::cout << i << " Failed curl VectorX by: " << std::abs(actual_curl_vector[0] - GMLS_Curl_DivFree_VectorX) << std::endl;
401  }
402  if(std::abs(actual_curl_vector[1] - GMLS_Curl_DivFree_VectorY) > failure_tolerance) {
403  all_passed = false;
404  std::cout << i << " Failed curl VectorY by: " << std::abs(actual_curl_vector[1] - GMLS_Curl_DivFree_VectorY) << std::endl;
405  }
406  if(std::abs(actual_curl_vector[2] - GMLS_Curl_DivFree_VectorZ) > failure_tolerance) {
407  all_passed = false;
408  std::cout << i << " Failed curl VectorZ by: " << std::abs(actual_curl_vector[2] - GMLS_Curl_DivFree_VectorZ) << std::endl;
409  }
410  }
411 
412  // check curlcurl curlcurl of vector evaluation
413  if(std::abs(actual_curlcurl_vector[0] - GMLS_CurlCurl_DivFree_VectorX) > failure_tolerance) {
414  all_passed = false;
415  std::cout << i << " Failed curl curl VectorX by: " << std::abs(actual_curlcurl_vector[0] - GMLS_CurlCurl_DivFree_VectorX) << std::endl;
416  }
417  if (dimension>1) {
418  if(std::abs(actual_curlcurl_vector[1] - GMLS_CurlCurl_DivFree_VectorY) > failure_tolerance) {
419  all_passed = false;
420  std::cout << i << " Failed curl curl VectorY by: " << std::abs(actual_curlcurl_vector[1] - GMLS_CurlCurl_DivFree_VectorY) << std::endl;
421  }
422  }
423  if (dimension>2) {
424  if(std::abs(actual_curlcurl_vector[2] - GMLS_CurlCurl_DivFree_VectorZ) > failure_tolerance) {
425  all_passed = false;
426  std::cout << i << " Failed curl curl VectorZ by: " << std::abs(actual_curlcurl_vector[2] - GMLS_CurlCurl_DivFree_VectorZ) << std::endl;
427  }
428  }
429 
430  // check gradient of vector evaluation
431  if (dimension==3) {
432  if (std::abs(actual_gradient_vector[0] - GMLS_Gradient_DivFree_VectorXX) > failure_tolerance) {
433  all_passed = false;
434  std::cout << i << " Failed gradient_x VectorX by: " << std::abs(actual_gradient_vector[0] - GMLS_Gradient_DivFree_VectorXX) << std::endl;
435  }
436  if (std::abs(actual_gradient_vector[1] - GMLS_Gradient_DivFree_VectorXY) > failure_tolerance) {
437  all_passed = false;
438  std::cout << i << " Failed gradient_y VectorX by: " << std::abs(actual_gradient_vector[1] - GMLS_Gradient_DivFree_VectorXY) << std::endl;
439  }
440  if (std::abs(actual_gradient_vector[2] - GMLS_Gradient_DivFree_VectorXZ) > failure_tolerance) {
441  all_passed = false;
442  std::cout << i << " Failed gradient_z VectorX by: " << std::abs(actual_gradient_vector[2] - GMLS_Gradient_DivFree_VectorXZ) << std::endl;
443  }
444 
445  if (std::abs(actual_gradient_vector[3] - GMLS_Gradient_DivFree_VectorYX) > failure_tolerance) {
446  all_passed = false;
447  std::cout << i << " Failed gradient_x VectorY by: " << std::abs(actual_gradient_vector[3] - GMLS_Gradient_DivFree_VectorYX) << std::endl;
448  }
449  if (std::abs(actual_gradient_vector[4] - GMLS_Gradient_DivFree_VectorYY) > failure_tolerance) {
450  all_passed = false;
451  std::cout << i << " Failed gradient_y VectorY by: " << std::abs(actual_gradient_vector[4] - GMLS_Gradient_DivFree_VectorYY) << std::endl;
452  }
453  if (std::abs(actual_gradient_vector[5] - GMLS_Gradient_DivFree_VectorYZ) > failure_tolerance) {
454  all_passed = false;
455  std::cout << i << " Failed gradient_z VectorY by: " << std::abs(actual_gradient_vector[5] - GMLS_Gradient_DivFree_VectorYZ) << std::endl;
456  }
457 
458  if (std::abs(actual_gradient_vector[6] - GMLS_Gradient_DivFree_VectorZX) > failure_tolerance) {
459  all_passed = false;
460  std::cout << i << " Failed gradient_x VectorZ by: " << std::abs(actual_gradient_vector[6] - GMLS_Gradient_DivFree_VectorZX) << std::endl;
461  }
462  if (std::abs(actual_gradient_vector[7] - GMLS_Gradient_DivFree_VectorZY) > failure_tolerance) {
463  all_passed = false;
464  std::cout << i << " Failed gradient_y VectorZ by: " << std::abs(actual_gradient_vector[7] - GMLS_Gradient_DivFree_VectorZY) << std::endl;
465  }
466  if (std::abs(actual_gradient_vector[8] - GMLS_Gradient_DivFree_VectorZZ) > failure_tolerance) {
467  all_passed = false;
468  std::cout << i << " Failed gradient_z VectorZ by: " << std::abs(actual_gradient_vector[8] - GMLS_Gradient_DivFree_VectorZZ) << std::endl;
469  }
470  }
471 
472  if (dimension==2) {
473  if (std::abs(actual_gradient_vector[0] - GMLS_Gradient_DivFree_VectorXX) > failure_tolerance) {
474  all_passed = false;
475  std::cout << i << " Failed gradient_x VectorX by: " << std::abs(actual_gradient_vector[0] - GMLS_Gradient_DivFree_VectorXX) << std::endl;
476  }
477  if (std::abs(actual_gradient_vector[1] - GMLS_Gradient_DivFree_VectorXY) > failure_tolerance) {
478  all_passed = false;
479  std::cout << i << " Failed gradient_y VectorX by: " << std::abs(actual_gradient_vector[1] - GMLS_Gradient_DivFree_VectorXY) << std::endl;
480  }
481 
482  if (std::abs(actual_gradient_vector[2] - GMLS_Gradient_DivFree_VectorYX) > failure_tolerance) {
483  all_passed = false;
484  std::cout << i << " Failed gradient_x VectorY by: " << std::abs(actual_gradient_vector[2] - GMLS_Gradient_DivFree_VectorYX) << std::endl;
485  }
486  if (std::abs(actual_gradient_vector[3] - GMLS_Gradient_DivFree_VectorYY) > failure_tolerance) {
487  all_passed = false;
488  std::cout << i << " Failed gradient_y VectorY by: " << (std::abs(actual_gradient_vector[3] - GMLS_Gradient_DivFree_VectorYY) > failure_tolerance) << std::endl;
489  }
490  }
491  }
492 
493  //! [Check That Solutions Are Correct]
494  // popRegion hidden from tutorial
495  // stop timing comparison loop
496  Kokkos::Profiling::popRegion();
497  //! [Finalize Program]
498 
499 
500 } // end of code block to reduce scope, causing Kokkos View de-allocations
501 // otherwise, Views may be deallocating when we call Kokkos::finalize() later
502 
503 // finalize Kokkos and MPI (if available)
504 Kokkos::finalize();
505 #ifdef COMPADRE_USE_MPI
506 MPI_Finalize();
507 #endif
508 
509 // output to user that test passed or failed
510 if(all_passed) {
511  fprintf(stdout, "Passed test \n");
512  return 0;
513 } else {
514  fprintf(stdout, "Failed test \n");
515  return -1;
516 }
517 
518 } // main
519 
520 
521 //! [Finalize Program]
Divergence-free vector polynomial basis.
KOKKOS_INLINE_FUNCTION double divfreeTestSolution_span_basis(double x, double y, double z, int component, int dimension, int exact_order)
Point evaluation of the curl of a curl of a vector (results in a vector)
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...
KOKKOS_INLINE_FUNCTION double curldivfreeTestSolution_span_basis(double x, double y, double z, int component, int dimension, int exact_order)
int main(int argc, char *args[])
[Parse Command Line Arguments]
KOKKOS_INLINE_FUNCTION double curlcurldivfreeTestSolution_single_polynomial(double x, double y, double z, int component, int dimension)
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...
KOKKOS_INLINE_FUNCTION double divfreeTestSolution_single_polynomial(double x, double y, double z, int component, int dimension)
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 (allocates memory for output)
Generalized Moving Least Squares (GMLS)
KOKKOS_INLINE_FUNCTION double gradientdivfreeTestSolution_span_basis(double x, double y, double z, int input_component, int output_component, int dimension, int exact_order)
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) ...
Point evaluation of the gradient of a vector (results in a matrix, NOT CURRENTLY IMPLEMENTED) ...
constexpr SamplingFunctional VectorPointSample
Point evaluations of the entire vector source function.