9 #include <Compadre_Config.h> 17 #ifdef COMPADRE_USE_MPI 21 #include <Kokkos_Timer.hpp> 22 #include <Kokkos_Core.hpp> 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;
48 auto order = clp.
order;
57 const double failure_tolerance = 1e-9;
64 Kokkos::Profiling::pushRegion(
"Setup Point Data");
68 double h_spacing = 0.05;
69 int n_neg1_to_1 = 2*(1/h_spacing) + 1;
72 const int number_source_coords = std::pow(n_neg1_to_1, dimension);
75 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> source_coords_device(
"source coordinates",
76 number_source_coords, 3);
77 Kokkos::View<double**>::HostMirror source_coords = Kokkos::create_mirror_view(source_coords_device);
80 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> target_coords_device (
"target coordinates", number_target_coords, 3);
81 Kokkos::View<double**>::HostMirror target_coords = Kokkos::create_mirror_view(target_coords_device);
84 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> additional_target_coords_device (
"additional target coordinates", 2*number_target_coords , 3);
85 Kokkos::View<double**>::HostMirror additional_target_coords = Kokkos::create_mirror_view(additional_target_coords_device);
88 Kokkos::View<int**, Kokkos::DefaultExecutionSpace> additional_target_indices_device (
"additional target indices", number_target_coords, 4 );
89 Kokkos::View<int**>::HostMirror additional_target_indices = Kokkos::create_mirror_view(additional_target_indices_device);
94 double this_coord[3] = {0,0,0};
95 for (
int i=-n_neg1_to_1/2; i<n_neg1_to_1/2+1; ++i) {
96 this_coord[0] = i*h_spacing;
97 for (
int j=-n_neg1_to_1/2; j<n_neg1_to_1/2+1; ++j) {
98 this_coord[1] = j*h_spacing;
99 for (
int k=-n_neg1_to_1/2; k<n_neg1_to_1/2+1; ++k) {
100 this_coord[2] = k*h_spacing;
102 source_coords(source_index,0) = this_coord[0];
103 source_coords(source_index,1) = this_coord[1];
104 source_coords(source_index,2) = this_coord[2];
109 source_coords(source_index,0) = this_coord[0];
110 source_coords(source_index,1) = this_coord[1];
111 source_coords(source_index,2) = 0;
116 source_coords(source_index,0) = this_coord[0];
117 source_coords(source_index,1) = 0;
118 source_coords(source_index,2) = 0;
124 for(
int i=0; i<number_target_coords; i++){
127 double rand_dir[3] = {0,0,0};
129 for (
int j=0; j<dimension; ++j) {
131 rand_dir[j] = ((double)rand() / (double) RAND_MAX) - 0.5;
135 for (
int j=0; j<dimension; ++j) {
136 target_coords(i,j) = rand_dir[j];
144 int extra_evaluation_coordinates_count = 0;
145 for(
int i=0; i<number_target_coords; i++){
148 additional_target_indices(i,0) = (i%3)+1;
151 for (
int k=0; k<(i%3+1); ++k) {
152 for (
int j=0; j<dimension; ++j) {
153 additional_target_coords(extra_evaluation_coordinates_count,j) = target_coords(i,j) + (j==0)*1e-3 + (j==1)*1e-2 + (j==1)*(-1e-1);
155 additional_target_indices(i,k+1) = extra_evaluation_coordinates_count;
156 extra_evaluation_coordinates_count++;
163 Kokkos::Profiling::popRegion();
164 Kokkos::Profiling::pushRegion(
"Creating Data");
170 Kokkos::deep_copy(source_coords_device, source_coords);
173 Kokkos::deep_copy(target_coords_device, target_coords);
176 Kokkos::deep_copy(additional_target_coords_device, additional_target_coords);
179 Kokkos::deep_copy(additional_target_indices_device, additional_target_indices);
182 Kokkos::View<double*, Kokkos::DefaultExecutionSpace> sampling_data_device(
"samples of true solution",
183 source_coords_device.extent(0));
185 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> gradient_sampling_data_device(
"samples of true gradient",
186 source_coords_device.extent(0), dimension);
188 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> divergence_sampling_data_device
189 (
"samples of true solution for divergence test", source_coords_device.extent(0), dimension);
191 Kokkos::parallel_for(
"Sampling Manufactured Solutions", Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>
192 (0,source_coords.extent(0)), KOKKOS_LAMBDA(
const int i) {
195 double xval = source_coords_device(i,0);
196 double yval = (dimension>1) ? source_coords_device(i,1) : 0;
197 double zval = (dimension>2) ? source_coords_device(i,2) : 0;
200 sampling_data_device(i) =
trueSolution(xval, yval, zval, order, dimension);
203 double true_grad[3] = {0,0,0};
204 trueGradient(true_grad, xval, yval,zval, order, dimension);
206 for (
int j=0; j<dimension; ++j) {
207 gradient_sampling_data_device(i,j) = true_grad[j];
218 Kokkos::Profiling::popRegion();
219 Kokkos::Profiling::pushRegion(
"Neighbor Search");
230 double epsilon_multiplier = 1.5;
231 int estimated_upper_bound_number_neighbors =
232 point_cloud_search.getEstimatedNumberNeighborsUpperBound(min_neighbors, dimension, epsilon_multiplier);
234 Kokkos::View<int**, Kokkos::DefaultExecutionSpace> neighbor_lists_device(
"neighbor lists",
235 number_target_coords, estimated_upper_bound_number_neighbors);
236 Kokkos::View<int**>::HostMirror neighbor_lists = Kokkos::create_mirror_view(neighbor_lists_device);
239 Kokkos::View<double*, Kokkos::DefaultExecutionSpace> epsilon_device(
"h supports", number_target_coords);
240 Kokkos::View<double*>::HostMirror epsilon = Kokkos::create_mirror_view(epsilon_device);
245 point_cloud_search.generate2DNeighborListsFromKNNSearch(
false , target_coords, neighbor_lists,
246 epsilon, min_neighbors, epsilon_multiplier);
251 Kokkos::Profiling::popRegion();
262 Kokkos::deep_copy(neighbor_lists_device, neighbor_lists);
263 Kokkos::deep_copy(epsilon_device, epsilon);
268 solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
285 my_GMLS.
setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
288 my_GMLS.setAdditionalEvaluationSitesData(additional_target_indices_device, additional_target_coords_device);
291 std::vector<TargetOperation> lro(2);
296 my_GMLS.addTargets(lro);
302 my_GMLS.setWeightingPower(2);
305 my_GMLS.generateAlphas();
310 double instantiation_time = timer.seconds();
311 std::cout <<
"Took " << instantiation_time <<
"s to complete alphas generation." << std::endl;
313 Kokkos::Profiling::pushRegion(
"Apply Alphas to Data");
356 Kokkos::Profiling::popRegion();
358 Kokkos::Profiling::pushRegion(
"Comparison");
372 extra_evaluation_coordinates_count = 0;
373 for (
int i=0; i<number_target_coords; i++) {
375 for (
int k=0; k<(i%3)+1; ++k) {
378 GMLS_value = output_value1(i);
380 GMLS_GradX = output_gradient1(i,0);
382 GMLS_GradY = (dimension>1) ? output_gradient1(i,1) : 0;
384 GMLS_GradZ = (dimension>2) ? output_gradient1(i,2) : 0;
387 GMLS_value = output_value2(i);
389 GMLS_GradX = output_gradient2(i,0);
391 GMLS_GradY = (dimension>1) ? output_gradient2(i,1) : 0;
393 GMLS_GradZ = (dimension>2) ? output_gradient2(i,2) : 0;
396 GMLS_value = output_value3(i);
398 GMLS_GradX = output_gradient3(i,0);
400 GMLS_GradY = (dimension>1) ? output_gradient3(i,1) : 0;
402 GMLS_GradZ = (dimension>2) ? output_gradient3(i,2) : 0;
406 double xval = additional_target_coords(extra_evaluation_coordinates_count,0);
407 double yval = additional_target_coords(extra_evaluation_coordinates_count,1);
408 double zval = additional_target_coords(extra_evaluation_coordinates_count,2);
411 double actual_value =
trueSolution(xval, yval, zval, order, dimension);
413 double actual_Gradient[3] = {0,0,0};
414 trueGradient(actual_Gradient, xval, yval, zval, order, dimension);
417 if(GMLS_value!=GMLS_value || std::abs(actual_value - GMLS_value) > failure_tolerance) {
419 std::cout << i <<
" Failed Actual by: " << std::abs(actual_value - GMLS_value) <<
" for evaluation site: " << k << std::endl;
423 if(std::abs(actual_Gradient[0] - GMLS_GradX) > failure_tolerance) {
425 std::cout << i <<
" Failed GradX by: " << std::abs(actual_Gradient[0] - GMLS_GradX) <<
" for evaluation site: " << k << std::endl;
427 if(std::abs(actual_Gradient[1] - GMLS_GradY) > failure_tolerance) {
429 std::cout << i <<
" Failed GradY by: " << std::abs(actual_Gradient[1] - GMLS_GradY) <<
" for evaluation site: " << k << std::endl;
433 if(std::abs(actual_Gradient[2] - GMLS_GradZ) > failure_tolerance) {
435 std::cout << i <<
" Failed GradZ by: " << std::abs(actual_Gradient[2] - GMLS_GradZ) <<
" for evaluation site: " << k << std::endl;
439 extra_evaluation_coordinates_count++;
447 Kokkos::Profiling::popRegion();
456 #ifdef COMPADRE_USE_MPI 462 fprintf(stdout,
"Passed test \n");
465 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...
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...
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)
int main(int argc, char *args[])
[Parse Command Line Arguments]
Point evaluation of the gradient of a scalar.
Generalized Moving Least Squares (GMLS)
KOKKOS_INLINE_FUNCTION double divergenceTestSamples(double x, double y, double z, int component, int dimension)
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) ...
constexpr SamplingFunctional PointSample
Available sampling functionals.
std::string constraint_name
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.