9 #include <Compadre_Config.h>
17 #ifdef COMPADRE_USE_MPI
21 #include <Kokkos_Timer.hpp>
22 #include <Kokkos_Core.hpp>
24 using namespace Compadre;
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::ceil(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)(1.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::ceil(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);
154 Kokkos::Profiling::popRegion();
156 Kokkos::Profiling::pushRegion(
"Creating Data");
162 Kokkos::deep_copy(source_coords_device, source_coords);
165 Kokkos::deep_copy(target_coords_device, target_coords);
172 Kokkos::View<double*, Kokkos::DefaultExecutionSpace> sampling_data_device(
"samples of true solution",
173 source_coords_device.extent(0));
175 Kokkos::View<double*, Kokkos::DefaultExecutionSpace> ones_data_device(
"samples of ones",
176 source_coords_device.extent(0));
177 Kokkos::deep_copy(ones_data_device, 1.0);
180 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> sampling_vector_data_device(
"samples of vector true solution",
181 source_coords_device.extent(0), 3);
183 Kokkos::parallel_for(
"Sampling Manufactured Solutions", Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>
184 (0,source_coords.extent(0)), KOKKOS_LAMBDA(
const int i) {
187 double xval = source_coords_device(i,0);
188 double yval = (dimension>1) ? source_coords_device(i,1) : 0;
189 double zval = (dimension>2) ? source_coords_device(i,2) : 0;
194 for (
int j=0; j<3; ++j) {
195 double gradient[3] = {0,0,0};
197 sampling_vector_data_device(i,j) = gradient[j];
204 Kokkos::Profiling::popRegion();
205 Kokkos::Profiling::pushRegion(
"Neighbor Search");
209 double epsilon_multiplier = 1.9;
215 Kokkos::View<int*> neighbor_lists_device(
"neighbor lists",
217 Kokkos::View<int*>::HostMirror neighbor_lists = Kokkos::create_mirror_view(neighbor_lists_device);
220 Kokkos::View<int*> number_of_neighbors_list_device(
"number of neighbor lists",
221 number_target_coords);
222 Kokkos::View<int*>::HostMirror number_of_neighbors_list = Kokkos::create_mirror_view(number_of_neighbors_list_device);
225 Kokkos::View<double*, Kokkos::DefaultExecutionSpace> epsilon_device(
"h supports", number_target_coords);
226 Kokkos::View<double*>::HostMirror epsilon = Kokkos::create_mirror_view(epsilon_device);
228 size_t storage_size = point_cloud_search.generateCRNeighborListsFromKNNSearch(
true , target_coords, neighbor_lists,
229 number_of_neighbors_list, epsilon, min_neighbors, epsilon_multiplier);
232 Kokkos::resize(neighbor_lists_device, storage_size);
233 neighbor_lists = Kokkos::create_mirror_view(neighbor_lists_device);
236 point_cloud_search.generateCRNeighborListsFromKNNSearch(
false , target_coords, neighbor_lists,
237 number_of_neighbors_list, epsilon, min_neighbors, epsilon_multiplier);
251 Kokkos::Profiling::popRegion();
262 Kokkos::deep_copy(neighbor_lists_device, neighbor_lists);
263 Kokkos::deep_copy(number_of_neighbors_list_device, number_of_neighbors_list);
264 Kokkos::deep_copy(epsilon_device, epsilon);
268 GMLS my_GMLS_scalar(order, dimension,
269 solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
286 my_GMLS_scalar.
setProblemData(neighbor_lists_device, number_of_neighbors_list_device, source_coords_device, target_coords_device, epsilon_device);
291 my_GMLS_scalar.setReferenceOutwardNormalDirection(target_coords_device,
true );
294 std::vector<TargetOperation> lro_scalar(4);
301 my_GMLS_scalar.addTargets(lro_scalar);
307 my_GMLS_scalar.setCurvatureWeightingParameter(2);
313 my_GMLS_scalar.setWeightingParameter(2);
316 my_GMLS_scalar.generateAlphas();
318 Kokkos::Profiling::pushRegion(
"Full Polynomial Basis GMLS Solution");
325 solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
328 my_GMLS_vector.
setProblemData(neighbor_lists_device, number_of_neighbors_list_device, source_coords_device, target_coords_device, epsilon_device);
329 std::vector<TargetOperation> lro_vector(2);
332 my_GMLS_vector.addTargets(lro_vector);
334 my_GMLS_vector.setCurvatureWeightingParameter(2);
336 my_GMLS_vector.setWeightingParameter(2);
337 my_GMLS_vector.generateAlphas();
338 Kokkos::Profiling::popRegion();
340 Kokkos::Profiling::pushRegion(
"Scalar Polynomial Basis Repeated to Form a Vector GMLS Solution");
365 solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
368 my_GMLS_vector_of_scalar_clones.
setProblemData(neighbor_lists_device, number_of_neighbors_list_device, source_coords_device, target_coords_device, epsilon_device);
369 std::vector<TargetOperation> lro_vector_of_scalar_clones(2);
372 my_GMLS_vector_of_scalar_clones.addTargets(lro_vector_of_scalar_clones);
374 my_GMLS_vector_of_scalar_clones.setCurvatureWeightingParameter(2);
376 my_GMLS_vector_of_scalar_clones.setWeightingParameter(2);
377 my_GMLS_vector_of_scalar_clones.generateAlphas();
378 Kokkos::Profiling::popRegion();
383 double instantiation_time = timer.seconds();
384 std::cout <<
"Took " << instantiation_time <<
"s to complete alphas generation." << std::endl;
386 Kokkos::Profiling::pushRegion(
"Apply Alphas to Data");
401 Evaluator scalar_gmls_evaluator(&my_GMLS_scalar);
402 Evaluator vector_gmls_evaluator(&my_GMLS_vector);
403 Evaluator vector_gmls_evaluator_of_scalar_clones(&my_GMLS_vector_of_scalar_clones);
423 auto output_vector_of_scalar_clones =
427 auto output_divergence_of_scalar_clones =
461 Kokkos::Profiling::popRegion();
463 Kokkos::Profiling::pushRegion(
"Comparison");
467 double tangent_bundle_error = 0;
468 double tangent_bundle_norm = 0;
469 double values_error = 0;
470 double values_norm = 0;
471 double laplacian_error = 0;
472 double laplacian_norm = 0;
473 double gradient_ambient_error = 0;
474 double gradient_ambient_norm = 0;
477 double vector_ambient_error = 0;
478 double vector_ambient_norm = 0;
479 double divergence_ambient_error = 0;
480 double divergence_ambient_norm = 0;
481 double vector_of_scalar_clones_ambient_error = 0;
482 double vector_of_scalar_clones_ambient_norm = 0;
483 double divergence_of_scalar_clones_ambient_error = 0;
484 double divergence_of_scalar_clones_ambient_norm = 0;
487 for (
int i=0; i<number_target_coords; i++) {
490 double GMLS_value = output_value(i);
491 double GMLS_gc = output_gc(i);
494 double GMLS_Laplacian = output_laplacian(i);
497 double xval = target_coords(i,0);
498 double yval = (dimension>1) ? target_coords(i,1) : 0;
499 double zval = (dimension>2) ? target_coords(i,2) : 0;
500 double coord[3] = {xval, yval, zval};
503 for (
unsigned int j=0; j<static_cast<unsigned int>(dimension-1); ++j) {
506 double tangent_inner_prod = 0;
507 for (
int k=0; k<dimension; ++k) {
508 tangent_inner_prod += coord[k] * my_GMLS_scalar.getTangentBundle(i, j, k);
510 tangent_bundle_error += tangent_inner_prod * tangent_inner_prod;
512 double normal_inner_prod = 0;
513 for (
int k=0; k<dimension; ++k) {
514 normal_inner_prod += coord[k] * my_GMLS_scalar.getTangentBundle(i, dimension-1, k);
517 double normal_error_to_sum = (normal_inner_prod > 0) ? normal_inner_prod - 1 : normal_inner_prod + 1;
518 tangent_bundle_error += normal_error_to_sum * normal_error_to_sum;
519 tangent_bundle_norm += 1;
524 double actual_Gradient_ambient[3] = {0,0,0};
528 values_error += (GMLS_value - actual_value)*(GMLS_value - actual_value);
529 values_norm += actual_value*actual_value;
531 laplacian_error += (GMLS_Laplacian - actual_Laplacian)*(GMLS_Laplacian - actual_Laplacian);
532 laplacian_norm += actual_Laplacian*actual_Laplacian;
534 double actual_gc = 1.0;
535 gc_error += (GMLS_gc - actual_gc)*(GMLS_gc - actual_gc);
536 gc_norm += actual_gc*actual_gc;
538 for (
int j=0; j<dimension; ++j) {
539 gradient_ambient_error += (output_gradient(i,j) - actual_Gradient_ambient[j])*(output_gradient(i,j) - actual_Gradient_ambient[j]);
540 gradient_ambient_norm += actual_Gradient_ambient[j]*actual_Gradient_ambient[j];
543 for (
int j=0; j<dimension; ++j) {
544 vector_ambient_error += (output_vector(i,j) - actual_Gradient_ambient[j])*(output_vector(i,j) - actual_Gradient_ambient[j]);
545 vector_ambient_norm += actual_Gradient_ambient[j]*actual_Gradient_ambient[j];
548 divergence_ambient_error += (output_divergence(i) - actual_Laplacian)*(output_divergence(i) - actual_Laplacian);
549 divergence_ambient_norm += actual_Laplacian*actual_Laplacian;
551 for (
int j=0; j<dimension; ++j) {
552 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]);
553 vector_of_scalar_clones_ambient_norm += actual_Gradient_ambient[j]*actual_Gradient_ambient[j];
556 divergence_of_scalar_clones_ambient_error += (output_divergence_of_scalar_clones(i) - actual_Laplacian)*(output_divergence_of_scalar_clones(i) - actual_Laplacian);
557 divergence_of_scalar_clones_ambient_norm += actual_Laplacian*actual_Laplacian;
561 tangent_bundle_error /= number_target_coords;
562 tangent_bundle_error = std::sqrt(tangent_bundle_error);
563 tangent_bundle_norm /= number_target_coords;
564 tangent_bundle_norm = std::sqrt(tangent_bundle_norm);
566 values_error /= number_target_coords;
567 values_error = std::sqrt(values_error);
568 values_norm /= number_target_coords;
569 values_norm = std::sqrt(values_norm);
571 laplacian_error /= number_target_coords;
572 laplacian_error = std::sqrt(laplacian_error);
573 laplacian_norm /= number_target_coords;
574 laplacian_norm = std::sqrt(laplacian_norm);
576 gradient_ambient_error /= number_target_coords;
577 gradient_ambient_error = std::sqrt(gradient_ambient_error);
578 gradient_ambient_norm /= number_target_coords;
579 gradient_ambient_norm = std::sqrt(gradient_ambient_norm);
581 gc_error /= number_target_coords;
582 gc_error = std::sqrt(gc_error);
583 gc_norm /= number_target_coords;
584 gc_norm = std::sqrt(gc_norm);
586 vector_ambient_error /= number_target_coords;
587 vector_ambient_error = std::sqrt(vector_ambient_error);
588 vector_ambient_norm /= number_target_coords;
589 vector_ambient_norm = std::sqrt(vector_ambient_norm);
591 divergence_ambient_error /= number_target_coords;
592 divergence_ambient_error = std::sqrt(divergence_ambient_error);
593 divergence_ambient_norm /= number_target_coords;
594 divergence_ambient_norm = std::sqrt(divergence_ambient_norm);
596 vector_of_scalar_clones_ambient_error /= number_target_coords;
597 vector_of_scalar_clones_ambient_error = std::sqrt(vector_of_scalar_clones_ambient_error);
598 vector_of_scalar_clones_ambient_norm /= number_target_coords;
599 vector_of_scalar_clones_ambient_norm = std::sqrt(vector_of_scalar_clones_ambient_norm);
601 divergence_of_scalar_clones_ambient_error /= number_target_coords;
602 divergence_of_scalar_clones_ambient_error = std::sqrt(divergence_of_scalar_clones_ambient_error);
603 divergence_of_scalar_clones_ambient_norm /= number_target_coords;
604 divergence_of_scalar_clones_ambient_norm = std::sqrt(divergence_of_scalar_clones_ambient_norm);
606 printf(
"Tangent Bundle Error: %g\n", tangent_bundle_error / tangent_bundle_norm);
607 printf(
"Point Value Error: %g\n", values_error / values_norm);
608 printf(
"Laplace-Beltrami Error: %g\n", laplacian_error / laplacian_norm);
609 printf(
"Gaussian Curvature Error: %g\n", gc_error / gc_norm);
610 printf(
"Surface Gradient (Ambient) Error: %g\n", gradient_ambient_error / gradient_ambient_norm);
611 printf(
"Surface Vector (VectorBasis) Error: %g\n", vector_ambient_error / vector_ambient_norm);
612 printf(
"Surface Divergence (VectorBasis) Error: %g\n", divergence_ambient_error / divergence_ambient_norm);
613 printf(
"Surface Vector (ScalarClones) Error: %g\n",
614 vector_of_scalar_clones_ambient_error / vector_of_scalar_clones_ambient_norm);
615 printf(
"Surface Divergence (ScalarClones) Error: %g\n",
616 divergence_of_scalar_clones_ambient_error / divergence_of_scalar_clones_ambient_norm);
620 Kokkos::Profiling::popRegion();
629 #ifdef COMPADRE_USE_MPI
Point evaluation of Gaussian curvature.
KOKKOS_INLINE_FUNCTION double sphere_harmonic54(double x, double y, double z)
Point evaluation of a scalar.
Scalar basis reused as many times as there are components in the vector resulting in a much cheaper p...
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 the gradient of a scalar.
Point evaluation of the laplacian of a scalar (could be on a manifold or not)
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...
int main(int argc, char **argv)
Vector polynomial basis having # of components _dimensions, or (_dimensions-1) in the case of manifol...
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...
import subprocess import os import re import math import sys import argparse d d d
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)
Generalized Moving Least Squares (GMLS)
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)
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
Point evaluation of a vector (reconstructs entire vector at once, requiring a ReconstructionSpace hav...