9#include <Compadre_Config.h>
17#ifdef COMPADRE_USE_MPI
21#include <Kokkos_Timer.hpp>
22#include <Kokkos_Core.hpp>
29int main (
int argc,
char* args[]) {
32#ifdef COMPADRE_USE_MPI
33MPI_Init(&argc, &args);
37Kokkos::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);
222 GMLS my_GMLS_vector_1(ReconstructionSpace::VectorTaylorPolynomial,
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);
272 GMLS my_GMLS_vector_2(ReconstructionSpace::VectorTaylorPolynomial,
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);
295 GMLS my_GMLS_scalar(ReconstructionSpace::ScalarTaylorPolynomial,
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);
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)
KOKKOS_INLINE_FUNCTION void gradient_sphereHarmonic54_ambient(double *gradient, double x, double y, double z)
KOKKOS_INLINE_FUNCTION double laplace_beltrami_sphere_harmonic54(double x, double y, double z)
int main(int argc, char *args[])
[Parse Command Line Arguments]
Lightweight Evaluator Helper This class is a lightweight wrapper for extracting and applying all rele...
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)
Generalized Moving Least Squares (GMLS)
void addTargets(TargetOperation lro)
Adds a target to the vector of target functional to be applied to the reconstruction.
void setDimensionOfQuadraturePoints(int dim)
Dimensions of quadrature points to use.
void setCurvatureWeightingParameter(int wp, int index=0)
Parameter for weighting kernel for curvature index = 0 sets p paramater for weighting kernel index = ...
void setWeightingParameter(int wp, int index=0)
Parameter for weighting kernel for GMLS problem index = 0 sets p paramater for weighting kernel index...
void generateAlphas(const int number_of_batches=1, const bool keep_coefficients=false, const bool clear_cache=true)
Meant to calculate target operations and apply the evaluations to the previously constructed polynomi...
void setOrderOfQuadraturePoints(int order)
Number quadrature points to use.
void setQuadratureType(std::string quadrature_type)
Type of quadrature points.
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)
void setWeightingType(const std::string &wt)
Type for weighting kernel for GMLS problem.
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....
void setCurvatureWeightingType(const std::string &wt)
Type for weighting kernel for curvature.
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 StaggeredEdgeIntegralSample
Samples consist of the result of integrals of a vector dotted with the tangent along edges between ne...
@ ChainedStaggeredLaplacianOfScalarPointEvaluation
Point evaluation of the chained staggered Laplacian acting on VectorTaylorPolynomial basis + Staggere...
@ DivergenceOfVectorPointEvaluation
Point evaluation of the divergence of a vector (results in a scalar)
constexpr SamplingFunctional StaggeredEdgeAnalyticGradientIntegralSample
Analytical integral of a gradient source vector is just a difference of the scalar source at neighbor...
std::string constraint_name