9 #include <Compadre_Config.h>
17 #ifdef COMPADRE_USE_MPI
21 #include <Kokkos_Timer.hpp>
22 #include <Kokkos_Core.hpp>
24 using namespace Compadre;
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;
47 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> tangent_bundles_device (
"tangent bundles", number_target_coords, dimension, dimension);
85 Kokkos::View<double***>::HostMirror tangent_bundles = Kokkos::create_mirror_view(tangent_bundles_device);
89 double this_coord[3] = {0,0,0};
90 for (
int i=-n_neg1_to_1/2; i<n_neg1_to_1/2+1; ++i) {
91 this_coord[0] = i*h_spacing;
92 for (
int j=-n_neg1_to_1/2; j<n_neg1_to_1/2+1; ++j) {
93 this_coord[1] = j*h_spacing;
94 for (
int k=-n_neg1_to_1/2; k<n_neg1_to_1/2+1; ++k) {
95 this_coord[2] = k*h_spacing;
97 source_coords(source_index,0) = this_coord[0];
98 source_coords(source_index,1) = this_coord[1];
99 source_coords(source_index,2) = this_coord[2];
104 source_coords(source_index,0) = this_coord[0];
105 source_coords(source_index,1) = this_coord[1];
106 source_coords(source_index,2) = 0;
111 source_coords(source_index,0) = this_coord[0];
112 source_coords(source_index,1) = 0;
113 source_coords(source_index,2) = 0;
119 for(
int i=0; i<number_target_coords; i++){
122 double rand_dir[3] = {0,0,0};
124 for (
int j=0; j<dimension; ++j) {
126 rand_dir[j] = ((double)rand() / (double) RAND_MAX) - 0.5;
130 for (
int j=0; j<dimension; ++j) {
131 target_coords(i,j) = rand_dir[j];
136 if (dimension == 2) {
137 tangent_bundles(i, 0, 0) = 0.0;
138 tangent_bundles(i, 0, 1) = 0.0;
139 tangent_bundles(i, 1, 0) = 1.0/(sqrt(2.0));
140 tangent_bundles(i, 1, 1) = 1.0/(sqrt(2.0));
141 }
else if (dimension == 3) {
142 tangent_bundles(i, 0, 0) = 0.0;
143 tangent_bundles(i, 0, 1) = 0.0;
144 tangent_bundles(i, 0, 2) = 0.0;
145 tangent_bundles(i, 1, 0) = 0.0;
146 tangent_bundles(i, 1, 1) = 0.0;
147 tangent_bundles(i, 1, 2) = 0.0;
148 tangent_bundles(i, 2, 0) = 1.0/(sqrt(3.0));
149 tangent_bundles(i, 2, 1) = 1.0/(sqrt(3.0));
150 tangent_bundles(i, 2, 2) = 1.0/(sqrt(3.0));
156 Kokkos::Profiling::popRegion();
157 Kokkos::Profiling::pushRegion(
"Creating Data");
163 Kokkos::deep_copy(source_coords_device, source_coords);
166 Kokkos::deep_copy(target_coords_device, target_coords);
169 Kokkos::deep_copy(tangent_bundles_device, tangent_bundles);
172 Kokkos::View<double*, Kokkos::DefaultExecutionSpace> sampling_data_device(
"samples of true solution",
173 source_coords_device.extent(0));
175 Kokkos::parallel_for(
"Sampling Manufactured Solutions", Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>
176 (0,source_coords.extent(0)), KOKKOS_LAMBDA(
const int i) {
179 double xval = source_coords_device(i,0);
180 double yval = (dimension>1) ? source_coords_device(i,1) : 0;
181 double zval = (dimension>2) ? source_coords_device(i,2) : 0;
184 sampling_data_device(i) =
trueSolution(xval, yval, zval, order, dimension);
189 Kokkos::Profiling::popRegion();
190 Kokkos::Profiling::pushRegion(
"Neighbor Search");
201 double epsilon_multiplier = 1.8;
202 int estimated_upper_bound_number_neighbors =
203 point_cloud_search.getEstimatedNumberNeighborsUpperBound(min_neighbors, dimension, epsilon_multiplier);
205 Kokkos::View<int**, Kokkos::DefaultExecutionSpace> neighbor_lists_device(
"neighbor lists",
206 number_target_coords, estimated_upper_bound_number_neighbors);
207 Kokkos::View<int**>::HostMirror neighbor_lists = Kokkos::create_mirror_view(neighbor_lists_device);
210 Kokkos::View<double*, Kokkos::DefaultExecutionSpace> epsilon_device(
"h supports", number_target_coords);
211 Kokkos::View<double*>::HostMirror epsilon = Kokkos::create_mirror_view(epsilon_device);
216 point_cloud_search.generate2DNeighborListsFromKNNSearch(
false , target_coords, neighbor_lists,
217 epsilon, min_neighbors, epsilon_multiplier);
221 Kokkos::Profiling::popRegion();
232 Kokkos::deep_copy(neighbor_lists_device, neighbor_lists);
233 Kokkos::deep_copy(epsilon_device, epsilon);
239 solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
256 my_GMLS.
setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
257 my_GMLS.setTangentBundle(tangent_bundles_device);
264 my_GMLS.addTargets(lro);
270 my_GMLS.setWeightingParameter(2);
273 my_GMLS.generateAlphas(number_of_batches);
277 double instantiation_time = timer.seconds();
278 std::cout <<
"Took " << instantiation_time <<
"s to complete normal vectors generation." << std::endl;
280 Kokkos::Profiling::pushRegion(
"Apply Alphas to Data");
300 Kokkos::Profiling::popRegion();
302 Kokkos::Profiling::pushRegion(
"Comparison");
307 for (
int i=0; i<number_target_coords; i++) {
309 double xval = target_coords(i,0);
310 double yval = (dimension>1) ? target_coords(i,1) : 0;
311 double zval = (dimension>2) ? target_coords(i,2) : 0;
314 int num_neigh_i = neighbor_lists(i, 0);
315 double b_i = my_GMLS.getSolutionSetHost()->getAlpha0TensorTo0Tensor(lro, i, num_neigh_i);
318 double GMLS_value = output_value(i);
321 double actual_Laplacian =
trueLaplacian(xval, yval, zval, order, dimension);
324 double actual_Gradient[3] = {0,0,0};
325 trueGradient(actual_Gradient, xval, yval, zval, order, dimension);
326 double g = (dimension == 3) ? (1.0/sqrt(3.0))*(actual_Gradient[0] + actual_Gradient[1] + actual_Gradient[2])
327 : (1.0/sqrt(2.0))*(actual_Gradient[0] + actual_Gradient[1]);
328 double adjusted_value = GMLS_value + b_i*g;
331 if(GMLS_value!=GMLS_value || std::abs(actual_Laplacian - adjusted_value) > failure_tolerance) {
333 std::cout << i <<
" Failed Actual by: " << std::abs(actual_Laplacian - adjusted_value) << std::endl;
340 Kokkos::Profiling::popRegion();
347 #ifdef COMPADRE_USE_MPI
353 fprintf(stdout,
"Passed test \n");
356 fprintf(stdout,
"Failed test \n");
Lightweight Evaluator Helper This class is a lightweight wrapper for extracting and applying all rele...
Point evaluation of the laplacian of a scalar (could be on a manifold or not)
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...
int main(int argc, char **argv)
TargetOperation
Available target functionals.
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 polynomial basis centered at the target site and scaled by sum of basis powers e...
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 (allocates memory for output)
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)