9 #include <Compadre_Config.h> 17 #ifdef COMPADRE_USE_MPI 21 #include <Kokkos_Timer.hpp> 22 #include <Kokkos_Core.hpp> 42 double local_vec[3] = {0,0};
43 for (
int j=0; j<3; ++j) {
44 local_vec[0] += T(0,j) * u[j];
45 local_vec[1] += T(1,j) * u[j];
48 for (
int j=0; j<3; ++j) u[j] = 0;
49 for (
int j=0; j<3; ++j) {
50 u[j] += P(0, j) * local_vec[0];
51 u[j] += P(1, j) * local_vec[1];
59 int main (
int argc,
char* args[]) {
62 #ifdef COMPADRE_USE_MPI 63 MPI_Init(&argc, &args);
67 Kokkos::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);
300 my_GMLS_scalar.setAdditionalEvaluationSitesData(neighbor_lists_device, source_coords_device);
305 my_GMLS_scalar.setReferenceOutwardNormalDirection(target_coords_device,
true );
308 std::vector<TargetOperation> lro_scalar(3);
314 my_GMLS_scalar.addTargets(lro_scalar);
320 my_GMLS_scalar.setCurvatureWeightingPower(2);
326 my_GMLS_scalar.setWeightingPower(2);
329 my_GMLS_scalar.generateAlphas();
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);
342 my_GMLS_vector.setAdditionalEvaluationSitesData(neighbor_lists_device, source_coords_device);
343 std::vector<TargetOperation> lro_vector(2);
346 my_GMLS_vector.addTargets(lro_vector);
348 my_GMLS_vector.setCurvatureWeightingPower(2);
350 my_GMLS_vector.setWeightingPower(2);
351 my_GMLS_vector.generateAlphas();
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);
383 my_GMLS_vector_of_scalar_clones.setAdditionalEvaluationSitesData(neighbor_lists_device, source_coords_device);
384 std::vector<TargetOperation> lro_vector_of_scalar_clones(2);
387 my_GMLS_vector_of_scalar_clones.addTargets(lro_vector_of_scalar_clones);
389 my_GMLS_vector_of_scalar_clones.setCurvatureWeightingPower(2);
391 my_GMLS_vector_of_scalar_clones.setWeightingPower(2);
392 my_GMLS_vector_of_scalar_clones.generateAlphas();
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;
516 auto d_prestencil_weights = my_GMLS_vector_of_scalar_clones.getPrestencilWeights();
517 auto prestencil_weights = Kokkos::create_mirror_view(d_prestencil_weights);
518 Kokkos::deep_copy(prestencil_weights, d_prestencil_weights);
521 auto d_tangent_directions = my_GMLS_vector_of_scalar_clones.getTangentDirections();
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
KOKKOS_INLINE_FUNCTION double sphere_harmonic54(double x, double y, double z)
Point evaluation of the curl of a vector (results in a vector)
Lightweight Evaluator Helper This class is a lightweight wrapper for extracting and applying all rele...
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...
constexpr SamplingFunctional VaryingManifoldVectorPointSample
For integrating polynomial dotted with normal over an edge.
KOKKOS_INLINE_FUNCTION void curl_sphere_harmonic54(double *curl, double x, double y, double z)
Point evaluation of Gaussian curvature.
void AmbientLocalAmbient(XYZ &u, double *T_data, double *P_data)
[Ambient to Local Back To Ambient Helper Function]
#define TO_GLOBAL(variable)
Point evaluation of a scalar.
Kokkos::View< double **, layout_right, host_execution_space, Kokkos::MemoryTraits< Kokkos::Unmanaged > > host_scratch_matrix_right_type
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)
int main(int argc, char *args[])
[Ambient to Local Back To Ambient Helper Function]
KOKKOS_INLINE_FUNCTION void gradient_sphereHarmonic54_ambient(double *gradient, double x, double y, double z)
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) ...
constexpr SamplingFunctional PointSample
Available sampling functionals.
std::string constraint_name
Vector polynomial basis having # of components _dimensions, or (_dimensions-1) in the case of manifol...