Compadre 1.5.5
Loading...
Searching...
No Matches
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>
12#include "GMLS_Tutorial.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
22using namespace Compadre;
23
24int main (int argc, char* args[])
25{
26
27#ifdef COMPADRE_USE_MPI
28 MPI_Init(&argc, &args);
29#endif
30
31bool 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
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
378if(all_passed) {
379 fprintf(stdout, "Passed test \n");
380 return 0;
381} else {
382 fprintf(stdout, "Failed test \n");
383 return -1;
384}
385}
int main(int argc, char *args[])
Manifold GMLS Example.
Definition GMLS_Host.cpp:24
KOKKOS_INLINE_FUNCTION double trueSolution(double x, double y, double z, int order, int dimension)
KOKKOS_INLINE_FUNCTION void trueGradient(double *ans, double x, double y, double z, int order, int dimension)
KOKKOS_INLINE_FUNCTION double divergenceTestSamples(double x, double y, double z, int component, int dimension)
KOKKOS_INLINE_FUNCTION double curlTestSolution(double x, double y, double z, int component, int dimension)
KOKKOS_INLINE_FUNCTION double trueLaplacian(double x, double y, double z, int order, int dimension)
Lightweight Evaluator Helper This class is a lightweight wrapper for extracting and applying all rele...
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...
Generalized Moving Least Squares (GMLS)
void addTargets(TargetOperation lro)
Adds a target to the vector of target functional to be applied to the reconstruction.
void setWeightingParameter(int wp, int index=0)
Parameter for weighting kernel for GMLS problem index = 0 sets p paramater for weighting kernel index...
void generateAlphas(const int number_of_batches=1, const bool keep_coefficients=false, const bool clear_cache=true)
Meant to calculate target operations and apply the evaluations to the previously constructed polynomi...
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)
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....
@ LaplacianOfScalarPointEvaluation
Point evaluation of the laplacian of a scalar (could be on a manifold or not)
@ GradientOfScalarPointEvaluation
Point evaluation of the gradient of a scalar.
@ CurlOfVectorPointEvaluation
Point evaluation of the curl of a vector (results in a vector)
@ DivergenceOfVectorPointEvaluation
Point evaluation of the divergence of a vector (results in a scalar)
@ ScalarPointEvaluation
Point evaluation of a scalar.