59int main (
int argc,
char* args[]) {
62#ifdef COMPADRE_USE_MPI
63MPI_Init(&argc, &args);
67Kokkos::initialize(argc, args);
74 auto order = clp.
order;
88 Kokkos::Profiling::pushRegion(
"Setup Point Data");
93 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> source_coords_device(
"source coordinates",
94 1.25*N_pts_on_sphere, 3);
95 Kokkos::View<double**>::HostMirror source_coords = Kokkos::create_mirror_view(source_coords_device);
100 int number_source_coords;
106 double a = 4*
PI*r*r/N_pts_on_sphere;
107 double d = std::sqrt(a);
108 int M_theta = std::round(
PI/d);
109 double d_theta =
PI/M_theta;
110 double d_phi = a/d_theta;
111 for (
int i=0; i<M_theta; ++i) {
112 double theta =
PI*(i + 0.5)/M_theta;
113 int M_phi = std::round(2*
PI*std::sin(theta)/d_phi);
114 for (
int j=0; j<M_phi; ++j) {
115 double phi = 2*
PI*j/M_phi;
116 source_coords(N_count, 0) = theta;
117 source_coords(N_count, 1) = phi;
122 for (
int i=0; i<N_count; ++i) {
123 double theta = source_coords(i,0);
124 double phi = source_coords(i,1);
125 source_coords(i,0) = r*std::sin(theta)*std::cos(phi);
126 source_coords(i,1) = r*std::sin(theta)*std::sin(phi);
127 source_coords(i,2) = r*cos(theta);
130 number_source_coords = N_count;
134 Kokkos::resize(source_coords, number_source_coords, 3);
135 Kokkos::resize(source_coords_device, number_source_coords, 3);
138 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> target_coords_device (
"target coordinates",
139 number_target_coords, 3);
140 Kokkos::View<double**>::HostMirror target_coords = Kokkos::create_mirror_view(target_coords_device);
143 bool enough_pts_found =
false;
147 double a = 4.0*
PI*r*r/((double)(5.0*number_target_coords));
148 double d = std::sqrt(a);
149 int M_theta = std::round(
PI/d);
150 double d_theta =
PI/((double)M_theta);
151 double d_phi = a/d_theta;
152 for (
int i=0; i<M_theta; ++i) {
153 double theta =
PI*(i + 0.5)/M_theta;
154 int M_phi = std::round(2*
PI*std::sin(theta)/d_phi);
155 for (
int j=0; j<M_phi; ++j) {
156 double phi = 2*
PI*j/M_phi;
157 target_coords(N_count, 0) = theta;
158 target_coords(N_count, 1) = phi;
160 if (N_count == number_target_coords) {
161 enough_pts_found =
true;
165 if (enough_pts_found)
break;
168 for (
int i=0; i<N_count; ++i) {
169 double theta = target_coords(i,0);
170 double phi = target_coords(i,1);
171 target_coords(i,0) = r*std::sin(theta)*std::cos(phi);
172 target_coords(i,1) = r*std::sin(theta)*std::sin(phi);
173 target_coords(i,2) = r*cos(theta);
181 Kokkos::Profiling::popRegion();
183 Kokkos::Profiling::pushRegion(
"Creating Data");
189 Kokkos::deep_copy(source_coords_device, source_coords);
192 Kokkos::deep_copy(target_coords_device, target_coords);
199 Kokkos::View<double*, Kokkos::DefaultExecutionSpace> sampling_data_device(
"samples of true solution",
200 source_coords_device.extent(0));
201 Kokkos::View<double*, Kokkos::DefaultExecutionSpace> ones_data_device(
"samples of ones",
202 source_coords_device.extent(0));
203 Kokkos::deep_copy(ones_data_device, 1.0);
206 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> sampling_vector_data_device(
"samples of vector true solution",
207 source_coords_device.extent(0), 3);
209 Kokkos::parallel_for(
"Sampling Manufactured Solutions", Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>
210 (0,source_coords.extent(0)), KOKKOS_LAMBDA(
const int i) {
213 double xval = source_coords_device(i,0);
214 double yval = (dimension>1) ? source_coords_device(i,1) : 0;
215 double zval = (dimension>2) ? source_coords_device(i,2) : 0;
220 for (
int j=0; j<3; ++j) {
221 double gradient[3] = {0,0,0};
223 sampling_vector_data_device(i,j) = gradient[j];
230 Kokkos::Profiling::popRegion();
231 Kokkos::Profiling::pushRegion(
"Neighbor Search");
242 double epsilon_multiplier = 1.7;
243 int estimated_upper_bound_number_neighbors =
244 point_cloud_search.getEstimatedNumberNeighborsUpperBound(min_neighbors, dimension, epsilon_multiplier);
246 Kokkos::View<int**, Kokkos::DefaultExecutionSpace> neighbor_lists_device(
"neighbor lists",
247 number_target_coords, estimated_upper_bound_number_neighbors);
248 Kokkos::View<int**>::HostMirror neighbor_lists = Kokkos::create_mirror_view(neighbor_lists_device);
251 Kokkos::View<double*, Kokkos::DefaultExecutionSpace> epsilon_device(
"h supports", number_target_coords);
252 Kokkos::View<double*>::HostMirror epsilon = Kokkos::create_mirror_view(epsilon_device);
257 point_cloud_search.generate2DNeighborListsFromKNNSearch(
false , target_coords, neighbor_lists,
258 epsilon, min_neighbors, epsilon_multiplier);
262 Kokkos::Profiling::popRegion();
273 Kokkos::deep_copy(neighbor_lists_device, neighbor_lists);
274 Kokkos::deep_copy(epsilon_device, epsilon);
278 GMLS my_GMLS_scalar(order, dimension,
279 solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
296 my_GMLS_scalar.
setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
308 std::vector<TargetOperation> lro_scalar(3);
331 Kokkos::Profiling::pushRegion(
"Full Polynomial Basis GMLS Solution");
338 solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
341 my_GMLS_vector.
setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
343 std::vector<TargetOperation> lro_vector(2);
352 Kokkos::Profiling::popRegion();
354 Kokkos::Profiling::pushRegion(
"Scalar Polynomial Basis Repeated to Form a Vector GMLS Solution");
379 solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
382 my_GMLS_vector_of_scalar_clones.
setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
384 std::vector<TargetOperation> lro_vector_of_scalar_clones(2);
387 my_GMLS_vector_of_scalar_clones.
addTargets(lro_vector_of_scalar_clones);
390 my_GMLS_vector_of_scalar_clones.
setWeightingType(WeightingFunctionType::Power);
393 Kokkos::Profiling::popRegion();
398 double instantiation_time = timer.seconds();
399 std::cout <<
"Took " << instantiation_time <<
"s to complete alphas generation." << std::endl;
401 Kokkos::Profiling::pushRegion(
"Apply Alphas to Data");
416 Evaluator scalar_gmls_evaluator(&my_GMLS_scalar);
417 Evaluator vector_gmls_evaluator(&my_GMLS_vector);
418 Evaluator vector_gmls_evaluator_of_scalar_clones(&my_GMLS_vector_of_scalar_clones);
448 auto output_vector_of_scalar_clones =
488 Kokkos::Profiling::popRegion();
490 Kokkos::Profiling::pushRegion(
"Comparison");
494 double tangent_bundle_error = 0;
495 double tangent_bundle_norm = 0;
496 double values_error = 0;
497 double values_norm = 0;
500 double curl_ambient_error = 0;
501 double curl_ambient_norm = 0;
506 double vector_ambient_error = 0;
507 double vector_ambient_norm = 0;
510 double vector_of_scalar_clones_ambient_error = 0;
511 double vector_of_scalar_clones_ambient_norm = 0;
517 auto prestencil_weights = Kokkos::create_mirror_view(d_prestencil_weights);
518 Kokkos::deep_copy(prestencil_weights, d_prestencil_weights);
522 auto tangent_directions = Kokkos::create_mirror_view(d_tangent_directions);
523 Kokkos::deep_copy(tangent_directions, d_tangent_directions);
526 for (
int i=0; i<number_target_coords; i++) {
534 double GMLS_value = output_value(i);
537 double GMLS_gc = output_gaussian_curvature(i);
547 double xval = source_coords(neighbor_lists(i,1),0);
548 double yval = (dimension>1) ? source_coords(neighbor_lists(i,1),1) : 0;
549 double zval = (dimension>2) ? source_coords(neighbor_lists(i,1),2) : 0;
550 double coord[3] = {xval, yval, zval};
553 for (
int j=0; j<dimension-1; ++j) {
554 double tangent_inner_prod = 0;
555 for (
int k=0; k<std::min(dimension,3); ++k) {
556 tangent_inner_prod += coord[k] * prestencil_weights(0, i, 0 , j, k);
558 tangent_bundle_error += tangent_inner_prod * tangent_inner_prod;
560 double normal_inner_prod = 0;
561 for (
int k=0; k<dimension; ++k) {
562 normal_inner_prod += coord[k] * my_GMLS_scalar.
getTangentBundle(i, dimension-1, k);
565 double normal_error_to_sum = (normal_inner_prod > 0) ? normal_inner_prod - 1 : normal_inner_prod + 1;
566 tangent_bundle_error += normal_error_to_sum * normal_error_to_sum;
567 tangent_bundle_norm += 1;
571 double actual_gc = 1.0;
573 double actual_Gradient_ambient[3] = {0,0,0};
576 double actual_Curl_ambient[3] = {0,0,0};
579 values_error += (GMLS_value - actual_value)*(GMLS_value - actual_value);
580 values_norm += actual_value*actual_value;
582 gc_error += (GMLS_gc - actual_gc)*(GMLS_gc - actual_gc);
585 for (
int j=0; j<dimension; ++j) u[j] = output_curl(i,j);
587 for (
int j=0; j<dimension; ++j) output_curl(i,j) = u[j];
589 for (
int j=0; j<dimension; ++j) {
590 curl_ambient_error += (output_curl(i,j) - actual_Curl_ambient[j])*(output_curl(i,j) - actual_Curl_ambient[j]);
591 curl_ambient_norm += actual_Curl_ambient[j]*actual_Curl_ambient[j];
602 for (
int j=0; j<dimension; ++j) u[j] = output_vector(i,j);
604 for (
int j=0; j<dimension; ++j) output_vector(i,j) = u[j];
606 for (
int j=0; j<dimension; ++j) {
607 vector_ambient_error += (output_vector(i,j) - actual_Gradient_ambient[j])*(output_vector(i,j) - actual_Gradient_ambient[j]);
608 vector_ambient_norm += actual_Gradient_ambient[j]*actual_Gradient_ambient[j];
615 for (
int j=0; j<dimension; ++j) u[j] = output_vector_of_scalar_clones(i,j);
617 for (
int j=0; j<dimension; ++j) output_vector_of_scalar_clones(i,j) = u[j];
632 for (
int j=0; j<dimension; ++j) {
633 vector_of_scalar_clones_ambient_error += (output_vector_of_scalar_clones(i,j) - actual_Gradient_ambient[j])*(output_vector_of_scalar_clones(i,j) - actual_Gradient_ambient[j]);
634 vector_of_scalar_clones_ambient_norm += actual_Gradient_ambient[j]*actual_Gradient_ambient[j];
642 tangent_bundle_error /= number_target_coords;
643 tangent_bundle_error = std::sqrt(tangent_bundle_error);
644 tangent_bundle_norm /= number_target_coords;
645 tangent_bundle_norm = std::sqrt(tangent_bundle_norm);
647 values_error /= number_target_coords;
648 values_error = std::sqrt(values_error);
649 values_norm /= number_target_coords;
650 values_norm = std::sqrt(values_norm);
652 gc_error /= number_target_coords;
653 gc_error = std::sqrt(gc_error);
654 gc_norm /= number_target_coords;
655 gc_norm = std::sqrt(gc_norm);
657 curl_ambient_error /= number_target_coords;
658 curl_ambient_error = std::sqrt(curl_ambient_error);
659 curl_ambient_norm /= number_target_coords;
660 curl_ambient_norm = std::sqrt(curl_ambient_norm);
672 vector_ambient_error /= number_target_coords;
673 vector_ambient_error = std::sqrt(vector_ambient_error);
674 vector_ambient_norm /= number_target_coords;
675 vector_ambient_norm = std::sqrt(vector_ambient_norm);
682 vector_of_scalar_clones_ambient_error /= number_target_coords;
683 vector_of_scalar_clones_ambient_error = std::sqrt(vector_of_scalar_clones_ambient_error);
684 vector_of_scalar_clones_ambient_norm /= number_target_coords;
685 vector_of_scalar_clones_ambient_norm = std::sqrt(vector_of_scalar_clones_ambient_norm);
692 printf(
"Tangent Bundle Error: %g\n", tangent_bundle_error / tangent_bundle_norm);
693 printf(
"Point Value Error: %g\n", values_error / values_norm);
694 printf(
"Gaussian Curvature Error: %g\n", gc_error / gc_norm);
695 printf(
"Surface Curl (Ambient) Error: %g\n", curl_ambient_error / curl_ambient_norm);
698 printf(
"Surface Vector (VectorBasis) Error: %g\n", vector_ambient_error / vector_ambient_norm);
700 printf(
"Surface Vector (ScalarClones) Error: %g\n",
701 vector_of_scalar_clones_ambient_error / vector_of_scalar_clones_ambient_norm);
707 Kokkos::Profiling::popRegion();
716#ifdef COMPADRE_USE_MPI