17 #include <Compadre_Config.h>
23 #ifdef COMPADRE_USE_MPI
27 #include <Kokkos_Timer.hpp>
28 #include <Kokkos_Core.hpp>
30 using namespace Compadre;
32 int main (
int argc,
char* args[])
35 #ifdef COMPADRE_USE_MPI
36 MPI_Init(&argc, &args);
39 bool all_passed =
true;
44 auto order = clp.
order;
51 const double failure_tolerance = 1e-9;
53 const int offset = 15;
57 std::cout << min_neighbors <<
" " << max_neighbors << std::endl;
58 std::uniform_int_distribution<int> gen_num_neighbors(min_neighbors, max_neighbors);
61 Kokkos::initialize(argc, args);
63 Kokkos::Profiling::pushRegion(
"Setup");
67 std::uniform_int_distribution<int> gen_neighbor_number(offset, N);
70 Kokkos::View<int**, Kokkos::HostSpace> neighbor_lists(
"neighbor lists", number_target_coords, max_neighbors+1);
71 Kokkos::View<double**, Kokkos::HostSpace> source_coords(
"neighbor coordinates", N, dimension);
72 Kokkos::View<double*, Kokkos::HostSpace> epsilon(
"h supports", number_target_coords);
74 for (
int i=0; i<number_target_coords; i++) {
79 for(
int i = 0; i < offset; i++){
80 for(
int j = 0; j < dimension; j++){
81 source_coords(i,j) = 0.1;
86 for(
int i = offset; i < N; i++){
87 double randx = (2.0*(double)rand() / (double) RAND_MAX - 1.0)*epsilon(0)/2.0;
88 double randy = (2.0*(double)rand() / (double) RAND_MAX - 1.0)*epsilon(0)/2.0;
89 double randz = (2.0*(double)rand() / (double) RAND_MAX - 1.0)*epsilon(0)/2.0;
90 source_coords(i,0) = randx;
91 if (dimension>1) source_coords(i,1) = randy;
92 if (dimension>2) source_coords(i,2) = randz;
95 const double target_epsilon = 0.1;
97 Kokkos::View<double**, Kokkos::HostSpace> target_coords (
"target coordinates", number_target_coords, dimension);
98 for(
int i = 0; i < number_target_coords; i++){
99 double randx = (2.0*(double)rand() / (double) RAND_MAX - 1.0)*target_epsilon/2.0;
100 double randy = (2.0*(double)rand() / (double) RAND_MAX - 1.0)*target_epsilon/2.0;
101 double randz = (2.0*(double)rand() / (double) RAND_MAX - 1.0)*target_epsilon/2.0;
102 target_coords(i,0) = randx;
103 if (dimension>1) target_coords(i,1) = randy;
104 if (dimension>2) target_coords(i,2) = randz;
108 for (
int i=0; i<number_target_coords; i++) {
111 int r = max_neighbors;
112 neighbor_lists(i,0) = r;
114 for(
int j=0; j<r; j++){
115 neighbor_lists(i,j+1) = offset + j + 1;
134 Kokkos::Profiling::popRegion();
137 GMLS my_GMLS(order, dimension,
138 solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
140 my_GMLS.
setProblemData(neighbor_lists, source_coords, target_coords, epsilon);
141 my_GMLS.setWeightingParameter(10);
143 std::vector<TargetOperation> lro(3);
147 my_GMLS.addTargets(lro);
151 my_GMLS.generateAlphas();
153 double instantiation_time = timer.seconds();
154 std::cout <<
"Took " << instantiation_time <<
"s to complete instantiation." << std::endl;
156 Kokkos::Profiling::pushRegion(
"Creating Data");
160 Kokkos::View<double*, Kokkos::HostSpace> sampling_data(
"samples of true solution", source_coords.extent(0));
161 Kokkos::View<double**, Kokkos::HostSpace> gradient_sampling_data(
"samples of true gradient", source_coords.extent(0), dimension);
162 Kokkos::View<double**, Kokkos::LayoutLeft, Kokkos::HostSpace> divergence_sampling_data(
"samples of true solution for divergence test", source_coords.extent(0), dimension);
163 Kokkos::parallel_for(
"Sampling Manufactured Solutions", Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>(0,source_coords.extent(0)), KOKKOS_LAMBDA(
const int i) {
164 double xval = source_coords(i,0);
165 double yval = (dimension>1) ? source_coords(i,1) : 0;
166 double zval = (dimension>2) ? source_coords(i,2) : 0;
167 sampling_data(i) =
trueSolution(xval, yval, zval, order, dimension);
168 double true_grad[3] = {0,0,0};
169 trueGradient(true_grad, xval, yval,zval, order, dimension);
170 for (
int j=0; j<dimension; ++j) {
172 gradient_sampling_data(i,j) = true_grad[j];
175 Kokkos::Profiling::popRegion();
179 for (
int i=0; i<number_target_coords; i++) {
181 Kokkos::Profiling::pushRegion(
"Apply Alphas to Data");
255 double GMLS_CurlX = 0.0;
256 double GMLS_CurlY = 0.0;
257 double GMLS_CurlZ = 0.0;
259 for (
int j=0; j<dimension; ++j) {
266 for (
int j=0; j<dimension; ++j) {
272 Kokkos::Profiling::popRegion();
304 Kokkos::Profiling::pushRegion(
"Comparison");
306 double xval = target_coords(i,0);
307 double yval = (dimension>1) ? target_coords(i,1) : 0;
308 double zval = (dimension>2) ? target_coords(i,2) : 0;
310 double actual_value =
trueSolution(xval, yval, zval, order, dimension);
311 double actual_Laplacian =
trueLaplacian(xval, yval, zval, order, dimension);
312 double actual_Gradient[3] = {0,0,0};
313 trueGradient(actual_Gradient, xval, yval, zval, order, dimension);
314 double actual_Divergence;
315 actual_Divergence =
trueLaplacian(xval, yval, zval, order, dimension);
317 double actual_CurlX = 0;
318 double actual_CurlY = 0;
319 double actual_CurlZ = 0;
333 if(GMLS_value!=GMLS_value || std::abs(actual_value - GMLS_value) > failure_tolerance) {
335 std::cout <<
"Failed Actual by: " << std::abs(actual_value - GMLS_value) << std::endl;
338 if(std::abs(actual_Laplacian - GMLS_Laplacian) > failure_tolerance) {
340 std::cout <<
"Failed Laplacian by: " << std::abs(actual_Laplacian - GMLS_Laplacian) << std::endl;
343 if(std::abs(actual_Gradient[0] - GMLS_GradX) > failure_tolerance) {
345 std::cout <<
"Failed GradX by: " << std::abs(actual_Gradient[0] - GMLS_GradX) << std::endl;
349 if(std::abs(actual_Gradient[1] - GMLS_GradY) > failure_tolerance) {
351 std::cout <<
"Failed GradY by: " << std::abs(actual_Gradient[1] - GMLS_GradY) << std::endl;
356 if(std::abs(actual_Gradient[2] - GMLS_GradZ) > failure_tolerance) {
358 std::cout <<
"Failed GradZ by: " << std::abs(actual_Gradient[2] - GMLS_GradZ) << std::endl;
362 if(std::abs(actual_Divergence - GMLS_Divergence) > failure_tolerance) {
364 std::cout <<
"Failed Divergence by: " << std::abs(actual_Divergence - GMLS_Divergence) << std::endl;
369 tmp_diff += std::abs(actual_CurlX - GMLS_CurlX) + std::abs(actual_CurlY - GMLS_CurlY);
371 tmp_diff += std::abs(actual_CurlZ - GMLS_CurlZ);
372 if(std::abs(tmp_diff) > failure_tolerance) {
374 std::cout <<
"Failed Curl by: " << std::abs(tmp_diff) << std::endl;
376 Kokkos::Profiling::popRegion();
382 #ifdef COMPADRE_USE_MPI
387 fprintf(stdout,
"Passed test \n");
390 fprintf(stdout,
"Failed test \n");
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) ...
std::string constraint_name
KOKKOS_INLINE_FUNCTION double trueSolution(double x, double y, double z, int order, int dimension)