16 #include <Compadre_Config.h>
23 #ifdef COMPADRE_USE_MPI
27 #include <Kokkos_Timer.hpp>
28 #include <Kokkos_Core.hpp>
30 using namespace Compadre;
35 int main (
int argc,
char* args[]) {
38 #ifdef COMPADRE_USE_MPI
39 MPI_Init(&argc, &args);
43 Kokkos::initialize(argc, args);
46 bool all_passed =
true;
55 int constraint_type = 0;
57 int arg7toi = atoi(args[6]);
59 constraint_type = arg7toi;
69 int arg6toi = atoi(args[5]);
71 problem_type = arg6toi;
82 int arg5toi = atoi(args[4]);
84 solver_type = 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]);
120 const double failure_tolerance = 1e-9;
123 const double laplacian_failure_tolerance = 1e-9;
133 double h_spacing = 0.05;
134 int n_neg1_to_1 = 2*(1/h_spacing) + 1;
137 const int number_source_coords = std::pow(n_neg1_to_1, dimension);
140 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> source_coords_device(
"source coordinates",
141 number_source_coords, 3);
142 Kokkos::View<double**>::HostMirror source_coords = Kokkos::create_mirror_view(source_coords_device);
145 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> target_coords_device (
"target coordinates", number_target_coords, 3);
146 Kokkos::View<double**>::HostMirror target_coords = Kokkos::create_mirror_view(target_coords_device);
150 int source_index = 0;
151 double this_coord[3] = {0,0,0};
152 for (
int i=-n_neg1_to_1/2; i<n_neg1_to_1/2+1; ++i) {
153 this_coord[0] = i*h_spacing;
154 for (
int j=-n_neg1_to_1/2; j<n_neg1_to_1/2+1; ++j) {
155 this_coord[1] = j*h_spacing;
156 for (
int k=-n_neg1_to_1/2; k<n_neg1_to_1/2+1; ++k) {
157 this_coord[2] = k*h_spacing;
159 source_coords(source_index,0) = this_coord[0];
160 source_coords(source_index,1) = this_coord[1];
161 source_coords(source_index,2) = this_coord[2];
166 source_coords(source_index,0) = this_coord[0];
167 source_coords(source_index,1) = this_coord[1];
168 source_coords(source_index,2) = 0;
173 source_coords(source_index,0) = this_coord[0];
174 source_coords(source_index,1) = 0;
175 source_coords(source_index,2) = 0;
181 for(
int i=0; i<number_target_coords; i++){
184 double rand_dir[3] = {0,0,0};
186 for (
int j=0; j<dimension; ++j) {
188 rand_dir[j] = ((double)rand() / (double) RAND_MAX) - 0.5;
192 for (
int j=0; j<dimension; ++j) {
193 target_coords(i,j) = rand_dir[j];
206 Kokkos::deep_copy(source_coords_device, source_coords);
209 Kokkos::deep_copy(target_coords_device, target_coords);
212 Kokkos::View<double*, Kokkos::DefaultExecutionSpace> sampling_data_device(
"samples of true solution",
213 source_coords_device.extent(0));
215 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> gradient_sampling_data_device(
"samples of true gradient",
216 source_coords_device.extent(0), dimension);
218 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> divergence_sampling_data_device
219 (
"samples of true solution for divergence test", source_coords_device.extent(0), dimension);
221 Kokkos::parallel_for(
"Sampling Manufactured Solutions", Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>
222 (0,source_coords.extent(0)), KOKKOS_LAMBDA(
const int i) {
225 double xval = source_coords_device(i,0);
226 double yval = (dimension>1) ? source_coords_device(i,1) : 0;
227 double zval = (dimension>2) ? source_coords_device(i,2) : 0;
230 sampling_data_device(i) =
trueSolution(xval, yval, zval, order, dimension);
233 double true_grad[3] = {0,0,0};
234 trueGradient(true_grad, xval, yval,zval, order, dimension);
236 for (
int j=0; j<dimension; ++j) {
237 gradient_sampling_data_device(i,j) = true_grad[j];
252 std::string solver_name;
253 if (solver_type == 0) {
255 }
else if (solver_type == 1) {
257 }
else if (solver_type == 2) {
262 std::string problem_name;
263 if (problem_type == 0) {
264 problem_name =
"STANDARD";
265 }
else if (problem_type == 1) {
266 problem_name =
"MANIFOLD";
270 std::string constraint_name;
271 if (constraint_type == 0) {
272 constraint_name =
"NO_CONSTRAINT";
273 }
else if (constraint_type == 1) {
274 constraint_name =
"NEUMANN_GRAD_SCALAR";
280 solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
284 std::vector<TargetOperation> lro(5);
292 my_GMLS.addTargets(lro);
298 my_GMLS.setWeightingPower(2);
301 my_GMLS.setSourceSites(source_coords_device);
309 for (
int i=0; i<number_target_coords; i++) {
328 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> single_target_coords_device (
"single target coordinates", 1, 3);
329 Kokkos::View<double**>::HostMirror single_target_coords = Kokkos::create_mirror_view(single_target_coords_device);
330 for (
int j=0; j<3; ++j) {
331 single_target_coords(0,j) = target_coords(i,j);
334 Kokkos::deep_copy(single_target_coords_device, single_target_coords);
341 double epsilon_multiplier = 1.5;
342 int estimated_upper_bound_number_neighbors =
343 point_cloud_search.getEstimatedNumberNeighborsUpperBound(min_neighbors, dimension, epsilon_multiplier);
345 Kokkos::View<int**, Kokkos::DefaultExecutionSpace> single_neighbor_lists_device(
"neighbor lists",
346 1, estimated_upper_bound_number_neighbors);
347 Kokkos::View<int**>::HostMirror single_neighbor_lists = Kokkos::create_mirror_view(single_neighbor_lists_device);
350 Kokkos::View<double*, Kokkos::DefaultExecutionSpace> single_epsilon_device(
"h supports", 1);
351 Kokkos::View<double*>::HostMirror single_epsilon = Kokkos::create_mirror_view(single_epsilon_device);
356 point_cloud_search.generate2DNeighborListsFromKNNSearch(
false , single_target_coords,
357 single_neighbor_lists, single_epsilon, min_neighbors, epsilon_multiplier);
366 Kokkos::deep_copy(single_neighbor_lists_device, single_neighbor_lists);
367 Kokkos::deep_copy(single_epsilon_device, single_epsilon);
371 my_GMLS.setNeighborLists(single_neighbor_lists_device);
372 my_GMLS.setTargetSites(single_target_coords_device);
373 my_GMLS.setWindowSizes(single_epsilon_device);
377 my_GMLS.generateAlphas(1,
true );
382 double instantiation_time = timer.seconds();
383 std::cout <<
"Took " << instantiation_time <<
"s to complete alphas generation." << std::endl;
418 (sampling_data_device);
428 double GMLS_value = output_value(0);
431 double GMLS_Laplacian = output_laplacian(0);
437 double GMLS_GradX = scalar_coefficients(0,1)*1./single_epsilon(0);
441 double GMLS_GradY = (dimension>1) ? output_gradient(0,1) : 0;
444 double GMLS_GradZ = (dimension>2) ? output_gradient(0,2) : 0;
447 double GMLS_Divergence = output_divergence(0);
450 double GMLS_CurlX = (dimension>1) ? output_curl(0,0) : 0;
451 double GMLS_CurlY = (dimension>1) ? output_curl(0,1) : 0;
452 double GMLS_CurlZ = (dimension>2) ? output_curl(0,2) : 0;
456 double xval = target_coords(i,0);
457 double yval = (dimension>1) ? target_coords(i,1) : 0;
458 double zval = (dimension>2) ? target_coords(i,2) : 0;
461 double actual_value =
trueSolution(xval, yval, zval, order, dimension);
462 double actual_Laplacian =
trueLaplacian(xval, yval, zval, order, dimension);
464 double actual_Gradient[3] = {0,0,0};
465 trueGradient(actual_Gradient, xval, yval, zval, order, dimension);
467 double actual_Divergence;
468 actual_Divergence =
trueLaplacian(xval, yval, zval, order, dimension);
470 double actual_Curl[3] = {0,0,0};
481 if(GMLS_value!=GMLS_value || std::abs(actual_value - GMLS_value) > failure_tolerance) {
483 std::cout << i <<
" Failed Actual by: " << std::abs(actual_value - GMLS_value) << std::endl;
487 if(std::abs(actual_Laplacian - GMLS_Laplacian) > laplacian_failure_tolerance) {
489 std::cout << i <<
" Failed Laplacian by: " << std::abs(actual_Laplacian - GMLS_Laplacian) << std::endl;
493 if(std::abs(actual_Gradient[0] - GMLS_GradX) > failure_tolerance) {
495 std::cout << i <<
" Failed GradX by: " << std::abs(actual_Gradient[0] - GMLS_GradX) << std::endl;
497 if(std::abs(actual_Gradient[1] - GMLS_GradY) > failure_tolerance) {
499 std::cout << i <<
" Failed GradY by: " << std::abs(actual_Gradient[1] - GMLS_GradY) << std::endl;
503 if(std::abs(actual_Gradient[2] - GMLS_GradZ) > failure_tolerance) {
505 std::cout << i <<
" Failed GradZ by: " << std::abs(actual_Gradient[2] - GMLS_GradZ) << std::endl;
511 if(std::abs(actual_Divergence - GMLS_Divergence) > failure_tolerance) {
513 std::cout << i <<
" Failed Divergence by: " << std::abs(actual_Divergence - GMLS_Divergence) << std::endl;
520 tmp_diff += std::abs(actual_Curl[0] - GMLS_CurlX) + std::abs(actual_Curl[1] - GMLS_CurlY);
522 tmp_diff += std::abs(actual_Curl[2] - GMLS_CurlZ);
523 if(std::abs(tmp_diff) > failure_tolerance) {
525 std::cout << i <<
" Failed Curl by: " << std::abs(tmp_diff) << std::endl;
541 #ifdef COMPADRE_USE_MPI
547 fprintf(stdout,
"Passed test \n");
550 fprintf(stdout,
"Failed test \n");
Point evaluation of the curl of a vector (results in a vector)
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::View< output_data_type, output_array_layout, output_memory_space > applyFullPolynomialCoefficientsBasisToDataAllComponents(view_type_input_data sampling_data, bool scalar_as_vector_if_needed=true) const
Generation of polynomial reconstruction coefficients by applying to data in GMLS. ...
int main(int argc, char *args[])
[Parse Command Line Arguments]
Point evaluation of a scalar.
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 void trueGradient(double *ans, double x, double y, double z, int order, int dimension)
Scalar basis reused as many times as there are components in the vector resulting in a much cheaper p...
Point evaluation of the laplacian of a scalar (could be on a manifold or not)
Point evaluation of the divergence of a vector (results in a scalar)
Point evaluation of the gradient of a scalar.
KOKKOS_INLINE_FUNCTION double trueLaplacian(double x, double y, double z, int order, int dimension)
Generalized Moving Least Squares (GMLS)
KOKKOS_INLINE_FUNCTION double curlTestSolution(double x, double y, double z, int component, int dimension)
KOKKOS_INLINE_FUNCTION double divergenceTestSamples(double x, double y, double z, int component, int dimension)
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.
KOKKOS_INLINE_FUNCTION double trueSolution(double x, double y, double z, int order, int dimension)
constexpr SamplingFunctional VectorPointSample
Point evaluations of the entire vector source function.