Compadre  1.5.5
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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.setWeightingParameter(10);
134 
135  std::vector<TargetOperation> lro(3);
136  lro[0] = ScalarPointEvaluation;
139  my_GMLS.addTargets(lro);
140  // add two more targets individually to test addTargets(...) function
141  my_GMLS.addTargets(DivergenceOfVectorPointEvaluation);
142  my_GMLS.addTargets(CurlOfVectorPointEvaluation);
143  my_GMLS.generateAlphas();
144 
145  double instantiation_time = timer.seconds();
146  std::cout << "Took " << instantiation_time << "s to complete instantiation." << std::endl;
147 
148  Kokkos::Profiling::pushRegion("Creating Data");
149 
150 
151  // need Kokkos View storing true solution
152  Kokkos::View<double*, Kokkos::HostSpace> sampling_data("samples of true solution", source_coords.extent(0));
153  Kokkos::View<double**, Kokkos::HostSpace> gradient_sampling_data("samples of true gradient", source_coords.extent(0), dimension);
154  Kokkos::View<double**, Kokkos::LayoutLeft, Kokkos::HostSpace> divergence_sampling_data("samples of true solution for divergence test", source_coords.extent(0), dimension);
155  Kokkos::parallel_for("Sampling Manufactured Solutions", Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>(0,source_coords.extent(0)), KOKKOS_LAMBDA(const int i) {
156  double xval = source_coords(i,0);
157  double yval = (dimension>1) ? source_coords(i,1) : 0;
158  double zval = (dimension>2) ? source_coords(i,2) : 0;
159  sampling_data(i) = trueSolution(xval, yval, zval, order, dimension);
160  double true_grad[3] = {0,0,0};
161  trueGradient(true_grad, xval, yval,zval, order, dimension);
162  for (int j=0; j<dimension; ++j) {
163  divergence_sampling_data(i,j) = divergenceTestSamples(xval, yval, zval, j, dimension);
164  gradient_sampling_data(i,j) = true_grad[j];
165  }
166  });
167  Kokkos::Profiling::popRegion();
168 
169  Evaluator gmls_evaluator(&my_GMLS);
170 
171  for (int i=0; i<number_target_coords; i++) {
172 
173  Kokkos::Profiling::pushRegion("Apply Alphas to Data");
174 
175  double GMLS_value = gmls_evaluator.applyAlphasToDataSingleComponentSingleTargetSite(sampling_data, 0, ScalarPointEvaluation, i, 0, 0, 0, 0, 0);
176  //for (int j = 0; j< neighbor_lists(i,0); j++){
177  // double xval = source_coords(neighbor_lists(i,j+1),0);
178  // double yval = (dimension>1) ? source_coords(neighbor_lists(i,j+1),1) : 0;
179  // double zval = (dimension>2) ? source_coords(neighbor_lists(i,j+1),2) : 0;
180  // GMLS_value += gmls_evaluator.getAlpha0TensorTo0Tensor(ScalarPointEvaluation, i, j)*trueSolution(xval, yval, zval, order, dimension);
181  //}
182 
183  double GMLS_Laplacian = gmls_evaluator.applyAlphasToDataSingleComponentSingleTargetSite(sampling_data, 0, LaplacianOfScalarPointEvaluation, i, 0, 0, 0, 0, 0);
184  //double GMLS_Laplacian = 0.0;
185  //for (int j = 0; j< neighbor_lists(i,0); j++){
186  // double xval = source_coords(neighbor_lists(i,j+1),0);
187  // double yval = (dimension>1) ? source_coords(neighbor_lists(i,j+1),1) : 0;
188  // double zval = (dimension>2) ? source_coords(neighbor_lists(i,j+1),2) : 0;
189  // GMLS_Laplacian += gmls_evaluator.getAlpha0TensorTo0Tensor(LaplacianOfScalarPointEvaluation, i, j)*trueSolution(xval, yval, zval, order, dimension);
190  //}
191 
192  double GMLS_GradX = gmls_evaluator.applyAlphasToDataSingleComponentSingleTargetSite(sampling_data, 0, GradientOfScalarPointEvaluation, i, 0, 0, 0, 0, 0);
193  //double GMLS_GradX = 0.0;
194  //for (int j = 0; j< neighbor_lists(i,0); j++){
195  // double xval = source_coords(neighbor_lists(i,j+1),0);
196  // double yval = (dimension>1) ? source_coords(neighbor_lists(i,j+1),1) : 0;
197  // double zval = (dimension>2) ? source_coords(neighbor_lists(i,j+1),2) : 0;
198  // GMLS_GradX += gmls_evaluator.getAlpha0TensorTo1Tensor(GradientOfScalarPointEvaluation, i, 0, j)*trueSolution(xval, yval, zval, order, dimension);
199  //}
200 
201  double GMLS_GradY = (dimension>1) ? gmls_evaluator.applyAlphasToDataSingleComponentSingleTargetSite(sampling_data, 0, GradientOfScalarPointEvaluation, i, 0, 1, 0, 0, 0) : 0;
202  //double GMLS_GradY = 0.0;
203  //if (dimension>1) {
204  // for (int j = 0; j< neighbor_lists(i,0); j++){
205  // double xval = source_coords(neighbor_lists(i,j+1),0);
206  // double yval = source_coords(neighbor_lists(i,j+1),1);
207  // double zval = (dimension>2) ? source_coords(neighbor_lists(i,j+1),2) : 0;
208  // GMLS_GradY += gmls_evaluator.getAlpha0TensorTo1Tensor(GradientOfScalarPointEvaluation, i, 1, j)*trueSolution(xval, yval, zval, order, dimension);
209  // }
210  //}
211 
212  double GMLS_GradZ = (dimension>2) ? gmls_evaluator.applyAlphasToDataSingleComponentSingleTargetSite(sampling_data, 0, GradientOfScalarPointEvaluation, i, 0, 2, 0, 0, 0) : 0;
213  //double GMLS_GradZ = 0.0;
214  //if (dimension>2) {
215  // for (int j = 0; j< neighbor_lists(i,0); j++){
216  // double xval = source_coords(neighbor_lists(i,j+1),0);
217  // double yval = source_coords(neighbor_lists(i,j+1),1);
218  // double zval = source_coords(neighbor_lists(i,j+1),2);
219  // GMLS_GradZ += gmls_evaluator.getAlpha0TensorTo1Tensor(GradientOfScalarPointEvaluation, i, 2, j)*trueSolution(xval, yval, zval, order, dimension);
220  // }
221  //}
222 
223  double GMLS_Divergence = gmls_evaluator.applyAlphasToDataSingleComponentSingleTargetSite(gradient_sampling_data, 0, DivergenceOfVectorPointEvaluation, i, 0, 0, 0, 0, 0);
224  if (dimension>1) GMLS_Divergence += gmls_evaluator.applyAlphasToDataSingleComponentSingleTargetSite(gradient_sampling_data, 1, DivergenceOfVectorPointEvaluation, i, 0, 0, 0, 1, 0);
225  if (dimension>2) GMLS_Divergence += gmls_evaluator.applyAlphasToDataSingleComponentSingleTargetSite(gradient_sampling_data, 2, DivergenceOfVectorPointEvaluation, i, 0, 0, 0, 2, 0);
226 
227  //double GMLS_Divergence = gmls_evaluator.applyAlphasToDataSingleComponentSingleTargetSite(Kokkos::subview(gradient_sampling_data,Kokkos::ALL,0), 0, DivergenceOfVectorPointEvaluation, i, 0, 0, 0, 0);
228  //if (dimension>1) GMLS_Divergence += gmls_evaluator.applyAlphasToDataSingleComponentSingleTargetSite(Kokkos::subview(gradient_sampling_data,Kokkos::ALL,1), 1, DivergenceOfVectorPointEvaluation, i, 0, 0, 1, 0);
229  //if (dimension>2) GMLS_Divergence += gmls_evaluator.applyAlphasToDataSingleComponentSingleTargetSite(Kokkos::subview(gradient_sampling_data,Kokkos::ALL,2), 2, DivergenceOfVectorPointEvaluation, i, 0, 0, 2, 0);
230  //double GMLS_Divergence = 0.0;
231  //for (int j = 0; j< neighbor_lists(i,0); j++){
232  // double xval = source_coords(neighbor_lists(i,j+1),0);
233  // double yval = (dimension>1) ? source_coords(neighbor_lists(i,j+1),1) : 0;
234  // double zval = (dimension>2) ? source_coords(neighbor_lists(i,j+1),2) : 0;
235  // // TODO: use different functions for the vector components
236  // if (use_arbitrary_order_divergence) {
237  // GMLS_Divergence += gmls_evaluator.getAlpha1TensorTo0Tensor(DivergenceOfVectorPointEvaluation, i, j, 0)*trueSolution(xval, yval, zval, order, dimension);
238  // if (dimension>1) GMLS_Divergence += gmls_evaluator.getAlpha1TensorTo0Tensor(DivergenceOfVectorPointEvaluation, i, j, 1)*trueSolution(xval, yval, zval, order, dimension);
239  // if (dimension>2) GMLS_Divergence += gmls_evaluator.getAlpha1TensorTo0Tensor(DivergenceOfVectorPointEvaluation, i, j, 2)*trueSolution(xval, yval, zval, order, dimension);
240  // } else {
241  // for (int k=0; k<dimension; ++k) {
242  // GMLS_Divergence += gmls_evaluator.getAlpha1TensorTo0Tensor(DivergenceOfVectorPointEvaluation, i, j, k)*divergenceTestSamples(xval, yval, zval, k, dimension);
243  // }
244  // }
245  //}
246 
247  double GMLS_CurlX = 0.0;
248  double GMLS_CurlY = 0.0;
249  double GMLS_CurlZ = 0.0;
250  if (dimension>1) {
251  for (int j=0; j<dimension; ++j) {
252  GMLS_CurlX += gmls_evaluator.applyAlphasToDataSingleComponentSingleTargetSite(divergence_sampling_data, j, CurlOfVectorPointEvaluation, i, 0, 0, 0, j, 0);
253  GMLS_CurlY += gmls_evaluator.applyAlphasToDataSingleComponentSingleTargetSite(divergence_sampling_data, j, CurlOfVectorPointEvaluation, i, 0, 1, 0, j, 0);
254  }
255  }
256 
257  if (dimension>2) {
258  for (int j=0; j<dimension; ++j) {
259  GMLS_CurlZ += gmls_evaluator.applyAlphasToDataSingleComponentSingleTargetSite(divergence_sampling_data, j, CurlOfVectorPointEvaluation, i, 0, 2, 0, j, 0);
260  }
261 
262  }
263 
264  Kokkos::Profiling::popRegion();
265  //if (dimension>1) {
266  // for (int j = 0; j< neighbor_lists(i,0); j++){
267  // double xval = source_coords(neighbor_lists(i,j+1),0);
268  // double yval = source_coords(neighbor_lists(i,j+1),1);
269  // double zval = (dimension>2) ? source_coords(neighbor_lists(i,j+1),2) : 0;
270  // for (int k=0; k<dimension; ++k) {
271  // GMLS_CurlX += my_GMLS.getAlpha1TensorTo1Tensor(CurlOfVectorPointEvaluation, i, 0, j, k)*divergenceTestSamples(xval, yval, zval, k, dimension);
272  // }
273  // }
274 
275  // for (int j = 0; j< neighbor_lists(i,0); j++){
276  // double xval = source_coords(neighbor_lists(i,j+1),0);
277  // double yval = source_coords(neighbor_lists(i,j+1),1);
278  // double zval = (dimension>2) ? source_coords(neighbor_lists(i,j+1),2) : 0;
279  // for (int k=0; k<dimension; ++k) {
280  // GMLS_CurlY += my_GMLS.getAlpha1TensorTo1Tensor(CurlOfVectorPointEvaluation, i, 1, j, k)*divergenceTestSamples(xval, yval, zval, k, dimension);
281  // }
282  // }
283  //}
284 
285  //if (dimension>2) {
286  // for (int j = 0; j< neighbor_lists(i,0); j++){
287  // double xval = source_coords(neighbor_lists(i,j+1),0);
288  // double yval = source_coords(neighbor_lists(i,j+1),1);
289  // double zval = source_coords(neighbor_lists(i,j+1),2);
290  // for (int k=0; k<dimension; ++k) {
291  // GMLS_CurlZ += my_GMLS.getAlpha1TensorTo1Tensor(CurlOfVectorPointEvaluation, i, 2, j, k)*divergenceTestSamples(xval, yval, zval, k, dimension);
292  // }
293  // }
294  //}
295  //
296  Kokkos::Profiling::pushRegion("Comparison");
297 
298  double xval = target_coords(i,0);
299  double yval = (dimension>1) ? target_coords(i,1) : 0;
300  double zval = (dimension>2) ? target_coords(i,2) : 0;
301 
302  double actual_value = trueSolution(xval, yval, zval, order, dimension);
303  double actual_Laplacian = trueLaplacian(xval, yval, zval, order, dimension);
304  double actual_Gradient[3] = {0,0,0};
305  trueGradient(actual_Gradient, xval, yval, zval, order, dimension);
306  double actual_Divergence;
307  actual_Divergence = trueLaplacian(xval, yval, zval, order, dimension);
308 
309  double actual_CurlX = 0;
310  double actual_CurlY = 0;
311  double actual_CurlZ = 0;
312  if (dimension>1) {
313  actual_CurlX = curlTestSolution(xval, yval, zval, 0, dimension);
314  actual_CurlY = curlTestSolution(xval, yval, zval, 1, dimension);
315  }
316  if (dimension>2) {
317  actual_CurlZ = curlTestSolution(xval, yval, zval, 2, dimension);
318  }
319 
320 // fprintf(stdout, "Reconstructed value: %f \n", GMLS_value);
321 // fprintf(stdout, "Actual value: %f \n", actual_value);
322 // fprintf(stdout, "Reconstructed Laplacian: %f \n", GMLS_Laplacian);
323 // fprintf(stdout, "Actual Laplacian: %f \n", actual_Laplacian);
324 
325  if(GMLS_value!=GMLS_value || std::abs(actual_value - GMLS_value) > failure_tolerance) {
326  all_passed = false;
327  std::cout << "Failed Actual by: " << std::abs(actual_value - GMLS_value) << std::endl;
328  }
329 
330  if(std::abs(actual_Laplacian - GMLS_Laplacian) > failure_tolerance) {
331  all_passed = false;
332  std::cout << "Failed Laplacian by: " << std::abs(actual_Laplacian - GMLS_Laplacian) << std::endl;
333  }
334 
335  if(std::abs(actual_Gradient[0] - GMLS_GradX) > failure_tolerance) {
336  all_passed = false;
337  std::cout << "Failed GradX by: " << std::abs(actual_Gradient[0] - GMLS_GradX) << std::endl;
338  }
339 
340  if (dimension>1) {
341  if(std::abs(actual_Gradient[1] - GMLS_GradY) > failure_tolerance) {
342  all_passed = false;
343  std::cout << "Failed GradY by: " << std::abs(actual_Gradient[1] - GMLS_GradY) << std::endl;
344  }
345  }
346 
347  if (dimension>2) {
348  if(std::abs(actual_Gradient[2] - GMLS_GradZ) > failure_tolerance) {
349  all_passed = false;
350  std::cout << "Failed GradZ by: " << std::abs(actual_Gradient[2] - GMLS_GradZ) << std::endl;
351  }
352  }
353 
354  if(std::abs(actual_Divergence - GMLS_Divergence) > failure_tolerance) {
355  all_passed = false;
356  std::cout << "Failed Divergence by: " << std::abs(actual_Divergence - GMLS_Divergence) << std::endl;
357  }
358 
359  double tmp_diff = 0;
360  if (dimension>1)
361  tmp_diff += std::abs(actual_CurlX - GMLS_CurlX) + std::abs(actual_CurlY - GMLS_CurlY);
362  if (dimension>2)
363  tmp_diff += std::abs(actual_CurlZ - GMLS_CurlZ);
364  if(std::abs(tmp_diff) > failure_tolerance) {
365  all_passed = false;
366  std::cout << "Failed Curl by: " << std::abs(tmp_diff) << std::endl;
367  }
368  Kokkos::Profiling::popRegion();
369  }
370 
371 }
372 
373  Kokkos::finalize();
374 #ifdef COMPADRE_USE_MPI
375  MPI_Finalize();
376 #endif
377 
378 if(all_passed) {
379  fprintf(stdout, "Passed test \n");
380  return 0;
381 } else {
382  fprintf(stdout, "Failed test \n");
383  return -1;
384 }
385 }
Point evaluation of a scalar.
Lightweight Evaluator Helper This class is a lightweight wrapper for extracting and applying all rele...
Point evaluation of the gradient of a scalar.
Point evaluation of the laplacian of a scalar (could be on a manifold or not)
int main(int argc, char **argv)
Point evaluation of the curl of a vector (results in a vector)
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...
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 divergence of a vector (results in 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)
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)