Compadre 1.5.5
Loading...
Searching...
No Matches
GMLS_NeumannGradScalar.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>
13
14#include "GMLS_Tutorial.hpp"
16
17#ifdef COMPADRE_USE_MPI
18#include <mpi.h>
19#endif
20
21#include <Kokkos_Timer.hpp>
22#include <Kokkos_Core.hpp>
23
24using namespace Compadre;
25
26//! [Parse Command Line Arguments]
27
28// called from command line
29int main (int argc, char* args[]) {
30
31// initializes MPI (if available) with command line arguments given
32#ifdef COMPADRE_USE_MPI
33MPI_Init(&argc, &args);
34#endif
35
36// initializes Kokkos with command line arguments given
37Kokkos::initialize(argc, args);
38
39// becomes false if the computed solution not within the failure_threshold of the actual solution
40bool all_passed = true;
41
42// code block to reduce scope for all Kokkos View allocations
43// otherwise, Views may be deallocating when we call Kokkos finalize() later
44{
45
46 CommandLineProcessor clp(argc, args);
47 auto order = clp.order;
48 auto dimension = clp.dimension;
49 auto number_target_coords = clp.number_target_coords;
50 auto constraint_name = clp.constraint_name;
51 auto solver_name = clp.solver_name;
52 auto problem_name = clp.problem_name;
53 auto number_of_batches = clp.number_of_batches;
54
55 // the functions we will be seeking to reconstruct are in the span of the basis
56 // of the reconstruction space we choose for GMLS, so the error should be very small
57 const double failure_tolerance = 1e-9;
58
59 // minimum neighbors for unisolvency is the same as the size of the polynomial basis
60 const int min_neighbors = Compadre::GMLS::getNP(order, dimension);
61
62 //! [Parse Command Line Arguments]
63 Kokkos::Timer timer;
64 Kokkos::Profiling::pushRegion("Setup Point Data");
65 //! [Setting Up The Point Cloud]
66
67 // approximate spacing of source sites
68 double h_spacing = 0.05;
69 int n_neg1_to_1 = 2*(1/h_spacing) + 1; // always odd
70
71 // number of source coordinate sites that will fill a box of [-1,1]x[-1,1]x[-1,1] with a spacing approximately h
72 const int number_source_coords = std::pow(n_neg1_to_1, dimension);
73
74 // coordinates of source sites
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);
78
79 // coordinates of target sites
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);
82
83 // tangent bundle for each target sites
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);
86
87 // fill source coordinates with a uniform grid
88 int source_index = 0;
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;
96 if (dimension==3) {
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];
100 source_index++;
101 }
102 }
103 if (dimension==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;
107 source_index++;
108 }
109 }
110 if (dimension==1) {
111 source_coords(source_index,0) = this_coord[0];
112 source_coords(source_index,1) = 0;
113 source_coords(source_index,2) = 0;
114 source_index++;
115 }
116 }
117
118 // fill target coords somewhere inside of [-0.5,0.5]x[-0.5,0.5]x[-0.5,0.5]
119 for(int i=0; i<number_target_coords; i++){
120
121 // first, we get a uniformly random distributed direction
122 double rand_dir[3] = {0,0,0};
123
124 for (int j=0; j<dimension; ++j) {
125 // rand_dir[j] is in [-0.5, 0.5]
126 rand_dir[j] = ((double)rand() / (double) RAND_MAX) - 0.5;
127 }
128
129 // then we get a uniformly random radius
130 for (int j=0; j<dimension; ++j) {
131 target_coords(i,j) = rand_dir[j];
132 }
133 // target_coords(i, 2) = 1.0;
134
135 // Set tangent bundles
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));
151 }
152 }
153
154 //! [Setting Up The Point Cloud]
155
156 Kokkos::Profiling::popRegion();
157 Kokkos::Profiling::pushRegion("Creating Data");
158
159 //! [Creating The Data]
160
161
162 // source coordinates need copied to device before using to construct sampling data
163 Kokkos::deep_copy(source_coords_device, source_coords);
164
165 // target coordinates copied next, because it is a convenient time to send them to device
166 Kokkos::deep_copy(target_coords_device, target_coords);
167
168 // tangent bundles copied next, because it is a convenient time to send them to device
169 Kokkos::deep_copy(tangent_bundles_device, tangent_bundles);
170
171 // need Kokkos View storing true solution
172 Kokkos::View<double*, Kokkos::DefaultExecutionSpace> sampling_data_device("samples of true solution",
173 source_coords_device.extent(0));
174
175 Kokkos::parallel_for("Sampling Manufactured Solutions", Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>
176 (0,source_coords.extent(0)), KOKKOS_LAMBDA(const int i) {
177
178 // coordinates of source site 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;
182
183 // data for targets with scalar input
184 sampling_data_device(i) = trueSolution(xval, yval, zval, order, dimension);
185 });
186
187 //! [Creating The Data]
188
189 Kokkos::Profiling::popRegion();
190 Kokkos::Profiling::pushRegion("Neighbor Search");
191
192 //! [Performing Neighbor Search]
193
194
195 // Point cloud construction for neighbor search
196 // CreatePointCloudSearch constructs an object of type PointCloudSearch, but deduces the templates for you
197 auto point_cloud_search(CreatePointCloudSearch(source_coords, dimension));
198
199 // each row is a neighbor list for a target site, with the first column of each row containing
200 // the number of neighbors for that rows corresponding target site
201 double epsilon_multiplier = 1.8;
202 int estimated_upper_bound_number_neighbors =
203 point_cloud_search.getEstimatedNumberNeighborsUpperBound(min_neighbors, dimension, epsilon_multiplier);
204
205 Kokkos::View<int**, Kokkos::DefaultExecutionSpace> neighbor_lists_device("neighbor lists",
206 number_target_coords, estimated_upper_bound_number_neighbors); // first column is # of neighbors
207 Kokkos::View<int**>::HostMirror neighbor_lists = Kokkos::create_mirror_view(neighbor_lists_device);
208
209 // each target site has a window size
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);
212
213 // query the point cloud to generate the neighbor lists using a kdtree to produce the n nearest neighbor
214 // to each target site, adding (epsilon_multiplier-1)*100% to whatever the distance away the further neighbor used is from
215 // each target to the view for epsilon
216 point_cloud_search.generate2DNeighborListsFromKNNSearch(false /*not dry run*/, target_coords, neighbor_lists,
217 epsilon, min_neighbors, epsilon_multiplier);
218
219 //! [Performing Neighbor Search]
220
221 Kokkos::Profiling::popRegion();
222 Kokkos::fence(); // let call to build neighbor lists complete before copying back to device
223 timer.reset();
224
225 //! [Setting Up The GMLS Object]
226
227
228 // Copy data back to device (they were filled on the host)
229 // We could have filled Kokkos Views with memory space on the host
230 // and used these instead, and then the copying of data to the device
231 // would be performed in the GMLS class
232 Kokkos::deep_copy(neighbor_lists_device, neighbor_lists);
233 Kokkos::deep_copy(epsilon_device, epsilon);
234
235 // initialize an instance of the GMLS class
238 order, dimension,
239 solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
240 0 /*manifold order*/);
241
242 // pass in neighbor lists, source coordinates, target coordinates, and window sizes
243 //
244 // neighbor lists have the format:
245 // dimensions: (# number of target sites) X (# maximum number of neighbors for any given target + 1)
246 // the first column contains the number of neighbors for that rows corresponding target index
247 //
248 // source coordinates have the format:
249 // dimensions: (# number of source sites) X (dimension)
250 // entries in the neighbor lists (integers) correspond to rows of this 2D array
251 //
252 // target coordinates have the format:
253 // dimensions: (# number of target sites) X (dimension)
254 // # of target sites is same as # of rows of neighbor lists
255 //
256 my_GMLS.setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
257 my_GMLS.setTangentBundle(tangent_bundles_device);
258
259 // create a vector of target operations
260 TargetOperation lro;
262
263 // and then pass them to the GMLS class
264 my_GMLS.addTargets(lro);
265
266 // sets the weighting kernel function from WeightingFunctionType
267 my_GMLS.setWeightingType(WeightingFunctionType::Power);
268
269 // power to use in that weighting kernel function
270 my_GMLS.setWeightingParameter(2);
271
272 // generate the alphas that to be combined with data for each target operation requested in lro
273 my_GMLS.generateAlphas(number_of_batches);
274
275 //! [Setting Up The GMLS Object]
276
277 double instantiation_time = timer.seconds();
278 std::cout << "Took " << instantiation_time << "s to complete normal vectors generation." << std::endl;
279 Kokkos::fence(); // let generateNormalVectors finish up before using alphas
280 Kokkos::Profiling::pushRegion("Apply Alphas to Data");
281
282 //! [Apply GMLS Alphas To Data]
283
284 // it is important to note that if you expect to use the data as a 1D view, then you should use double*
285 // however, if you know that the target operation will result in a 2D view (vector or matrix output),
286 // then you should template with double** as this is something that can not be infered from the input data
287 // or the target operator at compile time. Additionally, a template argument is required indicating either
288 // Kokkos::HostSpace or Kokkos::DefaultExecutionSpace::memory_space()
289
290 // The Evaluator class takes care of handling input data views as well as the output data views.
291 // It uses information from the GMLS class to determine how many components are in the input
292 // as well as output for any choice of target functionals and then performs the contactions
293 // on the data using the alpha coefficients generated by the GMLS class, all on the device.
294 Evaluator gmls_evaluator(&my_GMLS);
295
296 auto output_value = gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double*, Kokkos::HostSpace>
297 (sampling_data_device, LaplacianOfScalarPointEvaluation);
298
299 Kokkos::fence(); // let application of alphas to data finish before using results
300 Kokkos::Profiling::popRegion();
301 // times the Comparison in Kokkos
302 Kokkos::Profiling::pushRegion("Comparison");
303
304 //! [Check That Solutions Are Correct]
305
306 // loop through the target sites
307 for (int i=0; i<number_target_coords; i++) {
308 // target site i's coordinate
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;
312
313 // 0th entry is # of neighbors, which is the index beyond the last neighbor
314 int num_neigh_i = neighbor_lists(i, 0);
315 double b_i = my_GMLS.getSolutionSetHost()->getAlpha0TensorTo0Tensor(lro, i, num_neigh_i);
316
317 // load value from output
318 double GMLS_value = output_value(i);
319
320 // obtain the real Laplacian
321 double actual_Laplacian = trueLaplacian(xval, yval, zval, order, dimension);
322
323 // calculate value of g to reconstruct the computed Laplacian
324 double actual_Gradient[3] = {0,0,0}; // initialized for 3, but only filled up to dimension
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;
329
330 // check actual function value
331 if(GMLS_value!=GMLS_value || std::abs(actual_Laplacian - adjusted_value) > failure_tolerance) {
332 all_passed = false;
333 std::cout << i << " Failed Actual by: " << std::abs(actual_Laplacian - adjusted_value) << std::endl;
334 }
335 }
336
337 //! [Check That Solutions Are Correct]
338 // popRegion hidden from tutorial
339 // stop timing comparison loop
340 Kokkos::Profiling::popRegion();
341 //! [Finalize Program]
342} // end of code block to reduce scope, causing Kokkos View de-allocations
343// otherwise, Views may be deallocating when we call Kokkos finalize() later
344
345// finalize Kokkos and MPI (if available)
346Kokkos::finalize();
347#ifdef COMPADRE_USE_MPI
348MPI_Finalize();
349#endif
350
351// output to user that test passed or failed
352if(all_passed) {
353 fprintf(stdout, "Passed test \n");
354 return 0;
355} else {
356 fprintf(stdout, "Failed test \n");
357 return -1;
358}
359
360} // main
361
362
363//! [Finalize Program]
int main(int argc, char *args[])
[Parse Command Line Arguments]
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 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...
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)
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 setTangentBundle(view_type tangent_directions)
(OPTIONAL) Sets orthonormal tangent directions for reconstruction on a manifold.
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...
decltype(_h_ss) * getSolutionSetHost(bool alpha_validity_check=true)
Get solution set on host.
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)
void setWeightingType(const std::string &wt)
Type for weighting kernel for GMLS problem.
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....
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...
constexpr SamplingFunctional PointSample
Available sampling functionals.
TargetOperation
Available target functionals.
@ LaplacianOfScalarPointEvaluation
Point evaluation of the laplacian of a scalar (could be on a manifold or not)
@ ScalarTaylorPolynomial
Scalar polynomial basis centered at the target site and scaled by sum of basis powers e....