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;
56 const double failure_tolerance = 9e-8;
63 Kokkos::Profiling::pushRegion(
"Setup Point Data");
67 double h_spacing = 0.05;
68 int n_neg1_to_1 = 2*(1/h_spacing) + 1;
71 const int number_source_coords = std::pow(n_neg1_to_1, dimension);
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);
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);
83 Kokkos::View<double***, Kokkos::DefaultExecutionSpace> tangent_bundles_device (
"tangent bundles", number_target_coords, 3, 3);
84 Kokkos::View<double***>::HostMirror tangent_bundles = Kokkos::create_mirror_view(tangent_bundles_device);
88 double this_coord[3] = {0,0,0};
89 for (
int i=-n_neg1_to_1/2; i<n_neg1_to_1/2+1; ++i) {
90 this_coord[0] = i*h_spacing;
91 for (
int j=-n_neg1_to_1/2; j<n_neg1_to_1/2+1; ++j) {
92 this_coord[1] = j*h_spacing;
93 for (
int k=-n_neg1_to_1/2; k<n_neg1_to_1/2+1; ++k) {
94 this_coord[2] = k*h_spacing;
96 source_coords(source_index,0) = this_coord[0];
97 source_coords(source_index,1) = this_coord[1];
98 source_coords(source_index,2) = this_coord[2];
103 source_coords(source_index,0) = this_coord[0];
104 source_coords(source_index,1) = this_coord[1];
105 source_coords(source_index,2) = 0;
110 source_coords(source_index,0) = this_coord[0];
111 source_coords(source_index,1) = 0;
112 source_coords(source_index,2) = 0;
121 std::mt19937 rng(50);
123 std::uniform_int_distribution<int> gen_num_neighbours(0, source_index);
125 for (
int i=0; i<number_target_coords; ++i) {
126 const int source_site_to_copy = gen_num_neighbours(rng);
127 for (
int j=0; j<3; j++) {
128 target_coords(i, j) = source_coords(source_site_to_copy, j);
130 if (constraint_name ==
"NEUMANN_GRAD_SCALAR") {
131 if (dimension == 3) {
132 tangent_bundles(i, 0, 0) = 0.0;
133 tangent_bundles(i, 0, 1) = 0.0;
134 tangent_bundles(i, 0, 2) = 0.0;
135 tangent_bundles(i, 1, 0) = 0.0;
136 tangent_bundles(i, 1, 1) = 0.0;
137 tangent_bundles(i, 1, 2) = 0.0;
138 tangent_bundles(i, 2, 0) = 1.0/(sqrt(3.0));
139 tangent_bundles(i, 2, 1) = 1.0/(sqrt(3.0));
140 tangent_bundles(i, 2, 2) = 1.0/(sqrt(3.0));
142 if (dimension == 2) {
143 tangent_bundles(i, 0, 0) = 0.0;
144 tangent_bundles(i, 0, 1) = 0.0;
145 tangent_bundles(i, 1, 0) = 1.0/(sqrt(2.0));
146 tangent_bundles(i, 1, 1) = 1.0/(sqrt(2.0));
153 Kokkos::Profiling::popRegion();
154 Kokkos::Profiling::pushRegion(
"Creating Data");
160 Kokkos::deep_copy(source_coords_device, source_coords);
163 Kokkos::deep_copy(target_coords_device, target_coords);
166 Kokkos::View<double*, Kokkos::DefaultExecutionSpace> sampling_data_device(
"samples of true solution",
167 source_coords_device.extent(0));
169 Kokkos::parallel_for(
"Sampling Manufactured Solutions", Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>
170 (0,source_coords.extent(0)), KOKKOS_LAMBDA(
const int i) {
172 double xval = source_coords_device(i,0);
173 double yval = (dimension>1) ? source_coords_device(i,1) : 0;
174 double zval = (dimension>2) ? source_coords_device(i,2) : 0;
177 sampling_data_device(i) =
trueSolution(xval, yval, zval, order, dimension);
183 Kokkos::Profiling::popRegion();
184 Kokkos::Profiling::pushRegion(
"Neighbor Search");
196 double epsilon_multiplier = 2.2;
197 int estimated_upper_bound_number_neighbors =
198 point_cloud_search.getEstimatedNumberNeighborsUpperBound(min_neighbors, dimension, epsilon_multiplier);
200 Kokkos::View<int**, Kokkos::DefaultExecutionSpace> neighbor_lists_device(
"neighbor lists",
201 number_target_coords, estimated_upper_bound_number_neighbors);
202 Kokkos::View<int**>::HostMirror neighbor_lists = Kokkos::create_mirror_view(neighbor_lists_device);
205 Kokkos::View<double*, Kokkos::DefaultExecutionSpace> epsilon_device(
"h supports", number_target_coords);
206 Kokkos::View<double*>::HostMirror epsilon = Kokkos::create_mirror_view(epsilon_device);
211 point_cloud_search.generate2DNeighborListsFromKNNSearch(
false , target_coords, neighbor_lists,
212 epsilon, min_neighbors, epsilon_multiplier);
216 Kokkos::Profiling::popRegion();
227 Kokkos::deep_copy(neighbor_lists_device, neighbor_lists);
228 Kokkos::deep_copy(epsilon_device, epsilon);
236 solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
244 solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
261 scalar_basis_gmls.
setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
262 vector_basis_gmls.
setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
263 if (constraint_name ==
"NEUMANN_GRAD_SCALAR") {
269 std::vector<TargetOperation> lro(2);
296 double instantiation_time = timer.seconds();
297 std::cout <<
"Took " << instantiation_time <<
"s to complete alphas generation." << std::endl;
299 Kokkos::Profiling::pushRegion(
"Apply Alphas to Data");
313 Evaluator gmls_evaluator_scalar(&scalar_basis_gmls);
314 Evaluator gmls_evaluator_vector(&vector_basis_gmls);
332 Kokkos::Profiling::popRegion();
334 Kokkos::Profiling::pushRegion(
"Comparison");
339 for (
int i=0; i<number_target_coords; i++) {
342 double xval = target_coords(i,0);
343 double yval = (dimension>1) ? target_coords(i,1) : 0;
344 double zval = (dimension>2) ? target_coords(i,2) : 0;
347 double actual_Divergence;
348 actual_Divergence =
trueLaplacian(xval, yval, zval, order, dimension);
349 double actual_Gradient[3] = {0,0,0};
350 trueGradient(actual_Gradient, xval, yval, zval, order, dimension);
353 double g = (dimension == 3) ? (1.0/sqrt(3.0))*(actual_Gradient[0] + actual_Gradient[1] + actual_Gradient[2])
354 : (1.0/sqrt(2.0))*(actual_Gradient[0] + actual_Gradient[1]);
357 int num_neigh_i = neighbor_lists(i, 0);
360 double GMLS_Divergence_Scalar = output_divergence_scalar(i);
361 double GMLS_Divergence_Vector = output_divergence_vector(i);
364 if (constraint_name ==
"NEUMANN_GRAD_SCALAR") {
366 GMLS_Divergence_Scalar = GMLS_Divergence_Scalar + b_i_scalar*g;
369 GMLS_Divergence_Vector = GMLS_Divergence_Vector + b_i_vector*g;
373 if(std::abs(actual_Divergence - GMLS_Divergence_Scalar) > failure_tolerance) {
375 std::cout << i <<
" Failed Divergence on SCALAR basis by: " << std::abs(actual_Divergence - GMLS_Divergence_Scalar) << std::endl;
376 std::cout << i <<
" GMLS " << GMLS_Divergence_Scalar <<
" adjusted " << GMLS_Divergence_Scalar <<
" actual " << actual_Divergence << std::endl;
380 if(std::abs(actual_Divergence - GMLS_Divergence_Vector) > failure_tolerance) {
382 std::cout << i <<
" Failed Divergence on VECTOR basis by: " << std::abs(actual_Divergence - GMLS_Divergence_Vector) << std::endl;
383 std::cout << i <<
" GMLS " << GMLS_Divergence_Vector <<
" adjusted " << GMLS_Divergence_Vector <<
" actual " << actual_Divergence << std::endl;
387 double Scalar_GMLS_GradX = (dimension>1) ? output_gradient_scalar(i,0) : 0;
388 double Scalar_GMLS_GradY = (dimension>1) ? output_gradient_scalar(i,1) : 0;
389 double Scalar_GMLS_GradZ = (dimension>2) ? output_gradient_scalar(i,2) : 0;
391 double Vector_GMLS_GradX = (dimension>1) ? output_gradient_vector(i,0) : 0;
392 double Vector_GMLS_GradY = (dimension>1) ? output_gradient_vector(i,1) : 0;
393 double Vector_GMLS_GradZ = (dimension>2) ? output_gradient_vector(i,2) : 0;
396 if (constraint_name ==
"NEUMANN_GRAD_SCALAR") {
400 Scalar_GMLS_GradX = Scalar_GMLS_GradX + bx_i_scalar*g;
401 Scalar_GMLS_GradY = Scalar_GMLS_GradY + by_i_scalar*g;
402 Scalar_GMLS_GradZ = Scalar_GMLS_GradZ + bz_i_scalar*g;
407 Vector_GMLS_GradX = Vector_GMLS_GradX + bx_i_vector*g;
408 Vector_GMLS_GradY = Vector_GMLS_GradY + by_i_vector*g;
409 Vector_GMLS_GradZ = Vector_GMLS_GradZ + bz_i_vector*g;
413 if(std::abs(actual_Gradient[0] - Scalar_GMLS_GradX) > failure_tolerance) {
415 std::cout << i <<
" Failed Scalar GradX by: " << std::abs(actual_Gradient[0] - Scalar_GMLS_GradX) << std::endl;
417 if(std::abs(actual_Gradient[1] - Scalar_GMLS_GradY) > failure_tolerance) {
419 std::cout << i <<
" Failed Scalar GradY by: " << std::abs(actual_Gradient[1] - Scalar_GMLS_GradY) << std::endl;
423 if(std::abs(actual_Gradient[2] - Scalar_GMLS_GradZ) > failure_tolerance) {
425 std::cout << i <<
" Failed Scalar GradZ by: " << std::abs(actual_Gradient[2] - Scalar_GMLS_GradZ) << std::endl;
431 if(std::abs(actual_Gradient[0] - Vector_GMLS_GradX) > failure_tolerance) {
433 std::cout << i <<
" Failed Vector GradX by: " << std::abs(actual_Gradient[0] - Vector_GMLS_GradX) << std::endl;
435 if(std::abs(actual_Gradient[1] - Vector_GMLS_GradY) > failure_tolerance) {
437 std::cout << i <<
" Failed Vector GradY by: " << std::abs(actual_Gradient[1] - Vector_GMLS_GradY) << std::endl;
441 if(std::abs(actual_Gradient[2] - Vector_GMLS_GradZ) > failure_tolerance) {
443 std::cout << i <<
" Failed Vector GradZ by: " << std::abs(actual_Gradient[2] - Vector_GMLS_GradZ) << std::endl;
453 Kokkos::Profiling::popRegion();
461#ifdef COMPADRE_USE_MPI
467 fprintf(stdout,
"Passed test \n");
470 fprintf(stdout,
"Failed test \n");