Compadre 1.5.5
Loading...
Searching...
No Matches
GMLS_DivergenceFree.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// 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 CommandLineProcessor clp(argc, args);
47 auto order = clp.order;
48 auto dimension = clp.dimension;
49 auto number_target_coords = clp.number_target_coords;
50 auto constraint_name = clp.constraint_name;
51 auto solver_name = clp.solver_name;
52 auto problem_name = clp.problem_name;
53
54 // the functions we will be seeking to reconstruct are in the span of the basis
55 // of the reconstruction space we choose for GMLS, so the error should be very small
56 const double failure_tolerance = 1e-9;
57
58 // minimum neighbors for unisolvency is the same as the size of the polynomial basis
59 const int min_neighbors = Compadre::GMLS::getNP(order, dimension, DivergenceFreeVectorTaylorPolynomial);
60
61 //! [Parse Command Line Arguments]
62 Kokkos::Timer timer;
63 Kokkos::Profiling::pushRegion("Setup Point Data");
64 //! [Setting Up The Point Cloud]
65
66 // approximate spacing of source sites
67 double h_spacing = 0.05;
68 int n_neg1_to_1 = 2*(1/h_spacing) + 1; // always odd
69
70 // number of source coordinate sites that will fill a box of [-1,1]x[-1,1]x[-1,1] with a spacing approximately h
71 const int number_source_coords = std::pow(n_neg1_to_1, dimension);
72
73 // coordinates of source sites
74 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> source_coords_device("source coordinates",
75 number_source_coords, 3);
76 Kokkos::View<double**>::HostMirror source_coords = Kokkos::create_mirror_view(source_coords_device);
77
78 // coordinates of target sites
79 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> target_coords_device ("target coordinates", number_target_coords, 3);
80 Kokkos::View<double**>::HostMirror target_coords = Kokkos::create_mirror_view(target_coords_device);
81
82
83 // fill source coordinates with a uniform grid
84 int source_index = 0;
85 double this_coord[3] = {0,0,0};
86 for (int i=-n_neg1_to_1/2; i<n_neg1_to_1/2+1; ++i) {
87 this_coord[0] = i*h_spacing;
88 for (int j=-n_neg1_to_1/2; j<n_neg1_to_1/2+1; ++j) {
89 this_coord[1] = j*h_spacing;
90 for (int k=-n_neg1_to_1/2; k<n_neg1_to_1/2+1; ++k) {
91 this_coord[2] = k*h_spacing;
92 if (dimension==3) {
93 source_coords(source_index,0) = this_coord[0];
94 source_coords(source_index,1) = this_coord[1];
95 source_coords(source_index,2) = this_coord[2];
96 source_index++;
97 }
98 }
99 if (dimension==2) {
100 source_coords(source_index,0) = this_coord[0];
101 source_coords(source_index,1) = this_coord[1];
102 source_coords(source_index,2) = 0;
103 source_index++;
104 }
105 }
106 if (dimension==1) {
107 source_coords(source_index,0) = this_coord[0];
108 source_coords(source_index,1) = 0;
109 source_coords(source_index,2) = 0;
110 source_index++;
111 }
112 }
113
114 // Generate target points - these are random permutation from available source points
115 // Note that this is assuming that the number of targets in this test will not exceed
116 // the number of source coords, which is 41^3 = 68921
117 // seed random number generator
118 std::mt19937 rng(50);
119 // generate random integers in [0..number_source_coords-1] (used to pick target sites)
120 std::uniform_int_distribution<int> gen_num_neighbours(0, source_index);
121 // fill target sites with random selections from source sites
122 for (int i=0; i<number_target_coords; ++i) {
123 const int source_site_to_copy = gen_num_neighbours(rng);
124 for (int j=0; j<3; j++) {
125 target_coords(i, j) = source_coords(source_site_to_copy, j);
126 }
127 }
128
129 //! [Setting Up The Point Cloud]
130
131 Kokkos::Profiling::popRegion();
132 Kokkos::Profiling::pushRegion("Creating Data");
133
134 //! [Creating The Data]
135
136
137 // source coordinates need copied to device before using to construct sampling data
138 Kokkos::deep_copy(source_coords_device, source_coords);
139
140 // target coordinates copied next, because it is a convenient time to send them to device
141 Kokkos::deep_copy(target_coords_device, target_coords);
142
143 // need Kokkos View storing true solution
144 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> vector_sampling_span_basis_data_device("samples of true vector solution",
145 source_coords_device.extent(0), dimension);
146 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> vector_sampling_single_polynomial_data_device("samples of true vector solution",
147 source_coords_device.extent(0), dimension);
148
149 Kokkos::parallel_for("Sampling Manufactured Solutions", Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>
150 (0,source_coords.extent(0)), KOKKOS_LAMBDA(const int i) {
151 // coordinates of source site i
152 double xval = source_coords_device(i,0);
153 double yval = (dimension>1) ? source_coords_device(i,1) : 0;
154 double zval = (dimension>2) ? source_coords_device(i,2) : 0;
155
156 // data for targets with vector input
157 for (int j=0; j<dimension; ++j) {
158 vector_sampling_span_basis_data_device(i, j) = divfreeTestSolution_span_basis(xval, yval, zval, j, dimension, order);
159 vector_sampling_single_polynomial_data_device(i, j) = divfreeTestSolution_single_polynomial(xval, yval, zval, j, dimension);
160 }
161 });
162
163
164 //! [Creating The Data]
165
166 Kokkos::Profiling::popRegion();
167 Kokkos::Profiling::pushRegion("Neighbor Search");
168
169 //! [Performing Neighbor Search]
170
171
172 // Point cloud construction for neighbor search
173 // CreatePointCloudSearch constructs an object of type PointCloudSearch, but deduces the templates for you
174 auto point_cloud_search(CreatePointCloudSearch(source_coords, dimension));
175
176 // each row is a neighbor list for a target site, with the first column of each row containing
177 // the number of neighbors for that rows corresponding target site
178 // for the default values in this test, the multiplier is suggested to be 2.2
179 double epsilon_multiplier = 2.2;
180 int estimated_upper_bound_number_neighbors =
181 point_cloud_search.getEstimatedNumberNeighborsUpperBound(min_neighbors, dimension, epsilon_multiplier);
182
183 Kokkos::View<int**, Kokkos::DefaultExecutionSpace> neighbor_lists_device("neighbor lists",
184 number_target_coords, estimated_upper_bound_number_neighbors); // first column is # of neighbors
185 Kokkos::View<int**>::HostMirror neighbor_lists = Kokkos::create_mirror_view(neighbor_lists_device);
186
187 // each target site has a window size
188 Kokkos::View<double*, Kokkos::DefaultExecutionSpace> epsilon_device("h supports", number_target_coords);
189 Kokkos::View<double*>::HostMirror epsilon = Kokkos::create_mirror_view(epsilon_device);
190
191 // query the point cloud to generate the neighbor lists using a kdtree to produce the n nearest neighbor
192 // to each target site, adding (epsilon_multiplier-1)*100% to whatever the distance away the further neighbor used is from
193 // each target to the view for epsilon
194 point_cloud_search.generate2DNeighborListsFromKNNSearch(false /*not dry run*/, target_coords, neighbor_lists,
195 epsilon, min_neighbors, epsilon_multiplier);
196
197 //! [Performing Neighbor Search]
198
199 Kokkos::Profiling::popRegion();
200 Kokkos::fence(); // let call to build neighbor lists complete before copying back to device
201 timer.reset();
202
203 //! [Setting Up The GMLS Object]
204
205
206 // Copy data back to device (they were filled on the host)
207 // We could have filled Kokkos Views with memory space on the host
208 // and used these instead, and then the copying of data to the device
209 // would be performed in the GMLS class
210 Kokkos::deep_copy(neighbor_lists_device, neighbor_lists);
211 Kokkos::deep_copy(epsilon_device, epsilon);
212
213 // initialize an instance of the GMLS class
214 // NULL in manifold order indicates non-manifold case
215 // Vector basis to perform GMLS on divergence free basis
216 GMLS vector_divfree_basis_gmls(DivergenceFreeVectorTaylorPolynomial,
218 order, dimension,
219 solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
220 0 /*manifold order*/);
221
222 // pass in neighbor lists, source coordinates, target coordinates, and window sizes
223 //
224 // neighbor lists have the format:
225 // dimensions: (# number of target sites) X (# maximum number of neighbors for any given target + 1)
226 // the first column contains the number of neighbors for that rows corresponding target index
227 //
228 // source coordinates have the format:
229 // dimensions: (# number of source sites) X (dimension)
230 // entries in the neighbor lists (integers) correspond to rows of this 2D array
231 //
232 // target coordinates have the format:
233 // dimensions: (# number of target sites) X (dimension)
234 // # of target sites is same as # of rows of neighbor lists
235 //
236 vector_divfree_basis_gmls.setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
237
238 // create a vector of target operations
239 std::vector<TargetOperation> lro(4);
240 lro[0] = VectorPointEvaluation;
244
245 // and then pass them to the GMLS class
246 vector_divfree_basis_gmls.addTargets(lro);
247
248 // sets the weighting kernel function from WeightingFunctionType
249 vector_divfree_basis_gmls.setWeightingType(WeightingFunctionType::Power);
250
251 // power to use in that weighting kernel function
252 vector_divfree_basis_gmls.setWeightingParameter(2);
253
254 // generate the alphas that to be combined with data for each target operation requested in lro
255 vector_divfree_basis_gmls.generateAlphas(15 /* # batches */);
256
257 //! [Setting Up The GMLS Object]
258
259 double instantiation_time = timer.seconds();
260 std::cout << "Took " << instantiation_time << "s to complete alphas generation." << std::endl;
261 Kokkos::fence(); // let generateAlphas finish up before using alphas
262 Kokkos::Profiling::pushRegion("Apply Alphas to Data");
263
264 //! [Apply GMLS Alphas To Data]
265
266 // it is important to note that if you expect to use the data as a 1D view, then you should use double*
267 // however, if you know that the target operation will result in a 2D view (vector or matrix output),
268 // then you should template with double** as this is something that can not be infered from the input data
269 // or the target operator at compile time. Additionally, a template argument is required indicating either
270 // Kokkos::HostSpace or Kokkos::DefaultExecutionSpace::memory_space()
271
272 // The Evaluator class takes care of handling input data views as well as the output data views.
273 // It uses information from the GMLS class to determine how many components are in the input
274 // as well as output for any choice of target functionals and then performs the contactions
275 // on the data using the alpha coefficients generated by the GMLS class, all on the device.
276 Evaluator gmls_evaluator_vector(&vector_divfree_basis_gmls);
277
278 auto output_vector_evaluation = gmls_evaluator_vector.applyAlphasToDataAllComponentsAllTargetSites<double**, Kokkos::HostSpace>
279 (vector_sampling_span_basis_data_device, VectorPointEvaluation, VectorPointSample);
280 auto output_curl_vector_evaluation = gmls_evaluator_vector.applyAlphasToDataAllComponentsAllTargetSites<double**, Kokkos::HostSpace>
281 (vector_sampling_span_basis_data_device, CurlOfVectorPointEvaluation, VectorPointSample);
282 auto output_curlcurl_vector_evaluation = gmls_evaluator_vector.applyAlphasToDataAllComponentsAllTargetSites<double**, Kokkos::HostSpace>
283 (vector_sampling_single_polynomial_data_device, CurlCurlOfVectorPointEvaluation, VectorPointSample);
284 auto output_gradient_vector_evaluation = gmls_evaluator_vector.applyAlphasToDataAllComponentsAllTargetSites<double**, Kokkos::HostSpace>
285 (vector_sampling_span_basis_data_device, GradientOfVectorPointEvaluation, VectorPointSample);
286
287 //! [Apply GMLS Alphas To Data]
288
289 Kokkos::fence(); // let application of alphas to data finish before using results
290 Kokkos::Profiling::popRegion();
291 // times the Comparison in Kokkos
292 Kokkos::Profiling::pushRegion("Comparison");
293
294 //! [Check That Solutions Are Correct]
295
296 // loop through the target sites
297 for (int i=0; i<number_target_coords; i++) {
298 // load vector components from output
299 double GMLS_DivFree_VectorX = output_vector_evaluation(i, 0);
300 double GMLS_DivFree_VectorY = (dimension>1) ? output_vector_evaluation(i, 1) : 0;
301 double GMLS_DivFree_VectorZ = (dimension>2) ? output_vector_evaluation(i, 2) : 0;
302
303 // load curl of vector components from output
304 double GMLS_Curl_DivFree_VectorX = output_curl_vector_evaluation(i, 0);
305 double GMLS_Curl_DivFree_VectorY = (dimension>2) ? output_curl_vector_evaluation(i, 1) : 0;
306 double GMLS_Curl_DivFree_VectorZ = (dimension>2) ? output_curl_vector_evaluation(i, 2) : 0;
307
308 // load curl curl of vector components from output
309 double GMLS_CurlCurl_DivFree_VectorX = output_curlcurl_vector_evaluation(i, 0);
310 double GMLS_CurlCurl_DivFree_VectorY = (dimension>1) ? output_curlcurl_vector_evaluation(i, 1) : 0;
311 double GMLS_CurlCurl_DivFree_VectorZ = (dimension>2) ? output_curlcurl_vector_evaluation(i, 2) : 0;
312
313 // load gradient of vector components from output
314 double GMLS_Gradient_DivFree_VectorXX = output_gradient_vector_evaluation(i, 0);
315 double GMLS_Gradient_DivFree_VectorXY = output_gradient_vector_evaluation(i, 1);
316 double GMLS_Gradient_DivFree_VectorXZ = (dimension==3) ? output_gradient_vector_evaluation(i, 2) : 0.0;
317 double GMLS_Gradient_DivFree_VectorYX = (dimension==3) ? output_gradient_vector_evaluation(i, 3) : output_gradient_vector_evaluation(i, 2);
318 double GMLS_Gradient_DivFree_VectorYY = (dimension==3) ? output_gradient_vector_evaluation(i, 4) : output_gradient_vector_evaluation(i, 3);
319 double GMLS_Gradient_DivFree_VectorYZ = (dimension==3) ? output_gradient_vector_evaluation(i, 5) : 0.0;
320 double GMLS_Gradient_DivFree_VectorZX = (dimension==3) ? output_gradient_vector_evaluation(i, 6) : 0.0;
321 double GMLS_Gradient_DivFree_VectorZY = (dimension==3) ? output_gradient_vector_evaluation(i, 7) : 0.0;
322 double GMLS_Gradient_DivFree_VectorZZ = (dimension==3) ? output_gradient_vector_evaluation(i, 8) : 0.0;
323
324 // target site i's coordinate
325 double xval = target_coords(i,0);
326 double yval = (dimension>1) ? target_coords(i,1) : 0;
327 double zval = (dimension>2) ? target_coords(i,2) : 0;
328
329 // evaluation of vector exact solutions
330 double actual_vector[3] = {0, 0, 0};
331 if (dimension>=1) {
332 actual_vector[0] = divfreeTestSolution_span_basis(xval, yval, zval, 0, dimension, order);
333 if (dimension>=2) {
334 actual_vector[1] = divfreeTestSolution_span_basis(xval, yval, zval, 1, dimension, order);
335 if (dimension==3) {
336 actual_vector[2] = divfreeTestSolution_span_basis(xval, yval, zval, 2, dimension, order);
337 }
338 }
339 }
340
341 // evaluation of curl of vector exact solutions
342 double actual_curl_vector[3] = {0, 0, 0};
343 if (dimension>=1) {
344 actual_curl_vector[0] = curldivfreeTestSolution_span_basis(xval, yval, zval, 0, dimension, order);
345 if (dimension==3) {
346 actual_curl_vector[1] = curldivfreeTestSolution_span_basis(xval, yval, zval, 1, dimension, order);
347 actual_curl_vector[2] = curldivfreeTestSolution_span_basis(xval, yval, zval, 2, dimension, order);
348 }
349 }
350
351 // evaluation of curl of curl of vector exact solutions
352 double actual_curlcurl_vector[3] = {0, 0, 0};
353 if (dimension>=1) {
354 actual_curlcurl_vector[0] = curlcurldivfreeTestSolution_single_polynomial(xval, yval, zval, 0, dimension);
355 if (dimension>=2) {
356 actual_curlcurl_vector[1] = curlcurldivfreeTestSolution_single_polynomial(xval, yval, zval, 1, dimension);
357 if (dimension==3) {
358 actual_curlcurl_vector[2] = curlcurldivfreeTestSolution_single_polynomial(xval, yval, zval, 2, dimension);
359 }
360 }
361 }
362
363 double actual_gradient_vector[9] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
364 if (dimension==3) {
365 for (int axes = 0; axes < 9; ++axes)
366 actual_gradient_vector[axes] = gradientdivfreeTestSolution_span_basis(xval, yval, zval, axes/dimension, axes%dimension, dimension, order);
367 }
368 if (dimension==2) {
369 for (int axes = 0; axes < 4; ++axes)
370 actual_gradient_vector[axes] = gradientdivfreeTestSolution_span_basis(xval, yval, zval, axes/dimension, axes%dimension, dimension, order);
371 }
372
373 // check vector evaluation
374 if(std::abs(actual_vector[0] - GMLS_DivFree_VectorX) > failure_tolerance) {
375 all_passed = false;
376 std::cout << i << " Failed VectorX by: " << std::abs(actual_vector[0] - GMLS_DivFree_VectorX) << std::endl;
377 if (dimension>1) {
378 if(std::abs(actual_vector[1] - GMLS_DivFree_VectorY) > failure_tolerance) {
379 all_passed = false;
380 std::cout << i << " Failed VectorY by: " << std::abs(actual_vector[1] - GMLS_DivFree_VectorY) << std::endl;
381 }
382 }
383 if (dimension>2) {
384 if(std::abs(actual_vector[2] - GMLS_DivFree_VectorZ) > failure_tolerance) {
385 all_passed = false;
386 std::cout << i << " Failed VectorZ by: " << std::abs(actual_vector[2] - GMLS_DivFree_VectorZ) << std::endl;
387 }
388 }
389 }
390
391 // check curl of vector evaluation
392 if (dimension==2) {
393 if(std::abs(actual_curl_vector[0] - GMLS_Curl_DivFree_VectorX) > failure_tolerance) {
394 all_passed = false;
395 std::cout << i << " Failed curl VectorX by: " << std::abs(actual_curl_vector[2] - GMLS_Curl_DivFree_VectorX) << std::endl;
396 }
397 } else if (dimension>2) {
398 if(std::abs(actual_curl_vector[0] - GMLS_Curl_DivFree_VectorX) > failure_tolerance) {
399 all_passed = false;
400 std::cout << i << " Failed curl VectorX by: " << std::abs(actual_curl_vector[0] - GMLS_Curl_DivFree_VectorX) << std::endl;
401 }
402 if(std::abs(actual_curl_vector[1] - GMLS_Curl_DivFree_VectorY) > failure_tolerance) {
403 all_passed = false;
404 std::cout << i << " Failed curl VectorY by: " << std::abs(actual_curl_vector[1] - GMLS_Curl_DivFree_VectorY) << std::endl;
405 }
406 if(std::abs(actual_curl_vector[2] - GMLS_Curl_DivFree_VectorZ) > failure_tolerance) {
407 all_passed = false;
408 std::cout << i << " Failed curl VectorZ by: " << std::abs(actual_curl_vector[2] - GMLS_Curl_DivFree_VectorZ) << std::endl;
409 }
410 }
411
412 // check curlcurl curlcurl of vector evaluation
413 if(std::abs(actual_curlcurl_vector[0] - GMLS_CurlCurl_DivFree_VectorX) > failure_tolerance) {
414 all_passed = false;
415 std::cout << i << " Failed curl curl VectorX by: " << std::abs(actual_curlcurl_vector[0] - GMLS_CurlCurl_DivFree_VectorX) << std::endl;
416 }
417 if (dimension>1) {
418 if(std::abs(actual_curlcurl_vector[1] - GMLS_CurlCurl_DivFree_VectorY) > failure_tolerance) {
419 all_passed = false;
420 std::cout << i << " Failed curl curl VectorY by: " << std::abs(actual_curlcurl_vector[1] - GMLS_CurlCurl_DivFree_VectorY) << std::endl;
421 }
422 }
423 if (dimension>2) {
424 if(std::abs(actual_curlcurl_vector[2] - GMLS_CurlCurl_DivFree_VectorZ) > failure_tolerance) {
425 all_passed = false;
426 std::cout << i << " Failed curl curl VectorZ by: " << std::abs(actual_curlcurl_vector[2] - GMLS_CurlCurl_DivFree_VectorZ) << std::endl;
427 }
428 }
429
430 // check gradient of vector evaluation
431 if (dimension==3) {
432 if (std::abs(actual_gradient_vector[0] - GMLS_Gradient_DivFree_VectorXX) > failure_tolerance) {
433 all_passed = false;
434 std::cout << i << " Failed gradient_x VectorX by: " << std::abs(actual_gradient_vector[0] - GMLS_Gradient_DivFree_VectorXX) << std::endl;
435 }
436 if (std::abs(actual_gradient_vector[1] - GMLS_Gradient_DivFree_VectorXY) > failure_tolerance) {
437 all_passed = false;
438 std::cout << i << " Failed gradient_y VectorX by: " << std::abs(actual_gradient_vector[1] - GMLS_Gradient_DivFree_VectorXY) << std::endl;
439 }
440 if (std::abs(actual_gradient_vector[2] - GMLS_Gradient_DivFree_VectorXZ) > failure_tolerance) {
441 all_passed = false;
442 std::cout << i << " Failed gradient_z VectorX by: " << std::abs(actual_gradient_vector[2] - GMLS_Gradient_DivFree_VectorXZ) << std::endl;
443 }
444
445 if (std::abs(actual_gradient_vector[3] - GMLS_Gradient_DivFree_VectorYX) > failure_tolerance) {
446 all_passed = false;
447 std::cout << i << " Failed gradient_x VectorY by: " << std::abs(actual_gradient_vector[3] - GMLS_Gradient_DivFree_VectorYX) << std::endl;
448 }
449 if (std::abs(actual_gradient_vector[4] - GMLS_Gradient_DivFree_VectorYY) > failure_tolerance) {
450 all_passed = false;
451 std::cout << i << " Failed gradient_y VectorY by: " << std::abs(actual_gradient_vector[4] - GMLS_Gradient_DivFree_VectorYY) << std::endl;
452 }
453 if (std::abs(actual_gradient_vector[5] - GMLS_Gradient_DivFree_VectorYZ) > failure_tolerance) {
454 all_passed = false;
455 std::cout << i << " Failed gradient_z VectorY by: " << std::abs(actual_gradient_vector[5] - GMLS_Gradient_DivFree_VectorYZ) << std::endl;
456 }
457
458 if (std::abs(actual_gradient_vector[6] - GMLS_Gradient_DivFree_VectorZX) > failure_tolerance) {
459 all_passed = false;
460 std::cout << i << " Failed gradient_x VectorZ by: " << std::abs(actual_gradient_vector[6] - GMLS_Gradient_DivFree_VectorZX) << std::endl;
461 }
462 if (std::abs(actual_gradient_vector[7] - GMLS_Gradient_DivFree_VectorZY) > failure_tolerance) {
463 all_passed = false;
464 std::cout << i << " Failed gradient_y VectorZ by: " << std::abs(actual_gradient_vector[7] - GMLS_Gradient_DivFree_VectorZY) << std::endl;
465 }
466 if (std::abs(actual_gradient_vector[8] - GMLS_Gradient_DivFree_VectorZZ) > failure_tolerance) {
467 all_passed = false;
468 std::cout << i << " Failed gradient_z VectorZ by: " << std::abs(actual_gradient_vector[8] - GMLS_Gradient_DivFree_VectorZZ) << std::endl;
469 }
470 }
471
472 if (dimension==2) {
473 if (std::abs(actual_gradient_vector[0] - GMLS_Gradient_DivFree_VectorXX) > failure_tolerance) {
474 all_passed = false;
475 std::cout << i << " Failed gradient_x VectorX by: " << std::abs(actual_gradient_vector[0] - GMLS_Gradient_DivFree_VectorXX) << std::endl;
476 }
477 if (std::abs(actual_gradient_vector[1] - GMLS_Gradient_DivFree_VectorXY) > failure_tolerance) {
478 all_passed = false;
479 std::cout << i << " Failed gradient_y VectorX by: " << std::abs(actual_gradient_vector[1] - GMLS_Gradient_DivFree_VectorXY) << std::endl;
480 }
481
482 if (std::abs(actual_gradient_vector[2] - GMLS_Gradient_DivFree_VectorYX) > failure_tolerance) {
483 all_passed = false;
484 std::cout << i << " Failed gradient_x VectorY by: " << std::abs(actual_gradient_vector[2] - GMLS_Gradient_DivFree_VectorYX) << std::endl;
485 }
486 if (std::abs(actual_gradient_vector[3] - GMLS_Gradient_DivFree_VectorYY) > failure_tolerance) {
487 all_passed = false;
488 std::cout << i << " Failed gradient_y VectorY by: " << (std::abs(actual_gradient_vector[3] - GMLS_Gradient_DivFree_VectorYY) > failure_tolerance) << std::endl;
489 }
490 }
491 }
492
493 //! [Check That Solutions Are Correct]
494 // popRegion hidden from tutorial
495 // stop timing comparison loop
496 Kokkos::Profiling::popRegion();
497 //! [Finalize Program]
498
499
500} // end of code block to reduce scope, causing Kokkos View de-allocations
501// otherwise, Views may be deallocating when we call Kokkos::finalize() later
502
503// finalize Kokkos and MPI (if available)
504Kokkos::finalize();
505#ifdef COMPADRE_USE_MPI
506MPI_Finalize();
507#endif
508
509// output to user that test passed or failed
510if(all_passed) {
511 fprintf(stdout, "Passed test \n");
512 return 0;
513} else {
514 fprintf(stdout, "Failed test \n");
515 return -1;
516}
517
518} // main
519
520
521//! [Finalize Program]
int main(int argc, char *args[])
[Parse Command Line Arguments]
KOKKOS_INLINE_FUNCTION double divfreeTestSolution_single_polynomial(double x, double y, double z, int component, int dimension)
KOKKOS_INLINE_FUNCTION double gradientdivfreeTestSolution_span_basis(double x, double y, double z, int input_component, int output_component, int dimension, int exact_order)
KOKKOS_INLINE_FUNCTION double curldivfreeTestSolution_span_basis(double x, double y, double z, int component, int dimension, int exact_order)
KOKKOS_INLINE_FUNCTION double curlcurldivfreeTestSolution_single_polynomial(double x, double y, double z, int component, int dimension)
KOKKOS_INLINE_FUNCTION double divfreeTestSolution_span_basis(double x, double y, double z, int component, int dimension, int exact_order)
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 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...
@ CurlOfVectorPointEvaluation
Point evaluation of the curl of a vector (results in a vector)
@ CurlCurlOfVectorPointEvaluation
Point evaluation of the curl of a curl of a vector (results in a vector)
@ GradientOfVectorPointEvaluation
Point evaluation of the gradient of a vector (results in a matrix, NOT CURRENTLY IMPLEMENTED)
@ VectorPointEvaluation
Point evaluation of a vector (reconstructs entire vector at once, requiring a ReconstructionSpace hav...
constexpr SamplingFunctional VectorPointSample
Point evaluations of the entire vector source function.
@ DivergenceFreeVectorTaylorPolynomial
Divergence-free vector polynomial basis.