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;
75 double a = 4*
PI*r*r/N_pts_on_sphere;
76 double d = std::sqrt(a);
77 int M_theta = std::round(
PI/d);
78 double d_theta =
PI/M_theta;
79 double d_phi = a/d_theta;
80 for (
int i=0; i<M_theta; ++i) {
81 double theta =
PI*(i + 0.5)/M_theta;
82 int M_phi = std::round(2*
PI*std::sin(theta)/d_phi);
83 for (
int j=0; j<M_phi; ++j) {
84 double phi = 2*
PI*j/M_phi;
85 source_coords(N_count, 0) = theta;
86 source_coords(N_count, 1) = phi;
91 for (
int i=0; i<N_count; ++i) {
92 double theta = source_coords(i,0);
93 double phi = source_coords(i,1);
94 source_coords(i,0) = r*std::sin(theta)*std::cos(phi);
95 source_coords(i,1) = r*std::sin(theta)*std::sin(phi);
96 source_coords(i,2) = r*cos(theta);
99 number_source_coords = N_count;
103 Kokkos::resize(source_coords, number_source_coords, 3);
104 Kokkos::resize(source_coords_device, number_source_coords, 3);
107 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> target_coords_device(
"target coordinates",
108 number_target_coords, 3);
109 Kokkos::View<double**>::HostMirror target_coords = Kokkos::create_mirror_view(target_coords_device);
112 std::mt19937 rng(50);
115 std::uniform_int_distribution<int> gen_num_neighbors(0, number_source_coords-1);
118 for (
int i=0; i<number_target_coords; ++i) {
119 const int source_site_to_copy = gen_num_neighbors(rng);
120 for (
int j=0; j<3; ++j) {
121 target_coords(i,j) = source_coords(source_site_to_copy,j);
128 Kokkos::Profiling::popRegion();
130 Kokkos::Profiling::pushRegion(
"Creating Data");
136 Kokkos::deep_copy(source_coords_device, source_coords);
137 Kokkos::deep_copy(target_coords_device, target_coords);
144 Kokkos::View<double*, Kokkos::DefaultExecutionSpace> sampling_data_device(
"samples of true solution",
145 source_coords_device.extent(0));
148 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> sampling_vector_data_device(
"samples of vector true solution",
149 source_coords_device.extent(0), 3);
151 Kokkos::parallel_for(
"Sampling Manufactured Solutions", Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>
152 (0,source_coords.extent(0)), KOKKOS_LAMBDA(
const int i) {
155 double xval = source_coords_device(i,0);
156 double yval = (dimension>1) ? source_coords_device(i,1) : 0;
157 double zval = (dimension>2) ? source_coords_device(i,2) : 0;
163 for (
int j=0; j<3; ++j) {
164 double gradient[3] = {0,0,0};
166 sampling_vector_data_device(i,j) = gradient[j];
174 Kokkos::Profiling::popRegion();
175 Kokkos::Profiling::pushRegion(
"Neighbor Search");
186 double epsilon_multiplier = 1.5;
187 int estimated_upper_bound_number_neighbors =
188 point_cloud_search.getEstimatedNumberNeighborsUpperBound(min_neighbors, dimension, epsilon_multiplier);
190 Kokkos::View<int**, Kokkos::DefaultExecutionSpace> neighbor_lists_device(
"neighbor lists",
191 number_target_coords, estimated_upper_bound_number_neighbors);
192 Kokkos::View<int**>::HostMirror neighbor_lists = Kokkos::create_mirror_view(neighbor_lists_device);
195 Kokkos::View<double*, Kokkos::DefaultExecutionSpace> epsilon_device(
"h supports", number_target_coords);
196 Kokkos::View<double*>::HostMirror epsilon = Kokkos::create_mirror_view(epsilon_device);
201 point_cloud_search.generate2DNeighborListsFromKNNSearch(
false , target_coords, neighbor_lists,
202 epsilon, min_neighbors, epsilon_multiplier);
207 Kokkos::Profiling::popRegion();
218 Kokkos::deep_copy(neighbor_lists_device, neighbor_lists);
219 Kokkos::deep_copy(epsilon_device, epsilon);
225 solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
242 my_GMLS_vector_1.
setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
245 std::vector<TargetOperation> lro_vector_1(1);
249 my_GMLS_vector_1.addTargets(lro_vector_1);
255 my_GMLS_vector_1.setCurvatureWeightingParameter(2);
261 my_GMLS_vector_1.setWeightingParameter(2);
264 my_GMLS_vector_1.setOrderOfQuadraturePoints(2);
265 my_GMLS_vector_1.setDimensionOfQuadraturePoints(1);
266 my_GMLS_vector_1.setQuadratureType(
"LINE");
269 my_GMLS_vector_1.generateAlphas();
276 solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
279 my_GMLS_vector_2.
setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
280 std::vector<TargetOperation> lro_vector_2(2);
284 my_GMLS_vector_2.addTargets(lro_vector_2);
286 my_GMLS_vector_2.setCurvatureWeightingParameter(2);
288 my_GMLS_vector_2.setWeightingParameter(2);
289 my_GMLS_vector_2.setOrderOfQuadraturePoints(2);
290 my_GMLS_vector_2.setDimensionOfQuadraturePoints(1);
291 my_GMLS_vector_2.setQuadratureType(
"LINE");
292 my_GMLS_vector_2.generateAlphas();
298 solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
301 my_GMLS_scalar.
setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
303 std::vector<TargetOperation> lro_scalar(1);
306 my_GMLS_scalar.addTargets(lro_scalar);
308 my_GMLS_scalar.setCurvatureWeightingParameter(2);
310 my_GMLS_scalar.setWeightingParameter(2);
311 my_GMLS_scalar.generateAlphas();
316 double instantiation_time = timer.seconds();
317 std::cout <<
"Took " << instantiation_time <<
"s to complete alphas generation." << std::endl;
319 Kokkos::Profiling::pushRegion(
"Apply Alphas to Data");
334 Evaluator vector_1_gmls_evaluator(&my_GMLS_vector_1);
335 Evaluator vector_2_gmls_evaluator(&my_GMLS_vector_2);
336 Evaluator scalar_gmls_evaluator(&my_GMLS_scalar);
347 auto output_divergence_vectorsamples =
351 auto output_divergence_scalarsamples =
355 auto output_laplacian_vectorbasis =
359 auto output_laplacian_scalarbasis =
367 Kokkos::Profiling::popRegion();
369 Kokkos::Profiling::pushRegion(
"Comparison");
373 double laplacian_vectorbasis_error = 0;
374 double laplacian_vectorbasis_norm = 0;
376 double laplacian_scalarbasis_error = 0;
377 double laplacian_scalarbasis_norm = 0;
379 double gradient_vectorbasis_ambient_error = 0;
380 double gradient_vectorbasis_ambient_norm = 0;
382 double gradient_scalarbasis_ambient_error = 0;
383 double gradient_scalarbasis_ambient_norm = 0;
385 double divergence_vectorsamples_ambient_error = 0;
386 double divergence_vectorsamples_ambient_norm = 0;
388 double divergence_scalarsamples_ambient_error = 0;
389 double divergence_scalarsamples_ambient_norm = 0;
392 for (
int i=0; i<number_target_coords; i++) {
395 double xval = target_coords(i,0);
396 double yval = (dimension>1) ? target_coords(i,1) : 0;
397 double zval = (dimension>2) ? target_coords(i,2) : 0;
401 double actual_Gradient_ambient[3] = {0,0,0};
404 laplacian_vectorbasis_error += (output_laplacian_vectorbasis(i) - actual_Laplacian)*(output_laplacian_vectorbasis(i) - actual_Laplacian);
405 laplacian_vectorbasis_norm += actual_Laplacian*actual_Laplacian;
408 laplacian_scalarbasis_error += (output_laplacian_scalarbasis(i) - actual_Laplacian)*(output_laplacian_scalarbasis(i) - actual_Laplacian);
409 laplacian_scalarbasis_norm += actual_Laplacian*actual_Laplacian;
424 divergence_vectorsamples_ambient_error += (output_divergence_vectorsamples(i) - actual_Laplacian)*(output_divergence_vectorsamples(i) - actual_Laplacian);
425 divergence_vectorsamples_ambient_norm += actual_Laplacian*actual_Laplacian;
427 divergence_scalarsamples_ambient_error += (output_divergence_scalarsamples(i) - actual_Laplacian)*(output_divergence_scalarsamples(i) - actual_Laplacian);
428 divergence_scalarsamples_ambient_norm += actual_Laplacian*actual_Laplacian;
432 laplacian_vectorbasis_error /= number_target_coords;
433 laplacian_vectorbasis_error = std::sqrt(laplacian_vectorbasis_error);
434 laplacian_vectorbasis_norm /= number_target_coords;
435 laplacian_vectorbasis_norm = std::sqrt(laplacian_vectorbasis_norm);
437 laplacian_scalarbasis_error /= number_target_coords;
438 laplacian_scalarbasis_error = std::sqrt(laplacian_scalarbasis_error);
439 laplacian_scalarbasis_norm /= number_target_coords;
440 laplacian_scalarbasis_norm = std::sqrt(laplacian_scalarbasis_norm);
442 gradient_vectorbasis_ambient_error /= number_target_coords;
443 gradient_vectorbasis_ambient_error = std::sqrt(gradient_vectorbasis_ambient_error);
444 gradient_vectorbasis_ambient_norm /= number_target_coords;
445 gradient_vectorbasis_ambient_norm = std::sqrt(gradient_vectorbasis_ambient_norm);
447 gradient_scalarbasis_ambient_error /= number_target_coords;
448 gradient_scalarbasis_ambient_error = std::sqrt(gradient_scalarbasis_ambient_error);
449 gradient_scalarbasis_ambient_norm /= number_target_coords;
450 gradient_scalarbasis_ambient_norm = std::sqrt(gradient_scalarbasis_ambient_norm);
452 divergence_vectorsamples_ambient_error /= number_target_coords;
453 divergence_vectorsamples_ambient_error = std::sqrt(divergence_vectorsamples_ambient_error);
454 divergence_vectorsamples_ambient_norm /= number_target_coords;
455 divergence_vectorsamples_ambient_norm = std::sqrt(divergence_vectorsamples_ambient_norm);
457 divergence_scalarsamples_ambient_error /= number_target_coords;
458 divergence_scalarsamples_ambient_error = std::sqrt(divergence_scalarsamples_ambient_error);
459 divergence_scalarsamples_ambient_norm /= number_target_coords;
460 divergence_scalarsamples_ambient_norm = std::sqrt(divergence_scalarsamples_ambient_norm);
462 printf(
"Staggered Laplace-Beltrami (VectorBasis) Error: %g\n", laplacian_vectorbasis_error / laplacian_vectorbasis_norm);
463 printf(
"Staggered Laplace-Beltrami (ScalarBasis) Error: %g\n", laplacian_scalarbasis_error / laplacian_scalarbasis_norm);
464 printf(
"Surface Staggered Gradient (VectorBasis) Error: %g\n", gradient_vectorbasis_ambient_error / gradient_vectorbasis_ambient_norm);
465 printf(
"Surface Staggered Gradient (ScalarBasis) Error: %g\n", gradient_scalarbasis_ambient_error / gradient_scalarbasis_ambient_norm);
466 printf(
"Surface Staggered Divergence (VectorSamples) Error: %g\n", divergence_vectorsamples_ambient_error / divergence_vectorsamples_ambient_norm);
467 printf(
"Surface Staggered Divergence (ScalarSamples) Error: %g\n", divergence_scalarsamples_ambient_error / divergence_scalarsamples_ambient_norm);
471 Kokkos::Profiling::popRegion();
480 #ifdef COMPADRE_USE_MPI
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...
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...
constexpr SamplingFunctional StaggeredEdgeAnalyticGradientIntegralSample
Analytical integral of a gradient source vector is just a difference of the scalar source at neighbor...
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...
constexpr SamplingFunctional StaggeredEdgeIntegralSample
Samples consist of the result of integrals of a vector dotted with the tangent along edges between ne...
Scalar polynomial basis centered at the target site and scaled by sum of basis powers e...
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 the chained staggered Laplacian acting on VectorTaylorPolynomial basis + Staggere...