9 #include <Compadre_Config.h> 17 #ifdef COMPADRE_USE_MPI 21 #include <Kokkos_Timer.hpp> 22 #include <Kokkos_Core.hpp> 29 int main (
int argc,
char* args[]) {
32 #ifdef COMPADRE_USE_MPI 33 MPI_Init(&argc, &args);
37 Kokkos::initialize(argc, args);
44 auto order = clp.
order;
58 Kokkos::Profiling::pushRegion(
"Setup Point Data");
63 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> source_coords_device(
"source coordinates",
64 1.25*N_pts_on_sphere, 3);
65 Kokkos::View<double**>::HostMirror source_coords = Kokkos::create_mirror_view(source_coords_device);
70 int number_source_coords;
76 double a = 4*
PI*r*r/N_pts_on_sphere;
77 double d = std::sqrt(a);
78 int M_theta = std::round(
PI/d);
79 double d_theta =
PI/M_theta;
80 double d_phi = a/d_theta;
81 for (
int i=0; i<M_theta; ++i) {
82 double theta =
PI*(i + 0.5)/M_theta;
83 int M_phi = std::round(2*
PI*std::sin(theta)/d_phi);
84 for (
int j=0; j<M_phi; ++j) {
85 double phi = 2*
PI*j/M_phi;
86 source_coords(N_count, 0) = theta;
87 source_coords(N_count, 1) = phi;
92 for (
int i=0; i<N_count; ++i) {
93 double theta = source_coords(i,0);
94 double phi = source_coords(i,1);
95 source_coords(i,0) = r*std::sin(theta)*std::cos(phi);
96 source_coords(i,1) = r*std::sin(theta)*std::sin(phi);
97 source_coords(i,2) = r*cos(theta);
100 number_source_coords = N_count;
104 Kokkos::resize(source_coords, number_source_coords, 3);
105 Kokkos::resize(source_coords_device, number_source_coords, 3);
108 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> target_coords_device (
"target coordinates",
109 number_target_coords, 3);
110 Kokkos::View<double**>::HostMirror target_coords = Kokkos::create_mirror_view(target_coords_device);
113 bool enough_pts_found =
false;
117 double a = 4.0*
PI*r*r/((double)(5.0*number_target_coords));
118 double d = std::sqrt(a);
119 int M_theta = std::round(
PI/d);
120 double d_theta =
PI/((double)M_theta);
121 double d_phi = a/d_theta;
122 for (
int i=0; i<M_theta; ++i) {
123 double theta =
PI*(i + 0.5)/M_theta;
124 int M_phi = std::round(2*
PI*std::sin(theta)/d_phi);
125 for (
int j=0; j<M_phi; ++j) {
126 double phi = 2*
PI*j/M_phi;
127 target_coords(N_count, 0) = theta;
128 target_coords(N_count, 1) = phi;
130 if (N_count == number_target_coords) {
131 enough_pts_found =
true;
135 if (enough_pts_found)
break;
138 for (
int i=0; i<N_count; ++i) {
139 double theta = target_coords(i,0);
140 double phi = target_coords(i,1);
141 target_coords(i,0) = r*std::sin(theta)*std::cos(phi);
142 target_coords(i,1) = r*std::sin(theta)*std::sin(phi);
143 target_coords(i,2) = r*cos(theta);
151 Kokkos::Profiling::popRegion();
153 Kokkos::Profiling::pushRegion(
"Creating Data");
159 Kokkos::deep_copy(source_coords_device, source_coords);
162 Kokkos::deep_copy(target_coords_device, target_coords);
169 Kokkos::View<double*, Kokkos::DefaultExecutionSpace> sampling_data_device(
"samples of true solution",
170 source_coords_device.extent(0));
172 Kokkos::View<double*, Kokkos::DefaultExecutionSpace> ones_data_device(
"samples of ones",
173 source_coords_device.extent(0));
174 Kokkos::deep_copy(ones_data_device, 1.0);
177 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> sampling_vector_data_device(
"samples of vector true solution",
178 source_coords_device.extent(0), 3);
180 Kokkos::parallel_for(
"Sampling Manufactured Solutions", Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>
181 (0,source_coords.extent(0)), KOKKOS_LAMBDA(
const int i) {
184 double xval = source_coords_device(i,0);
185 double yval = (dimension>1) ? source_coords_device(i,1) : 0;
186 double zval = (dimension>2) ? source_coords_device(i,2) : 0;
191 for (
int j=0; j<3; ++j) {
192 double gradient[3] = {0,0,0};
194 sampling_vector_data_device(i,j) = gradient[j];
201 Kokkos::Profiling::popRegion();
202 Kokkos::Profiling::pushRegion(
"Neighbor Search");
213 double epsilon_multiplier = 1.9;
214 int estimated_upper_bound_number_neighbors =
215 point_cloud_search.getEstimatedNumberNeighborsUpperBound(min_neighbors, dimension, epsilon_multiplier);
217 Kokkos::View<int**, Kokkos::DefaultExecutionSpace> neighbor_lists_device(
"neighbor lists",
218 number_target_coords, estimated_upper_bound_number_neighbors);
219 Kokkos::View<int**>::HostMirror neighbor_lists = Kokkos::create_mirror_view(neighbor_lists_device);
222 Kokkos::View<double*, Kokkos::DefaultExecutionSpace> epsilon_device(
"h supports", number_target_coords);
223 Kokkos::View<double*>::HostMirror epsilon = Kokkos::create_mirror_view(epsilon_device);
228 point_cloud_search.generate2DNeighborListsFromKNNSearch(
false , target_coords, neighbor_lists,
229 epsilon, min_neighbors, epsilon_multiplier);
233 Kokkos::Profiling::popRegion();
244 Kokkos::deep_copy(neighbor_lists_device, neighbor_lists);
245 Kokkos::deep_copy(epsilon_device, epsilon);
249 GMLS my_GMLS_scalar(order, dimension,
250 solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
267 my_GMLS_scalar.
setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
272 my_GMLS_scalar.setReferenceOutwardNormalDirection(target_coords_device,
true );
275 std::vector<TargetOperation> lro_scalar(4);
282 my_GMLS_scalar.addTargets(lro_scalar);
288 my_GMLS_scalar.setCurvatureWeightingPower(2);
294 my_GMLS_scalar.setWeightingPower(2);
297 my_GMLS_scalar.generateAlphas();
299 Kokkos::Profiling::pushRegion(
"Full Polynomial Basis GMLS Solution");
306 solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
309 my_GMLS_vector.
setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
310 std::vector<TargetOperation> lro_vector(2);
313 my_GMLS_vector.addTargets(lro_vector);
315 my_GMLS_vector.setCurvatureWeightingPower(2);
317 my_GMLS_vector.setWeightingPower(2);
318 my_GMLS_vector.generateAlphas();
319 Kokkos::Profiling::popRegion();
321 Kokkos::Profiling::pushRegion(
"Scalar Polynomial Basis Repeated to Form a Vector GMLS Solution");
346 solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
349 my_GMLS_vector_of_scalar_clones.
setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
350 std::vector<TargetOperation> lro_vector_of_scalar_clones(2);
353 my_GMLS_vector_of_scalar_clones.addTargets(lro_vector_of_scalar_clones);
355 my_GMLS_vector_of_scalar_clones.setCurvatureWeightingPower(2);
357 my_GMLS_vector_of_scalar_clones.setWeightingPower(2);
358 my_GMLS_vector_of_scalar_clones.generateAlphas();
359 Kokkos::Profiling::popRegion();
364 double instantiation_time = timer.seconds();
365 std::cout <<
"Took " << instantiation_time <<
"s to complete alphas generation." << std::endl;
367 Kokkos::Profiling::pushRegion(
"Apply Alphas to Data");
382 Evaluator scalar_gmls_evaluator(&my_GMLS_scalar);
383 Evaluator vector_gmls_evaluator(&my_GMLS_vector);
384 Evaluator vector_gmls_evaluator_of_scalar_clones(&my_GMLS_vector_of_scalar_clones);
404 auto output_vector_of_scalar_clones =
408 auto output_divergence_of_scalar_clones =
442 Kokkos::Profiling::popRegion();
444 Kokkos::Profiling::pushRegion(
"Comparison");
448 double tangent_bundle_error = 0;
449 double tangent_bundle_norm = 0;
450 double values_error = 0;
451 double values_norm = 0;
452 double laplacian_error = 0;
453 double laplacian_norm = 0;
454 double gradient_ambient_error = 0;
455 double gradient_ambient_norm = 0;
458 double vector_ambient_error = 0;
459 double vector_ambient_norm = 0;
460 double divergence_ambient_error = 0;
461 double divergence_ambient_norm = 0;
462 double vector_of_scalar_clones_ambient_error = 0;
463 double vector_of_scalar_clones_ambient_norm = 0;
464 double divergence_of_scalar_clones_ambient_error = 0;
465 double divergence_of_scalar_clones_ambient_norm = 0;
468 for (
int i=0; i<number_target_coords; i++) {
471 double GMLS_value = output_value(i);
472 double GMLS_gc = output_gc(i);
475 double GMLS_Laplacian = output_laplacian(i);
478 double xval = target_coords(i,0);
479 double yval = (dimension>1) ? target_coords(i,1) : 0;
480 double zval = (dimension>2) ? target_coords(i,2) : 0;
481 double coord[3] = {xval, yval, zval};
484 for (
int j=0; j<dimension-1; ++j) {
485 double tangent_inner_prod = 0;
486 for (
int k=0; k<dimension; ++k) {
487 tangent_inner_prod += coord[k] * my_GMLS_scalar.getTangentBundle(i, j, k);
489 tangent_bundle_error += tangent_inner_prod * tangent_inner_prod;
491 double normal_inner_prod = 0;
492 for (
int k=0; k<dimension; ++k) {
493 normal_inner_prod += coord[k] * my_GMLS_scalar.getTangentBundle(i, dimension-1, k);
496 double normal_error_to_sum = (normal_inner_prod > 0) ? normal_inner_prod - 1 : normal_inner_prod + 1;
497 tangent_bundle_error += normal_error_to_sum * normal_error_to_sum;
498 tangent_bundle_norm += 1;
503 double actual_Gradient_ambient[3] = {0,0,0};
507 values_error += (GMLS_value - actual_value)*(GMLS_value - actual_value);
508 values_norm += actual_value*actual_value;
510 laplacian_error += (GMLS_Laplacian - actual_Laplacian)*(GMLS_Laplacian - actual_Laplacian);
511 laplacian_norm += actual_Laplacian*actual_Laplacian;
513 double actual_gc = 1.0;
514 gc_error += (GMLS_gc - actual_gc)*(GMLS_gc - actual_gc);
515 gc_norm += actual_gc*actual_gc;
517 for (
int j=0; j<dimension; ++j) {
518 gradient_ambient_error += (output_gradient(i,j) - actual_Gradient_ambient[j])*(output_gradient(i,j) - actual_Gradient_ambient[j]);
519 gradient_ambient_norm += actual_Gradient_ambient[j]*actual_Gradient_ambient[j];
522 for (
int j=0; j<dimension; ++j) {
523 vector_ambient_error += (output_vector(i,j) - actual_Gradient_ambient[j])*(output_vector(i,j) - actual_Gradient_ambient[j]);
524 vector_ambient_norm += actual_Gradient_ambient[j]*actual_Gradient_ambient[j];
527 divergence_ambient_error += (output_divergence(i) - actual_Laplacian)*(output_divergence(i) - actual_Laplacian);
528 divergence_ambient_norm += actual_Laplacian*actual_Laplacian;
530 for (
int j=0; j<dimension; ++j) {
531 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]);
532 vector_of_scalar_clones_ambient_norm += actual_Gradient_ambient[j]*actual_Gradient_ambient[j];
535 divergence_of_scalar_clones_ambient_error += (output_divergence_of_scalar_clones(i) - actual_Laplacian)*(output_divergence_of_scalar_clones(i) - actual_Laplacian);
536 divergence_of_scalar_clones_ambient_norm += actual_Laplacian*actual_Laplacian;
540 tangent_bundle_error /= number_target_coords;
541 tangent_bundle_error = std::sqrt(tangent_bundle_error);
542 tangent_bundle_norm /= number_target_coords;
543 tangent_bundle_norm = std::sqrt(tangent_bundle_norm);
545 values_error /= number_target_coords;
546 values_error = std::sqrt(values_error);
547 values_norm /= number_target_coords;
548 values_norm = std::sqrt(values_norm);
550 laplacian_error /= number_target_coords;
551 laplacian_error = std::sqrt(laplacian_error);
552 laplacian_norm /= number_target_coords;
553 laplacian_norm = std::sqrt(laplacian_norm);
555 gradient_ambient_error /= number_target_coords;
556 gradient_ambient_error = std::sqrt(gradient_ambient_error);
557 gradient_ambient_norm /= number_target_coords;
558 gradient_ambient_norm = std::sqrt(gradient_ambient_norm);
560 gc_error /= number_target_coords;
561 gc_error = std::sqrt(gc_error);
562 gc_norm /= number_target_coords;
563 gc_norm = std::sqrt(gc_norm);
565 vector_ambient_error /= number_target_coords;
566 vector_ambient_error = std::sqrt(vector_ambient_error);
567 vector_ambient_norm /= number_target_coords;
568 vector_ambient_norm = std::sqrt(vector_ambient_norm);
570 divergence_ambient_error /= number_target_coords;
571 divergence_ambient_error = std::sqrt(divergence_ambient_error);
572 divergence_ambient_norm /= number_target_coords;
573 divergence_ambient_norm = std::sqrt(divergence_ambient_norm);
575 vector_of_scalar_clones_ambient_error /= number_target_coords;
576 vector_of_scalar_clones_ambient_error = std::sqrt(vector_of_scalar_clones_ambient_error);
577 vector_of_scalar_clones_ambient_norm /= number_target_coords;
578 vector_of_scalar_clones_ambient_norm = std::sqrt(vector_of_scalar_clones_ambient_norm);
580 divergence_of_scalar_clones_ambient_error /= number_target_coords;
581 divergence_of_scalar_clones_ambient_error = std::sqrt(divergence_of_scalar_clones_ambient_error);
582 divergence_of_scalar_clones_ambient_norm /= number_target_coords;
583 divergence_of_scalar_clones_ambient_norm = std::sqrt(divergence_of_scalar_clones_ambient_norm);
585 printf(
"Tangent Bundle Error: %g\n", tangent_bundle_error / tangent_bundle_norm);
586 printf(
"Point Value Error: %g\n", values_error / values_norm);
587 printf(
"Laplace-Beltrami Error: %g\n", laplacian_error / laplacian_norm);
588 printf(
"Gaussian Curvature Error: %g\n", gc_error / gc_norm);
589 printf(
"Surface Gradient (Ambient) Error: %g\n", gradient_ambient_error / gradient_ambient_norm);
590 printf(
"Surface Vector (VectorBasis) Error: %g\n", vector_ambient_error / vector_ambient_norm);
591 printf(
"Surface Divergence (VectorBasis) Error: %g\n", divergence_ambient_error / divergence_ambient_norm);
592 printf(
"Surface Vector (ScalarClones) Error: %g\n",
593 vector_of_scalar_clones_ambient_error / vector_of_scalar_clones_ambient_norm);
594 printf(
"Surface Divergence (ScalarClones) Error: %g\n",
595 divergence_of_scalar_clones_ambient_error / divergence_of_scalar_clones_ambient_norm);
599 Kokkos::Profiling::popRegion();
608 #ifdef COMPADRE_USE_MPI
int main(int argc, char *args[])
[Parse Command Line Arguments]
KOKKOS_INLINE_FUNCTION double sphere_harmonic54(double x, double y, double z)
Lightweight Evaluator Helper This class is a lightweight wrapper for extracting and applying all rele...
constexpr SamplingFunctional ManifoldVectorPointSample
Point evaluations of the entire vector source function (but on a manifold, so it includes a transform...
Point evaluation of a vector (reconstructs entire vector at once, requiring a ReconstructionSpace hav...
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...
Point evaluation of Gaussian curvature.
Point evaluation of a scalar.
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...
Scalar basis reused as many times as there are components in the vector resulting in a much cheaper p...
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)
Point evaluation of the laplacian of a scalar (could be on a manifold or not)
Point evaluation of the divergence of a vector (results in a scalar)
KOKKOS_INLINE_FUNCTION void gradient_sphereHarmonic54_ambient(double *gradient, double x, double y, double z)
Point evaluation of the gradient of a scalar.
Generalized Moving Least Squares (GMLS)
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) ...
KOKKOS_INLINE_FUNCTION double laplace_beltrami_sphere_harmonic54(double x, double y, double z)
std::string constraint_name
Vector polynomial basis having # of components _dimensions, or (_dimensions-1) in the case of manifol...