9 #include <Compadre_Config.h> 16 #ifdef COMPADRE_USE_MPI 20 #include <Kokkos_Timer.hpp> 21 #include <Kokkos_Core.hpp> 28 int main (
int argc,
char* args[]) {
31 #ifdef COMPADRE_USE_MPI 32 MPI_Init(&argc, &args);
36 Kokkos::initialize(argc, args);
45 int constraint_type = 0;
47 int arg8toi = atoi(args[7]);
49 constraint_type = arg8toi;
59 int arg7toi = atoi(args[6]);
61 problem_type = arg7toi;
72 int arg6toi = atoi(args[5]);
74 solver_type = arg6toi;
80 int N_pts_on_sphere = 1000;
82 int arg5toi = atoi(args[4]);
84 N_pts_on_sphere = arg5toi;
92 int arg4toi = atoi(args[3]);
100 int number_target_coords = 200;
102 int arg3toi = atoi(args[2]);
104 number_target_coords = arg3toi;
112 int arg2toi = atoi(args[1]);
124 Kokkos::Profiling::pushRegion(
"Setup Point Data");
129 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> source_coords_device(
"source coordinates",
130 1.25*N_pts_on_sphere, 3);
131 Kokkos::View<double**>::HostMirror source_coords = Kokkos::create_mirror_view(source_coords_device);
136 int number_source_coords;
141 double a = 4*
PI*r*r/N_pts_on_sphere;
142 double d = std::sqrt(a);
143 int M_theta = std::round(
PI/d);
144 double d_theta =
PI/M_theta;
145 double d_phi = a/d_theta;
146 for (
int i=0; i<M_theta; ++i) {
147 double theta =
PI*(i + 0.5)/M_theta;
148 int M_phi = std::round(2*
PI*std::sin(theta)/d_phi);
149 for (
int j=0; j<M_phi; ++j) {
150 double phi = 2*
PI*j/M_phi;
151 source_coords(N_count, 0) = theta;
152 source_coords(N_count, 1) = phi;
157 for (
int i=0; i<N_count; ++i) {
158 double theta = source_coords(i,0);
159 double phi = source_coords(i,1);
160 source_coords(i,0) = r*std::sin(theta)*std::cos(phi);
161 source_coords(i,1) = r*std::sin(theta)*std::sin(phi);
162 source_coords(i,2) = r*cos(theta);
165 number_source_coords = N_count;
169 Kokkos::resize(source_coords, number_source_coords, 3);
170 Kokkos::resize(source_coords_device, number_source_coords, 3);
173 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> target_coords_device(
"target coordinates",
174 number_target_coords, 3);
175 Kokkos::View<double**>::HostMirror target_coords = Kokkos::create_mirror_view(target_coords_device);
178 std::mt19937 rng(50);
181 std::uniform_int_distribution<int> gen_num_neighbors(0, number_source_coords-1);
184 for (
int i=0; i<number_target_coords; ++i) {
185 const int source_site_to_copy = gen_num_neighbors(rng);
186 for (
int j=0; j<3; ++j) {
187 target_coords(i,j) = source_coords(source_site_to_copy,j);
194 Kokkos::Profiling::popRegion();
196 Kokkos::Profiling::pushRegion(
"Creating Data");
202 Kokkos::deep_copy(source_coords_device, source_coords);
203 Kokkos::deep_copy(target_coords_device, target_coords);
210 Kokkos::View<double*, Kokkos::DefaultExecutionSpace> sampling_data_device(
"samples of true solution",
211 source_coords_device.extent(0));
214 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> sampling_vector_data_device(
"samples of vector true solution",
215 source_coords_device.extent(0), 3);
217 Kokkos::parallel_for(
"Sampling Manufactured Solutions", Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>
218 (0,source_coords.extent(0)), KOKKOS_LAMBDA(
const int i) {
221 double xval = source_coords_device(i,0);
222 double yval = (dimension>1) ? source_coords_device(i,1) : 0;
223 double zval = (dimension>2) ? source_coords_device(i,2) : 0;
229 for (
int j=0; j<3; ++j) {
230 double gradient[3] = {0,0,0};
232 sampling_vector_data_device(i,j) = gradient[j];
240 Kokkos::Profiling::popRegion();
241 Kokkos::Profiling::pushRegion(
"Neighbor Search");
252 double epsilon_multiplier = 1.5;
253 int estimated_upper_bound_number_neighbors =
254 point_cloud_search.getEstimatedNumberNeighborsUpperBound(min_neighbors, dimension, epsilon_multiplier);
256 Kokkos::View<int**, Kokkos::DefaultExecutionSpace> neighbor_lists_device(
"neighbor lists",
257 number_target_coords, estimated_upper_bound_number_neighbors);
258 Kokkos::View<int**>::HostMirror neighbor_lists = Kokkos::create_mirror_view(neighbor_lists_device);
261 Kokkos::View<double*, Kokkos::DefaultExecutionSpace> epsilon_device(
"h supports", number_target_coords);
262 Kokkos::View<double*>::HostMirror epsilon = Kokkos::create_mirror_view(epsilon_device);
267 point_cloud_search.generate2DNeighborListsFromKNNSearch(
false , target_coords, neighbor_lists,
268 epsilon, min_neighbors, epsilon_multiplier);
273 Kokkos::Profiling::popRegion();
284 Kokkos::deep_copy(neighbor_lists_device, neighbor_lists);
285 Kokkos::deep_copy(epsilon_device, epsilon);
288 std::string solver_name;
289 if (solver_type == 0) {
291 }
else if (solver_type == 1) {
293 }
else if (solver_type == 2) {
298 std::string problem_name;
299 if (problem_type == 0) {
300 problem_name =
"STANDARD";
301 }
else if (problem_type == 1) {
302 problem_name =
"MANIFOLD";
306 std::string constraint_name;
307 if (constraint_type == 0) {
308 constraint_name =
"NO_CONSTRAINT";
309 }
else if (constraint_type == 1) {
310 constraint_name =
"NEUMANN_GRAD_SCALAR";
317 solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
334 my_GMLS_vector_1.
setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
337 std::vector<TargetOperation> lro_vector_1(1);
341 my_GMLS_vector_1.addTargets(lro_vector_1);
347 my_GMLS_vector_1.setCurvatureWeightingPower(2);
353 my_GMLS_vector_1.setWeightingPower(2);
356 my_GMLS_vector_1.setOrderOfQuadraturePoints(2);
357 my_GMLS_vector_1.setDimensionOfQuadraturePoints(1);
358 my_GMLS_vector_1.setQuadratureType(
"LINE");
361 my_GMLS_vector_1.generateAlphas();
368 solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
371 my_GMLS_vector_2.
setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
372 std::vector<TargetOperation> lro_vector_2(2);
376 my_GMLS_vector_2.addTargets(lro_vector_2);
378 my_GMLS_vector_2.setCurvatureWeightingPower(2);
380 my_GMLS_vector_2.setWeightingPower(2);
381 my_GMLS_vector_2.setOrderOfQuadraturePoints(2);
382 my_GMLS_vector_2.setDimensionOfQuadraturePoints(1);
383 my_GMLS_vector_2.setQuadratureType(
"LINE");
384 my_GMLS_vector_2.generateAlphas();
390 solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
393 my_GMLS_scalar.
setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
395 std::vector<TargetOperation> lro_scalar(1);
398 my_GMLS_scalar.addTargets(lro_scalar);
400 my_GMLS_scalar.setCurvatureWeightingPower(2);
402 my_GMLS_scalar.setWeightingPower(2);
403 my_GMLS_scalar.generateAlphas();
408 double instantiation_time = timer.seconds();
409 std::cout <<
"Took " << instantiation_time <<
"s to complete alphas generation." << std::endl;
411 Kokkos::Profiling::pushRegion(
"Apply Alphas to Data");
426 Evaluator vector_1_gmls_evaluator(&my_GMLS_vector_1);
427 Evaluator vector_2_gmls_evaluator(&my_GMLS_vector_2);
428 Evaluator scalar_gmls_evaluator(&my_GMLS_scalar);
439 auto output_divergence_vectorsamples =
443 auto output_divergence_scalarsamples =
447 auto output_laplacian_vectorbasis =
451 auto output_laplacian_scalarbasis =
459 Kokkos::Profiling::popRegion();
461 Kokkos::Profiling::pushRegion(
"Comparison");
465 double laplacian_vectorbasis_error = 0;
466 double laplacian_vectorbasis_norm = 0;
468 double laplacian_scalarbasis_error = 0;
469 double laplacian_scalarbasis_norm = 0;
471 double gradient_vectorbasis_ambient_error = 0;
472 double gradient_vectorbasis_ambient_norm = 0;
474 double gradient_scalarbasis_ambient_error = 0;
475 double gradient_scalarbasis_ambient_norm = 0;
477 double divergence_vectorsamples_ambient_error = 0;
478 double divergence_vectorsamples_ambient_norm = 0;
480 double divergence_scalarsamples_ambient_error = 0;
481 double divergence_scalarsamples_ambient_norm = 0;
484 for (
int i=0; i<number_target_coords; i++) {
487 double xval = target_coords(i,0);
488 double yval = (dimension>1) ? target_coords(i,1) : 0;
489 double zval = (dimension>2) ? target_coords(i,2) : 0;
493 double actual_Gradient_ambient[3] = {0,0,0};
496 laplacian_vectorbasis_error += (output_laplacian_vectorbasis(i) - actual_Laplacian)*(output_laplacian_vectorbasis(i) - actual_Laplacian);
497 laplacian_vectorbasis_norm += actual_Laplacian*actual_Laplacian;
500 laplacian_scalarbasis_error += (output_laplacian_scalarbasis(i) - actual_Laplacian)*(output_laplacian_scalarbasis(i) - actual_Laplacian);
501 laplacian_scalarbasis_norm += actual_Laplacian*actual_Laplacian;
516 divergence_vectorsamples_ambient_error += (output_divergence_vectorsamples(i) - actual_Laplacian)*(output_divergence_vectorsamples(i) - actual_Laplacian);
517 divergence_vectorsamples_ambient_norm += actual_Laplacian*actual_Laplacian;
519 divergence_scalarsamples_ambient_error += (output_divergence_scalarsamples(i) - actual_Laplacian)*(output_divergence_scalarsamples(i) - actual_Laplacian);
520 divergence_scalarsamples_ambient_norm += actual_Laplacian*actual_Laplacian;
524 laplacian_vectorbasis_error /= number_target_coords;
525 laplacian_vectorbasis_error = std::sqrt(laplacian_vectorbasis_error);
526 laplacian_vectorbasis_norm /= number_target_coords;
527 laplacian_vectorbasis_norm = std::sqrt(laplacian_vectorbasis_norm);
529 laplacian_scalarbasis_error /= number_target_coords;
530 laplacian_scalarbasis_error = std::sqrt(laplacian_scalarbasis_error);
531 laplacian_scalarbasis_norm /= number_target_coords;
532 laplacian_scalarbasis_norm = std::sqrt(laplacian_scalarbasis_norm);
534 gradient_vectorbasis_ambient_error /= number_target_coords;
535 gradient_vectorbasis_ambient_error = std::sqrt(gradient_vectorbasis_ambient_error);
536 gradient_vectorbasis_ambient_norm /= number_target_coords;
537 gradient_vectorbasis_ambient_norm = std::sqrt(gradient_vectorbasis_ambient_norm);
539 gradient_scalarbasis_ambient_error /= number_target_coords;
540 gradient_scalarbasis_ambient_error = std::sqrt(gradient_scalarbasis_ambient_error);
541 gradient_scalarbasis_ambient_norm /= number_target_coords;
542 gradient_scalarbasis_ambient_norm = std::sqrt(gradient_scalarbasis_ambient_norm);
544 divergence_vectorsamples_ambient_error /= number_target_coords;
545 divergence_vectorsamples_ambient_error = std::sqrt(divergence_vectorsamples_ambient_error);
546 divergence_vectorsamples_ambient_norm /= number_target_coords;
547 divergence_vectorsamples_ambient_norm = std::sqrt(divergence_vectorsamples_ambient_norm);
549 divergence_scalarsamples_ambient_error /= number_target_coords;
550 divergence_scalarsamples_ambient_error = std::sqrt(divergence_scalarsamples_ambient_error);
551 divergence_scalarsamples_ambient_norm /= number_target_coords;
552 divergence_scalarsamples_ambient_norm = std::sqrt(divergence_scalarsamples_ambient_norm);
554 printf(
"Staggered Laplace-Beltrami (VectorBasis) Error: %g\n", laplacian_vectorbasis_error / laplacian_vectorbasis_norm);
555 printf(
"Staggered Laplace-Beltrami (ScalarBasis) Error: %g\n", laplacian_scalarbasis_error / laplacian_scalarbasis_norm);
556 printf(
"Surface Staggered Gradient (VectorBasis) Error: %g\n", gradient_vectorbasis_ambient_error / gradient_vectorbasis_ambient_norm);
557 printf(
"Surface Staggered Gradient (ScalarBasis) Error: %g\n", gradient_scalarbasis_ambient_error / gradient_scalarbasis_ambient_norm);
558 printf(
"Surface Staggered Divergence (VectorSamples) Error: %g\n", divergence_vectorsamples_ambient_error / divergence_vectorsamples_ambient_norm);
559 printf(
"Surface Staggered Divergence (ScalarSamples) Error: %g\n", divergence_scalarsamples_ambient_error / divergence_scalarsamples_ambient_norm);
563 Kokkos::Profiling::popRegion();
572 #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.
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)
Vector polynomial basis having # of components _dimensions, or (_dimensions-1) in the case of manifol...