9 #include <Compadre_Config.h> 
   14 #ifdef COMPADRE_USE_MPI 
   18 #include <Kokkos_Timer.hpp> 
   19 #include <Kokkos_Core.hpp> 
   21 using namespace Compadre;
 
   23 int main (
int argc, 
char* args[])
 
   26 #ifdef COMPADRE_USE_MPI 
   27     MPI_Init(&argc, &args);
 
   30 bool all_passed = 
true;
 
   37     int constraint_type = 0; 
 
   39         int arg7toi = atoi(args[6]);
 
   41             constraint_type = arg7toi;
 
   51         int arg6toi = atoi(args[5]);
 
   53             problem_type = arg6toi;
 
   64         int arg5toi = atoi(args[4]);
 
   66             solver_type = arg5toi;
 
   72         int arg4toi = atoi(args[3]);
 
   78     const double failure_tolerance = 1e-9;
 
   80     const int order = atoi(args[1]); 
 
   82     const int offset = 15;
 
   86     std::cout << min_neighbors << 
" " << max_neighbors << std::endl;
 
   87     std::uniform_int_distribution<int> gen_num_neighbors(min_neighbors, max_neighbors); 
 
   90     Kokkos::initialize(argc, args);
 
   92     Kokkos::Profiling::pushRegion(
"Setup");
 
   95     const int number_target_coords = atoi(args[2]); 
 
   97     std::uniform_int_distribution<int> gen_neighbor_number(offset, N); 
 
  100     Kokkos::View<int**, Kokkos::HostSpace>    neighbor_lists(
"neighbor lists", number_target_coords, max_neighbors+1); 
 
  101     Kokkos::View<double**, Kokkos::HostSpace> source_coords(
"neighbor coordinates", N, dimension);
 
  102     Kokkos::View<double*, Kokkos::HostSpace> epsilon(
"h supports", number_target_coords);
 
  104     for (
int i=0; i<number_target_coords; i++) {
 
  109     for(
int i = 0; i < offset; i++){
 
  110         for(
int j = 0; j < dimension; j++){
 
  111             source_coords(i,j) = 0.1;
 
  116     for(
int i = offset; i < N; i++){ 
 
  117         double randx = (2.0*(double)rand() / (double) RAND_MAX - 1.0)*epsilon(0)/2.0;
 
  118         double randy = (2.0*(double)rand() / (double) RAND_MAX - 1.0)*epsilon(0)/2.0;
 
  119         double randz = (2.0*(double)rand() / (double) RAND_MAX - 1.0)*epsilon(0)/2.0;
 
  120         source_coords(i,0) = randx;
 
  121         if (dimension>1) source_coords(i,1) = randy;
 
  122         if (dimension>2) source_coords(i,2) = randz;
 
  125     const double target_epsilon = 0.1;
 
  127     Kokkos::View<double**, Kokkos::HostSpace> target_coords (
"target coordinates", number_target_coords, dimension);
 
  128     for(
int i = 0; i < number_target_coords; i++){ 
 
  129         double randx = (2.0*(double)rand() / (double) RAND_MAX - 1.0)*target_epsilon/2.0;
 
  130         double randy = (2.0*(double)rand() / (double) RAND_MAX - 1.0)*target_epsilon/2.0;
 
  131         double randz = (2.0*(double)rand() / (double) RAND_MAX - 1.0)*target_epsilon/2.0;
 
  132         target_coords(i,0) = randx;
 
  133         if (dimension>1) target_coords(i,1) = randy;
 
  134         if (dimension>2) target_coords(i,2) = randz;
 
  138     for (
int i=0; i<number_target_coords; i++) {
 
  141         int r = max_neighbors;
 
  142         neighbor_lists(i,0) = r; 
 
  144         for(
int j=0; j<r; j++){
 
  145             neighbor_lists(i,j+1) = offset + j + 1;
 
  164     Kokkos::Profiling::popRegion();
 
  168     std::string solver_name;
 
  169     if (solver_type == 0) { 
 
  171     } 
else if (solver_type == 1) { 
 
  173     } 
else if (solver_type == 2) { 
 
  178     std::string problem_name;
 
  179     if (problem_type == 0) { 
 
  180         problem_name = 
"STANDARD";
 
  181     } 
else if (problem_type == 1) { 
 
  182         problem_name = 
"MANIFOLD";
 
  186     std::string constraint_name;
 
  187     if (constraint_type == 0) { 
 
  188         constraint_name = 
"NO_CONSTRAINT";
 
  189     } 
else if (constraint_type == 1) { 
 
  190         constraint_name = 
"NEUMANN_GRAD_SCALAR";
 
  193     GMLS my_GMLS(order, dimension,
 
  194                  solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
 
  196     my_GMLS.
setProblemData(neighbor_lists, source_coords, target_coords, epsilon);
 
  197     my_GMLS.setWeightingPower(10);
 
  199     std::vector<TargetOperation> lro(5);
 
  205     my_GMLS.addTargets(lro);
 
  206     my_GMLS.generateAlphas();
 
  208     double instantiation_time = timer.seconds();
 
  209     std::cout << 
"Took " << instantiation_time << 
"s to complete instantiation." << std::endl;
 
  211     Kokkos::Profiling::pushRegion(
"Creating Data");
 
  215     Kokkos::View<double*, Kokkos::HostSpace> sampling_data(
"samples of true solution", source_coords.extent(0));
 
  216     Kokkos::View<double**, Kokkos::HostSpace> gradient_sampling_data(
"samples of true gradient", source_coords.extent(0), dimension);
 
  217     Kokkos::View<double**, Kokkos::LayoutLeft, Kokkos::HostSpace> divergence_sampling_data(
"samples of true solution for divergence test", source_coords.extent(0), dimension);
 
  218     Kokkos::parallel_for(
"Sampling Manufactured Solutions", Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>(0,source_coords.extent(0)), KOKKOS_LAMBDA(
const int i) {
 
  219         double xval = source_coords(i,0);
 
  220         double yval = (dimension>1) ? source_coords(i,1) : 0;
 
  221         double zval = (dimension>2) ? source_coords(i,2) : 0;
 
  222         sampling_data(i) = 
trueSolution(xval, yval, zval, order, dimension);
 
  223         double true_grad[3] = {0,0,0};
 
  224         trueGradient(true_grad, xval, yval,zval, order, dimension);
 
  225         for (
int j=0; j<dimension; ++j) {
 
  227             gradient_sampling_data(i,j) = true_grad[j];
 
  230     Kokkos::Profiling::popRegion();
 
  234     for (
int i=0; i<number_target_coords; i++) {
 
  236         Kokkos::Profiling::pushRegion(
"Apply Alphas to Data");
 
  310         double GMLS_CurlX = 0.0;
 
  311         double GMLS_CurlY = 0.0;
 
  312         double GMLS_CurlZ = 0.0;
 
  314             for (
int j=0; j<dimension; ++j) {
 
  321             for (
int j=0; j<dimension; ++j) {
 
  327         Kokkos::Profiling::popRegion();
 
  359         Kokkos::Profiling::pushRegion(
"Comparison");
 
  361         double xval = target_coords(i,0);
 
  362         double yval = (dimension>1) ? target_coords(i,1) : 0;
 
  363         double zval = (dimension>2) ? target_coords(i,2) : 0;
 
  365         double actual_value = 
trueSolution(xval, yval, zval, order, dimension);
 
  366         double actual_Laplacian = 
trueLaplacian(xval, yval, zval, order, dimension);
 
  367         double actual_Gradient[3] = {0,0,0};
 
  368         trueGradient(actual_Gradient, xval, yval, zval, order, dimension);
 
  369         double actual_Divergence;
 
  370         actual_Divergence = 
trueLaplacian(xval, yval, zval, order, dimension);
 
  372         double actual_CurlX = 0;
 
  373         double actual_CurlY = 0;
 
  374         double actual_CurlZ = 0;
 
  388         if(GMLS_value!=GMLS_value || std::abs(actual_value - GMLS_value) > failure_tolerance) {
 
  390             std::cout << 
"Failed Actual by: " << std::abs(actual_value - GMLS_value) << std::endl;
 
  393         if(std::abs(actual_Laplacian - GMLS_Laplacian) > failure_tolerance) {
 
  395             std::cout << 
"Failed Laplacian by: " << std::abs(actual_Laplacian - GMLS_Laplacian) << std::endl;
 
  398         if(std::abs(actual_Gradient[0] - GMLS_GradX) > failure_tolerance) {
 
  400             std::cout << 
"Failed GradX by: " << std::abs(actual_Gradient[0] - GMLS_GradX) << std::endl;
 
  404             if(std::abs(actual_Gradient[1] - GMLS_GradY) > failure_tolerance) {
 
  406                 std::cout << 
"Failed GradY by: " << std::abs(actual_Gradient[1] - GMLS_GradY) << std::endl;
 
  411             if(std::abs(actual_Gradient[2] - GMLS_GradZ) > failure_tolerance) {
 
  413                 std::cout << 
"Failed GradZ by: " << std::abs(actual_Gradient[2] - GMLS_GradZ) << std::endl;
 
  417         if(std::abs(actual_Divergence - GMLS_Divergence) > failure_tolerance) {
 
  419             std::cout << 
"Failed Divergence by: " << std::abs(actual_Divergence - GMLS_Divergence) << std::endl;
 
  424             tmp_diff += std::abs(actual_CurlX - GMLS_CurlX) + std::abs(actual_CurlY - GMLS_CurlY);
 
  426             tmp_diff += std::abs(actual_CurlZ - GMLS_CurlZ);
 
  427         if(std::abs(tmp_diff) > failure_tolerance) {
 
  429             std::cout << 
"Failed Curl by: " << std::abs(tmp_diff) << std::endl;
 
  431         Kokkos::Profiling::popRegion();
 
  437 #ifdef COMPADRE_USE_MPI 
  442     fprintf(stdout, 
"Passed test \n");
 
  445     fprintf(stdout, 
"Failed test \n");
 
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...
 
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...
 
int main(int argc, char *args[])
[Parse Command Line Arguments] 
 
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)
 
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)