9 #include <Compadre_Config.h> 16 #ifdef COMPADRE_USE_MPI 20 #include <Kokkos_Timer.hpp> 21 #include <Kokkos_Core.hpp> 41 double local_vec[3] = {0,0};
42 for (
int j=0; j<3; ++j) {
43 local_vec[0] += T(0,j) * u[j];
44 local_vec[1] += T(1,j) * u[j];
47 for (
int j=0; j<3; ++j) u[j] = 0;
48 for (
int j=0; j<3; ++j) {
49 u[j] += P(0, j) * local_vec[0];
50 u[j] += P(1, j) * local_vec[1];
58 int main (
int argc,
char* args[]) {
61 #ifdef COMPADRE_USE_MPI 62 MPI_Init(&argc, &args);
66 Kokkos::initialize(argc, args);
75 int constraint_type = 0;
77 int arg8toi = atoi(args[7]);
79 constraint_type = arg8toi;
89 int arg7toi = atoi(args[6]);
91 problem_type = arg7toi;
102 int arg6toi = atoi(args[5]);
104 solver_type = arg6toi;
110 int N_pts_on_sphere = 1000;
112 int arg5toi = atoi(args[4]);
114 N_pts_on_sphere = arg5toi;
122 int arg4toi = atoi(args[3]);
130 int number_target_coords = 200;
132 int arg3toi = atoi(args[2]);
134 number_target_coords = arg3toi;
142 int arg2toi = atoi(args[1]);
154 Kokkos::Profiling::pushRegion(
"Setup Point Data");
159 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> source_coords_device(
"source coordinates",
160 1.25*N_pts_on_sphere, 3);
161 Kokkos::View<double**>::HostMirror source_coords = Kokkos::create_mirror_view(source_coords_device);
166 int number_source_coords;
172 double a = 4*
PI*r*r/N_pts_on_sphere;
173 double d = std::sqrt(a);
174 int M_theta = std::round(
PI/d);
175 double d_theta =
PI/M_theta;
176 double d_phi = a/d_theta;
177 for (
int i=0; i<M_theta; ++i) {
178 double theta =
PI*(i + 0.5)/M_theta;
179 int M_phi = std::round(2*
PI*std::sin(theta)/d_phi);
180 for (
int j=0; j<M_phi; ++j) {
181 double phi = 2*
PI*j/M_phi;
182 source_coords(N_count, 0) = theta;
183 source_coords(N_count, 1) = phi;
188 for (
int i=0; i<N_count; ++i) {
189 double theta = source_coords(i,0);
190 double phi = source_coords(i,1);
191 source_coords(i,0) = r*std::sin(theta)*std::cos(phi);
192 source_coords(i,1) = r*std::sin(theta)*std::sin(phi);
193 source_coords(i,2) = r*cos(theta);
196 number_source_coords = N_count;
200 Kokkos::resize(source_coords, number_source_coords, 3);
201 Kokkos::resize(source_coords_device, number_source_coords, 3);
204 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> target_coords_device (
"target coordinates",
205 number_target_coords, 3);
206 Kokkos::View<double**>::HostMirror target_coords = Kokkos::create_mirror_view(target_coords_device);
209 bool enough_pts_found =
false;
213 double a = 4.0*
PI*r*r/((double)(5.0*number_target_coords));
214 double d = std::sqrt(a);
215 int M_theta = std::round(
PI/d);
216 double d_theta =
PI/((double)M_theta);
217 double d_phi = a/d_theta;
218 for (
int i=0; i<M_theta; ++i) {
219 double theta =
PI*(i + 0.5)/M_theta;
220 int M_phi = std::round(2*
PI*std::sin(theta)/d_phi);
221 for (
int j=0; j<M_phi; ++j) {
222 double phi = 2*
PI*j/M_phi;
223 target_coords(N_count, 0) = theta;
224 target_coords(N_count, 1) = phi;
226 if (N_count == number_target_coords) {
227 enough_pts_found =
true;
231 if (enough_pts_found)
break;
234 for (
int i=0; i<N_count; ++i) {
235 double theta = target_coords(i,0);
236 double phi = target_coords(i,1);
237 target_coords(i,0) = r*std::sin(theta)*std::cos(phi);
238 target_coords(i,1) = r*std::sin(theta)*std::sin(phi);
239 target_coords(i,2) = r*cos(theta);
247 Kokkos::Profiling::popRegion();
249 Kokkos::Profiling::pushRegion(
"Creating Data");
255 Kokkos::deep_copy(source_coords_device, source_coords);
258 Kokkos::deep_copy(target_coords_device, target_coords);
265 Kokkos::View<double*, Kokkos::DefaultExecutionSpace> sampling_data_device(
"samples of true solution",
266 source_coords_device.extent(0));
267 Kokkos::View<double*, Kokkos::DefaultExecutionSpace> ones_data_device(
"samples of ones",
268 source_coords_device.extent(0));
269 Kokkos::deep_copy(ones_data_device, 1.0);
272 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> sampling_vector_data_device(
"samples of vector true solution",
273 source_coords_device.extent(0), 3);
275 Kokkos::parallel_for(
"Sampling Manufactured Solutions", Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>
276 (0,source_coords.extent(0)), KOKKOS_LAMBDA(
const int i) {
279 double xval = source_coords_device(i,0);
280 double yval = (dimension>1) ? source_coords_device(i,1) : 0;
281 double zval = (dimension>2) ? source_coords_device(i,2) : 0;
286 for (
int j=0; j<3; ++j) {
287 double gradient[3] = {0,0,0};
289 sampling_vector_data_device(i,j) = gradient[j];
296 Kokkos::Profiling::popRegion();
297 Kokkos::Profiling::pushRegion(
"Neighbor Search");
308 double epsilon_multiplier = 1.7;
309 int estimated_upper_bound_number_neighbors =
310 point_cloud_search.getEstimatedNumberNeighborsUpperBound(min_neighbors, dimension, epsilon_multiplier);
312 Kokkos::View<int**, Kokkos::DefaultExecutionSpace> neighbor_lists_device(
"neighbor lists",
313 number_target_coords, estimated_upper_bound_number_neighbors);
314 Kokkos::View<int**>::HostMirror neighbor_lists = Kokkos::create_mirror_view(neighbor_lists_device);
317 Kokkos::View<double*, Kokkos::DefaultExecutionSpace> epsilon_device(
"h supports", number_target_coords);
318 Kokkos::View<double*>::HostMirror epsilon = Kokkos::create_mirror_view(epsilon_device);
323 point_cloud_search.generate2DNeighborListsFromKNNSearch(
false , target_coords, neighbor_lists,
324 epsilon, min_neighbors, epsilon_multiplier);
328 Kokkos::Profiling::popRegion();
339 Kokkos::deep_copy(neighbor_lists_device, neighbor_lists);
340 Kokkos::deep_copy(epsilon_device, epsilon);
343 std::string solver_name;
344 if (solver_type == 0) {
346 }
else if (solver_type == 1) {
348 }
else if (solver_type == 2) {
353 std::string problem_name;
354 if (problem_type == 0) {
355 problem_name =
"STANDARD";
356 }
else if (problem_type == 1) {
357 problem_name =
"MANIFOLD";
361 std::string constraint_name;
362 if (constraint_type == 0) {
363 constraint_name =
"NO_CONSTRAINT";
364 }
else if (constraint_type == 1) {
365 constraint_name =
"NEUMANN_GRAD_SCALAR";
370 GMLS my_GMLS_scalar(order, dimension,
371 solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
388 my_GMLS_scalar.
setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
392 my_GMLS_scalar.setAdditionalEvaluationSitesData(neighbor_lists_device, source_coords_device);
397 my_GMLS_scalar.setReferenceOutwardNormalDirection(target_coords_device,
true );
400 std::vector<TargetOperation> lro_scalar(3);
406 my_GMLS_scalar.addTargets(lro_scalar);
412 my_GMLS_scalar.setCurvatureWeightingPower(2);
418 my_GMLS_scalar.setWeightingPower(2);
421 my_GMLS_scalar.generateAlphas();
423 Kokkos::Profiling::pushRegion(
"Full Polynomial Basis GMLS Solution");
430 solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
433 my_GMLS_vector.
setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
434 my_GMLS_vector.setAdditionalEvaluationSitesData(neighbor_lists_device, source_coords_device);
435 std::vector<TargetOperation> lro_vector(2);
438 my_GMLS_vector.addTargets(lro_vector);
440 my_GMLS_vector.setCurvatureWeightingPower(2);
442 my_GMLS_vector.setWeightingPower(2);
443 my_GMLS_vector.generateAlphas();
444 Kokkos::Profiling::popRegion();
446 Kokkos::Profiling::pushRegion(
"Scalar Polynomial Basis Repeated to Form a Vector GMLS Solution");
471 solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
474 my_GMLS_vector_of_scalar_clones.
setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
475 my_GMLS_vector_of_scalar_clones.setAdditionalEvaluationSitesData(neighbor_lists_device, source_coords_device);
476 std::vector<TargetOperation> lro_vector_of_scalar_clones(2);
479 my_GMLS_vector_of_scalar_clones.addTargets(lro_vector_of_scalar_clones);
481 my_GMLS_vector_of_scalar_clones.setCurvatureWeightingPower(2);
483 my_GMLS_vector_of_scalar_clones.setWeightingPower(2);
484 my_GMLS_vector_of_scalar_clones.generateAlphas();
485 Kokkos::Profiling::popRegion();
490 double instantiation_time = timer.seconds();
491 std::cout <<
"Took " << instantiation_time <<
"s to complete alphas generation." << std::endl;
493 Kokkos::Profiling::pushRegion(
"Apply Alphas to Data");
508 Evaluator scalar_gmls_evaluator(&my_GMLS_scalar);
509 Evaluator vector_gmls_evaluator(&my_GMLS_vector);
510 Evaluator vector_gmls_evaluator_of_scalar_clones(&my_GMLS_vector_of_scalar_clones);
540 auto output_vector_of_scalar_clones =
580 Kokkos::Profiling::popRegion();
582 Kokkos::Profiling::pushRegion(
"Comparison");
586 double tangent_bundle_error = 0;
587 double tangent_bundle_norm = 0;
588 double values_error = 0;
589 double values_norm = 0;
592 double curl_ambient_error = 0;
593 double curl_ambient_norm = 0;
598 double vector_ambient_error = 0;
599 double vector_ambient_norm = 0;
602 double vector_of_scalar_clones_ambient_error = 0;
603 double vector_of_scalar_clones_ambient_norm = 0;
608 auto d_prestencil_weights = my_GMLS_vector_of_scalar_clones.getPrestencilWeights();
609 auto prestencil_weights = Kokkos::create_mirror_view(d_prestencil_weights);
610 Kokkos::deep_copy(prestencil_weights, d_prestencil_weights);
613 auto d_tangent_directions = my_GMLS_vector_of_scalar_clones.getTangentDirections();
614 auto tangent_directions = Kokkos::create_mirror_view(d_tangent_directions);
615 Kokkos::deep_copy(tangent_directions, d_tangent_directions);
618 for (
int i=0; i<number_target_coords; i++) {
626 double GMLS_value = output_value(i);
629 double GMLS_gc = output_gaussian_curvature(i);
639 double xval = source_coords(neighbor_lists(i,1),0);
640 double yval = (dimension>1) ? source_coords(neighbor_lists(i,1),1) : 0;
641 double zval = (dimension>2) ? source_coords(neighbor_lists(i,1),2) : 0;
642 double coord[3] = {xval, yval, zval};
645 for (
int j=0; j<dimension-1; ++j) {
646 double tangent_inner_prod = 0;
647 for (
int k=0; k<dimension; ++k) {
648 tangent_inner_prod += coord[k] * prestencil_weights(0, i, 0 , j, k);
650 tangent_bundle_error += tangent_inner_prod * tangent_inner_prod;
652 double normal_inner_prod = 0;
653 for (
int k=0; k<dimension; ++k) {
654 normal_inner_prod += coord[k] * my_GMLS_scalar.getTangentBundle(i, dimension-1, k);
657 double normal_error_to_sum = (normal_inner_prod > 0) ? normal_inner_prod - 1 : normal_inner_prod + 1;
658 tangent_bundle_error += normal_error_to_sum * normal_error_to_sum;
659 tangent_bundle_norm += 1;
663 double actual_gc = 1.0;
665 double actual_Gradient_ambient[3] = {0,0,0};
668 double actual_Curl_ambient[3] = {0,0,0};
671 values_error += (GMLS_value - actual_value)*(GMLS_value - actual_value);
672 values_norm += actual_value*actual_value;
674 gc_error += (GMLS_gc - actual_gc)*(GMLS_gc - actual_gc);
677 for (
int j=0; j<dimension; ++j) u[j] = output_curl(i,j);
679 for (
int j=0; j<dimension; ++j) output_curl(i,j) = u[j];
681 for (
int j=0; j<dimension; ++j) {
682 curl_ambient_error += (output_curl(i,j) - actual_Curl_ambient[j])*(output_curl(i,j) - actual_Curl_ambient[j]);
683 curl_ambient_norm += actual_Curl_ambient[j]*actual_Curl_ambient[j];
694 for (
int j=0; j<dimension; ++j) u[j] = output_vector(i,j);
696 for (
int j=0; j<dimension; ++j) output_vector(i,j) = u[j];
698 for (
int j=0; j<dimension; ++j) {
699 vector_ambient_error += (output_vector(i,j) - actual_Gradient_ambient[j])*(output_vector(i,j) - actual_Gradient_ambient[j]);
700 vector_ambient_norm += actual_Gradient_ambient[j]*actual_Gradient_ambient[j];
707 for (
int j=0; j<dimension; ++j) u[j] = output_vector_of_scalar_clones(i,j);
709 for (
int j=0; j<dimension; ++j) output_vector_of_scalar_clones(i,j) = u[j];
724 for (
int j=0; j<dimension; ++j) {
725 vector_of_scalar_clones_ambient_error += (output_vector_of_scalar_clones(i,j) - actual_Gradient_ambient[j])*(output_vector_of_scalar_clones(i,j) - actual_Gradient_ambient[j]);
726 vector_of_scalar_clones_ambient_norm += actual_Gradient_ambient[j]*actual_Gradient_ambient[j];
734 tangent_bundle_error /= number_target_coords;
735 tangent_bundle_error = std::sqrt(tangent_bundle_error);
736 tangent_bundle_norm /= number_target_coords;
737 tangent_bundle_norm = std::sqrt(tangent_bundle_norm);
739 values_error /= number_target_coords;
740 values_error = std::sqrt(values_error);
741 values_norm /= number_target_coords;
742 values_norm = std::sqrt(values_norm);
744 gc_error /= number_target_coords;
745 gc_error = std::sqrt(gc_error);
746 gc_norm /= number_target_coords;
747 gc_norm = std::sqrt(gc_norm);
749 curl_ambient_error /= number_target_coords;
750 curl_ambient_error = std::sqrt(curl_ambient_error);
751 curl_ambient_norm /= number_target_coords;
752 curl_ambient_norm = std::sqrt(curl_ambient_norm);
764 vector_ambient_error /= number_target_coords;
765 vector_ambient_error = std::sqrt(vector_ambient_error);
766 vector_ambient_norm /= number_target_coords;
767 vector_ambient_norm = std::sqrt(vector_ambient_norm);
774 vector_of_scalar_clones_ambient_error /= number_target_coords;
775 vector_of_scalar_clones_ambient_error = std::sqrt(vector_of_scalar_clones_ambient_error);
776 vector_of_scalar_clones_ambient_norm /= number_target_coords;
777 vector_of_scalar_clones_ambient_norm = std::sqrt(vector_of_scalar_clones_ambient_norm);
784 printf(
"Tangent Bundle Error: %g\n", tangent_bundle_error / tangent_bundle_norm);
785 printf(
"Point Value Error: %g\n", values_error / values_norm);
786 printf(
"Gaussian Curvature Error: %g\n", gc_error / gc_norm);
787 printf(
"Surface Curl (Ambient) Error: %g\n", curl_ambient_error / curl_ambient_norm);
790 printf(
"Surface Vector (VectorBasis) Error: %g\n", vector_ambient_error / vector_ambient_norm);
792 printf(
"Surface Vector (ScalarClones) Error: %g\n",
793 vector_of_scalar_clones_ambient_error / vector_of_scalar_clones_ambient_norm);
799 Kokkos::Profiling::popRegion();
808 #ifdef COMPADRE_USE_MPI KOKKOS_INLINE_FUNCTION double sphere_harmonic54(double x, double y, double z)
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...
Point evaluation of a vector (reconstructs entire vector at once, requiring a ReconstructionSpace hav...
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 VaryingManifoldVectorPointSample
For integrating polynomial dotted with normal over an edge.
KOKKOS_INLINE_FUNCTION void curl_sphere_harmonic54(double *curl, double x, double y, double z)
Point evaluation of Gaussian curvature.
void AmbientLocalAmbient(XYZ &u, double *T_data, double *P_data)
[Ambient to Local Back To Ambient Helper Function]
#define TO_GLOBAL(variable)
Point evaluation of a scalar.
Kokkos::View< double **, layout_right, host_execution_space, Kokkos::MemoryTraits< Kokkos::Unmanaged > > host_scratch_matrix_right_type
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...
Scalar basis reused as many times as there are components in the vector resulting in a much cheaper p...
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.
int main(int argc, char *args[])
[Ambient to Local Back To Ambient Helper Function]
KOKKOS_INLINE_FUNCTION void gradient_sphereHarmonic54_ambient(double *gradient, double x, double y, double z)
Generalized Moving Least Squares (GMLS)
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) ...
constexpr SamplingFunctional PointSample
Available sampling functionals.
Vector polynomial basis having # of components _dimensions, or (_dimensions-1) in the case of manifol...