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);
40 bool all_passed =
true;
47 int number_of_batches = 1;
49 int arg8toi = atoi(args[7]);
51 number_of_batches = arg8toi;
59 int constraint_type = 1;
61 int arg7toi = atoi(args[6]);
63 constraint_type = arg7toi;
73 int arg6toi = atoi(args[5]);
75 problem_type = arg6toi;
86 int arg5toi = atoi(args[4]);
88 solver_type = arg5toi;
96 int arg4toi = atoi(args[3]);
105 int number_target_coords = 200;
107 int arg3toi = atoi(args[2]);
109 number_target_coords = arg3toi;
117 int arg2toi = atoi(args[1]);
125 const double failure_tolerance = 1e-9;
132 Kokkos::Profiling::pushRegion(
"Setup Point Data");
136 double h_spacing = 0.05;
137 int n_neg1_to_1 = 2*(1/h_spacing) + 1;
140 const int number_source_coords = std::pow(n_neg1_to_1, dimension);
143 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> source_coords_device(
"source coordinates",
144 number_source_coords, 3);
145 Kokkos::View<double**>::HostMirror source_coords = Kokkos::create_mirror_view(source_coords_device);
148 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> target_coords_device (
"target coordinates", number_target_coords, 3);
149 Kokkos::View<double**>::HostMirror target_coords = Kokkos::create_mirror_view(target_coords_device);
152 Kokkos::View<double***, Kokkos::DefaultExecutionSpace> tangent_bundles_device (
"tangent bundles", number_target_coords, dimension, dimension);
153 Kokkos::View<double***>::HostMirror tangent_bundles = Kokkos::create_mirror_view(tangent_bundles_device);
156 int source_index = 0;
157 double this_coord[3] = {0,0,0};
158 for (
int i=-n_neg1_to_1/2; i<n_neg1_to_1/2+1; ++i) {
159 this_coord[0] = i*h_spacing;
160 for (
int j=-n_neg1_to_1/2; j<n_neg1_to_1/2+1; ++j) {
161 this_coord[1] = j*h_spacing;
162 for (
int k=-n_neg1_to_1/2; k<n_neg1_to_1/2+1; ++k) {
163 this_coord[2] = k*h_spacing;
165 source_coords(source_index,0) = this_coord[0];
166 source_coords(source_index,1) = this_coord[1];
167 source_coords(source_index,2) = this_coord[2];
172 source_coords(source_index,0) = this_coord[0];
173 source_coords(source_index,1) = this_coord[1];
174 source_coords(source_index,2) = 0;
179 source_coords(source_index,0) = this_coord[0];
180 source_coords(source_index,1) = 0;
181 source_coords(source_index,2) = 0;
187 for(
int i=0; i<number_target_coords; i++){
190 double rand_dir[3] = {0,0,0};
192 for (
int j=0; j<dimension; ++j) {
194 rand_dir[j] = ((double)rand() / (double) RAND_MAX) - 0.5;
198 for (
int j=0; j<dimension; ++j) {
199 target_coords(i,j) = rand_dir[j];
204 if (dimension == 2) {
205 tangent_bundles(i, 0, 0) = 0.0;
206 tangent_bundles(i, 0, 1) = 0.0;
207 tangent_bundles(i, 1, 0) = 1.0/(sqrt(2.0));
208 tangent_bundles(i, 1, 1) = 1.0/(sqrt(2.0));
209 }
else if (dimension == 3) {
210 tangent_bundles(i, 0, 0) = 0.0;
211 tangent_bundles(i, 0, 1) = 0.0;
212 tangent_bundles(i, 0, 2) = 0.0;
213 tangent_bundles(i, 1, 0) = 0.0;
214 tangent_bundles(i, 1, 1) = 0.0;
215 tangent_bundles(i, 1, 2) = 0.0;
216 tangent_bundles(i, 2, 0) = 1.0/(sqrt(3.0));
217 tangent_bundles(i, 2, 1) = 1.0/(sqrt(3.0));
218 tangent_bundles(i, 2, 2) = 1.0/(sqrt(3.0));
224 Kokkos::Profiling::popRegion();
225 Kokkos::Profiling::pushRegion(
"Creating Data");
231 Kokkos::deep_copy(source_coords_device, source_coords);
234 Kokkos::deep_copy(target_coords_device, target_coords);
237 Kokkos::deep_copy(tangent_bundles_device, tangent_bundles);
240 Kokkos::View<double*, Kokkos::DefaultExecutionSpace> sampling_data_device(
"samples of true solution",
241 source_coords_device.extent(0));
243 Kokkos::parallel_for(
"Sampling Manufactured Solutions", Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>
244 (0,source_coords.extent(0)), KOKKOS_LAMBDA(
const int i) {
247 double xval = source_coords_device(i,0);
248 double yval = (dimension>1) ? source_coords_device(i,1) : 0;
249 double zval = (dimension>2) ? source_coords_device(i,2) : 0;
252 sampling_data_device(i) =
trueSolution(xval, yval, zval, order, dimension);
257 Kokkos::Profiling::popRegion();
258 Kokkos::Profiling::pushRegion(
"Neighbor Search");
269 double epsilon_multiplier = 1.8;
270 int estimated_upper_bound_number_neighbors =
271 point_cloud_search.getEstimatedNumberNeighborsUpperBound(min_neighbors, dimension, epsilon_multiplier);
273 Kokkos::View<int**, Kokkos::DefaultExecutionSpace> neighbor_lists_device(
"neighbor lists",
274 number_target_coords, estimated_upper_bound_number_neighbors);
275 Kokkos::View<int**>::HostMirror neighbor_lists = Kokkos::create_mirror_view(neighbor_lists_device);
278 Kokkos::View<double*, Kokkos::DefaultExecutionSpace> epsilon_device(
"h supports", number_target_coords);
279 Kokkos::View<double*>::HostMirror epsilon = Kokkos::create_mirror_view(epsilon_device);
284 point_cloud_search.generate2DNeighborListsFromKNNSearch(
false , target_coords, neighbor_lists,
285 epsilon, min_neighbors, epsilon_multiplier);
289 Kokkos::Profiling::popRegion();
300 Kokkos::deep_copy(neighbor_lists_device, neighbor_lists);
301 Kokkos::deep_copy(epsilon_device, epsilon);
304 std::string solver_name;
305 if (solver_type == 0) {
307 }
else if (solver_type == 1) {
309 }
else if (solver_type == 2) {
314 std::string problem_name;
315 if (problem_type == 0) {
316 problem_name =
"STANDARD";
317 }
else if (problem_type == 1) {
318 problem_name =
"MANIFOLD";
322 std::string constraint_name;
323 if (constraint_type == 0) {
324 constraint_name =
"NO_CONSTRAINT";
325 }
else if (constraint_type == 1) {
326 constraint_name =
"NEUMANN_GRAD_SCALAR";
333 solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
350 my_GMLS.
setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
351 my_GMLS.setTangentBundle(tangent_bundles_device);
358 my_GMLS.addTargets(lro);
364 my_GMLS.setWeightingPower(2);
367 my_GMLS.generateAlphas(number_of_batches);
371 double instantiation_time = timer.seconds();
372 std::cout <<
"Took " << instantiation_time <<
"s to complete normal vectors generation." << std::endl;
374 Kokkos::Profiling::pushRegion(
"Apply Alphas to Data");
394 Kokkos::Profiling::popRegion();
396 Kokkos::Profiling::pushRegion(
"Comparison");
401 for (
int i=0; i<number_target_coords; i++) {
403 double xval = target_coords(i,0);
404 double yval = (dimension>1) ? target_coords(i,1) : 0;
405 double zval = (dimension>2) ? target_coords(i,2) : 0;
408 int num_neigh_i = neighbor_lists(i, 0);
409 double b_i = my_GMLS.getAlpha0TensorTo0Tensor(lro, i, num_neigh_i);
412 double GMLS_value = output_value(i);
415 double actual_Laplacian =
trueLaplacian(xval, yval, zval, order, dimension);
418 double actual_Gradient[3] = {0,0,0};
419 trueGradient(actual_Gradient, xval, yval, zval, order, dimension);
420 double g = (dimension == 3) ? (1.0/sqrt(3.0))*(actual_Gradient[0] + actual_Gradient[1] + actual_Gradient[2])
421 : (1.0/sqrt(2.0))*(actual_Gradient[0] + actual_Gradient[1]);
422 double adjusted_value = GMLS_value + b_i*g;
425 if(GMLS_value!=GMLS_value || std::abs(actual_Laplacian - adjusted_value) > failure_tolerance) {
427 std::cout << i <<
" Failed Actual by: " << std::abs(actual_Laplacian - adjusted_value) << std::endl;
434 Kokkos::Profiling::popRegion();
441 #ifdef COMPADRE_USE_MPI 447 fprintf(stdout,
"Passed test \n");
450 fprintf(stdout,
"Failed test \n");
Class handling Kokkos command line arguments and returning parameters.
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...
Scalar polynomial basis centered at the target site and scaled by sum of basis powers e...
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)
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 laplacian of a scalar (could be on a manifold or not)
KOKKOS_INLINE_FUNCTION double trueLaplacian(double x, double y, double z, int order, int dimension)
Generalized Moving Least Squares (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) ...
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...
TargetOperation
Available target functionals.