Compadre  1.3.3
GMLS_Host.cpp
Go to the documentation of this file.
1 #include <iostream>
2 #include <string>
3 #include <vector>
4 #include <map>
5 #include <stdlib.h>
6 #include <cstdio>
7 #include <random>
8 
9 #include <Compadre_Config.h>
10 #include <Compadre_GMLS.hpp>
11 #include <Compadre_Evaluator.hpp>
12 #include "GMLS_Tutorial.hpp"
13 #include "CommandLineProcessor.hpp"
14 
15 #ifdef COMPADRE_USE_MPI
16 #include <mpi.h>
17 #endif
18 
19 #include <Kokkos_Timer.hpp>
20 #include <Kokkos_Core.hpp>
21 
22 using namespace Compadre;
23 
24 int main (int argc, char* args[])
25 {
26 
27 #ifdef COMPADRE_USE_MPI
28  MPI_Init(&argc, &args);
29 #endif
30 
31 bool all_passed = true;
32 
33 {
34 
35  CommandLineProcessor clp(argc, args);
36  auto order = clp.order;
37  auto dimension = clp.dimension;
38  auto number_target_coords = clp.number_target_coords;
39  auto constraint_name = clp.constraint_name;
40  auto solver_name = clp.solver_name;
41  auto problem_name = clp.problem_name;
42 
43  const double failure_tolerance = 1e-9;
44 
45  const int offset = 15;
46  std::mt19937 rng(50);
47  const int min_neighbors = 1*Compadre::GMLS::getNP(order);
48  const int max_neighbors = 1*Compadre::GMLS::getNP(order)*1.15;
49  std::cout << min_neighbors << " " << max_neighbors << std::endl;
50  std::uniform_int_distribution<int> gen_num_neighbors(min_neighbors, max_neighbors); // uniform, unbiased
51 
52 
53  Kokkos::initialize(argc, args);
54  Kokkos::Timer timer;
55  Kokkos::Profiling::pushRegion("Setup");
56 
57 
58  const int N = 40000;
59  std::uniform_int_distribution<int> gen_neighbor_number(offset, N); // 0 to 10 are junk (part of test)
60 
61 
62  Kokkos::View<int**, Kokkos::HostSpace> neighbor_lists("neighbor lists", number_target_coords, max_neighbors+1); // first column is # of neighbors
63  Kokkos::View<double**, Kokkos::HostSpace> source_coords("neighbor coordinates", N, dimension);
64  Kokkos::View<double*, Kokkos::HostSpace> epsilon("h supports", number_target_coords);
65 
66  for (int i=0; i<number_target_coords; i++) {
67  epsilon(i) = 0.5;
68  }
69 
70 // // fake coordinates not to be used
71  for(int i = 0; i < offset; i++){
72  for(int j = 0; j < dimension; j++){
73  source_coords(i,j) = 0.1;
74  }
75  }
76 
77  // filling others with random coordinates
78  for(int i = offset; i < N; i++){ //ignore first ten entries
79  double randx = (2.0*(double)rand() / (double) RAND_MAX - 1.0)*epsilon(0)/2.0;
80  double randy = (2.0*(double)rand() / (double) RAND_MAX - 1.0)*epsilon(0)/2.0;
81  double randz = (2.0*(double)rand() / (double) RAND_MAX - 1.0)*epsilon(0)/2.0;
82  source_coords(i,0) = randx;
83  if (dimension>1) source_coords(i,1) = randy;
84  if (dimension>2) source_coords(i,2) = randz;
85  }
86 
87  const double target_epsilon = 0.1;
88  // fill target coords
89  Kokkos::View<double**, Kokkos::HostSpace> target_coords ("target coordinates", number_target_coords, dimension);
90  for(int i = 0; i < number_target_coords; i++){ //ignore first ten entries
91  double randx = (2.0*(double)rand() / (double) RAND_MAX - 1.0)*target_epsilon/2.0;
92  double randy = (2.0*(double)rand() / (double) RAND_MAX - 1.0)*target_epsilon/2.0;
93  double randz = (2.0*(double)rand() / (double) RAND_MAX - 1.0)*target_epsilon/2.0;
94  target_coords(i,0) = randx;
95  if (dimension>1) target_coords(i,1) = randy;
96  if (dimension>2) target_coords(i,2) = randz;
97  }
98 
99  // randomly fill neighbor lists
100  for (int i=0; i<number_target_coords; i++) {
101 // int r = gen_num_neighbors(rng);
102 // assert(r<source_coords.extent(0)-offset);
103  int r = max_neighbors;
104  neighbor_lists(i,0) = r; // number of neighbors is random between max and min
105 
106  for(int j=0; j<r; j++){
107  neighbor_lists(i,j+1) = offset + j + 1;
108 // bool placed = false;
109 // while (!placed) {
110 // int ind = gen_neighbor_number(rng);
111 // bool located = false;
112 // for (int k=1; k<j+1; k++) {
113 // if (neighbor_lists(i,k) == ind) {
114 // located = true;
115 // break;
116 // }
117 // }
118 // if (!located) {
119 // neighbor_lists(i,j+1) = ind;
120 // placed = true;
121 // } // neighbors can be from 10,...,N-1
122 // }
123  }
124  }
125 
126  Kokkos::Profiling::popRegion();
127  timer.reset();
128 
129  GMLS my_GMLS(order, dimension,
130  solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
131  2 /*manifold order*/);
132  my_GMLS.setProblemData(neighbor_lists, source_coords, target_coords, epsilon);
133  my_GMLS.setWeightingPower(10);
134 
135  std::vector<TargetOperation> lro(5);
136  lro[0] = ScalarPointEvaluation;
141  my_GMLS.addTargets(lro);
142  my_GMLS.generateAlphas();
143 
144  double instantiation_time = timer.seconds();
145  std::cout << "Took " << instantiation_time << "s to complete instantiation." << std::endl;
146 
147  Kokkos::Profiling::pushRegion("Creating Data");
148 
149 
150  // need Kokkos View storing true solution
151  Kokkos::View<double*, Kokkos::HostSpace> sampling_data("samples of true solution", source_coords.extent(0));
152  Kokkos::View<double**, Kokkos::HostSpace> gradient_sampling_data("samples of true gradient", source_coords.extent(0), dimension);
153  Kokkos::View<double**, Kokkos::LayoutLeft, Kokkos::HostSpace> divergence_sampling_data("samples of true solution for divergence test", source_coords.extent(0), dimension);
154  Kokkos::parallel_for("Sampling Manufactured Solutions", Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>(0,source_coords.extent(0)), KOKKOS_LAMBDA(const int i) {
155  double xval = source_coords(i,0);
156  double yval = (dimension>1) ? source_coords(i,1) : 0;
157  double zval = (dimension>2) ? source_coords(i,2) : 0;
158  sampling_data(i) = trueSolution(xval, yval, zval, order, dimension);
159  double true_grad[3] = {0,0,0};
160  trueGradient(true_grad, xval, yval,zval, order, dimension);
161  for (int j=0; j<dimension; ++j) {
162  divergence_sampling_data(i,j) = divergenceTestSamples(xval, yval, zval, j, dimension);
163  gradient_sampling_data(i,j) = true_grad[j];
164  }
165  });
166  Kokkos::Profiling::popRegion();
167 
168  Evaluator gmls_evaluator(&my_GMLS);
169 
170  for (int i=0; i<number_target_coords; i++) {
171 
172  Kokkos::Profiling::pushRegion("Apply Alphas to Data");
173 
174  double GMLS_value = gmls_evaluator.applyAlphasToDataSingleComponentSingleTargetSite(sampling_data, 0, ScalarPointEvaluation, i, 0, 0, 0, 0, 0);
175  //for (int j = 0; j< neighbor_lists(i,0); j++){
176  // double xval = source_coords(neighbor_lists(i,j+1),0);
177  // double yval = (dimension>1) ? source_coords(neighbor_lists(i,j+1),1) : 0;
178  // double zval = (dimension>2) ? source_coords(neighbor_lists(i,j+1),2) : 0;
179  // GMLS_value += gmls_evaluator.getAlpha0TensorTo0Tensor(ScalarPointEvaluation, i, j)*trueSolution(xval, yval, zval, order, dimension);
180  //}
181 
182  double GMLS_Laplacian = gmls_evaluator.applyAlphasToDataSingleComponentSingleTargetSite(sampling_data, 0, LaplacianOfScalarPointEvaluation, i, 0, 0, 0, 0, 0);
183  //double GMLS_Laplacian = 0.0;
184  //for (int j = 0; j< neighbor_lists(i,0); j++){
185  // double xval = source_coords(neighbor_lists(i,j+1),0);
186  // double yval = (dimension>1) ? source_coords(neighbor_lists(i,j+1),1) : 0;
187  // double zval = (dimension>2) ? source_coords(neighbor_lists(i,j+1),2) : 0;
188  // GMLS_Laplacian += gmls_evaluator.getAlpha0TensorTo0Tensor(LaplacianOfScalarPointEvaluation, i, j)*trueSolution(xval, yval, zval, order, dimension);
189  //}
190 
191  double GMLS_GradX = gmls_evaluator.applyAlphasToDataSingleComponentSingleTargetSite(sampling_data, 0, GradientOfScalarPointEvaluation, i, 0, 0, 0, 0, 0);
192  //double GMLS_GradX = 0.0;
193  //for (int j = 0; j< neighbor_lists(i,0); j++){
194  // double xval = source_coords(neighbor_lists(i,j+1),0);
195  // double yval = (dimension>1) ? source_coords(neighbor_lists(i,j+1),1) : 0;
196  // double zval = (dimension>2) ? source_coords(neighbor_lists(i,j+1),2) : 0;
197  // GMLS_GradX += gmls_evaluator.getAlpha0TensorTo1Tensor(GradientOfScalarPointEvaluation, i, 0, j)*trueSolution(xval, yval, zval, order, dimension);
198  //}
199 
200  double GMLS_GradY = (dimension>1) ? gmls_evaluator.applyAlphasToDataSingleComponentSingleTargetSite(sampling_data, 0, GradientOfScalarPointEvaluation, i, 0, 1, 0, 0, 0) : 0;
201  //double GMLS_GradY = 0.0;
202  //if (dimension>1) {
203  // for (int j = 0; j< neighbor_lists(i,0); j++){
204  // double xval = source_coords(neighbor_lists(i,j+1),0);
205  // double yval = source_coords(neighbor_lists(i,j+1),1);
206  // double zval = (dimension>2) ? source_coords(neighbor_lists(i,j+1),2) : 0;
207  // GMLS_GradY += gmls_evaluator.getAlpha0TensorTo1Tensor(GradientOfScalarPointEvaluation, i, 1, j)*trueSolution(xval, yval, zval, order, dimension);
208  // }
209  //}
210 
211  double GMLS_GradZ = (dimension>2) ? gmls_evaluator.applyAlphasToDataSingleComponentSingleTargetSite(sampling_data, 0, GradientOfScalarPointEvaluation, i, 0, 2, 0, 0, 0) : 0;
212  //double GMLS_GradZ = 0.0;
213  //if (dimension>2) {
214  // for (int j = 0; j< neighbor_lists(i,0); j++){
215  // double xval = source_coords(neighbor_lists(i,j+1),0);
216  // double yval = source_coords(neighbor_lists(i,j+1),1);
217  // double zval = source_coords(neighbor_lists(i,j+1),2);
218  // GMLS_GradZ += gmls_evaluator.getAlpha0TensorTo1Tensor(GradientOfScalarPointEvaluation, i, 2, j)*trueSolution(xval, yval, zval, order, dimension);
219  // }
220  //}
221 
222  double GMLS_Divergence = gmls_evaluator.applyAlphasToDataSingleComponentSingleTargetSite(gradient_sampling_data, 0, DivergenceOfVectorPointEvaluation, i, 0, 0, 0, 0, 0);
223  if (dimension>1) GMLS_Divergence += gmls_evaluator.applyAlphasToDataSingleComponentSingleTargetSite(gradient_sampling_data, 1, DivergenceOfVectorPointEvaluation, i, 0, 0, 0, 1, 0);
224  if (dimension>2) GMLS_Divergence += gmls_evaluator.applyAlphasToDataSingleComponentSingleTargetSite(gradient_sampling_data, 2, DivergenceOfVectorPointEvaluation, i, 0, 0, 0, 2, 0);
225 
226  //double GMLS_Divergence = gmls_evaluator.applyAlphasToDataSingleComponentSingleTargetSite(Kokkos::subview(gradient_sampling_data,Kokkos::ALL,0), 0, DivergenceOfVectorPointEvaluation, i, 0, 0, 0, 0);
227  //if (dimension>1) GMLS_Divergence += gmls_evaluator.applyAlphasToDataSingleComponentSingleTargetSite(Kokkos::subview(gradient_sampling_data,Kokkos::ALL,1), 1, DivergenceOfVectorPointEvaluation, i, 0, 0, 1, 0);
228  //if (dimension>2) GMLS_Divergence += gmls_evaluator.applyAlphasToDataSingleComponentSingleTargetSite(Kokkos::subview(gradient_sampling_data,Kokkos::ALL,2), 2, DivergenceOfVectorPointEvaluation, i, 0, 0, 2, 0);
229  //double GMLS_Divergence = 0.0;
230  //for (int j = 0; j< neighbor_lists(i,0); j++){
231  // double xval = source_coords(neighbor_lists(i,j+1),0);
232  // double yval = (dimension>1) ? source_coords(neighbor_lists(i,j+1),1) : 0;
233  // double zval = (dimension>2) ? source_coords(neighbor_lists(i,j+1),2) : 0;
234  // // TODO: use different functions for the vector components
235  // if (use_arbitrary_order_divergence) {
236  // GMLS_Divergence += gmls_evaluator.getAlpha1TensorTo0Tensor(DivergenceOfVectorPointEvaluation, i, j, 0)*trueSolution(xval, yval, zval, order, dimension);
237  // if (dimension>1) GMLS_Divergence += gmls_evaluator.getAlpha1TensorTo0Tensor(DivergenceOfVectorPointEvaluation, i, j, 1)*trueSolution(xval, yval, zval, order, dimension);
238  // if (dimension>2) GMLS_Divergence += gmls_evaluator.getAlpha1TensorTo0Tensor(DivergenceOfVectorPointEvaluation, i, j, 2)*trueSolution(xval, yval, zval, order, dimension);
239  // } else {
240  // for (int k=0; k<dimension; ++k) {
241  // GMLS_Divergence += gmls_evaluator.getAlpha1TensorTo0Tensor(DivergenceOfVectorPointEvaluation, i, j, k)*divergenceTestSamples(xval, yval, zval, k, dimension);
242  // }
243  // }
244  //}
245 
246  double GMLS_CurlX = 0.0;
247  double GMLS_CurlY = 0.0;
248  double GMLS_CurlZ = 0.0;
249  if (dimension>1) {
250  for (int j=0; j<dimension; ++j) {
251  GMLS_CurlX += gmls_evaluator.applyAlphasToDataSingleComponentSingleTargetSite(divergence_sampling_data, j, CurlOfVectorPointEvaluation, i, 0, 0, 0, j, 0);
252  GMLS_CurlY += gmls_evaluator.applyAlphasToDataSingleComponentSingleTargetSite(divergence_sampling_data, j, CurlOfVectorPointEvaluation, i, 0, 1, 0, j, 0);
253  }
254  }
255 
256  if (dimension>2) {
257  for (int j=0; j<dimension; ++j) {
258  GMLS_CurlZ += gmls_evaluator.applyAlphasToDataSingleComponentSingleTargetSite(divergence_sampling_data, j, CurlOfVectorPointEvaluation, i, 0, 2, 0, j, 0);
259  }
260 
261  }
262 
263  Kokkos::Profiling::popRegion();
264  //if (dimension>1) {
265  // for (int j = 0; j< neighbor_lists(i,0); j++){
266  // double xval = source_coords(neighbor_lists(i,j+1),0);
267  // double yval = source_coords(neighbor_lists(i,j+1),1);
268  // double zval = (dimension>2) ? source_coords(neighbor_lists(i,j+1),2) : 0;
269  // for (int k=0; k<dimension; ++k) {
270  // GMLS_CurlX += my_GMLS.getAlpha1TensorTo1Tensor(CurlOfVectorPointEvaluation, i, 0, j, k)*divergenceTestSamples(xval, yval, zval, k, dimension);
271  // }
272  // }
273 
274  // for (int j = 0; j< neighbor_lists(i,0); j++){
275  // double xval = source_coords(neighbor_lists(i,j+1),0);
276  // double yval = source_coords(neighbor_lists(i,j+1),1);
277  // double zval = (dimension>2) ? source_coords(neighbor_lists(i,j+1),2) : 0;
278  // for (int k=0; k<dimension; ++k) {
279  // GMLS_CurlY += my_GMLS.getAlpha1TensorTo1Tensor(CurlOfVectorPointEvaluation, i, 1, j, k)*divergenceTestSamples(xval, yval, zval, k, dimension);
280  // }
281  // }
282  //}
283 
284  //if (dimension>2) {
285  // for (int j = 0; j< neighbor_lists(i,0); j++){
286  // double xval = source_coords(neighbor_lists(i,j+1),0);
287  // double yval = source_coords(neighbor_lists(i,j+1),1);
288  // double zval = source_coords(neighbor_lists(i,j+1),2);
289  // for (int k=0; k<dimension; ++k) {
290  // GMLS_CurlZ += my_GMLS.getAlpha1TensorTo1Tensor(CurlOfVectorPointEvaluation, i, 2, j, k)*divergenceTestSamples(xval, yval, zval, k, dimension);
291  // }
292  // }
293  //}
294  //
295  Kokkos::Profiling::pushRegion("Comparison");
296 
297  double xval = target_coords(i,0);
298  double yval = (dimension>1) ? target_coords(i,1) : 0;
299  double zval = (dimension>2) ? target_coords(i,2) : 0;
300 
301  double actual_value = trueSolution(xval, yval, zval, order, dimension);
302  double actual_Laplacian = trueLaplacian(xval, yval, zval, order, dimension);
303  double actual_Gradient[3] = {0,0,0};
304  trueGradient(actual_Gradient, xval, yval, zval, order, dimension);
305  double actual_Divergence;
306  actual_Divergence = trueLaplacian(xval, yval, zval, order, dimension);
307 
308  double actual_CurlX = 0;
309  double actual_CurlY = 0;
310  double actual_CurlZ = 0;
311  if (dimension>1) {
312  actual_CurlX = curlTestSolution(xval, yval, zval, 0, dimension);
313  actual_CurlY = curlTestSolution(xval, yval, zval, 1, dimension);
314  }
315  if (dimension>2) {
316  actual_CurlZ = curlTestSolution(xval, yval, zval, 2, dimension);
317  }
318 
319 // fprintf(stdout, "Reconstructed value: %f \n", GMLS_value);
320 // fprintf(stdout, "Actual value: %f \n", actual_value);
321 // fprintf(stdout, "Reconstructed Laplacian: %f \n", GMLS_Laplacian);
322 // fprintf(stdout, "Actual Laplacian: %f \n", actual_Laplacian);
323 
324  if(GMLS_value!=GMLS_value || std::abs(actual_value - GMLS_value) > failure_tolerance) {
325  all_passed = false;
326  std::cout << "Failed Actual by: " << std::abs(actual_value - GMLS_value) << std::endl;
327  }
328 
329  if(std::abs(actual_Laplacian - GMLS_Laplacian) > failure_tolerance) {
330  all_passed = false;
331  std::cout << "Failed Laplacian by: " << std::abs(actual_Laplacian - GMLS_Laplacian) << std::endl;
332  }
333 
334  if(std::abs(actual_Gradient[0] - GMLS_GradX) > failure_tolerance) {
335  all_passed = false;
336  std::cout << "Failed GradX by: " << std::abs(actual_Gradient[0] - GMLS_GradX) << std::endl;
337  }
338 
339  if (dimension>1) {
340  if(std::abs(actual_Gradient[1] - GMLS_GradY) > failure_tolerance) {
341  all_passed = false;
342  std::cout << "Failed GradY by: " << std::abs(actual_Gradient[1] - GMLS_GradY) << std::endl;
343  }
344  }
345 
346  if (dimension>2) {
347  if(std::abs(actual_Gradient[2] - GMLS_GradZ) > failure_tolerance) {
348  all_passed = false;
349  std::cout << "Failed GradZ by: " << std::abs(actual_Gradient[2] - GMLS_GradZ) << std::endl;
350  }
351  }
352 
353  if(std::abs(actual_Divergence - GMLS_Divergence) > failure_tolerance) {
354  all_passed = false;
355  std::cout << "Failed Divergence by: " << std::abs(actual_Divergence - GMLS_Divergence) << std::endl;
356  }
357 
358  double tmp_diff = 0;
359  if (dimension>1)
360  tmp_diff += std::abs(actual_CurlX - GMLS_CurlX) + std::abs(actual_CurlY - GMLS_CurlY);
361  if (dimension>2)
362  tmp_diff += std::abs(actual_CurlZ - GMLS_CurlZ);
363  if(std::abs(tmp_diff) > failure_tolerance) {
364  all_passed = false;
365  std::cout << "Failed Curl by: " << std::abs(tmp_diff) << std::endl;
366  }
367  Kokkos::Profiling::popRegion();
368  }
369 
370 }
371 
372  Kokkos::finalize();
373 #ifdef COMPADRE_USE_MPI
374  MPI_Finalize();
375 #endif
376 
377 if(all_passed) {
378  fprintf(stdout, "Passed test \n");
379  return 0;
380 } else {
381  fprintf(stdout, "Failed test \n");
382  return -1;
383 }
384 }
int main(int argc, char *args[])
Manifold GMLS Example.
Definition: GMLS_Host.cpp:24
Point evaluation of the curl of a vector (results in a vector)
Lightweight Evaluator Helper This class is a lightweight wrapper for extracting and applying all rele...
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)
Point evaluation of the laplacian of a scalar (could be on a manifold or not)
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_INLINE_FUNCTION double curlTestSolution(double x, double y, double z, int component, int dimension)
KOKKOS_INLINE_FUNCTION double divergenceTestSamples(double x, double y, double z, int component, int dimension)
double applyAlphasToDataSingleComponentSingleTargetSite(view_type_data sampling_input_data, const int column_of_input, TargetOperation lro, const int target_index, const int evaluation_site_local_index, const int output_component_axis_1, const int output_component_axis_2, const int input_component_axis_1, const int input_component_axis_2, bool scalar_as_vector_if_needed=true) const
Dot product of alphas with sampling data, FOR A SINGLE target_index, where sampling data is in a 1D/2...
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)