Compadre  1.3.3
GMLS_Device.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>
14 
15 #include "GMLS_Tutorial.hpp"
16 #include "CommandLineProcessor.hpp"
17 
18 #ifdef COMPADRE_USE_MPI
19 #include <mpi.h>
20 #endif
21 
22 #include <Kokkos_Timer.hpp>
23 #include <Kokkos_Core.hpp>
24 
25 using namespace Compadre;
26 
27 //! [Parse Command Line Arguments]
28 
29 // called from command line
30 int main (int argc, char* args[]) {
31 
32 // initializes MPI (if available) with command line arguments given
33 #ifdef COMPADRE_USE_MPI
34 MPI_Init(&argc, &args);
35 #endif
36 
37 // initializes Kokkos with command line arguments given
38 auto kp = KokkosParser(argc, args, true);
39 
40 // becomes false if the computed solution not within the failure_threshold of the actual solution
41 bool all_passed = true;
42 
43 // code block to reduce scope for all Kokkos View allocations
44 // otherwise, Views may be deallocating when we call Kokkos finalize() later
45 {
46 
47  CommandLineProcessor clp(argc, args);
48  auto order = clp.order;
49  auto dimension = clp.dimension;
50  auto number_target_coords = clp.number_target_coords;
51  auto constraint_name = clp.constraint_name;
52  auto solver_name = clp.solver_name;
53  auto problem_name = clp.problem_name;
54  auto number_of_batches = clp.number_of_batches;
55  bool keep_coefficients = number_of_batches==1;
56 
57  // the functions we will be seeking to reconstruct are in the span of the basis
58  // of the reconstruction space we choose for GMLS, so the error should be very small
59  const double failure_tolerance = 1e-9;
60 
61  // Laplacian is a second order differential operator, which we expect to be slightly less accurate
62  const double laplacian_failure_tolerance = 1e-9;
63 
64  // minimum neighbors for unisolvency is the same as the size of the polynomial basis
65  const int min_neighbors = Compadre::GMLS::getNP(order, dimension);
66 
67  //! [Parse Command Line Arguments]
68  Kokkos::Timer timer;
69  Kokkos::Profiling::pushRegion("Setup Point Data");
70  //! [Setting Up The Point Cloud]
71 
72  // approximate spacing of source sites
73  double h_spacing = 0.05;
74  int n_neg1_to_1 = 2*(1/h_spacing) + 1; // always odd
75 
76  // number of source coordinate sites that will fill a box of [-1,1]x[-1,1]x[-1,1] with a spacing approximately h
77  const int number_source_coords = std::pow(n_neg1_to_1, dimension);
78 
79  // coordinates of source sites
80  Kokkos::View<double**, Kokkos::DefaultExecutionSpace> source_coords_device("source coordinates",
81  number_source_coords, 3);
82  Kokkos::View<double**>::HostMirror source_coords = Kokkos::create_mirror_view(source_coords_device);
83 
84  // coordinates of target sites
85  Kokkos::View<double**, Kokkos::DefaultExecutionSpace> target_coords_device ("target coordinates", number_target_coords, 3);
86  Kokkos::View<double**>::HostMirror target_coords = Kokkos::create_mirror_view(target_coords_device);
87 
88 
89  // fill source coordinates with a uniform grid
90  int source_index = 0;
91  double this_coord[3] = {0,0,0};
92  for (int i=-n_neg1_to_1/2; i<n_neg1_to_1/2+1; ++i) {
93  this_coord[0] = i*h_spacing;
94  for (int j=-n_neg1_to_1/2; j<n_neg1_to_1/2+1; ++j) {
95  this_coord[1] = j*h_spacing;
96  for (int k=-n_neg1_to_1/2; k<n_neg1_to_1/2+1; ++k) {
97  this_coord[2] = k*h_spacing;
98  if (dimension==3) {
99  source_coords(source_index,0) = this_coord[0];
100  source_coords(source_index,1) = this_coord[1];
101  source_coords(source_index,2) = this_coord[2];
102  source_index++;
103  }
104  }
105  if (dimension==2) {
106  source_coords(source_index,0) = this_coord[0];
107  source_coords(source_index,1) = this_coord[1];
108  source_coords(source_index,2) = 0;
109  source_index++;
110  }
111  }
112  if (dimension==1) {
113  source_coords(source_index,0) = this_coord[0];
114  source_coords(source_index,1) = 0;
115  source_coords(source_index,2) = 0;
116  source_index++;
117  }
118  }
119 
120  // fill target coords somewhere inside of [-0.5,0.5]x[-0.5,0.5]x[-0.5,0.5]
121  for(int i=0; i<number_target_coords; i++){
122 
123  // first, we get a uniformly random distributed direction
124  double rand_dir[3] = {0,0,0};
125 
126  for (int j=0; j<dimension; ++j) {
127  // rand_dir[j] is in [-0.5, 0.5]
128  rand_dir[j] = ((double)rand() / (double) RAND_MAX) - 0.5;
129  }
130 
131  // then we get a uniformly random radius
132  for (int j=0; j<dimension; ++j) {
133  target_coords(i,j) = rand_dir[j];
134  }
135 
136  }
137 
138 
139  //! [Setting Up The Point Cloud]
140 
141  Kokkos::Profiling::popRegion();
142  Kokkos::Profiling::pushRegion("Creating Data");
143 
144  //! [Creating The Data]
145 
146 
147  // source coordinates need copied to device before using to construct sampling data
148  Kokkos::deep_copy(source_coords_device, source_coords);
149 
150  // target coordinates copied next, because it is a convenient time to send them to device
151  Kokkos::deep_copy(target_coords_device, target_coords);
152 
153  // need Kokkos View storing true solution
154  Kokkos::View<double*, Kokkos::DefaultExecutionSpace> sampling_data_device("samples of true solution",
155  source_coords_device.extent(0));
156 
157  Kokkos::View<double**, Kokkos::DefaultExecutionSpace> gradient_sampling_data_device("samples of true gradient",
158  source_coords_device.extent(0), dimension);
159 
160  Kokkos::View<double**, Kokkos::DefaultExecutionSpace> divergence_sampling_data_device
161  ("samples of true solution for divergence test", source_coords_device.extent(0), dimension);
162 
163  Kokkos::parallel_for("Sampling Manufactured Solutions", Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>
164  (0,source_coords.extent(0)), KOKKOS_LAMBDA(const int i) {
165 
166  // coordinates of source site i
167  double xval = source_coords_device(i,0);
168  double yval = (dimension>1) ? source_coords_device(i,1) : 0;
169  double zval = (dimension>2) ? source_coords_device(i,2) : 0;
170 
171  // data for targets with scalar input
172  sampling_data_device(i) = trueSolution(xval, yval, zval, order, dimension);
173 
174  // data for targets with vector input (divergence)
175  double true_grad[3] = {0,0,0};
176  trueGradient(true_grad, xval, yval,zval, order, dimension);
177 
178  for (int j=0; j<dimension; ++j) {
179  gradient_sampling_data_device(i,j) = true_grad[j];
180 
181  // data for target with vector input (curl)
182  divergence_sampling_data_device(i,j) = divergenceTestSamples(xval, yval, zval, j, dimension);
183  }
184 
185  });
186 
187 
188  //! [Creating The Data]
189 
190  Kokkos::Profiling::popRegion();
191  Kokkos::Profiling::pushRegion("Neighbor Search");
192 
193  //! [Performing Neighbor Search]
194 
195 
196  // Point cloud construction for neighbor search
197  // CreatePointCloudSearch constructs an object of type PointCloudSearch, but deduces the templates for you
198  auto point_cloud_search(CreatePointCloudSearch(source_coords, dimension));
199 
200  double epsilon_multiplier = 1.4;
201 
202  // neighbor_lists_device will contain all neighbor lists (for each target site) in a compressed row format
203  // Initially, we do a dry-run to calculate neighborhood sizes before actually storing the result. This is
204  // why we can start with a neighbor_lists_device size of 0.
205  Kokkos::View<int*> neighbor_lists_device("neighbor lists",
206  0); // first column is # of neighbors
207  Kokkos::View<int*>::HostMirror neighbor_lists = Kokkos::create_mirror_view(neighbor_lists_device);
208  // number_of_neighbors_list must be the same size as the number of target sites so that it can be populated
209  // with the number of neighbors for each target site.
210  Kokkos::View<int*> number_of_neighbors_list_device("number of neighbor lists",
211  number_target_coords); // first column is # of neighbors
212  Kokkos::View<int*>::HostMirror number_of_neighbors_list = Kokkos::create_mirror_view(number_of_neighbors_list_device);
213 
214  // each target site has a window size
215  Kokkos::View<double*, Kokkos::DefaultExecutionSpace> epsilon_device("h supports", number_target_coords);
216  Kokkos::View<double*>::HostMirror epsilon = Kokkos::create_mirror_view(epsilon_device);
217 
218  // query the point cloud to generate the neighbor lists using a kdtree to produce the n nearest neighbor
219  // to each target site, adding (epsilon_multiplier-1)*100% to whatever the distance away the further neighbor used is from
220  // each target to the view for epsilon
221  //
222  // This dry run populates number_of_neighbors_list with neighborhood sizes
223  size_t storage_size = point_cloud_search.generateCRNeighborListsFromKNNSearch(true /*dry run*/, target_coords, neighbor_lists,
224  number_of_neighbors_list, epsilon, min_neighbors, epsilon_multiplier);
225 
226  // resize neighbor_lists_device so as to be large enough to contain all neighborhoods
227  Kokkos::resize(neighbor_lists_device, storage_size);
228  neighbor_lists = Kokkos::create_mirror_view(neighbor_lists_device);
229 
230  // query the point cloud a second time, but this time storing results into neighbor_lists
231  point_cloud_search.generateCRNeighborListsFromKNNSearch(false /*not dry run*/, target_coords, neighbor_lists,
232  number_of_neighbors_list, epsilon, min_neighbors, epsilon_multiplier);
233 
234  //! [Performing Neighbor Search]
235 
236  Kokkos::Profiling::popRegion();
237  Kokkos::fence(); // let call to build neighbor lists complete before copying back to device
238  timer.reset();
239 
240  //! [Setting Up The GMLS Object]
241 
242 
243  // Copy data back to device (they were filled on the host)
244  // We could have filled Kokkos Views with memory space on the host
245  // and used these instead, and then the copying of data to the device
246  // would be performed in the GMLS class
247  Kokkos::deep_copy(neighbor_lists_device, neighbor_lists);
248  Kokkos::deep_copy(number_of_neighbors_list_device, number_of_neighbors_list);
249  Kokkos::deep_copy(epsilon_device, epsilon);
250 
251  // initialize an instance of the GMLS class
253  order, dimension,
254  solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
255  2 /*manifold order*/);
256 
257  // pass in neighbor lists, number of neighbor lists, source coordinates, target coordinates, and window sizes
258  //
259  // neighbor lists has a compressed row format and is a 1D view:
260  // dimensions: ? (determined by neighbor search, but total size of the sum of the number of neighbors over all target sites)
261  //
262  // number of neighbors list is a 1D view:
263  // dimensions: (# number of target sites)
264  // each entry contains the number of neighbors for a target site
265  //
266  // source coordinates have the format:
267  // dimensions: (# number of source sites) X (dimension)
268  // entries in the neighbor lists (integers) correspond to rows of this 2D array
269  //
270  // target coordinates have the format:
271  // dimensions: (# number of target sites) X (dimension)
272  // # of target sites is same as # of rows of neighbor lists
273  //
274  my_GMLS.setProblemData(neighbor_lists_device, number_of_neighbors_list_device, source_coords_device, target_coords_device, epsilon_device);
275 
276  // create a vector of target operations
277  std::vector<TargetOperation> lro(5);
278  lro[0] = ScalarPointEvaluation;
283 
284  // and then pass them to the GMLS class
285  my_GMLS.addTargets(lro);
286 
287  // sets the weighting kernel function from WeightingFunctionType
288  my_GMLS.setWeightingType(WeightingFunctionType::Power);
289 
290  // power to use in that weighting kernel function
291  my_GMLS.setWeightingPower(2);
292 
293  // generate the alphas that to be combined with data for each target operation requested in lro
294  my_GMLS.generateAlphas(number_of_batches, keep_coefficients /* keep polynomial coefficients, only needed for a test later in this program */);
295 
296 
297  //! [Setting Up The GMLS Object]
298 
299  double instantiation_time = timer.seconds();
300  std::cout << "Took " << instantiation_time << "s to complete alphas generation." << std::endl;
301  Kokkos::fence(); // let generateAlphas finish up before using alphas
302  Kokkos::Profiling::pushRegion("Apply Alphas to Data");
303 
304  //! [Apply GMLS Alphas To Data]
305 
306  // it is important to note that if you expect to use the data as a 1D view, then you should use double*
307  // however, if you know that the target operation will result in a 2D view (vector or matrix output),
308  // then you should template with double** as this is something that can not be infered from the input data
309  // or the target operator at compile time. Additionally, a template argument is required indicating either
310  // Kokkos::HostSpace or Kokkos::DefaultExecutionSpace::memory_space()
311 
312  // The Evaluator class takes care of handling input data views as well as the output data views.
313  // It uses information from the GMLS class to determine how many components are in the input
314  // as well as output for any choice of target functionals and then performs the contactions
315  // on the data using the alpha coefficients generated by the GMLS class, all on the device.
316  Evaluator gmls_evaluator(&my_GMLS);
317 
318  auto output_value = gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double*, Kokkos::HostSpace>
319  (sampling_data_device, ScalarPointEvaluation);
320 
321  auto output_laplacian = gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double*, Kokkos::HostSpace>
322  (sampling_data_device, LaplacianOfScalarPointEvaluation);
323 
324  auto output_gradient = gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double**, Kokkos::HostSpace>
325  (sampling_data_device, GradientOfScalarPointEvaluation);
326 
327  auto output_divergence = gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double*, Kokkos::HostSpace>
328  (gradient_sampling_data_device, DivergenceOfVectorPointEvaluation, VectorPointSample);
329 
330  auto output_curl = gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double**, Kokkos::HostSpace>
331  (divergence_sampling_data_device, CurlOfVectorPointEvaluation);
332 
333  // retrieves polynomial coefficients instead of remapped field
334  decltype(output_curl) scalar_coefficients;
335  if (number_of_batches==1)
336  scalar_coefficients =
337  gmls_evaluator.applyFullPolynomialCoefficientsBasisToDataAllComponents<double**, Kokkos::HostSpace>
338  (sampling_data_device);
339 
340  //! [Apply GMLS Alphas To Data]
341 
342  Kokkos::fence(); // let application of alphas to data finish before using results
343  Kokkos::Profiling::popRegion();
344  // times the Comparison in Kokkos
345  Kokkos::Profiling::pushRegion("Comparison");
346 
347  //! [Check That Solutions Are Correct]
348 
349 
350  // loop through the target sites
351  for (int i=0; i<number_target_coords; i++) {
352 
353  // load value from output
354  double GMLS_value = output_value(i);
355 
356  // load laplacian from output
357  double GMLS_Laplacian = output_laplacian(i);
358 
359  // load partial x from gradient
360  // this is a test that the scalar_coefficients 2d array returned hold valid entries
361  // scalar_coefficients(i,1)*1./epsilon(i) is equivalent to the target operation acting
362  // on the polynomials applied to the polynomial coefficients
363  double GMLS_GradX = (number_of_batches==1) ? scalar_coefficients(i,1)*1./epsilon(i) : output_gradient(i,0);
364 
365  // load partial y from gradient
366  double GMLS_GradY = (dimension>1) ? output_gradient(i,1) : 0;
367 
368  // load partial z from gradient
369  double GMLS_GradZ = (dimension>2) ? output_gradient(i,2) : 0;
370 
371  // load divergence from output
372  double GMLS_Divergence = output_divergence(i);
373 
374  // load curl from output
375  double GMLS_CurlX = (dimension>1) ? output_curl(i,0) : 0;
376  double GMLS_CurlY = (dimension>1) ? output_curl(i,1) : 0;
377  double GMLS_CurlZ = (dimension>2) ? output_curl(i,2) : 0;
378 
379 
380  // target site i's coordinate
381  double xval = target_coords(i,0);
382  double yval = (dimension>1) ? target_coords(i,1) : 0;
383  double zval = (dimension>2) ? target_coords(i,2) : 0;
384 
385  // evaluation of various exact solutions
386  double actual_value = trueSolution(xval, yval, zval, order, dimension);
387  double actual_Laplacian = trueLaplacian(xval, yval, zval, order, dimension);
388 
389  double actual_Gradient[3] = {0,0,0}; // initialized for 3, but only filled up to dimension
390  trueGradient(actual_Gradient, xval, yval, zval, order, dimension);
391 
392  double actual_Divergence;
393  actual_Divergence = trueLaplacian(xval, yval, zval, order, dimension);
394 
395  double actual_Curl[3] = {0,0,0}; // initialized for 3, but only filled up to dimension
396  // (and not at all for dimimension = 1)
397  if (dimension>1) {
398  actual_Curl[0] = curlTestSolution(xval, yval, zval, 0, dimension);
399  actual_Curl[1] = curlTestSolution(xval, yval, zval, 1, dimension);
400  if (dimension>2) {
401  actual_Curl[2] = curlTestSolution(xval, yval, zval, 2, dimension);
402  }
403  }
404 
405  // check actual function value
406  if(GMLS_value!=GMLS_value || std::abs(actual_value - GMLS_value) > failure_tolerance) {
407  all_passed = false;
408  std::cout << i << " Failed Actual by: " << std::abs(actual_value - GMLS_value) << std::endl;
409  }
410 
411  // check Laplacian
412  if(std::abs(actual_Laplacian - GMLS_Laplacian) > laplacian_failure_tolerance) {
413  all_passed = false;
414  std::cout << i <<" Failed Laplacian by: " << std::abs(actual_Laplacian - GMLS_Laplacian) << std::endl;
415  }
416 
417  // check gradient
418  if(std::abs(actual_Gradient[0] - GMLS_GradX) > failure_tolerance) {
419  all_passed = false;
420  std::cout << i << " Failed GradX by: " << std::abs(actual_Gradient[0] - GMLS_GradX) << std::endl;
421  if (dimension>1) {
422  if(std::abs(actual_Gradient[1] - GMLS_GradY) > failure_tolerance) {
423  all_passed = false;
424  std::cout << i << " Failed GradY by: " << std::abs(actual_Gradient[1] - GMLS_GradY) << std::endl;
425  }
426  }
427  if (dimension>2) {
428  if(std::abs(actual_Gradient[2] - GMLS_GradZ) > failure_tolerance) {
429  all_passed = false;
430  std::cout << i << " Failed GradZ by: " << std::abs(actual_Gradient[2] - GMLS_GradZ) << std::endl;
431  }
432  }
433  }
434 
435  // check divergence
436  if(std::abs(actual_Divergence - GMLS_Divergence) > failure_tolerance) {
437  all_passed = false;
438  std::cout << i << " Failed Divergence by: " << std::abs(actual_Divergence - GMLS_Divergence) << std::endl;
439  }
440 
441  // check curl
442  if (order > 2) { // reconstructed solution not in basis unless order greater than 2 used
443  double tmp_diff = 0;
444  if (dimension>1)
445  tmp_diff += std::abs(actual_Curl[0] - GMLS_CurlX) + std::abs(actual_Curl[1] - GMLS_CurlY);
446  if (dimension>2)
447  tmp_diff += std::abs(actual_Curl[2] - GMLS_CurlZ);
448  if(std::abs(tmp_diff) > failure_tolerance) {
449  all_passed = false;
450  std::cout << i << " Failed Curl by: " << std::abs(tmp_diff) << std::endl;
451  }
452  }
453  }
454 
455 
456  //! [Check That Solutions Are Correct]
457  // popRegion hidden from tutorial
458  // stop timing comparison loop
459  Kokkos::Profiling::popRegion();
460  //! [Finalize Program]
461 
462 
463 } // end of code block to reduce scope, causing Kokkos View de-allocations
464 // otherwise, Views may be deallocating when we call Kokkos finalize() later
465 
466 // finalize Kokkos and MPI (if available)
467 kp.finalize();
468 #ifdef COMPADRE_USE_MPI
469 MPI_Finalize();
470 #endif
471 
472 // output to user that test passed or failed
473 if(all_passed) {
474  fprintf(stdout, "Passed test \n");
475  return 0;
476 } else {
477  fprintf(stdout, "Failed test \n");
478  return -1;
479 }
480 
481 } // main
482 
483 
484 //! [Finalize Program]
Point evaluation of the curl of a vector (results in a vector)
Class handling Kokkos command line arguments and returning parameters.
Lightweight Evaluator Helper This class is a lightweight wrapper for extracting and applying all rele...
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...
int main(int argc, char *args[])
[Parse Command Line Arguments]
Definition: GMLS_Device.cpp:30
Point evaluation of a scalar.
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 void trueGradient(double *ans, double x, double y, double z, int order, int dimension)
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 (allocates memory for output)
Kokkos::View< output_data_type, output_array_layout, output_memory_space > applyFullPolynomialCoefficientsBasisToDataAllComponents(view_type_input_data sampling_data, bool scalar_as_vector_if_needed=true) const
Generation of polynomial reconstruction coefficients by applying to data in GMLS (allocates memory fo...
Point evaluation of the laplacian of a scalar (could be on a manifold or not)
Point evaluation of the divergence of a vector (results in a scalar)
Point evaluation of the gradient of a scalar.
KOKKOS_INLINE_FUNCTION double trueLaplacian(double x, double y, double z, int order, int dimension)
Generalized Moving Least Squares (GMLS)
KOKKOS_INLINE_FUNCTION double curlTestSolution(double x, double y, double z, int component, int dimension)
KOKKOS_INLINE_FUNCTION double divergenceTestSamples(double x, double y, double z, int component, int dimension)
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) ...
KOKKOS_INLINE_FUNCTION double trueSolution(double x, double y, double z, int order, int dimension)
constexpr SamplingFunctional VectorPointSample
Point evaluations of the entire vector source function.