Compadre  1.3.3
GMLS_SmallBatchReuse_Device.cpp
Go to the documentation of this file.
1 /*
2  *
3  * This examples tests the ability to reuse a GMLS class instance,
4  * changing the target and neighbor list for the new target.
5  *
6  */
7 
8 #include <iostream>
9 #include <string>
10 #include <vector>
11 #include <map>
12 #include <stdlib.h>
13 #include <cstdio>
14 #include <random>
15 
16 #include <Compadre_Config.h>
17 #include <Compadre_GMLS.hpp>
18 #include <Compadre_Evaluator.hpp>
20 
21 #include "GMLS_Tutorial.hpp"
22 #include "CommandLineProcessor.hpp"
23 
24 #ifdef COMPADRE_USE_MPI
25 #include <mpi.h>
26 #endif
27 
28 #include <Kokkos_Timer.hpp>
29 #include <Kokkos_Core.hpp>
30 
31 using namespace Compadre;
32 
33 //! [Parse Command Line Arguments]
34 
35 // called from command line
36 int main (int argc, char* args[]) {
37 
38 // initializes MPI (if available) with command line arguments given
39 #ifdef COMPADRE_USE_MPI
40 MPI_Init(&argc, &args);
41 #endif
42 
43 // initializes Kokkos with command line arguments given
44 Kokkos::initialize(argc, args);
45 
46 // becomes false if the computed solution not within the failure_threshold of the actual solution
47 bool all_passed = true;
48 
49 // code block to reduce scope for all Kokkos View allocations
50 // otherwise, Views may be deallocating when we call Kokkos::finalize() later
51 {
52 
53  CommandLineProcessor clp(argc, args);
54  auto order = clp.order;
55  auto dimension = clp.dimension;
56  auto number_target_coords = clp.number_target_coords;
57  auto constraint_name = clp.constraint_name;
58  auto solver_name = clp.solver_name;
59  auto problem_name = clp.problem_name;
60 
61  // the functions we will be seeking to reconstruct are in the span of the basis
62  // of the reconstruction space we choose for GMLS, so the error should be very small
63  const double failure_tolerance = 1e-9;
64 
65  // Laplacian is a second order differential operator, which we expect to be slightly less accurate
66  const double laplacian_failure_tolerance = 1e-9;
67 
68  // minimum neighbors for unisolvency is the same as the size of the polynomial basis
69  const int min_neighbors = Compadre::GMLS::getNP(order, dimension);
70 
71  //! [Parse Command Line Arguments]
72  Kokkos::Timer timer;
73  //! [Setting Up The Point Cloud]
74 
75  // approximate spacing of source sites
76  double h_spacing = 0.05;
77  int n_neg1_to_1 = 2*(1/h_spacing) + 1; // always odd
78 
79  // number of source coordinate sites that will fill a box of [-1,1]x[-1,1]x[-1,1] with a spacing approximately h
80  const int number_source_coords = std::pow(n_neg1_to_1, dimension);
81 
82  // coordinates of source sites
83  Kokkos::View<double**, Kokkos::DefaultExecutionSpace> source_coords_device("source coordinates",
84  number_source_coords, 3);
85  Kokkos::View<double**>::HostMirror source_coords = Kokkos::create_mirror_view(source_coords_device);
86 
87  // coordinates of target sites
88  Kokkos::View<double**, Kokkos::DefaultExecutionSpace> target_coords_device ("target coordinates", number_target_coords, 3);
89  Kokkos::View<double**>::HostMirror target_coords = Kokkos::create_mirror_view(target_coords_device);
90 
91 
92  // fill source coordinates with a uniform grid
93  int source_index = 0;
94  double this_coord[3] = {0,0,0};
95  for (int i=-n_neg1_to_1/2; i<n_neg1_to_1/2+1; ++i) {
96  this_coord[0] = i*h_spacing;
97  for (int j=-n_neg1_to_1/2; j<n_neg1_to_1/2+1; ++j) {
98  this_coord[1] = j*h_spacing;
99  for (int k=-n_neg1_to_1/2; k<n_neg1_to_1/2+1; ++k) {
100  this_coord[2] = k*h_spacing;
101  if (dimension==3) {
102  source_coords(source_index,0) = this_coord[0];
103  source_coords(source_index,1) = this_coord[1];
104  source_coords(source_index,2) = this_coord[2];
105  source_index++;
106  }
107  }
108  if (dimension==2) {
109  source_coords(source_index,0) = this_coord[0];
110  source_coords(source_index,1) = this_coord[1];
111  source_coords(source_index,2) = 0;
112  source_index++;
113  }
114  }
115  if (dimension==1) {
116  source_coords(source_index,0) = this_coord[0];
117  source_coords(source_index,1) = 0;
118  source_coords(source_index,2) = 0;
119  source_index++;
120  }
121  }
122 
123  // fill target coords somewhere inside of [-0.5,0.5]x[-0.5,0.5]x[-0.5,0.5]
124  for(int i=0; i<number_target_coords; i++){
125 
126  // first, we get a uniformly random distributed direction
127  double rand_dir[3] = {0,0,0};
128 
129  for (int j=0; j<dimension; ++j) {
130  // rand_dir[j] is in [-0.5, 0.5]
131  rand_dir[j] = ((double)rand() / (double) RAND_MAX) - 0.5;
132  }
133 
134  // then we get a uniformly random radius
135  for (int j=0; j<dimension; ++j) {
136  target_coords(i,j) = rand_dir[j];
137  }
138 
139  }
140 
141 
142  //! [Setting Up The Point Cloud]
143 
144 
145  //! [Creating The Data]
146 
147 
148  // source coordinates need copied to device before using to construct sampling data
149  Kokkos::deep_copy(source_coords_device, source_coords);
150 
151  // target coordinates copied next, because it is a convenient time to send them to device
152  Kokkos::deep_copy(target_coords_device, target_coords);
153 
154  // need Kokkos View storing true solution
155  Kokkos::View<double*, Kokkos::DefaultExecutionSpace> sampling_data_device("samples of true solution",
156  source_coords_device.extent(0));
157 
158  Kokkos::View<double**, Kokkos::DefaultExecutionSpace> gradient_sampling_data_device("samples of true gradient",
159  source_coords_device.extent(0), dimension);
160 
161  Kokkos::View<double**, Kokkos::DefaultExecutionSpace> divergence_sampling_data_device
162  ("samples of true solution for divergence test", source_coords_device.extent(0), dimension);
163 
164  Kokkos::parallel_for("Sampling Manufactured Solutions", Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>
165  (0,source_coords.extent(0)), KOKKOS_LAMBDA(const int i) {
166 
167  // coordinates of source site i
168  double xval = source_coords_device(i,0);
169  double yval = (dimension>1) ? source_coords_device(i,1) : 0;
170  double zval = (dimension>2) ? source_coords_device(i,2) : 0;
171 
172  // data for targets with scalar input
173  sampling_data_device(i) = trueSolution(xval, yval, zval, order, dimension);
174 
175  // data for targets with vector input (divergence)
176  double true_grad[3] = {0,0,0};
177  trueGradient(true_grad, xval, yval,zval, order, dimension);
178 
179  for (int j=0; j<dimension; ++j) {
180  gradient_sampling_data_device(i,j) = true_grad[j];
181 
182  // data for target with vector input (curl)
183  divergence_sampling_data_device(i,j) = divergenceTestSamples(xval, yval, zval, j, dimension);
184  }
185 
186  });
187 
188 
189  //! [Creating The Data]
190 
191 
192  //! [Setting Up The GMLS Object]
193 
194  // initialize an instance of the GMLS class
196  order, dimension,
197  solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
198  2 /*manifold order*/);
199 
200  // create a vector of target operations
201  std::vector<TargetOperation> lro(5);
202  lro[0] = ScalarPointEvaluation;
207 
208  // and then pass them to the GMLS class
209  my_GMLS.addTargets(lro);
210 
211  // sets the weighting kernel function from WeightingFunctionType
212  my_GMLS.setWeightingType(WeightingFunctionType::Power);
213 
214  // power to use in that weighting kernel function
215  my_GMLS.setWeightingPower(2);
216 
217  // set source sites once for all targets
218  my_GMLS.setSourceSites(source_coords_device);
219 
220 
221  // Point cloud construction for neighbor search
222  // CreatePointCloudSearch constructs an object of type PointCloudSearch, but deduces the templates for you
223  auto point_cloud_search(CreatePointCloudSearch(source_coords, dimension));
224 
225  // loop through the target sites
226  for (int i=0; i<number_target_coords; i++) {
227  timer.reset();
228  // pass in neighbor lists, source coordinates, target coordinates, and window sizes
229  //
230  // single neighbor lists have the format:
231  // dimensions: (1 single neighbor list for one target) X (# maximum number of neighbors for any given target + 1)
232  // the first column contains the number of neighbors for that rows corresponding target index
233  //
234  // source coordinates have the format:
235  // dimensions: (# number of source sites) X (dimension)
236  // entries in the neighbor lists (integers) correspond to rows of this 2D array
237  //
238  // single target coordinates have the format:
239  // dimensions: (1 single target site) X (dimension)
240  // # of target sites is same as # of rows of neighbor lists
241  //
242 
243 
244  // coordinates of single target sites
245  Kokkos::View<double**, Kokkos::DefaultExecutionSpace> single_target_coords_device ("single target coordinates", 1, 3);
246  Kokkos::View<double**>::HostMirror single_target_coords = Kokkos::create_mirror_view(single_target_coords_device);
247  for (int j=0; j<3; ++j) {
248  single_target_coords(0,j) = target_coords(i,j);
249  }
250  // target coordinates copied next, because it is a convenient time to send them to device
251  Kokkos::deep_copy(single_target_coords_device, single_target_coords);
252  Kokkos::fence();
253 
254  //! [Performing Neighbor Search]
255 
256  // each row is a neighbor list for a target site, with the first column of each row containing
257  // the number of neighbors for that rows corresponding target site
258  double epsilon_multiplier = 1.5;
259  int estimated_upper_bound_number_neighbors =
260  point_cloud_search.getEstimatedNumberNeighborsUpperBound(min_neighbors, dimension, epsilon_multiplier);
261 
262  Kokkos::View<int**, Kokkos::DefaultExecutionSpace> single_neighbor_lists_device("neighbor lists",
263  1, estimated_upper_bound_number_neighbors); // first column is # of neighbors
264  Kokkos::View<int**>::HostMirror single_neighbor_lists = Kokkos::create_mirror_view(single_neighbor_lists_device);
265 
266  // each target site has a window size
267  Kokkos::View<double*, Kokkos::DefaultExecutionSpace> single_epsilon_device("h supports", 1);
268  Kokkos::View<double*>::HostMirror single_epsilon = Kokkos::create_mirror_view(single_epsilon_device);
269 
270  // query the point cloud to generate the neighbor lists using a kdtree to produce the n nearest neighbor
271  // to each target site, adding (epsilon_multiplier-1)*100% to whatever the distance away the further neighbor used is from
272  // each target to the view for epsilon
273  point_cloud_search.generate2DNeighborListsFromKNNSearch(false /*not dry run*/, single_target_coords,
274  single_neighbor_lists, single_epsilon, min_neighbors, epsilon_multiplier);
275 
276  //! [Performing Neighbor Search]
277 
278 
279  // Copy data back to device (they were filled on the host)
280  // We could have filled Kokkos Views with memory space on the host
281  // and used these instead, and then the copying of data to the device
282  // would be performed in the GMLS class
283  Kokkos::deep_copy(single_neighbor_lists_device, single_neighbor_lists);
284  Kokkos::deep_copy(single_epsilon_device, single_epsilon);
285  Kokkos::fence(); // let call to build neighbor lists complete before copying back to device
286 
287  // Create temporary small views to hold just one target's information at a time
288  my_GMLS.setNeighborLists(single_neighbor_lists_device);
289  my_GMLS.setTargetSites(single_target_coords_device);
290  my_GMLS.setWindowSizes(single_epsilon_device);
291 
292 
293  // generate the alphas that to be combined with data for each target operation requested in lro
294  my_GMLS.generateAlphas(1, true /* 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 
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  auto scalar_coefficients = gmls_evaluator.applyFullPolynomialCoefficientsBasisToDataAllComponents<double**, Kokkos::HostSpace>
335  (sampling_data_device);
336 
337  //! [Apply GMLS Alphas To Data]
338 
339  Kokkos::fence(); // let application of alphas to data finish before using results
340  // times the Comparison in Kokkos
341 
342  //! [Check That Solutions Are Correct]
343 
344  // load value from output
345  double GMLS_value = output_value(0);
346 
347  // load laplacian from output
348  double GMLS_Laplacian = output_laplacian(0);
349 
350  // load partial x from gradient
351  // this is a test that the scalar_coefficients 2d array returned hold valid entries
352  // scalar_coefficients(i,1)*1./epsilon(i) is equivalent to the target operation acting
353  // on the polynomials applied to the polynomial coefficients
354  double GMLS_GradX = scalar_coefficients(0,1)*1./single_epsilon(0);
355  //output_gradient(i,0);
356 
357  // load partial y from gradient
358  double GMLS_GradY = (dimension>1) ? output_gradient(0,1) : 0;
359 
360  // load partial z from gradient
361  double GMLS_GradZ = (dimension>2) ? output_gradient(0,2) : 0;
362 
363  // load divergence from output
364  double GMLS_Divergence = output_divergence(0);
365 
366  // load curl from output
367  double GMLS_CurlX = (dimension>1) ? output_curl(0,0) : 0;
368  double GMLS_CurlY = (dimension>1) ? output_curl(0,1) : 0;
369  double GMLS_CurlZ = (dimension>2) ? output_curl(0,2) : 0;
370 
371 
372  // target site i's coordinate
373  double xval = target_coords(i,0);
374  double yval = (dimension>1) ? target_coords(i,1) : 0;
375  double zval = (dimension>2) ? target_coords(i,2) : 0;
376 
377  // evaluation of various exact solutions
378  double actual_value = trueSolution(xval, yval, zval, order, dimension);
379  double actual_Laplacian = trueLaplacian(xval, yval, zval, order, dimension);
380 
381  double actual_Gradient[3] = {0,0,0}; // initialized for 3, but only filled up to dimension
382  trueGradient(actual_Gradient, xval, yval, zval, order, dimension);
383 
384  double actual_Divergence;
385  actual_Divergence = trueLaplacian(xval, yval, zval, order, dimension);
386 
387  double actual_Curl[3] = {0,0,0}; // initialized for 3, but only filled up to dimension
388  // (and not at all for dimimension = 1)
389  if (dimension>1) {
390  actual_Curl[0] = curlTestSolution(xval, yval, zval, 0, dimension);
391  actual_Curl[1] = curlTestSolution(xval, yval, zval, 1, dimension);
392  if (dimension>2) {
393  actual_Curl[2] = curlTestSolution(xval, yval, zval, 2, dimension);
394  }
395  }
396 
397  // check actual function value
398  if(GMLS_value!=GMLS_value || std::abs(actual_value - GMLS_value) > failure_tolerance) {
399  all_passed = false;
400  std::cout << i << " Failed Actual by: " << std::abs(actual_value - GMLS_value) << std::endl;
401  }
402 
403  // check Laplacian
404  if(std::abs(actual_Laplacian - GMLS_Laplacian) > laplacian_failure_tolerance) {
405  all_passed = false;
406  std::cout << i <<" Failed Laplacian by: " << std::abs(actual_Laplacian - GMLS_Laplacian) << std::endl;
407  }
408 
409  // check gradient
410  if(std::abs(actual_Gradient[0] - GMLS_GradX) > failure_tolerance) {
411  all_passed = false;
412  std::cout << i << " Failed GradX by: " << std::abs(actual_Gradient[0] - GMLS_GradX) << std::endl;
413  if (dimension>1) {
414  if(std::abs(actual_Gradient[1] - GMLS_GradY) > failure_tolerance) {
415  all_passed = false;
416  std::cout << i << " Failed GradY by: " << std::abs(actual_Gradient[1] - GMLS_GradY) << std::endl;
417  }
418  }
419  if (dimension>2) {
420  if(std::abs(actual_Gradient[2] - GMLS_GradZ) > failure_tolerance) {
421  all_passed = false;
422  std::cout << i << " Failed GradZ by: " << std::abs(actual_Gradient[2] - GMLS_GradZ) << std::endl;
423  }
424  }
425  }
426 
427  // check divergence
428  if(std::abs(actual_Divergence - GMLS_Divergence) > failure_tolerance) {
429  all_passed = false;
430  std::cout << i << " Failed Divergence by: " << std::abs(actual_Divergence - GMLS_Divergence) << std::endl;
431  }
432 
433  // check curl
434  if (order > 2) { // reconstructed solution not in basis unless order greater than 2 used
435  double tmp_diff = 0;
436  if (dimension>1)
437  tmp_diff += std::abs(actual_Curl[0] - GMLS_CurlX) + std::abs(actual_Curl[1] - GMLS_CurlY);
438  if (dimension>2)
439  tmp_diff += std::abs(actual_Curl[2] - GMLS_CurlZ);
440  if(std::abs(tmp_diff) > failure_tolerance) {
441  all_passed = false;
442  std::cout << i << " Failed Curl by: " << std::abs(tmp_diff) << std::endl;
443  }
444  }
445  }
446 
447 
448  //! [Check That Solutions Are Correct]
449  // stop timing comparison loop
450  //! [Finalize Program]
451 
452 
453 } // end of code block to reduce scope, causing Kokkos View de-allocations
454 // otherwise, Views may be deallocating when we call Kokkos::finalize() later
455 
456 // finalize Kokkos and MPI (if available)
457 Kokkos::finalize();
458 #ifdef COMPADRE_USE_MPI
459 MPI_Finalize();
460 #endif
461 
462 // output to user that test passed or failed
463 if(all_passed) {
464  fprintf(stdout, "Passed test \n");
465  return 0;
466 } else {
467  fprintf(stdout, "Failed test \n");
468  return -1;
469 }
470 
471 } // main
472 
473 
474 //! [Finalize Program]
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...
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...
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.
int main(int argc, char *args[])
[Parse Command Line Arguments]
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)
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.