29int main (
int argc,
char* args[]) {
32#ifdef COMPADRE_USE_MPI
33MPI_Init(&argc, &args);
37Kokkos::initialize(argc, args);
40bool all_passed =
true;
47 auto order = clp.
order;
54 bool keep_coefficients = number_of_batches==1;
58 const double failure_tolerance = 1e-9;
61 const double laplacian_failure_tolerance = 1e-9;
68 Kokkos::Profiling::pushRegion(
"Setup Point Data");
72 double h_spacing = 0.05;
73 int n_neg1_to_1 = 2*(1/h_spacing) + 1;
76 const int number_source_coords = std::pow(n_neg1_to_1, dimension);
79 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> source_coords_device(
"source coordinates",
80 number_source_coords, 3);
81 Kokkos::View<double**>::HostMirror source_coords = Kokkos::create_mirror_view(source_coords_device);
84 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> target_coords_device (
"target coordinates", number_target_coords, 3);
85 Kokkos::View<double**>::HostMirror target_coords = Kokkos::create_mirror_view(target_coords_device);
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;
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];
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;
112 source_coords(source_index,0) = this_coord[0];
113 source_coords(source_index,1) = 0;
114 source_coords(source_index,2) = 0;
120 for(
int i=0; i<number_target_coords; i++){
123 double rand_dir[3] = {0,0,0};
125 for (
int j=0; j<dimension; ++j) {
127 rand_dir[j] = ((double)rand() / (double) RAND_MAX) - 0.5;
131 for (
int j=0; j<dimension; ++j) {
132 target_coords(i,j) = rand_dir[j];
140 Kokkos::Profiling::popRegion();
141 Kokkos::Profiling::pushRegion(
"Creating Data");
147 Kokkos::deep_copy(source_coords_device, source_coords);
150 Kokkos::deep_copy(target_coords_device, target_coords);
153 Kokkos::View<double*, Kokkos::DefaultExecutionSpace> sampling_data_device(
"samples of true solution",
154 source_coords_device.extent(0));
156 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> gradient_sampling_data_device(
"samples of true gradient",
157 source_coords_device.extent(0), dimension);
159 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> divergence_sampling_data_device
160 (
"samples of true solution for divergence test", source_coords_device.extent(0), dimension);
162 Kokkos::parallel_for(
"Sampling Manufactured Solutions", Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>
163 (0,source_coords.extent(0)), KOKKOS_LAMBDA(
const int i) {
166 double xval = source_coords_device(i,0);
167 double yval = (dimension>1) ? source_coords_device(i,1) : 0;
168 double zval = (dimension>2) ? source_coords_device(i,2) : 0;
171 sampling_data_device(i) =
trueSolution(xval, yval, zval, order, dimension);
174 double true_grad[3] = {0,0,0};
175 trueGradient(true_grad, xval, yval,zval, order, dimension);
177 for (
int j=0; j<dimension; ++j) {
178 gradient_sampling_data_device(i,j) = true_grad[j];
189 Kokkos::Profiling::popRegion();
190 Kokkos::Profiling::pushRegion(
"Neighbor Search");
199 double epsilon_multiplier = 1.4;
204 Kokkos::View<int*> neighbor_lists_device(
"neighbor lists",
206 Kokkos::View<int*>::HostMirror neighbor_lists = Kokkos::create_mirror_view(neighbor_lists_device);
209 Kokkos::View<int*> number_of_neighbors_list_device(
"number of neighbor lists",
210 number_target_coords);
211 Kokkos::View<int*>::HostMirror number_of_neighbors_list = Kokkos::create_mirror_view(number_of_neighbors_list_device);
214 Kokkos::View<double*, Kokkos::DefaultExecutionSpace> epsilon_device(
"h supports", number_target_coords);
215 Kokkos::View<double*>::HostMirror epsilon = Kokkos::create_mirror_view(epsilon_device);
222 size_t storage_size = point_cloud_search.generateCRNeighborListsFromKNNSearch(
true , target_coords, neighbor_lists,
223 number_of_neighbors_list, epsilon, min_neighbors, epsilon_multiplier);
226 Kokkos::resize(neighbor_lists_device, storage_size);
227 neighbor_lists = Kokkos::create_mirror_view(neighbor_lists_device);
230 point_cloud_search.generateCRNeighborListsFromKNNSearch(
false , target_coords, neighbor_lists,
231 number_of_neighbors_list, epsilon, min_neighbors, epsilon_multiplier);
235 Kokkos::Profiling::popRegion();
246 Kokkos::deep_copy(neighbor_lists_device, neighbor_lists);
247 Kokkos::deep_copy(number_of_neighbors_list_device, number_of_neighbors_list);
248 Kokkos::deep_copy(epsilon_device, epsilon);
253 solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
273 my_GMLS.
setProblemData(neighbor_lists_device, number_of_neighbors_list_device, source_coords_device, target_coords_device, epsilon_device);
276 std::vector<TargetOperation> lro(5);
298 double instantiation_time = timer.seconds();
299 std::cout <<
"Took " << instantiation_time <<
"s to complete alphas generation." << std::endl;
301 Kokkos::Profiling::pushRegion(
"Apply Alphas to Data");
333 decltype(output_curl) scalar_coefficients;
334 if (number_of_batches==1)
335 scalar_coefficients =
337 (sampling_data_device);
342 Kokkos::Profiling::popRegion();
344 Kokkos::Profiling::pushRegion(
"Comparison");
350 for (
int i=0; i<number_target_coords; i++) {
353 double GMLS_value = output_value(i);
356 double GMLS_Laplacian = output_laplacian(i);
362 double GMLS_GradX = (number_of_batches==1) ? scalar_coefficients(i,1)*1./epsilon(i) : output_gradient(i,0);
365 double GMLS_GradY = (dimension>1) ? output_gradient(i,1) : 0;
368 double GMLS_GradZ = (dimension>2) ? output_gradient(i,2) : 0;
371 double GMLS_Divergence = output_divergence(i);
374 double GMLS_CurlX = (dimension>1) ? output_curl(i,0) : 0;
375 double GMLS_CurlY = (dimension>1) ? output_curl(i,1) : 0;
376 double GMLS_CurlZ = (dimension>2) ? output_curl(i,2) : 0;
380 double xval = target_coords(i,0);
381 double yval = (dimension>1) ? target_coords(i,1) : 0;
382 double zval = (dimension>2) ? target_coords(i,2) : 0;
385 double actual_value =
trueSolution(xval, yval, zval, order, dimension);
386 double actual_Laplacian =
trueLaplacian(xval, yval, zval, order, dimension);
388 double actual_Gradient[3] = {0,0,0};
389 trueGradient(actual_Gradient, xval, yval, zval, order, dimension);
391 double actual_Divergence;
392 actual_Divergence =
trueLaplacian(xval, yval, zval, order, dimension);
394 double actual_Curl[3] = {0,0,0};
405 if(GMLS_value!=GMLS_value || std::abs(actual_value - GMLS_value) > failure_tolerance) {
407 std::cout << i <<
" Failed Actual by: " << std::abs(actual_value - GMLS_value) << std::endl;
411 if(std::abs(actual_Laplacian - GMLS_Laplacian) > laplacian_failure_tolerance) {
413 std::cout << i <<
" Failed Laplacian by: " << std::abs(actual_Laplacian - GMLS_Laplacian) << std::endl;
417 if(std::abs(actual_Gradient[0] - GMLS_GradX) > failure_tolerance) {
419 std::cout << i <<
" Failed GradX by: " << std::abs(actual_Gradient[0] - GMLS_GradX) << std::endl;
421 if(std::abs(actual_Gradient[1] - GMLS_GradY) > failure_tolerance) {
423 std::cout << i <<
" Failed GradY by: " << std::abs(actual_Gradient[1] - GMLS_GradY) << std::endl;
427 if(std::abs(actual_Gradient[2] - GMLS_GradZ) > failure_tolerance) {
429 std::cout << i <<
" Failed GradZ by: " << std::abs(actual_Gradient[2] - GMLS_GradZ) << std::endl;
435 if(std::abs(actual_Divergence - GMLS_Divergence) > failure_tolerance) {
437 std::cout << i <<
" Failed Divergence by: " << std::abs(actual_Divergence - GMLS_Divergence) << std::endl;
444 tmp_diff += std::abs(actual_Curl[0] - GMLS_CurlX) + std::abs(actual_Curl[1] - GMLS_CurlY);
446 tmp_diff += std::abs(actual_Curl[2] - GMLS_CurlZ);
447 if(std::abs(tmp_diff) > failure_tolerance) {
449 std::cout << i <<
" Failed Curl by: " << std::abs(tmp_diff) << std::endl;
458 Kokkos::Profiling::popRegion();
467#ifdef COMPADRE_USE_MPI
473 fprintf(stdout,
"Passed test \n");
476 fprintf(stdout,
"Failed test \n");