Compadre  1.2.0
GMLS_NeumannGradScalar.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 
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 auto kp = KokkosParser(argc, args, true);
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  // check if 8 arguments are given from the command line, the first being the program name
47  int number_of_batches = 1; // 1 batch by default
48  if (argc >= 8) {
49  int arg8toi = atoi(args[7]);
50  if (arg8toi > 0) {
51  number_of_batches = arg8toi;
52  }
53  }
54 
55  // check if 7 arguments are given from the command line, the first being the program name
56  // constraint_type used in solving each GMLS problem:
57  // 0 - No constraints used in solving each GMLS problem
58  // 1 - Neumann Gradient Scalar used in solving each GMLS problem
59  int constraint_type = 1; // Neumann Gradient Scalar by default
60  if (argc >= 7) {
61  int arg7toi = atoi(args[6]);
62  if (arg7toi > 0) {
63  constraint_type = arg7toi;
64  }
65  }
66 
67  // check if 6 arguments are given from the command line, the first being the program name
68  // problem_type used in solving each GMLS problem:
69  // 0 - Standard GMLS problem
70  // 1 - Manifold GMLS problem
71  int problem_type = 0; // Standard by default
72  if (argc >= 6) {
73  int arg6toi = atoi(args[5]);
74  if (arg6toi > 0) {
75  problem_type = arg6toi;
76  }
77  }
78 
79  // check if 5 arguments are given from the command line, the first being the program name
80  // solver_type used for factorization in solving each GMLS problem:
81  // 0 - SVD used for factorization in solving each GMLS problem
82  // 1 - QR used for factorization in solving each GMLS problem
83  // 2 - LU used for factorization in solving each GMLS problem
84  int solver_type = 2; // LU by default
85  if (argc >= 5) {
86  int arg5toi = atoi(args[4]);
87  if (arg5toi >= 0) {
88  solver_type = arg5toi;
89  }
90  }
91 
92  // check if 4 arguments are given from the command line
93  // dimension for the coordinates and the solution
94  int dimension = 3; // dimension 3 by default
95  if (argc >= 4) {
96  int arg4toi = atoi(args[3]);
97  if (arg4toi > 0) {
98  dimension = arg4toi;
99  }
100  }
101  compadre_assert_release((dimension==2 || dimension==3) && "Only supports 2D or 3D.");
102 
103  // check if 3 arguments are given from the command line
104  // set the number of target sites where we will reconstruct the target functionals at
105  int number_target_coords = 200;
106  if (argc >= 3) {
107  int arg3toi = atoi(args[2]);
108  if (arg3toi > 0) {
109  number_target_coords = arg3toi;
110  }
111  }
112 
113  // check if 2 arguments are given from the command line
114  // set the number of target sites where we will reconstruct the target functionals at
115  int order = 3;
116  if (argc >= 2) {
117  int arg2toi = atoi(args[1]);
118  if (arg2toi > 0) {
119  order = arg2toi;
120  }
121  }
122 
123  // the functions we will be seeking to reconstruct are in the span of the basis
124  // of the reconstruction space we choose for GMLS, so the error should be very small
125  const double failure_tolerance = 1e-9;
126 
127  // minimum neighbors for unisolvency is the same as the size of the polynomial basis
128  const int min_neighbors = Compadre::GMLS::getNP(order, dimension);
129 
130  //! [Parse Command Line Arguments]
131  Kokkos::Timer timer;
132  Kokkos::Profiling::pushRegion("Setup Point Data");
133  //! [Setting Up The Point Cloud]
134 
135  // approximate spacing of source sites
136  double h_spacing = 0.05;
137  int n_neg1_to_1 = 2*(1/h_spacing) + 1; // always odd
138 
139  // number of source coordinate sites that will fill a box of [-1,1]x[-1,1]x[-1,1] with a spacing approximately h
140  const int number_source_coords = std::pow(n_neg1_to_1, dimension);
141 
142  // coordinates of source sites
143  Kokkos::View<double**, Kokkos::DefaultExecutionSpace> source_coords_device("source coordinates",
144  number_source_coords, 3);
145  Kokkos::View<double**>::HostMirror source_coords = Kokkos::create_mirror_view(source_coords_device);
146 
147  // coordinates of target sites
148  Kokkos::View<double**, Kokkos::DefaultExecutionSpace> target_coords_device ("target coordinates", number_target_coords, 3);
149  Kokkos::View<double**>::HostMirror target_coords = Kokkos::create_mirror_view(target_coords_device);
150 
151  // tangent bundle for each target sites
152  Kokkos::View<double***, Kokkos::DefaultExecutionSpace> tangent_bundles_device ("tangent bundles", number_target_coords, dimension, dimension);
153  Kokkos::View<double***>::HostMirror tangent_bundles = Kokkos::create_mirror_view(tangent_bundles_device);
154 
155  // fill source coordinates with a uniform grid
156  int source_index = 0;
157  double this_coord[3] = {0,0,0};
158  for (int i=-n_neg1_to_1/2; i<n_neg1_to_1/2+1; ++i) {
159  this_coord[0] = i*h_spacing;
160  for (int j=-n_neg1_to_1/2; j<n_neg1_to_1/2+1; ++j) {
161  this_coord[1] = j*h_spacing;
162  for (int k=-n_neg1_to_1/2; k<n_neg1_to_1/2+1; ++k) {
163  this_coord[2] = k*h_spacing;
164  if (dimension==3) {
165  source_coords(source_index,0) = this_coord[0];
166  source_coords(source_index,1) = this_coord[1];
167  source_coords(source_index,2) = this_coord[2];
168  source_index++;
169  }
170  }
171  if (dimension==2) {
172  source_coords(source_index,0) = this_coord[0];
173  source_coords(source_index,1) = this_coord[1];
174  source_coords(source_index,2) = 0;
175  source_index++;
176  }
177  }
178  if (dimension==1) {
179  source_coords(source_index,0) = this_coord[0];
180  source_coords(source_index,1) = 0;
181  source_coords(source_index,2) = 0;
182  source_index++;
183  }
184  }
185 
186  // fill target coords somewhere inside of [-0.5,0.5]x[-0.5,0.5]x[-0.5,0.5]
187  for(int i=0; i<number_target_coords; i++){
188 
189  // first, we get a uniformly random distributed direction
190  double rand_dir[3] = {0,0,0};
191 
192  for (int j=0; j<dimension; ++j) {
193  // rand_dir[j] is in [-0.5, 0.5]
194  rand_dir[j] = ((double)rand() / (double) RAND_MAX) - 0.5;
195  }
196 
197  // then we get a uniformly random radius
198  for (int j=0; j<dimension; ++j) {
199  target_coords(i,j) = rand_dir[j];
200  }
201  // target_coords(i, 2) = 1.0;
202 
203  // Set tangent bundles
204  if (dimension == 2) {
205  tangent_bundles(i, 0, 0) = 0.0;
206  tangent_bundles(i, 0, 1) = 0.0;
207  tangent_bundles(i, 1, 0) = 1.0/(sqrt(2.0));
208  tangent_bundles(i, 1, 1) = 1.0/(sqrt(2.0));
209  } else if (dimension == 3) {
210  tangent_bundles(i, 0, 0) = 0.0;
211  tangent_bundles(i, 0, 1) = 0.0;
212  tangent_bundles(i, 0, 2) = 0.0;
213  tangent_bundles(i, 1, 0) = 0.0;
214  tangent_bundles(i, 1, 1) = 0.0;
215  tangent_bundles(i, 1, 2) = 0.0;
216  tangent_bundles(i, 2, 0) = 1.0/(sqrt(3.0));
217  tangent_bundles(i, 2, 1) = 1.0/(sqrt(3.0));
218  tangent_bundles(i, 2, 2) = 1.0/(sqrt(3.0));
219  }
220  }
221 
222  //! [Setting Up The Point Cloud]
223 
224  Kokkos::Profiling::popRegion();
225  Kokkos::Profiling::pushRegion("Creating Data");
226 
227  //! [Creating The Data]
228 
229 
230  // source coordinates need copied to device before using to construct sampling data
231  Kokkos::deep_copy(source_coords_device, source_coords);
232 
233  // target coordinates copied next, because it is a convenient time to send them to device
234  Kokkos::deep_copy(target_coords_device, target_coords);
235 
236  // tangent bundles copied next, because it is a convenient time to send them to device
237  Kokkos::deep_copy(tangent_bundles_device, tangent_bundles);
238 
239  // need Kokkos View storing true solution
240  Kokkos::View<double*, Kokkos::DefaultExecutionSpace> sampling_data_device("samples of true solution",
241  source_coords_device.extent(0));
242 
243  Kokkos::parallel_for("Sampling Manufactured Solutions", Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>
244  (0,source_coords.extent(0)), KOKKOS_LAMBDA(const int i) {
245 
246  // coordinates of source site i
247  double xval = source_coords_device(i,0);
248  double yval = (dimension>1) ? source_coords_device(i,1) : 0;
249  double zval = (dimension>2) ? source_coords_device(i,2) : 0;
250 
251  // data for targets with scalar input
252  sampling_data_device(i) = trueSolution(xval, yval, zval, order, dimension);
253  });
254 
255  //! [Creating The Data]
256 
257  Kokkos::Profiling::popRegion();
258  Kokkos::Profiling::pushRegion("Neighbor Search");
259 
260  //! [Performing Neighbor Search]
261 
262 
263  // Point cloud construction for neighbor search
264  // CreatePointCloudSearch constructs an object of type PointCloudSearch, but deduces the templates for you
265  auto point_cloud_search(CreatePointCloudSearch(source_coords, dimension));
266 
267  // each row is a neighbor list for a target site, with the first column of each row containing
268  // the number of neighbors for that rows corresponding target site
269  double epsilon_multiplier = 1.8;
270  int estimated_upper_bound_number_neighbors =
271  point_cloud_search.getEstimatedNumberNeighborsUpperBound(min_neighbors, dimension, epsilon_multiplier);
272 
273  Kokkos::View<int**, Kokkos::DefaultExecutionSpace> neighbor_lists_device("neighbor lists",
274  number_target_coords, estimated_upper_bound_number_neighbors); // first column is # of neighbors
275  Kokkos::View<int**>::HostMirror neighbor_lists = Kokkos::create_mirror_view(neighbor_lists_device);
276 
277  // each target site has a window size
278  Kokkos::View<double*, Kokkos::DefaultExecutionSpace> epsilon_device("h supports", number_target_coords);
279  Kokkos::View<double*>::HostMirror epsilon = Kokkos::create_mirror_view(epsilon_device);
280 
281  // query the point cloud to generate the neighbor lists using a kdtree to produce the n nearest neighbor
282  // to each target site, adding (epsilon_multiplier-1)*100% to whatever the distance away the further neighbor used is from
283  // each target to the view for epsilon
284  point_cloud_search.generate2DNeighborListsFromKNNSearch(false /*not dry run*/, target_coords, neighbor_lists,
285  epsilon, min_neighbors, epsilon_multiplier);
286 
287  //! [Performing Neighbor Search]
288 
289  Kokkos::Profiling::popRegion();
290  Kokkos::fence(); // let call to build neighbor lists complete before copying back to device
291  timer.reset();
292 
293  //! [Setting Up The GMLS Object]
294 
295 
296  // Copy data back to device (they were filled on the host)
297  // We could have filled Kokkos Views with memory space on the host
298  // and used these instead, and then the copying of data to the device
299  // would be performed in the GMLS class
300  Kokkos::deep_copy(neighbor_lists_device, neighbor_lists);
301  Kokkos::deep_copy(epsilon_device, epsilon);
302 
303  // solver name for passing into the GMLS class
304  std::string solver_name;
305  if (solver_type == 0) { // SVD
306  solver_name = "SVD";
307  } else if (solver_type == 1) { // QR
308  solver_name = "QR";
309  } else if (solver_type == 2) { // LU
310  solver_name = "LU";
311  }
312 
313  // problem name for passing into the GMLS class
314  std::string problem_name;
315  if (problem_type == 0) { // Standard
316  problem_name = "STANDARD";
317  } else if (problem_type == 1) { // Manifold
318  problem_name = "MANIFOLD";
319  }
320 
321  // boundary name for passing into the GMLS class
322  std::string constraint_name;
323  if (constraint_type == 0) { // No constraints
324  constraint_name = "NO_CONSTRAINT";
325  } else if (constraint_type == 1) { // Neumann Gradient Scalar
326  constraint_name = "NEUMANN_GRAD_SCALAR";
327  }
328 
329  // initialize an instance of the GMLS class
331  PointSample,
332  order, dimension,
333  solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
334  0 /*manifold order*/);
335 
336  // pass in neighbor lists, source coordinates, target coordinates, and window sizes
337  //
338  // neighbor lists have the format:
339  // dimensions: (# number of target sites) X (# maximum number of neighbors for any given target + 1)
340  // the first column contains the number of neighbors for that rows corresponding target index
341  //
342  // source coordinates have the format:
343  // dimensions: (# number of source sites) X (dimension)
344  // entries in the neighbor lists (integers) correspond to rows of this 2D array
345  //
346  // target coordinates have the format:
347  // dimensions: (# number of target sites) X (dimension)
348  // # of target sites is same as # of rows of neighbor lists
349  //
350  my_GMLS.setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
351  my_GMLS.setTangentBundle(tangent_bundles_device);
352 
353  // create a vector of target operations
354  TargetOperation lro;
356 
357  // and then pass them to the GMLS class
358  my_GMLS.addTargets(lro);
359 
360  // sets the weighting kernel function from WeightingFunctionType
361  my_GMLS.setWeightingType(WeightingFunctionType::Power);
362 
363  // power to use in that weighting kernel function
364  my_GMLS.setWeightingPower(2);
365 
366  // generate the alphas that to be combined with data for each target operation requested in lro
367  my_GMLS.generateAlphas(number_of_batches);
368 
369  //! [Setting Up The GMLS Object]
370 
371  double instantiation_time = timer.seconds();
372  std::cout << "Took " << instantiation_time << "s to complete normal vectors generation." << std::endl;
373  Kokkos::fence(); // let generateNormalVectors finish up before using alphas
374  Kokkos::Profiling::pushRegion("Apply Alphas to Data");
375 
376  //! [Apply GMLS Alphas To Data]
377 
378  // it is important to note that if you expect to use the data as a 1D view, then you should use double*
379  // however, if you know that the target operation will result in a 2D view (vector or matrix output),
380  // then you should template with double** as this is something that can not be infered from the input data
381  // or the target operator at compile time. Additionally, a template argument is required indicating either
382  // Kokkos::HostSpace or Kokkos::DefaultExecutionSpace::memory_space()
383 
384  // The Evaluator class takes care of handling input data views as well as the output data views.
385  // It uses information from the GMLS class to determine how many components are in the input
386  // as well as output for any choice of target functionals and then performs the contactions
387  // on the data using the alpha coefficients generated by the GMLS class, all on the device.
388  Evaluator gmls_evaluator(&my_GMLS);
389 
390  auto output_value = gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double*, Kokkos::HostSpace>
391  (sampling_data_device, LaplacianOfScalarPointEvaluation);
392 
393  Kokkos::fence(); // let application of alphas to data finish before using results
394  Kokkos::Profiling::popRegion();
395  // times the Comparison in Kokkos
396  Kokkos::Profiling::pushRegion("Comparison");
397 
398  //! [Check That Solutions Are Correct]
399 
400  // loop through the target sites
401  for (int i=0; i<number_target_coords; i++) {
402  // target site i's coordinate
403  double xval = target_coords(i,0);
404  double yval = (dimension>1) ? target_coords(i,1) : 0;
405  double zval = (dimension>2) ? target_coords(i,2) : 0;
406 
407  // 0th entry is # of neighbors, which is the index beyond the last neighbor
408  int num_neigh_i = neighbor_lists(i, 0);
409  double b_i = my_GMLS.getAlpha0TensorTo0Tensor(lro, i, num_neigh_i);
410 
411  // load value from output
412  double GMLS_value = output_value(i);
413 
414  // obtain the real Laplacian
415  double actual_Laplacian = trueLaplacian(xval, yval, zval, order, dimension);
416 
417  // calculate value of g to reconstruct the computed Laplacian
418  double actual_Gradient[3] = {0,0,0}; // initialized for 3, but only filled up to dimension
419  trueGradient(actual_Gradient, xval, yval, zval, order, dimension);
420  double g = (dimension == 3) ? (1.0/sqrt(3.0))*(actual_Gradient[0] + actual_Gradient[1] + actual_Gradient[2])
421  : (1.0/sqrt(2.0))*(actual_Gradient[0] + actual_Gradient[1]);
422  double adjusted_value = GMLS_value + b_i*g;
423 
424  // check actual function value
425  if(GMLS_value!=GMLS_value || std::abs(actual_Laplacian - adjusted_value) > failure_tolerance) {
426  all_passed = false;
427  std::cout << i << " Failed Actual by: " << std::abs(actual_Laplacian - adjusted_value) << std::endl;
428  }
429  }
430 
431  //! [Check That Solutions Are Correct]
432  // popRegion hidden from tutorial
433  // stop timing comparison loop
434  Kokkos::Profiling::popRegion();
435  //! [Finalize Program]
436 } // end of code block to reduce scope, causing Kokkos View de-allocations
437 // otherwise, Views may be deallocating when we call Kokkos finalize() later
438 
439 // finalize Kokkos and MPI (if available)
440 kp.finalize();
441 #ifdef COMPADRE_USE_MPI
442 MPI_Finalize();
443 #endif
444 
445 // output to user that test passed or failed
446 if(all_passed) {
447  fprintf(stdout, "Passed test \n");
448  return 0;
449 } else {
450  fprintf(stdout, "Failed test \n");
451  return -1;
452 }
453 
454 } // main
455 
456 
457 //! [Finalize Program]
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...
Scalar polynomial basis centered at the target site and scaled by sum of basis powers e...
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)
Kokkos::View< output_data_type, output_array_layout, output_memory_space > applyAlphasToDataAllComponentsAllTargetSites(view_type_input_data sampling_data, TargetOperation lro, const SamplingFunctional sro_in=PointSample, bool scalar_as_vector_if_needed=true, const int evaluation_site_local_index=0) const
Transformation of data under GMLS.
int main(int argc, char *args[])
[Parse Command Line Arguments]
Point evaluation of the laplacian of a scalar (could be on a manifold or not)
KOKKOS_INLINE_FUNCTION double trueLaplacian(double x, double y, double z, int order, int dimension)
Generalized Moving Least Squares (GMLS)
void setProblemData(view_type_1 neighbor_lists, view_type_2 source_coordinates, view_type_3 target_coordinates, view_type_4 epsilons)
Sets basic problem data (neighbor lists, source coordinates, and target coordinates) ...
constexpr SamplingFunctional PointSample
Available sampling functionals.
KOKKOS_INLINE_FUNCTION double trueSolution(double x, double y, double z, int order, int dimension)
#define compadre_assert_release(condition)
compadre_assert_release is used for assertions that should always be checked, but generally are not e...
TargetOperation
Available target functionals.