Compadre  1.3.3
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 #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 
56  // the functions we will be seeking to reconstruct are in the span of the basis
57  // of the reconstruction space we choose for GMLS, so the error should be very small
58  const double failure_tolerance = 1e-9;
59 
60  // minimum neighbors for unisolvency is the same as the size of the polynomial basis
61  const int min_neighbors = Compadre::GMLS::getNP(order, dimension);
62 
63  //! [Parse Command Line Arguments]
64  Kokkos::Timer timer;
65  Kokkos::Profiling::pushRegion("Setup Point Data");
66  //! [Setting Up The Point Cloud]
67 
68  // approximate spacing of source sites
69  double h_spacing = 0.05;
70  int n_neg1_to_1 = 2*(1/h_spacing) + 1; // always odd
71 
72  // number of source coordinate sites that will fill a box of [-1,1]x[-1,1]x[-1,1] with a spacing approximately h
73  const int number_source_coords = std::pow(n_neg1_to_1, dimension);
74 
75  // coordinates of source sites
76  Kokkos::View<double**, Kokkos::DefaultExecutionSpace> source_coords_device("source coordinates",
77  number_source_coords, 3);
78  Kokkos::View<double**>::HostMirror source_coords = Kokkos::create_mirror_view(source_coords_device);
79 
80  // coordinates of target sites
81  Kokkos::View<double**, Kokkos::DefaultExecutionSpace> target_coords_device ("target coordinates", number_target_coords, 3);
82  Kokkos::View<double**>::HostMirror target_coords = Kokkos::create_mirror_view(target_coords_device);
83 
84  // tangent bundle for each target sites
85  Kokkos::View<double***, Kokkos::DefaultExecutionSpace> tangent_bundles_device ("tangent bundles", number_target_coords, dimension, dimension);
86  Kokkos::View<double***>::HostMirror tangent_bundles = Kokkos::create_mirror_view(tangent_bundles_device);
87 
88  // fill source coordinates with a uniform grid
89  int source_index = 0;
90  double this_coord[3] = {0,0,0};
91  for (int i=-n_neg1_to_1/2; i<n_neg1_to_1/2+1; ++i) {
92  this_coord[0] = i*h_spacing;
93  for (int j=-n_neg1_to_1/2; j<n_neg1_to_1/2+1; ++j) {
94  this_coord[1] = j*h_spacing;
95  for (int k=-n_neg1_to_1/2; k<n_neg1_to_1/2+1; ++k) {
96  this_coord[2] = k*h_spacing;
97  if (dimension==3) {
98  source_coords(source_index,0) = this_coord[0];
99  source_coords(source_index,1) = this_coord[1];
100  source_coords(source_index,2) = this_coord[2];
101  source_index++;
102  }
103  }
104  if (dimension==2) {
105  source_coords(source_index,0) = this_coord[0];
106  source_coords(source_index,1) = this_coord[1];
107  source_coords(source_index,2) = 0;
108  source_index++;
109  }
110  }
111  if (dimension==1) {
112  source_coords(source_index,0) = this_coord[0];
113  source_coords(source_index,1) = 0;
114  source_coords(source_index,2) = 0;
115  source_index++;
116  }
117  }
118 
119  // fill target coords somewhere inside of [-0.5,0.5]x[-0.5,0.5]x[-0.5,0.5]
120  for(int i=0; i<number_target_coords; i++){
121 
122  // first, we get a uniformly random distributed direction
123  double rand_dir[3] = {0,0,0};
124 
125  for (int j=0; j<dimension; ++j) {
126  // rand_dir[j] is in [-0.5, 0.5]
127  rand_dir[j] = ((double)rand() / (double) RAND_MAX) - 0.5;
128  }
129 
130  // then we get a uniformly random radius
131  for (int j=0; j<dimension; ++j) {
132  target_coords(i,j) = rand_dir[j];
133  }
134  // target_coords(i, 2) = 1.0;
135 
136  // Set tangent bundles
137  if (dimension == 2) {
138  tangent_bundles(i, 0, 0) = 0.0;
139  tangent_bundles(i, 0, 1) = 0.0;
140  tangent_bundles(i, 1, 0) = 1.0/(sqrt(2.0));
141  tangent_bundles(i, 1, 1) = 1.0/(sqrt(2.0));
142  } else if (dimension == 3) {
143  tangent_bundles(i, 0, 0) = 0.0;
144  tangent_bundles(i, 0, 1) = 0.0;
145  tangent_bundles(i, 0, 2) = 0.0;
146  tangent_bundles(i, 1, 0) = 0.0;
147  tangent_bundles(i, 1, 1) = 0.0;
148  tangent_bundles(i, 1, 2) = 0.0;
149  tangent_bundles(i, 2, 0) = 1.0/(sqrt(3.0));
150  tangent_bundles(i, 2, 1) = 1.0/(sqrt(3.0));
151  tangent_bundles(i, 2, 2) = 1.0/(sqrt(3.0));
152  }
153  }
154 
155  //! [Setting Up The Point Cloud]
156 
157  Kokkos::Profiling::popRegion();
158  Kokkos::Profiling::pushRegion("Creating Data");
159 
160  //! [Creating The Data]
161 
162 
163  // source coordinates need copied to device before using to construct sampling data
164  Kokkos::deep_copy(source_coords_device, source_coords);
165 
166  // target coordinates copied next, because it is a convenient time to send them to device
167  Kokkos::deep_copy(target_coords_device, target_coords);
168 
169  // tangent bundles copied next, because it is a convenient time to send them to device
170  Kokkos::deep_copy(tangent_bundles_device, tangent_bundles);
171 
172  // need Kokkos View storing true solution
173  Kokkos::View<double*, Kokkos::DefaultExecutionSpace> sampling_data_device("samples of true solution",
174  source_coords_device.extent(0));
175 
176  Kokkos::parallel_for("Sampling Manufactured Solutions", Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>
177  (0,source_coords.extent(0)), KOKKOS_LAMBDA(const int i) {
178 
179  // coordinates of source site i
180  double xval = source_coords_device(i,0);
181  double yval = (dimension>1) ? source_coords_device(i,1) : 0;
182  double zval = (dimension>2) ? source_coords_device(i,2) : 0;
183 
184  // data for targets with scalar input
185  sampling_data_device(i) = trueSolution(xval, yval, zval, order, dimension);
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  // each row is a neighbor list for a target site, with the first column of each row containing
201  // the number of neighbors for that rows corresponding target site
202  double epsilon_multiplier = 1.8;
203  int estimated_upper_bound_number_neighbors =
204  point_cloud_search.getEstimatedNumberNeighborsUpperBound(min_neighbors, dimension, epsilon_multiplier);
205 
206  Kokkos::View<int**, Kokkos::DefaultExecutionSpace> neighbor_lists_device("neighbor lists",
207  number_target_coords, estimated_upper_bound_number_neighbors); // first column is # of neighbors
208  Kokkos::View<int**>::HostMirror neighbor_lists = Kokkos::create_mirror_view(neighbor_lists_device);
209 
210  // each target site has a window size
211  Kokkos::View<double*, Kokkos::DefaultExecutionSpace> epsilon_device("h supports", number_target_coords);
212  Kokkos::View<double*>::HostMirror epsilon = Kokkos::create_mirror_view(epsilon_device);
213 
214  // query the point cloud to generate the neighbor lists using a kdtree to produce the n nearest neighbor
215  // to each target site, adding (epsilon_multiplier-1)*100% to whatever the distance away the further neighbor used is from
216  // each target to the view for epsilon
217  point_cloud_search.generate2DNeighborListsFromKNNSearch(false /*not dry run*/, target_coords, neighbor_lists,
218  epsilon, min_neighbors, epsilon_multiplier);
219 
220  //! [Performing Neighbor Search]
221 
222  Kokkos::Profiling::popRegion();
223  Kokkos::fence(); // let call to build neighbor lists complete before copying back to device
224  timer.reset();
225 
226  //! [Setting Up The GMLS Object]
227 
228 
229  // Copy data back to device (they were filled on the host)
230  // We could have filled Kokkos Views with memory space on the host
231  // and used these instead, and then the copying of data to the device
232  // would be performed in the GMLS class
233  Kokkos::deep_copy(neighbor_lists_device, neighbor_lists);
234  Kokkos::deep_copy(epsilon_device, epsilon);
235 
236  // initialize an instance of the GMLS class
238  PointSample,
239  order, dimension,
240  solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
241  0 /*manifold order*/);
242 
243  // pass in neighbor lists, source coordinates, target coordinates, and window sizes
244  //
245  // neighbor lists have the format:
246  // dimensions: (# number of target sites) X (# maximum number of neighbors for any given target + 1)
247  // the first column contains the number of neighbors for that rows corresponding target index
248  //
249  // source coordinates have the format:
250  // dimensions: (# number of source sites) X (dimension)
251  // entries in the neighbor lists (integers) correspond to rows of this 2D array
252  //
253  // target coordinates have the format:
254  // dimensions: (# number of target sites) X (dimension)
255  // # of target sites is same as # of rows of neighbor lists
256  //
257  my_GMLS.setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
258  my_GMLS.setTangentBundle(tangent_bundles_device);
259 
260  // create a vector of target operations
261  TargetOperation lro;
263 
264  // and then pass them to the GMLS class
265  my_GMLS.addTargets(lro);
266 
267  // sets the weighting kernel function from WeightingFunctionType
268  my_GMLS.setWeightingType(WeightingFunctionType::Power);
269 
270  // power to use in that weighting kernel function
271  my_GMLS.setWeightingPower(2);
272 
273  // generate the alphas that to be combined with data for each target operation requested in lro
274  my_GMLS.generateAlphas(number_of_batches);
275 
276  //! [Setting Up The GMLS Object]
277 
278  double instantiation_time = timer.seconds();
279  std::cout << "Took " << instantiation_time << "s to complete normal vectors generation." << std::endl;
280  Kokkos::fence(); // let generateNormalVectors finish up before using alphas
281  Kokkos::Profiling::pushRegion("Apply Alphas to Data");
282 
283  //! [Apply GMLS Alphas To Data]
284 
285  // it is important to note that if you expect to use the data as a 1D view, then you should use double*
286  // however, if you know that the target operation will result in a 2D view (vector or matrix output),
287  // then you should template with double** as this is something that can not be infered from the input data
288  // or the target operator at compile time. Additionally, a template argument is required indicating either
289  // Kokkos::HostSpace or Kokkos::DefaultExecutionSpace::memory_space()
290 
291  // The Evaluator class takes care of handling input data views as well as the output data views.
292  // It uses information from the GMLS class to determine how many components are in the input
293  // as well as output for any choice of target functionals and then performs the contactions
294  // on the data using the alpha coefficients generated by the GMLS class, all on the device.
295  Evaluator gmls_evaluator(&my_GMLS);
296 
297  auto output_value = gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double*, Kokkos::HostSpace>
298  (sampling_data_device, LaplacianOfScalarPointEvaluation);
299 
300  Kokkos::fence(); // let application of alphas to data finish before using results
301  Kokkos::Profiling::popRegion();
302  // times the Comparison in Kokkos
303  Kokkos::Profiling::pushRegion("Comparison");
304 
305  //! [Check That Solutions Are Correct]
306 
307  // loop through the target sites
308  for (int i=0; i<number_target_coords; i++) {
309  // target site i's coordinate
310  double xval = target_coords(i,0);
311  double yval = (dimension>1) ? target_coords(i,1) : 0;
312  double zval = (dimension>2) ? target_coords(i,2) : 0;
313 
314  // 0th entry is # of neighbors, which is the index beyond the last neighbor
315  int num_neigh_i = neighbor_lists(i, 0);
316  double b_i = my_GMLS.getAlpha0TensorTo0Tensor(lro, i, num_neigh_i);
317 
318  // load value from output
319  double GMLS_value = output_value(i);
320 
321  // obtain the real Laplacian
322  double actual_Laplacian = trueLaplacian(xval, yval, zval, order, dimension);
323 
324  // calculate value of g to reconstruct the computed Laplacian
325  double actual_Gradient[3] = {0,0,0}; // initialized for 3, but only filled up to dimension
326  trueGradient(actual_Gradient, xval, yval, zval, order, dimension);
327  double g = (dimension == 3) ? (1.0/sqrt(3.0))*(actual_Gradient[0] + actual_Gradient[1] + actual_Gradient[2])
328  : (1.0/sqrt(2.0))*(actual_Gradient[0] + actual_Gradient[1]);
329  double adjusted_value = GMLS_value + b_i*g;
330 
331  // check actual function value
332  if(GMLS_value!=GMLS_value || std::abs(actual_Laplacian - adjusted_value) > failure_tolerance) {
333  all_passed = false;
334  std::cout << i << " Failed Actual by: " << std::abs(actual_Laplacian - adjusted_value) << std::endl;
335  }
336  }
337 
338  //! [Check That Solutions Are Correct]
339  // popRegion hidden from tutorial
340  // stop timing comparison loop
341  Kokkos::Profiling::popRegion();
342  //! [Finalize Program]
343 } // end of code block to reduce scope, causing Kokkos View de-allocations
344 // otherwise, Views may be deallocating when we call Kokkos finalize() later
345 
346 // finalize Kokkos and MPI (if available)
347 kp.finalize();
348 #ifdef COMPADRE_USE_MPI
349 MPI_Finalize();
350 #endif
351 
352 // output to user that test passed or failed
353 if(all_passed) {
354  fprintf(stdout, "Passed test \n");
355  return 0;
356 } else {
357  fprintf(stdout, "Failed test \n");
358  return -1;
359 }
360 
361 } // main
362 
363 
364 //! [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 (allocates memory for output)
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)
TargetOperation
Available target functionals.