Compadre 1.5.5
Loading...
Searching...
No Matches
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>
13
14#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
24using namespace Compadre;
25
26//! [Parse Command Line Arguments]
27
28// called from command line
29int main (int argc, char* args[]) {
30
31// initializes MPI (if available) with command line arguments given
32#ifdef COMPADRE_USE_MPI
33MPI_Init(&argc, &args);
34#endif
35
36// initializes Kokkos with command line arguments given
37Kokkos::initialize(argc, args);
38
39// becomes false if the computed solution not within the failure_threshold of the actual solution
40bool 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.setWeightingParameter(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)
455Kokkos::finalize();
456#ifdef COMPADRE_USE_MPI
457MPI_Finalize();
458#endif
459
460// output to user that test passed or failed
461if(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]
int main(int argc, char *args[])
[Parse Command Line Arguments]
KOKKOS_INLINE_FUNCTION double trueSolution(double x, double y, double z, int order, int dimension)
KOKKOS_INLINE_FUNCTION void trueGradient(double *ans, double x, double y, double z, int order, int dimension)
KOKKOS_INLINE_FUNCTION double divergenceTestSamples(double x, double y, double z, int component, int dimension)
Lightweight Evaluator Helper This class is a lightweight wrapper for extracting and applying all rele...
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)
Generalized Moving Least Squares (GMLS)
void addTargets(TargetOperation lro)
Adds a target to the vector of target functional to be applied to the reconstruction.
void setAdditionalEvaluationSitesData(view_type_1 additional_evaluation_indices, view_type_2 additional_evaluation_coordinates)
(OPTIONAL) Sets additional evaluation sites for each target site
void setWeightingParameter(int wp, int index=0)
Parameter for weighting kernel for GMLS problem index = 0 sets p paramater for weighting kernel index...
void generateAlphas(const int number_of_batches=1, const bool keep_coefficients=false, const bool clear_cache=true)
Meant to calculate target operations and apply the evaluations to the previously constructed polynomi...
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)
void setWeightingType(const std::string &wt)
Type for weighting kernel for GMLS problem.
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....
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...
constexpr SamplingFunctional PointSample
Available sampling functionals.
@ GradientOfScalarPointEvaluation
Point evaluation of the gradient of a scalar.
@ ScalarPointEvaluation
Point evaluation of a scalar.
constexpr SamplingFunctional VectorPointSample
Point evaluations of the entire vector source function.
@ VectorOfScalarClonesTaylorPolynomial
Scalar basis reused as many times as there are components in the vector resulting in a much cheaper p...