9 #include <Compadre_Config.h>
16 #ifdef COMPADRE_USE_MPI
20 #include <Kokkos_Timer.hpp>
21 #include <Kokkos_Core.hpp>
23 using namespace Compadre;
28 int main (
int argc,
char* args[]) {
31 #ifdef COMPADRE_USE_MPI
32 MPI_Init(&argc, &args);
36 Kokkos::initialize(argc, args);
39 bool all_passed =
true;
48 int constraint_type = 0;
50 int arg7toi = atoi(args[6]);
52 constraint_type = arg7toi;
62 int arg6toi = atoi(args[5]);
64 problem_type = arg6toi;
75 int arg5toi = atoi(args[4]);
77 solver_type = arg5toi;
85 int arg4toi = atoi(args[3]);
93 int number_target_coords = 200;
95 int arg3toi = atoi(args[2]);
97 number_target_coords = arg3toi;
105 int arg2toi = atoi(args[1]);
113 const double failure_tolerance = 9e-8;
120 Kokkos::Profiling::pushRegion(
"Setup Point Data");
124 double h_spacing = 0.05;
125 int n_neg1_to_1 = 2*(1/h_spacing) + 1;
128 const int number_source_coords = std::pow(n_neg1_to_1, dimension);
131 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> source_coords_device(
"source coordinates",
132 number_source_coords, 3);
133 Kokkos::View<double**>::HostMirror source_coords = Kokkos::create_mirror_view(source_coords_device);
136 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> target_coords_device (
"target coordinates", number_target_coords, 3);
137 Kokkos::View<double**>::HostMirror target_coords = Kokkos::create_mirror_view(target_coords_device);
140 Kokkos::View<double***, Kokkos::DefaultExecutionSpace> tangent_bundles_device (
"tangent bundles", number_target_coords, 3, 3);
141 Kokkos::View<double***>::HostMirror tangent_bundles = Kokkos::create_mirror_view(tangent_bundles_device);
144 int source_index = 0;
145 double this_coord[3] = {0,0,0};
146 for (
int i=-n_neg1_to_1/2; i<n_neg1_to_1/2+1; ++i) {
147 this_coord[0] = i*h_spacing;
148 for (
int j=-n_neg1_to_1/2; j<n_neg1_to_1/2+1; ++j) {
149 this_coord[1] = j*h_spacing;
150 for (
int k=-n_neg1_to_1/2; k<n_neg1_to_1/2+1; ++k) {
151 this_coord[2] = k*h_spacing;
153 source_coords(source_index,0) = this_coord[0];
154 source_coords(source_index,1) = this_coord[1];
155 source_coords(source_index,2) = this_coord[2];
160 source_coords(source_index,0) = this_coord[0];
161 source_coords(source_index,1) = this_coord[1];
162 source_coords(source_index,2) = 0;
167 source_coords(source_index,0) = this_coord[0];
168 source_coords(source_index,1) = 0;
169 source_coords(source_index,2) = 0;
178 std::mt19937 rng(50);
180 std::uniform_int_distribution<int> gen_num_neighbours(0, source_index);
182 for (
int i=0; i<number_target_coords; ++i) {
183 const int source_site_to_copy = gen_num_neighbours(rng);
184 for (
int j=0; j<3; j++) {
185 target_coords(i, j) = source_coords(source_site_to_copy, j);
187 if (constraint_type == 1) {
188 if (dimension == 3) {
189 tangent_bundles(i, 0, 0) = 0.0;
190 tangent_bundles(i, 0, 1) = 0.0;
191 tangent_bundles(i, 0, 2) = 0.0;
192 tangent_bundles(i, 1, 0) = 0.0;
193 tangent_bundles(i, 1, 1) = 0.0;
194 tangent_bundles(i, 1, 2) = 0.0;
195 tangent_bundles(i, 2, 0) = 1.0/(sqrt(3.0));
196 tangent_bundles(i, 2, 1) = 1.0/(sqrt(3.0));
197 tangent_bundles(i, 2, 2) = 1.0/(sqrt(3.0));
199 if (dimension == 2) {
200 tangent_bundles(i, 0, 0) = 0.0;
201 tangent_bundles(i, 0, 1) = 0.0;
202 tangent_bundles(i, 1, 0) = 1.0/(sqrt(2.0));
203 tangent_bundles(i, 1, 1) = 1.0/(sqrt(2.0));
210 Kokkos::Profiling::popRegion();
211 Kokkos::Profiling::pushRegion(
"Creating Data");
217 Kokkos::deep_copy(source_coords_device, source_coords);
220 Kokkos::deep_copy(target_coords_device, target_coords);
223 Kokkos::View<double*, Kokkos::DefaultExecutionSpace> sampling_data_device(
"samples of true solution",
224 source_coords_device.extent(0));
226 Kokkos::parallel_for(
"Sampling Manufactured Solutions", Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>
227 (0,source_coords.extent(0)), KOKKOS_LAMBDA(
const int i) {
229 double xval = source_coords_device(i,0);
230 double yval = (dimension>1) ? source_coords_device(i,1) : 0;
231 double zval = (dimension>2) ? source_coords_device(i,2) : 0;
234 sampling_data_device(i) =
trueSolution(xval, yval, zval, order, dimension);
240 Kokkos::Profiling::popRegion();
241 Kokkos::Profiling::pushRegion(
"Neighbor Search");
253 double epsilon_multiplier = 2.2;
254 int estimated_upper_bound_number_neighbors =
255 point_cloud_search.getEstimatedNumberNeighborsUpperBound(min_neighbors, dimension, epsilon_multiplier);
257 Kokkos::View<int**, Kokkos::DefaultExecutionSpace> neighbor_lists_device(
"neighbor lists",
258 number_target_coords, estimated_upper_bound_number_neighbors);
259 Kokkos::View<int**>::HostMirror neighbor_lists = Kokkos::create_mirror_view(neighbor_lists_device);
262 Kokkos::View<double*, Kokkos::DefaultExecutionSpace> epsilon_device(
"h supports", number_target_coords);
263 Kokkos::View<double*>::HostMirror epsilon = Kokkos::create_mirror_view(epsilon_device);
268 point_cloud_search.generate2DNeighborListsFromKNNSearch(
false , target_coords, neighbor_lists,
269 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";
319 solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
327 solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
344 scalar_basis_gmls.
setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
345 vector_basis_gmls.setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
346 if (constraint_name ==
"NEUMANN_GRAD_SCALAR") {
347 scalar_basis_gmls.setTangentBundle(tangent_bundles_device);
348 vector_basis_gmls.setTangentBundle(tangent_bundles_device);
352 std::vector<TargetOperation> lro(2);
357 scalar_basis_gmls.addTargets(lro);
358 vector_basis_gmls.addTargets(lro);
365 scalar_basis_gmls.setWeightingPower(2);
366 vector_basis_gmls.setWeightingPower(2);
369 vector_basis_gmls.setOrderOfQuadraturePoints(order);
370 vector_basis_gmls.setDimensionOfQuadraturePoints(1);
371 vector_basis_gmls.setQuadratureType(
"LINE");
374 scalar_basis_gmls.generateAlphas();
375 vector_basis_gmls.generateAlphas();
379 double instantiation_time = timer.seconds();
380 std::cout <<
"Took " << instantiation_time <<
"s to complete alphas generation." << std::endl;
382 Kokkos::Profiling::pushRegion(
"Apply Alphas to Data");
396 Evaluator gmls_evaluator_scalar(&scalar_basis_gmls);
397 Evaluator gmls_evaluator_vector(&vector_basis_gmls);
415 Kokkos::Profiling::popRegion();
417 Kokkos::Profiling::pushRegion(
"Comparison");
422 for (
int i=0; i<number_target_coords; i++) {
425 double xval = target_coords(i,0);
426 double yval = (dimension>1) ? target_coords(i,1) : 0;
427 double zval = (dimension>2) ? target_coords(i,2) : 0;
430 double actual_Divergence;
431 actual_Divergence =
trueLaplacian(xval, yval, zval, order, dimension);
432 double actual_Gradient[3] = {0,0,0};
433 trueGradient(actual_Gradient, xval, yval, zval, order, dimension);
436 double g = (dimension == 3) ? (1.0/sqrt(3.0))*(actual_Gradient[0] + actual_Gradient[1] + actual_Gradient[2])
437 : (1.0/sqrt(2.0))*(actual_Gradient[0] + actual_Gradient[1]);
440 int num_neigh_i = neighbor_lists(i, 0);
443 double GMLS_Divergence_Scalar = output_divergence_scalar(i);
444 double GMLS_Divergence_Vector = output_divergence_vector(i);
447 if (constraint_name ==
"NEUMANN_GRAD_SCALAR") {
448 double b_i_scalar = scalar_basis_gmls.getAlpha0TensorTo0Tensor(DivergenceOfVectorPointEvaluation, i, num_neigh_i);
449 GMLS_Divergence_Scalar = GMLS_Divergence_Scalar + b_i_scalar*g;
451 double b_i_vector = vector_basis_gmls.getAlpha0TensorTo0Tensor(DivergenceOfVectorPointEvaluation, i, num_neigh_i);
452 GMLS_Divergence_Vector = GMLS_Divergence_Vector + b_i_vector*g;
456 if(std::abs(actual_Divergence - GMLS_Divergence_Scalar) > failure_tolerance) {
458 std::cout << i <<
" Failed Divergence on SCALAR basis by: " << std::abs(actual_Divergence - GMLS_Divergence_Scalar) << std::endl;
459 std::cout << i <<
" GMLS " << GMLS_Divergence_Scalar <<
" adjusted " << GMLS_Divergence_Scalar <<
" actual " << actual_Divergence << std::endl;
463 if(std::abs(actual_Divergence - GMLS_Divergence_Vector) > failure_tolerance) {
465 std::cout << i <<
" Failed Divergence on VECTOR basis by: " << std::abs(actual_Divergence - GMLS_Divergence_Vector) << std::endl;
466 std::cout << i <<
" GMLS " << GMLS_Divergence_Vector <<
" adjusted " << GMLS_Divergence_Vector <<
" actual " << actual_Divergence << std::endl;
470 double Scalar_GMLS_GradX = (dimension>1) ? output_gradient_scalar(i,0) : 0;
471 double Scalar_GMLS_GradY = (dimension>1) ? output_gradient_scalar(i,1) : 0;
472 double Scalar_GMLS_GradZ = (dimension>2) ? output_gradient_scalar(i,2) : 0;
474 double Vector_GMLS_GradX = (dimension>1) ? output_gradient_vector(i,0) : 0;
475 double Vector_GMLS_GradY = (dimension>1) ? output_gradient_vector(i,1) : 0;
476 double Vector_GMLS_GradZ = (dimension>2) ? output_gradient_vector(i,2) : 0;
479 if (constraint_name ==
"NEUMANN_GRAD_SCALAR") {
480 double bx_i_scalar = scalar_basis_gmls.getAlpha0TensorTo1Tensor(GradientOfScalarPointEvaluation, i, 0, num_neigh_i);
481 double by_i_scalar = scalar_basis_gmls.getAlpha0TensorTo1Tensor(GradientOfScalarPointEvaluation, i, 1, num_neigh_i);
482 double bz_i_scalar = scalar_basis_gmls.getAlpha0TensorTo1Tensor(GradientOfScalarPointEvaluation, i, 2, num_neigh_i);
483 Scalar_GMLS_GradX = Scalar_GMLS_GradX + bx_i_scalar*g;
484 Scalar_GMLS_GradY = Scalar_GMLS_GradY + by_i_scalar*g;
485 Scalar_GMLS_GradZ = Scalar_GMLS_GradZ + bz_i_scalar*g;
487 double bx_i_vector = vector_basis_gmls.getAlpha0TensorTo1Tensor(GradientOfScalarPointEvaluation, i, 0, num_neigh_i);
488 double by_i_vector = vector_basis_gmls.getAlpha0TensorTo1Tensor(GradientOfScalarPointEvaluation, i, 1, num_neigh_i);
489 double bz_i_vector = vector_basis_gmls.getAlpha0TensorTo1Tensor(GradientOfScalarPointEvaluation, i, 2, num_neigh_i);
490 Vector_GMLS_GradX = Vector_GMLS_GradX + bx_i_vector*g;
491 Vector_GMLS_GradY = Vector_GMLS_GradY + by_i_vector*g;
492 Vector_GMLS_GradZ = Vector_GMLS_GradZ + bz_i_vector*g;
496 if(std::abs(actual_Gradient[0] - Scalar_GMLS_GradX) > failure_tolerance) {
498 std::cout << i <<
" Failed Scalar GradX by: " << std::abs(actual_Gradient[0] - Scalar_GMLS_GradX) << std::endl;
500 if(std::abs(actual_Gradient[1] - Scalar_GMLS_GradY) > failure_tolerance) {
502 std::cout << i <<
" Failed Scalar GradY by: " << std::abs(actual_Gradient[1] - Scalar_GMLS_GradY) << std::endl;
506 if(std::abs(actual_Gradient[2] - Scalar_GMLS_GradZ) > failure_tolerance) {
508 std::cout << i <<
" Failed Scalar GradZ by: " << std::abs(actual_Gradient[2] - Scalar_GMLS_GradZ) << std::endl;
514 if(std::abs(actual_Gradient[0] - Vector_GMLS_GradX) > failure_tolerance) {
516 std::cout << i <<
" Failed Vector GradX by: " << std::abs(actual_Gradient[0] - Vector_GMLS_GradX) << std::endl;
518 if(std::abs(actual_Gradient[1] - Vector_GMLS_GradY) > failure_tolerance) {
520 std::cout << i <<
" Failed Vector GradY by: " << std::abs(actual_Gradient[1] - Vector_GMLS_GradY) << std::endl;
524 if(std::abs(actual_Gradient[2] - Vector_GMLS_GradZ) > failure_tolerance) {
526 std::cout << i <<
" Failed Vector GradZ by: " << std::abs(actual_Gradient[2] - Vector_GMLS_GradZ) << std::endl;
536 Kokkos::Profiling::popRegion();
544 #ifdef COMPADRE_USE_MPI
550 fprintf(stdout,
"Passed test \n");
553 fprintf(stdout,
"Failed test \n");
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...
KOKKOS_INLINE_FUNCTION void trueGradient(double *ans, double x, double y, double z, int order, int dimension)
constexpr SamplingFunctional StaggeredEdgeIntegralSample
Samples consist of the result of integrals of a vector dotted with the tangent along edges between ne...
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::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.
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 trueSolution(double x, double y, double z, int order, int dimension)
Vector polynomial basis having # of components _dimensions, or (_dimensions-1) in the case of manifol...