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);
40 bool all_passed =
true;
47 auto order = clp.
order;
56 const double failure_tolerance = 1e-9;
63 Kokkos::Profiling::pushRegion(
"Setup Point Data");
67 double h_spacing = 0.05;
68 int n_neg1_to_1 = 2*(1/h_spacing) + 1;
71 const int number_source_coords = std::pow(n_neg1_to_1, dimension);
74 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> source_coords_device(
"source coordinates",
75 number_source_coords, 3);
76 Kokkos::View<double**>::HostMirror source_coords = Kokkos::create_mirror_view(source_coords_device);
79 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> target_coords_device (
"target coordinates", number_target_coords, 3);
80 Kokkos::View<double**>::HostMirror target_coords = Kokkos::create_mirror_view(target_coords_device);
85 double this_coord[3] = {0,0,0};
86 for (
int i=-n_neg1_to_1/2; i<n_neg1_to_1/2+1; ++i) {
87 this_coord[0] = i*h_spacing;
88 for (
int j=-n_neg1_to_1/2; j<n_neg1_to_1/2+1; ++j) {
89 this_coord[1] = j*h_spacing;
90 for (
int k=-n_neg1_to_1/2; k<n_neg1_to_1/2+1; ++k) {
91 this_coord[2] = k*h_spacing;
93 source_coords(source_index,0) = this_coord[0];
94 source_coords(source_index,1) = this_coord[1];
95 source_coords(source_index,2) = this_coord[2];
100 source_coords(source_index,0) = this_coord[0];
101 source_coords(source_index,1) = this_coord[1];
102 source_coords(source_index,2) = 0;
107 source_coords(source_index,0) = this_coord[0];
108 source_coords(source_index,1) = 0;
109 source_coords(source_index,2) = 0;
118 std::mt19937 rng(50);
120 std::uniform_int_distribution<int> gen_num_neighbours(0, source_index);
122 for (
int i=0; i<number_target_coords; ++i) {
123 const int source_site_to_copy = gen_num_neighbours(rng);
124 for (
int j=0; j<3; j++) {
125 target_coords(i, j) = source_coords(source_site_to_copy, j);
131 Kokkos::Profiling::popRegion();
132 Kokkos::Profiling::pushRegion(
"Creating Data");
138 Kokkos::deep_copy(source_coords_device, source_coords);
141 Kokkos::deep_copy(target_coords_device, target_coords);
144 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> vector_sampling_span_basis_data_device(
"samples of true vector solution",
145 source_coords_device.extent(0), dimension);
146 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> vector_sampling_single_polynomial_data_device(
"samples of true vector solution",
147 source_coords_device.extent(0), dimension);
149 Kokkos::parallel_for(
"Sampling Manufactured Solutions", Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>
150 (0,source_coords.extent(0)), KOKKOS_LAMBDA(
const int i) {
152 double xval = source_coords_device(i,0);
153 double yval = (dimension>1) ? source_coords_device(i,1) : 0;
154 double zval = (dimension>2) ? source_coords_device(i,2) : 0;
157 for (
int j=0; j<dimension; ++j) {
166 Kokkos::Profiling::popRegion();
167 Kokkos::Profiling::pushRegion(
"Neighbor Search");
179 double epsilon_multiplier = 2.2;
180 int estimated_upper_bound_number_neighbors =
181 point_cloud_search.getEstimatedNumberNeighborsUpperBound(min_neighbors, dimension, epsilon_multiplier);
183 Kokkos::View<int**, Kokkos::DefaultExecutionSpace> neighbor_lists_device(
"neighbor lists",
184 number_target_coords, estimated_upper_bound_number_neighbors);
185 Kokkos::View<int**>::HostMirror neighbor_lists = Kokkos::create_mirror_view(neighbor_lists_device);
188 Kokkos::View<double*, Kokkos::DefaultExecutionSpace> epsilon_device(
"h supports", number_target_coords);
189 Kokkos::View<double*>::HostMirror epsilon = Kokkos::create_mirror_view(epsilon_device);
194 point_cloud_search.generate2DNeighborListsFromKNNSearch(
false , target_coords, neighbor_lists,
195 epsilon, min_neighbors, epsilon_multiplier);
199 Kokkos::Profiling::popRegion();
210 Kokkos::deep_copy(neighbor_lists_device, neighbor_lists);
211 Kokkos::deep_copy(epsilon_device, epsilon);
219 solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
236 vector_divfree_basis_gmls.
setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
239 std::vector<TargetOperation> lro(4);
246 vector_divfree_basis_gmls.addTargets(lro);
252 vector_divfree_basis_gmls.setWeightingParameter(2);
255 vector_divfree_basis_gmls.generateAlphas(15 );
259 double instantiation_time = timer.seconds();
260 std::cout <<
"Took " << instantiation_time <<
"s to complete alphas generation." << std::endl;
262 Kokkos::Profiling::pushRegion(
"Apply Alphas to Data");
276 Evaluator gmls_evaluator_vector(&vector_divfree_basis_gmls);
290 Kokkos::Profiling::popRegion();
292 Kokkos::Profiling::pushRegion(
"Comparison");
297 for (
int i=0; i<number_target_coords; i++) {
299 double GMLS_DivFree_VectorX = output_vector_evaluation(i, 0);
300 double GMLS_DivFree_VectorY = (dimension>1) ? output_vector_evaluation(i, 1) : 0;
301 double GMLS_DivFree_VectorZ = (dimension>2) ? output_vector_evaluation(i, 2) : 0;
304 double GMLS_Curl_DivFree_VectorX = output_curl_vector_evaluation(i, 0);
305 double GMLS_Curl_DivFree_VectorY = (dimension>2) ? output_curl_vector_evaluation(i, 1) : 0;
306 double GMLS_Curl_DivFree_VectorZ = (dimension>2) ? output_curl_vector_evaluation(i, 2) : 0;
309 double GMLS_CurlCurl_DivFree_VectorX = output_curlcurl_vector_evaluation(i, 0);
310 double GMLS_CurlCurl_DivFree_VectorY = (dimension>1) ? output_curlcurl_vector_evaluation(i, 1) : 0;
311 double GMLS_CurlCurl_DivFree_VectorZ = (dimension>2) ? output_curlcurl_vector_evaluation(i, 2) : 0;
314 double GMLS_Gradient_DivFree_VectorXX = output_gradient_vector_evaluation(i, 0);
315 double GMLS_Gradient_DivFree_VectorXY = output_gradient_vector_evaluation(i, 1);
316 double GMLS_Gradient_DivFree_VectorXZ = (dimension==3) ? output_gradient_vector_evaluation(i, 2) : 0.0;
317 double GMLS_Gradient_DivFree_VectorYX = (dimension==3) ? output_gradient_vector_evaluation(i, 3) : output_gradient_vector_evaluation(i, 2);
318 double GMLS_Gradient_DivFree_VectorYY = (dimension==3) ? output_gradient_vector_evaluation(i, 4) : output_gradient_vector_evaluation(i, 3);
319 double GMLS_Gradient_DivFree_VectorYZ = (dimension==3) ? output_gradient_vector_evaluation(i, 5) : 0.0;
320 double GMLS_Gradient_DivFree_VectorZX = (dimension==3) ? output_gradient_vector_evaluation(i, 6) : 0.0;
321 double GMLS_Gradient_DivFree_VectorZY = (dimension==3) ? output_gradient_vector_evaluation(i, 7) : 0.0;
322 double GMLS_Gradient_DivFree_VectorZZ = (dimension==3) ? output_gradient_vector_evaluation(i, 8) : 0.0;
325 double xval = target_coords(i,0);
326 double yval = (dimension>1) ? target_coords(i,1) : 0;
327 double zval = (dimension>2) ? target_coords(i,2) : 0;
330 double actual_vector[3] = {0, 0, 0};
342 double actual_curl_vector[3] = {0, 0, 0};
352 double actual_curlcurl_vector[3] = {0, 0, 0};
363 double actual_gradient_vector[9] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
365 for (
int axes = 0; axes < 9; ++axes)
369 for (
int axes = 0; axes < 4; ++axes)
374 if(std::abs(actual_vector[0] - GMLS_DivFree_VectorX) > failure_tolerance) {
376 std::cout << i <<
" Failed VectorX by: " << std::abs(actual_vector[0] - GMLS_DivFree_VectorX) << std::endl;
378 if(std::abs(actual_vector[1] - GMLS_DivFree_VectorY) > failure_tolerance) {
380 std::cout << i <<
" Failed VectorY by: " << std::abs(actual_vector[1] - GMLS_DivFree_VectorY) << std::endl;
384 if(std::abs(actual_vector[2] - GMLS_DivFree_VectorZ) > failure_tolerance) {
386 std::cout << i <<
" Failed VectorZ by: " << std::abs(actual_vector[2] - GMLS_DivFree_VectorZ) << std::endl;
393 if(std::abs(actual_curl_vector[0] - GMLS_Curl_DivFree_VectorX) > failure_tolerance) {
395 std::cout << i <<
" Failed curl VectorX by: " << std::abs(actual_curl_vector[2] - GMLS_Curl_DivFree_VectorX) << std::endl;
397 }
else if (dimension>2) {
398 if(std::abs(actual_curl_vector[0] - GMLS_Curl_DivFree_VectorX) > failure_tolerance) {
400 std::cout << i <<
" Failed curl VectorX by: " << std::abs(actual_curl_vector[0] - GMLS_Curl_DivFree_VectorX) << std::endl;
402 if(std::abs(actual_curl_vector[1] - GMLS_Curl_DivFree_VectorY) > failure_tolerance) {
404 std::cout << i <<
" Failed curl VectorY by: " << std::abs(actual_curl_vector[1] - GMLS_Curl_DivFree_VectorY) << std::endl;
406 if(std::abs(actual_curl_vector[2] - GMLS_Curl_DivFree_VectorZ) > failure_tolerance) {
408 std::cout << i <<
" Failed curl VectorZ by: " << std::abs(actual_curl_vector[2] - GMLS_Curl_DivFree_VectorZ) << std::endl;
413 if(std::abs(actual_curlcurl_vector[0] - GMLS_CurlCurl_DivFree_VectorX) > failure_tolerance) {
415 std::cout << i <<
" Failed curl curl VectorX by: " << std::abs(actual_curlcurl_vector[0] - GMLS_CurlCurl_DivFree_VectorX) << std::endl;
418 if(std::abs(actual_curlcurl_vector[1] - GMLS_CurlCurl_DivFree_VectorY) > failure_tolerance) {
420 std::cout << i <<
" Failed curl curl VectorY by: " << std::abs(actual_curlcurl_vector[1] - GMLS_CurlCurl_DivFree_VectorY) << std::endl;
424 if(std::abs(actual_curlcurl_vector[2] - GMLS_CurlCurl_DivFree_VectorZ) > failure_tolerance) {
426 std::cout << i <<
" Failed curl curl VectorZ by: " << std::abs(actual_curlcurl_vector[2] - GMLS_CurlCurl_DivFree_VectorZ) << std::endl;
432 if (std::abs(actual_gradient_vector[0] - GMLS_Gradient_DivFree_VectorXX) > failure_tolerance) {
434 std::cout << i <<
" Failed gradient_x VectorX by: " << std::abs(actual_gradient_vector[0] - GMLS_Gradient_DivFree_VectorXX) << std::endl;
436 if (std::abs(actual_gradient_vector[1] - GMLS_Gradient_DivFree_VectorXY) > failure_tolerance) {
438 std::cout << i <<
" Failed gradient_y VectorX by: " << std::abs(actual_gradient_vector[1] - GMLS_Gradient_DivFree_VectorXY) << std::endl;
440 if (std::abs(actual_gradient_vector[2] - GMLS_Gradient_DivFree_VectorXZ) > failure_tolerance) {
442 std::cout << i <<
" Failed gradient_z VectorX by: " << std::abs(actual_gradient_vector[2] - GMLS_Gradient_DivFree_VectorXZ) << std::endl;
445 if (std::abs(actual_gradient_vector[3] - GMLS_Gradient_DivFree_VectorYX) > failure_tolerance) {
447 std::cout << i <<
" Failed gradient_x VectorY by: " << std::abs(actual_gradient_vector[3] - GMLS_Gradient_DivFree_VectorYX) << std::endl;
449 if (std::abs(actual_gradient_vector[4] - GMLS_Gradient_DivFree_VectorYY) > failure_tolerance) {
451 std::cout << i <<
" Failed gradient_y VectorY by: " << std::abs(actual_gradient_vector[4] - GMLS_Gradient_DivFree_VectorYY) << std::endl;
453 if (std::abs(actual_gradient_vector[5] - GMLS_Gradient_DivFree_VectorYZ) > failure_tolerance) {
455 std::cout << i <<
" Failed gradient_z VectorY by: " << std::abs(actual_gradient_vector[5] - GMLS_Gradient_DivFree_VectorYZ) << std::endl;
458 if (std::abs(actual_gradient_vector[6] - GMLS_Gradient_DivFree_VectorZX) > failure_tolerance) {
460 std::cout << i <<
" Failed gradient_x VectorZ by: " << std::abs(actual_gradient_vector[6] - GMLS_Gradient_DivFree_VectorZX) << std::endl;
462 if (std::abs(actual_gradient_vector[7] - GMLS_Gradient_DivFree_VectorZY) > failure_tolerance) {
464 std::cout << i <<
" Failed gradient_y VectorZ by: " << std::abs(actual_gradient_vector[7] - GMLS_Gradient_DivFree_VectorZY) << std::endl;
466 if (std::abs(actual_gradient_vector[8] - GMLS_Gradient_DivFree_VectorZZ) > failure_tolerance) {
468 std::cout << i <<
" Failed gradient_z VectorZ by: " << std::abs(actual_gradient_vector[8] - GMLS_Gradient_DivFree_VectorZZ) << std::endl;
473 if (std::abs(actual_gradient_vector[0] - GMLS_Gradient_DivFree_VectorXX) > failure_tolerance) {
475 std::cout << i <<
" Failed gradient_x VectorX by: " << std::abs(actual_gradient_vector[0] - GMLS_Gradient_DivFree_VectorXX) << std::endl;
477 if (std::abs(actual_gradient_vector[1] - GMLS_Gradient_DivFree_VectorXY) > failure_tolerance) {
479 std::cout << i <<
" Failed gradient_y VectorX by: " << std::abs(actual_gradient_vector[1] - GMLS_Gradient_DivFree_VectorXY) << std::endl;
482 if (std::abs(actual_gradient_vector[2] - GMLS_Gradient_DivFree_VectorYX) > failure_tolerance) {
484 std::cout << i <<
" Failed gradient_x VectorY by: " << std::abs(actual_gradient_vector[2] - GMLS_Gradient_DivFree_VectorYX) << std::endl;
486 if (std::abs(actual_gradient_vector[3] - GMLS_Gradient_DivFree_VectorYY) > failure_tolerance) {
488 std::cout << i <<
" Failed gradient_y VectorY by: " << (std::abs(actual_gradient_vector[3] - GMLS_Gradient_DivFree_VectorYY) > failure_tolerance) << std::endl;
496 Kokkos::Profiling::popRegion();
505 #ifdef COMPADRE_USE_MPI
511 fprintf(stdout,
"Passed test \n");
514 fprintf(stdout,
"Failed test \n");
KOKKOS_INLINE_FUNCTION double divfreeTestSolution_span_basis(double x, double y, double z, int component, int dimension, int exact_order)
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...
KOKKOS_INLINE_FUNCTION double curldivfreeTestSolution_span_basis(double x, double y, double z, int component, int dimension, int exact_order)
int main(int argc, char **argv)
Point evaluation of the curl of a vector (results in a vector)
KOKKOS_INLINE_FUNCTION double curlcurldivfreeTestSolution_single_polynomial(double x, double y, double z, int component, int dimension)
Point evaluation of the curl of a curl of a vector (results in a vector)
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...
KOKKOS_INLINE_FUNCTION double divfreeTestSolution_single_polynomial(double x, double y, double z, int component, int dimension)
Generalized Moving Least Squares (GMLS)
Divergence-free vector polynomial basis.
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)
KOKKOS_INLINE_FUNCTION double gradientdivfreeTestSolution_span_basis(double x, double y, double z, int input_component, int output_component, int dimension, int exact_order)
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) ...
std::string constraint_name
constexpr SamplingFunctional VectorPointSample
Point evaluations of the entire vector source function.
Point evaluation of the gradient of a vector (results in a matrix, NOT CURRENTLY IMPLEMENTED) ...
Point evaluation of a vector (reconstructs entire vector at once, requiring a ReconstructionSpace hav...