Compadre  1.3.3
GMLS_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_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 
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 
55  // the functions we will be seeking to reconstruct are in the span of the basis
56  // of the reconstruction space we choose for GMLS, so the error should be very small
57  const double failure_tolerance = 1e-9;
58 
59  // minimum neighbors for unisolvency is the same as the size of the polynomial basis
60  const int min_neighbors = Compadre::GMLS::getNP(order, dimension);
61 
62  //! [Parse Command Line Arguments]
63  Kokkos::Timer timer;
64  Kokkos::Profiling::pushRegion("Setup Point Data");
65  //! [Setting Up The Point Cloud]
66 
67  // approximate spacing of source sites
68  double h_spacing = 0.05;
69  int n_neg1_to_1 = 2*(1/h_spacing) + 1; // always odd
70 
71  // number of source coordinate sites that will fill a box of [-1,1]x[-1,1]x[-1,1] with a spacing approximately h
72  const int number_source_coords = std::pow(n_neg1_to_1, dimension);
73 
74  // coordinates of source sites
75  Kokkos::View<double**, Kokkos::DefaultExecutionSpace> source_coords_device("source coordinates",
76  number_source_coords, 3);
77  Kokkos::View<double**>::HostMirror source_coords = Kokkos::create_mirror_view(source_coords_device);
78 
79  // coordinates of target sites
80  Kokkos::View<double**, Kokkos::DefaultExecutionSpace> target_coords_device ("target coordinates", number_target_coords, 3);
81  Kokkos::View<double**>::HostMirror target_coords = Kokkos::create_mirror_view(target_coords_device);
82 
83  // coordinates of additional evaluation sites
84  Kokkos::View<double**, Kokkos::DefaultExecutionSpace> additional_target_coords_device ("additional target coordinates", 2*number_target_coords /* multiple evaluation sites for each target index */, 3);
85  Kokkos::View<double**>::HostMirror additional_target_coords = Kokkos::create_mirror_view(additional_target_coords_device);
86 
87  // additional target site indices
88  Kokkos::View<int**, Kokkos::DefaultExecutionSpace> additional_target_indices_device ("additional target indices", number_target_coords, 4 /* # of extra evaluation sites plus index for each */);
89  Kokkos::View<int**>::HostMirror additional_target_indices = Kokkos::create_mirror_view(additional_target_indices_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  // generate coordinates to test multiple site evaluations
142  // strategy is to have a variable number of evaluation sites per target site
143  // so as to fully test the multi-site evaluation
144  int extra_evaluation_coordinates_count = 0;
145  for(int i=0; i<number_target_coords; i++){
146 
147  // set list of indices for extra evaluations
148  additional_target_indices(i,0) = (i%3)+1;
149 
150  // evaluation sites are same as target plus some perturbation
151  for (int k=0; k<(i%3+1); ++k) {
152  for (int j=0; j<dimension; ++j) {
153  additional_target_coords(extra_evaluation_coordinates_count,j) = target_coords(i,j) + (j==0)*1e-3 + (j==1)*1e-2 + (j==1)*(-1e-1);
154  }
155  additional_target_indices(i,k+1) = extra_evaluation_coordinates_count;
156  extra_evaluation_coordinates_count++;
157  }
158  }
159 
160 
161  //! [Setting Up The Point Cloud]
162 
163  Kokkos::Profiling::popRegion();
164  Kokkos::Profiling::pushRegion("Creating Data");
165 
166  //! [Creating The Data]
167 
168 
169  // source coordinates need copied to device before using to construct sampling data
170  Kokkos::deep_copy(source_coords_device, source_coords);
171 
172  // target coordinates copied next, because it is a convenient time to send them to device
173  Kokkos::deep_copy(target_coords_device, target_coords);
174 
175  // additional evaluation coordinates copied next, because it is a convenient time to send them to device
176  Kokkos::deep_copy(additional_target_coords_device, additional_target_coords);
177 
178  // additional evaluation indices copied next, because it is a convenient time to send them to device
179  Kokkos::deep_copy(additional_target_indices_device, additional_target_indices);
180 
181  // need Kokkos View storing true solution
182  Kokkos::View<double*, Kokkos::DefaultExecutionSpace> sampling_data_device("samples of true solution",
183  source_coords_device.extent(0));
184 
185  Kokkos::View<double**, Kokkos::DefaultExecutionSpace> gradient_sampling_data_device("samples of true gradient",
186  source_coords_device.extent(0), dimension);
187 
188  Kokkos::View<double**, Kokkos::DefaultExecutionSpace> divergence_sampling_data_device
189  ("samples of true solution for divergence test", source_coords_device.extent(0), dimension);
190 
191  Kokkos::parallel_for("Sampling Manufactured Solutions", Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>
192  (0,source_coords.extent(0)), KOKKOS_LAMBDA(const int i) {
193 
194  // coordinates of source site i
195  double xval = source_coords_device(i,0);
196  double yval = (dimension>1) ? source_coords_device(i,1) : 0;
197  double zval = (dimension>2) ? source_coords_device(i,2) : 0;
198 
199  // data for targets with scalar input
200  sampling_data_device(i) = trueSolution(xval, yval, zval, order, dimension);
201 
202  // data for targets with vector input (divergence)
203  double true_grad[3] = {0,0,0};
204  trueGradient(true_grad, xval, yval,zval, order, dimension);
205 
206  for (int j=0; j<dimension; ++j) {
207  gradient_sampling_data_device(i,j) = true_grad[j];
208 
209  // data for target with vector input (curl)
210  divergence_sampling_data_device(i,j) = divergenceTestSamples(xval, yval, zval, j, dimension);
211  }
212 
213  });
214 
215 
216  //! [Creating The Data]
217 
218  Kokkos::Profiling::popRegion();
219  Kokkos::Profiling::pushRegion("Neighbor Search");
220 
221  //! [Performing Neighbor Search]
222 
223 
224  // Point cloud construction for neighbor search
225  // CreatePointCloudSearch constructs an object of type PointCloudSearch, but deduces the templates for you
226  auto point_cloud_search(CreatePointCloudSearch(source_coords, dimension));
227 
228  // each row is a neighbor list for a target site, with the first column of each row containing
229  // the number of neighbors for that rows corresponding target site
230  double epsilon_multiplier = 1.5;
231  int estimated_upper_bound_number_neighbors =
232  point_cloud_search.getEstimatedNumberNeighborsUpperBound(min_neighbors, dimension, epsilon_multiplier);
233 
234  Kokkos::View<int**, Kokkos::DefaultExecutionSpace> neighbor_lists_device("neighbor lists",
235  number_target_coords, estimated_upper_bound_number_neighbors); // first column is # of neighbors
236  Kokkos::View<int**>::HostMirror neighbor_lists = Kokkos::create_mirror_view(neighbor_lists_device);
237 
238  // each target site has a window size
239  Kokkos::View<double*, Kokkos::DefaultExecutionSpace> epsilon_device("h supports", number_target_coords);
240  Kokkos::View<double*>::HostMirror epsilon = Kokkos::create_mirror_view(epsilon_device);
241 
242  // query the point cloud to generate the neighbor lists using a kdtree to produce the n nearest neighbor
243  // to each target site, adding (epsilon_multiplier-1)*100% to whatever the distance away the further neighbor used is from
244  // each target to the view for epsilon
245  point_cloud_search.generate2DNeighborListsFromKNNSearch(false /*not dry run*/, target_coords, neighbor_lists,
246  epsilon, min_neighbors, epsilon_multiplier);
247 
248 
249  //! [Performing Neighbor Search]
250 
251  Kokkos::Profiling::popRegion();
252  Kokkos::fence(); // let call to build neighbor lists complete before copying back to device
253  timer.reset();
254 
255  //! [Setting Up The GMLS Object]
256 
257 
258  // Copy data back to device (they were filled on the host)
259  // We could have filled Kokkos Views with memory space on the host
260  // and used these instead, and then the copying of data to the device
261  // would be performed in the GMLS class
262  Kokkos::deep_copy(neighbor_lists_device, neighbor_lists);
263  Kokkos::deep_copy(epsilon_device, epsilon);
264 
265  // initialize an instance of the GMLS class
267  order, dimension,
268  solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
269  2 /*manifold order*/);
270 
271  // pass in neighbor lists, source coordinates, target coordinates, and window sizes
272  //
273  // neighbor lists have the format:
274  // dimensions: (# number of target sites) X (# maximum number of neighbors for any given target + 1)
275  // the first column contains the number of neighbors for that rows corresponding target index
276  //
277  // source coordinates have the format:
278  // dimensions: (# number of source sites) X (dimension)
279  // entries in the neighbor lists (integers) correspond to rows of this 2D array
280  //
281  // target coordinates have the format:
282  // dimensions: (# number of target sites) X (dimension)
283  // # of target sites is same as # of rows of neighbor lists
284  //
285  my_GMLS.setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
286 
287  // set up additional sites to evaluate target operators
288  my_GMLS.setAdditionalEvaluationSitesData(additional_target_indices_device, additional_target_coords_device);
289 
290  // create a vector of target operations
291  std::vector<TargetOperation> lro(2);
292  lro[0] = ScalarPointEvaluation;
294 
295  // and then pass them to the GMLS class
296  my_GMLS.addTargets(lro);
297 
298  // sets the weighting kernel function from WeightingFunctionType
299  my_GMLS.setWeightingType(WeightingFunctionType::Power);
300 
301  // power to use in that weighting kernel function
302  my_GMLS.setWeightingPower(2);
303 
304  // generate the alphas that to be combined with data for each target operation requested in lro
305  my_GMLS.generateAlphas();
306 
307 
308  //! [Setting Up The GMLS Object]
309 
310  double instantiation_time = timer.seconds();
311  std::cout << "Took " << instantiation_time << "s to complete alphas generation." << std::endl;
312  Kokkos::fence(); // let generateAlphas finish up before using alphas
313  Kokkos::Profiling::pushRegion("Apply Alphas to Data");
314 
315  //! [Apply GMLS Alphas To Data]
316 
317  // it is important to note that if you expect to use the data as a 1D view, then you should use double*
318  // however, if you know that the target operation will result in a 2D view (vector or matrix output),
319  // then you should template with double** as this is something that can not be infered from the input data
320  // or the target operator at compile time. Additionally, a template argument is required indicating either
321  // Kokkos::HostSpace or Kokkos::DefaultExecutionSpace::memory_space()
322 
323  // The Evaluator class takes care of handling input data views as well as the output data views.
324  // It uses information from the GMLS class to determine how many components are in the input
325  // as well as output for any choice of target functionals and then performs the contactions
326  // on the data using the alpha coefficients generated by the GMLS class, all on the device.
327  Evaluator gmls_evaluator(&my_GMLS);
328 
329  auto output_value1 = gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double*, Kokkos::HostSpace>
330  (sampling_data_device, ScalarPointEvaluation, PointSample,
331  true /*scalar_as_vector_if_needed*/, 1 /*evaluation site index*/);
332 
333  auto output_gradient1 = gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double**, Kokkos::HostSpace>
334  (sampling_data_device, GradientOfScalarPointEvaluation, PointSample,
335  true /*scalar_as_vector_if_needed*/, 1 /*evaluation site index*/);
336 
337  auto output_value2 = gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double*, Kokkos::HostSpace>
338  (sampling_data_device, ScalarPointEvaluation, PointSample,
339  true /*scalar_as_vector_if_needed*/, 2 /*evaluation site index*/);
340 
341  auto output_gradient2 = gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double**, Kokkos::HostSpace>
342  (sampling_data_device, GradientOfScalarPointEvaluation, PointSample,
343  true /*scalar_as_vector_if_needed*/, 2 /*evaluation site index*/);
344 
345  auto output_value3 = gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double*, Kokkos::HostSpace>
346  (sampling_data_device, ScalarPointEvaluation, PointSample,
347  true /*scalar_as_vector_if_needed*/, 3 /*evaluation site index*/);
348 
349  auto output_gradient3 = gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double**, Kokkos::HostSpace>
350  (sampling_data_device, GradientOfScalarPointEvaluation, PointSample,
351  true /*scalar_as_vector_if_needed*/, 3 /*evaluation site index*/);
352 
353  //! [Apply GMLS Alphas To Data]
354 
355  Kokkos::fence(); // let application of alphas to data finish before using results
356  Kokkos::Profiling::popRegion();
357  // times the Comparison in Kokkos
358  Kokkos::Profiling::pushRegion("Comparison");
359 
360  //! [Check That Solutions Are Correct]
361 
362  // load value from output
363  double GMLS_value;
364  // load partial x from gradient
365  double GMLS_GradX;
366  // load partial y from gradient
367  double GMLS_GradY;
368  // load partial z from gradient
369  double GMLS_GradZ;
370 
371  // loop through the target sites
372  extra_evaluation_coordinates_count = 0;
373  for (int i=0; i<number_target_coords; i++) {
374 
375  for (int k=0; k<(i%3)+1; ++k) {
376  if (k==0) {
377  // load value from output
378  GMLS_value = output_value1(i);
379  // load partial x from gradient
380  GMLS_GradX = output_gradient1(i,0);
381  // load partial y from gradient
382  GMLS_GradY = (dimension>1) ? output_gradient1(i,1) : 0;
383  // load partial z from gradient
384  GMLS_GradZ = (dimension>2) ? output_gradient1(i,2) : 0;
385  } else if (k==1) {
386  // load value from output
387  GMLS_value = output_value2(i);
388  // load partial x from gradient
389  GMLS_GradX = output_gradient2(i,0);
390  // load partial y from gradient
391  GMLS_GradY = (dimension>1) ? output_gradient2(i,1) : 0;
392  // load partial z from gradient
393  GMLS_GradZ = (dimension>2) ? output_gradient2(i,2) : 0;
394  } else if (k==2) {
395  // load value from output
396  GMLS_value = output_value3(i);
397  // load partial x from gradient
398  GMLS_GradX = output_gradient3(i,0);
399  // load partial y from gradient
400  GMLS_GradY = (dimension>1) ? output_gradient3(i,1) : 0;
401  // load partial z from gradient
402  GMLS_GradZ = (dimension>2) ? output_gradient3(i,2) : 0;
403  }
404 
405  // target site i's coordinate
406  double xval = additional_target_coords(extra_evaluation_coordinates_count,0);
407  double yval = additional_target_coords(extra_evaluation_coordinates_count,1);
408  double zval = additional_target_coords(extra_evaluation_coordinates_count,2);
409 
410  // evaluation of various exact solutions
411  double actual_value = trueSolution(xval, yval, zval, order, dimension);
412 
413  double actual_Gradient[3] = {0,0,0}; // initialized for 3, but only filled up to dimension
414  trueGradient(actual_Gradient, xval, yval, zval, order, dimension);
415 
416  // check actual function value
417  if(GMLS_value!=GMLS_value || std::abs(actual_value - GMLS_value) > failure_tolerance) {
418  all_passed = false;
419  std::cout << i << " Failed Actual by: " << std::abs(actual_value - GMLS_value) << " for evaluation site: " << k << std::endl;
420  }
421 
422  // check gradient
423  if(std::abs(actual_Gradient[0] - GMLS_GradX) > failure_tolerance) {
424  all_passed = false;
425  std::cout << i << " Failed GradX by: " << std::abs(actual_Gradient[0] - GMLS_GradX) << " for evaluation site: " << k << std::endl;
426  if (dimension>1) {
427  if(std::abs(actual_Gradient[1] - GMLS_GradY) > failure_tolerance) {
428  all_passed = false;
429  std::cout << i << " Failed GradY by: " << std::abs(actual_Gradient[1] - GMLS_GradY) << " for evaluation site: " << k << std::endl;
430  }
431  }
432  if (dimension>2) {
433  if(std::abs(actual_Gradient[2] - GMLS_GradZ) > failure_tolerance) {
434  all_passed = false;
435  std::cout << i << " Failed GradZ by: " << std::abs(actual_Gradient[2] - GMLS_GradZ) << " for evaluation site: " << k << std::endl;
436  }
437  }
438  }
439  extra_evaluation_coordinates_count++;
440  }
441  }
442 
443 
444  //! [Check That Solutions Are Correct]
445  // popRegion hidden from tutorial
446  // stop timing comparison loop
447  Kokkos::Profiling::popRegion();
448  //! [Finalize Program]
449 
450 
451 } // end of code block to reduce scope, causing Kokkos View de-allocations
452 // otherwise, Views may be deallocating when we call Kokkos::finalize() later
453 
454 // finalize Kokkos and MPI (if available)
455 Kokkos::finalize();
456 #ifdef COMPADRE_USE_MPI
457 MPI_Finalize();
458 #endif
459 
460 // output to user that test passed or failed
461 if(all_passed) {
462  fprintf(stdout, "Passed test \n");
463  return 0;
464 } else {
465  fprintf(stdout, "Failed test \n");
466  return -1;
467 }
468 
469 } // main
470 
471 
472 //! [Finalize Program]
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)
int main(int argc, char *args[])
[Parse Command Line Arguments]
Point evaluation of the gradient of a scalar.
Generalized Moving Least Squares (GMLS)
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) ...
constexpr SamplingFunctional PointSample
Available sampling functionals.
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.