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);
39 bool all_passed =
true;
49 int constraint_type = 0;
51 int arg7toi = atoi(args[6]);
53 constraint_type = arg7toi;
63 int arg6toi = atoi(args[5]);
65 problem_type = arg6toi;
76 int arg5toi = atoi(args[4]);
78 solver_type = arg5toi;
86 int arg4toi = atoi(args[3]);
95 int number_target_coords = 200;
97 int arg3toi = atoi(args[2]);
99 number_target_coords = arg3toi;
107 int arg2toi = atoi(args[1]);
115 const double failure_tolerance = 1e-9;
122 Kokkos::Profiling::pushRegion(
"Setup Point Data");
126 double h_spacing = 0.05;
127 int n_neg1_to_1 = 2*(1/h_spacing) + 1;
130 const int number_source_coords = std::pow(n_neg1_to_1, dimension);
133 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> source_coords_device(
"source coordinates",
134 number_source_coords, 3);
135 Kokkos::View<double**>::HostMirror source_coords = Kokkos::create_mirror_view(source_coords_device);
138 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> target_coords_device (
"target coordinates", number_target_coords, 3);
139 Kokkos::View<double**>::HostMirror target_coords = Kokkos::create_mirror_view(target_coords_device);
142 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> additional_target_coords_device (
"additional target coordinates", 2*number_target_coords , 3);
143 Kokkos::View<double**>::HostMirror additional_target_coords = Kokkos::create_mirror_view(additional_target_coords_device);
146 Kokkos::View<int**, Kokkos::DefaultExecutionSpace> additional_target_indices_device (
"additional target indices", number_target_coords, 4 );
147 Kokkos::View<int**>::HostMirror additional_target_indices = Kokkos::create_mirror_view(additional_target_indices_device);
151 int source_index = 0;
152 double this_coord[3] = {0,0,0};
153 for (
int i=-n_neg1_to_1/2; i<n_neg1_to_1/2+1; ++i) {
154 this_coord[0] = i*h_spacing;
155 for (
int j=-n_neg1_to_1/2; j<n_neg1_to_1/2+1; ++j) {
156 this_coord[1] = j*h_spacing;
157 for (
int k=-n_neg1_to_1/2; k<n_neg1_to_1/2+1; ++k) {
158 this_coord[2] = k*h_spacing;
160 source_coords(source_index,0) = this_coord[0];
161 source_coords(source_index,1) = this_coord[1];
162 source_coords(source_index,2) = this_coord[2];
167 source_coords(source_index,0) = this_coord[0];
168 source_coords(source_index,1) = this_coord[1];
169 source_coords(source_index,2) = 0;
174 source_coords(source_index,0) = this_coord[0];
175 source_coords(source_index,1) = 0;
176 source_coords(source_index,2) = 0;
182 for(
int i=0; i<number_target_coords; i++){
185 double rand_dir[3] = {0,0,0};
187 for (
int j=0; j<dimension; ++j) {
189 rand_dir[j] = ((double)rand() / (double) RAND_MAX) - 0.5;
193 for (
int j=0; j<dimension; ++j) {
194 target_coords(i,j) = rand_dir[j];
202 int extra_evaluation_coordinates_count = 0;
203 for(
int i=0; i<number_target_coords; i++){
206 additional_target_indices(i,0) = (i%3)+1;
209 for (
int k=0; k<(i%3+1); ++k) {
210 for (
int j=0; j<dimension; ++j) {
211 additional_target_coords(extra_evaluation_coordinates_count,j) = target_coords(i,j) + (j==0)*1e-3 + (j==1)*1e-2 + (j==1)*(-1e-1);
213 additional_target_indices(i,k+1) = extra_evaluation_coordinates_count;
214 extra_evaluation_coordinates_count++;
221 Kokkos::Profiling::popRegion();
222 Kokkos::Profiling::pushRegion(
"Creating Data");
228 Kokkos::deep_copy(source_coords_device, source_coords);
231 Kokkos::deep_copy(target_coords_device, target_coords);
234 Kokkos::deep_copy(additional_target_coords_device, additional_target_coords);
237 Kokkos::deep_copy(additional_target_indices_device, additional_target_indices);
240 Kokkos::View<double*, Kokkos::DefaultExecutionSpace> sampling_data_device(
"samples of true solution",
241 source_coords_device.extent(0));
243 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> gradient_sampling_data_device(
"samples of true gradient",
244 source_coords_device.extent(0), dimension);
246 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> divergence_sampling_data_device
247 (
"samples of true solution for divergence test", source_coords_device.extent(0), dimension);
249 Kokkos::parallel_for(
"Sampling Manufactured Solutions", Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>
250 (0,source_coords.extent(0)), KOKKOS_LAMBDA(
const int i) {
253 double xval = source_coords_device(i,0);
254 double yval = (dimension>1) ? source_coords_device(i,1) : 0;
255 double zval = (dimension>2) ? source_coords_device(i,2) : 0;
258 sampling_data_device(i) =
trueSolution(xval, yval, zval, order, dimension);
261 double true_grad[3] = {0,0,0};
262 trueGradient(true_grad, xval, yval,zval, order, dimension);
264 for (
int j=0; j<dimension; ++j) {
265 gradient_sampling_data_device(i,j) = true_grad[j];
276 Kokkos::Profiling::popRegion();
277 Kokkos::Profiling::pushRegion(
"Neighbor Search");
288 double epsilon_multiplier = 1.5;
289 int estimated_upper_bound_number_neighbors =
290 point_cloud_search.getEstimatedNumberNeighborsUpperBound(min_neighbors, dimension, epsilon_multiplier);
292 Kokkos::View<int**, Kokkos::DefaultExecutionSpace> neighbor_lists_device(
"neighbor lists",
293 number_target_coords, estimated_upper_bound_number_neighbors);
294 Kokkos::View<int**>::HostMirror neighbor_lists = Kokkos::create_mirror_view(neighbor_lists_device);
297 Kokkos::View<double*, Kokkos::DefaultExecutionSpace> epsilon_device(
"h supports", number_target_coords);
298 Kokkos::View<double*>::HostMirror epsilon = Kokkos::create_mirror_view(epsilon_device);
303 point_cloud_search.generate2DNeighborListsFromKNNSearch(
false , target_coords, neighbor_lists,
304 epsilon, min_neighbors, epsilon_multiplier);
309 Kokkos::Profiling::popRegion();
320 Kokkos::deep_copy(neighbor_lists_device, neighbor_lists);
321 Kokkos::deep_copy(epsilon_device, epsilon);
324 std::string solver_name;
325 if (solver_type == 0) {
327 }
else if (solver_type == 1) {
329 }
else if (solver_type == 2) {
334 std::string problem_name;
335 if (problem_type == 0) {
336 problem_name =
"STANDARD";
337 }
else if (problem_type == 1) {
338 problem_name =
"MANIFOLD";
342 std::string constraint_name;
343 if (constraint_type == 0) {
344 constraint_name =
"NO_CONSTRAINT";
345 }
else if (constraint_type == 1) {
346 constraint_name =
"NEUMANN_GRAD_SCALAR";
352 solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
369 my_GMLS.
setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
372 my_GMLS.setAdditionalEvaluationSitesData(additional_target_indices_device, additional_target_coords_device);
375 std::vector<TargetOperation> lro(2);
380 my_GMLS.addTargets(lro);
386 my_GMLS.setWeightingPower(2);
389 my_GMLS.generateAlphas();
394 double instantiation_time = timer.seconds();
395 std::cout <<
"Took " << instantiation_time <<
"s to complete alphas generation." << std::endl;
397 Kokkos::Profiling::pushRegion(
"Apply Alphas to Data");
440 Kokkos::Profiling::popRegion();
442 Kokkos::Profiling::pushRegion(
"Comparison");
456 extra_evaluation_coordinates_count = 0;
457 for (
int i=0; i<number_target_coords; i++) {
459 for (
int k=0; k<(i%3)+1; ++k) {
462 GMLS_value = output_value1(i);
464 GMLS_GradX = output_gradient1(i,0);
466 GMLS_GradY = (dimension>1) ? output_gradient1(i,1) : 0;
468 GMLS_GradZ = (dimension>2) ? output_gradient1(i,2) : 0;
471 GMLS_value = output_value2(i);
473 GMLS_GradX = output_gradient2(i,0);
475 GMLS_GradY = (dimension>1) ? output_gradient2(i,1) : 0;
477 GMLS_GradZ = (dimension>2) ? output_gradient2(i,2) : 0;
480 GMLS_value = output_value3(i);
482 GMLS_GradX = output_gradient3(i,0);
484 GMLS_GradY = (dimension>1) ? output_gradient3(i,1) : 0;
486 GMLS_GradZ = (dimension>2) ? output_gradient3(i,2) : 0;
490 double xval = additional_target_coords(extra_evaluation_coordinates_count,0);
491 double yval = additional_target_coords(extra_evaluation_coordinates_count,1);
492 double zval = additional_target_coords(extra_evaluation_coordinates_count,2);
495 double actual_value =
trueSolution(xval, yval, zval, order, dimension);
497 double actual_Gradient[3] = {0,0,0};
498 trueGradient(actual_Gradient, xval, yval, zval, order, dimension);
501 if(GMLS_value!=GMLS_value || std::abs(actual_value - GMLS_value) > failure_tolerance) {
503 std::cout << i <<
" Failed Actual by: " << std::abs(actual_value - GMLS_value) <<
" for evaluation site: " << k << std::endl;
507 if(std::abs(actual_Gradient[0] - GMLS_GradX) > failure_tolerance) {
509 std::cout << i <<
" Failed GradX by: " << std::abs(actual_Gradient[0] - GMLS_GradX) <<
" for evaluation site: " << k << std::endl;
511 if(std::abs(actual_Gradient[1] - GMLS_GradY) > failure_tolerance) {
513 std::cout << i <<
" Failed GradY by: " << std::abs(actual_Gradient[1] - GMLS_GradY) <<
" for evaluation site: " << k << std::endl;
517 if(std::abs(actual_Gradient[2] - GMLS_GradZ) > failure_tolerance) {
519 std::cout << i <<
" Failed GradZ by: " << std::abs(actual_Gradient[2] - GMLS_GradZ) <<
" for evaluation site: " << k << std::endl;
523 extra_evaluation_coordinates_count++;
531 Kokkos::Profiling::popRegion();
540 #ifdef COMPADRE_USE_MPI 546 fprintf(stdout,
"Passed test \n");
549 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.
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.
KOKKOS_INLINE_FUNCTION double trueSolution(double x, double y, double z, int order, int dimension)
#define compadre_assert_release(condition)
compadre_assert_release is used for assertions that should always be checked, but generally are not e...
constexpr SamplingFunctional VectorPointSample
Point evaluations of the entire vector source function.