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;
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.setCurvatureWeightingPower(2);
261 my_GMLS_vector_1.setWeightingPower(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.setCurvatureWeightingPower(2);
288 my_GMLS_vector_2.setWeightingPower(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.setCurvatureWeightingPower(2);
310 my_GMLS_scalar.setWeightingPower(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...
constexpr SamplingFunctional StaggeredEdgeAnalyticGradientIntegralSample
Analytical integral of a gradient source vector is just a difference of the scalar source at neighbor...
Scalar polynomial basis centered at the target site and scaled by sum of basis powers e...
int main(int argc, char *args[])
[Parse Command Line Arguments]
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...
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 chained staggered Laplacian acting on VectorTaylorPolynomial basis + Staggere...
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)
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...