Merge pull request #182 from SNLComputation/gradient_of_vector

Added GradientOfVectorPointEvaluation for several bases

Author: kuberry

examples/GMLS_Tutorial.hpp

@@ -142,6 +142,46 @@ double trueDivergence(double x, double y, double z, int order, int dimension) {
     return ans;
 }
 
+KOKKOS_INLINE_FUNCTION
+void trueHessian(double* ans, double x, double y, double z, int order, int dimension) {
+    for (int i=0; i<order+1; i++) {
+        for (int j=0; j<order+1; j++) {
+            for (int k=0; k<order+1; k++) {
+                if (i+j+k <= order) {
+                    // XX
+                    ans[0] += device_max(0,i)*device_max(0,i-1)*
+                               std::pow(x,device_max(0,i-2))*std::pow(y,j)*std::pow(z,k);
+                    if (dimension>1) {
+                        // XY
+                        ans[1] += device_max(0,i)*device_max(0,j)*
+                               std::pow(x,device_max(0,i-1))*std::pow(y,device_max(0,j-1))*std::pow(z,k);
+                        // YX = XY
+                        ans[1*dimension+0] = ans[1];
+                        // YY
+                        ans[1*dimension+1] += device_max(0,j)*device_max(0,j-1)*
+                               std::pow(x,i)*std::pow(y,device_max(0,j-2))*std::pow(z,k);
+                    }
+                    if (dimension>2) {
+                        // XZ
+                        ans[2] += device_max(0,i)*device_max(0,k)*
+                               std::pow(x,device_max(0,i-1))*std::pow(y,j)*std::pow(z,device_max(0,k-1));
+                        // YZ
+                        ans[1*dimension+2] += device_max(0,j)*device_max(0,k)*
+                               std::pow(x,i)*std::pow(y,device_max(0,j-1))*std::pow(z,device_max(0,k-1));
+                        // ZX = XZ
+                        ans[2*dimension+0] = ans[2];
+                        // ZY = YZ
+                        ans[2*dimension+1] = ans[1*dimension+2];
+                        // ZZ
+                        ans[2*dimension+2] += device_max(0,k)*device_max(0,k-1)*
+                               std::pow(x,i)*std::pow(y,j)*std::pow(z,device_max(0,k-2));
+                    }
+                }
+            }
+        }
+    }
+}
+
 KOKKOS_INLINE_FUNCTION
 double divergenceTestSamples(double x, double y, double z, int component, int dimension) {
     // solution can be captured exactly by at least 2rd order

examples/GMLS_Vector.cpp

@@ -347,11 +347,12 @@ bool all_passed = true;
     my_GMLS.setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
 
     // create a vector of target operations
-    std::vector<TargetOperation> lro(4);
+    std::vector<TargetOperation> lro(5);
     lro[0] = ScalarPointEvaluation;
     lro[1] = DivergenceOfVectorPointEvaluation;
     lro[2] = CurlOfVectorPointEvaluation;
     lro[3] = VectorPointEvaluation;
+    lro[4] = GradientOfVectorPointEvaluation;
 
     // and then pass them to the GMLS class
     my_GMLS.addTargets(lro);
@@ -398,6 +399,9 @@ bool all_passed = true;
 
     auto output_gradient = gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double**, Kokkos::HostSpace>
             (gradient_sampling_data_device, VectorPointEvaluation);
+
+    auto output_hessian = gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double**, Kokkos::HostSpace>
+            (gradient_sampling_data_device, GradientOfVectorPointEvaluation);
 
     // retrieves polynomial coefficients instead of remapped field
     auto scalar_coefficients = gmls_evaluator.applyFullPolynomialCoefficientsBasisToDataAllComponents<double**, Kokkos::HostSpace>
@@ -439,6 +443,16 @@ bool all_passed = true;
         double GMLS_CurlY = (dimension>1) ? output_curl(i,1) : 0;
         double GMLS_CurlZ = (dimension>2) ? output_curl(i,2) : 0;
 
+        // load hessian
+        double GMLS_GradXX = output_hessian(i,0*dimension+0);
+        double GMLS_GradXY = (dimension>1) ? output_hessian(i,0*dimension+1) : 0;
+        double GMLS_GradXZ = (dimension>2) ? output_hessian(i,0*dimension+2) : 0;
+        double GMLS_GradYX = (dimension>1) ? output_hessian(i,1*dimension+0) : 0;
+        double GMLS_GradYY = (dimension>1) ? output_hessian(i,1*dimension+1) : 0;
+        double GMLS_GradYZ = (dimension>2) ? output_hessian(i,1*dimension+2) : 0;
+        double GMLS_GradZX = (dimension>2) ? output_hessian(i,2*dimension+0) : 0;
+        double GMLS_GradZY = (dimension>2) ? output_hessian(i,2*dimension+1) : 0;
+        double GMLS_GradZZ = (dimension>2) ? output_hessian(i,2*dimension+2) : 0;
 
         // target site i's coordinate
         double xval = target_coords(i,0);
@@ -453,6 +467,9 @@ bool all_passed = true;
 
         double actual_Divergence;
         actual_Divergence = trueLaplacian(xval, yval, zval, order, dimension);
+
+        double actual_Hessian[9] = {0,0,0,0,0,0,0,0,0}; // initialized for 3, but only filled up to dimension
+        trueHessian(actual_Hessian, xval, yval, zval, order, dimension);
 
         double actual_Curl[3] = {0,0,0}; // initialized for 3, but only filled up to dimension 
         // (and not at all for dimimension = 1)
@@ -493,6 +510,48 @@ bool all_passed = true;
             all_passed = false;
             std::cout << i << " Failed Divergence by: " << std::abs(actual_Divergence - GMLS_Divergence) << std::endl;
         }
+
+        // check matrix (which is the hessian)
+        if(std::abs(actual_Hessian[0] - GMLS_GradXX) > failure_tolerance) {
+            all_passed = false;
+            std::cout << i << " Failed GradXX by: " << std::abs(actual_Hessian[0] - GMLS_GradXX) << std::endl;
+            if (dimension>1) {
+                if(std::abs(actual_Hessian[1] - GMLS_GradXY) > failure_tolerance) {
+                    all_passed = false;
+                    std::cout << i << " Failed GradXY by: " << std::abs(actual_Hessian[1] - GMLS_GradXY) << std::endl;
+                }
+                if(std::abs(actual_Hessian[1*dimension+1] - GMLS_GradYY) > failure_tolerance) {
+                    all_passed = false;
+                    std::cout << i << " Failed GradYY by: " << std::abs(actual_Hessian[1*dimension+1] - GMLS_GradYY) << std::endl;
+                }
+                if(std::abs(actual_Hessian[1*dimension+0] - GMLS_GradYX) > failure_tolerance) {
+                    all_passed = false;
+                    std::cout << i << " Failed GradYX by: " << std::abs(actual_Hessian[1*dimension+0] - GMLS_GradYX) << std::endl;
+                }
+            }
+            if (dimension>2) {
+                if(std::abs(actual_Hessian[2] - GMLS_GradXZ) > failure_tolerance) {
+                    all_passed = false;
+                    std::cout << i << " Failed GradXZ by: " << std::abs(actual_Hessian[2] - GMLS_GradXZ) << std::endl;
+                }
+                if(std::abs(actual_Hessian[1*dimension+2] - GMLS_GradYZ) > failure_tolerance) {
+                    all_passed = false;
+                    std::cout << i << " Failed GradYZ by: " << std::abs(actual_Hessian[1*dimension+2] - GMLS_GradYZ) << std::endl;
+                }
+                if(std::abs(actual_Hessian[2*dimension+0] - GMLS_GradZX) > failure_tolerance) {
+                    all_passed = false;
+                    std::cout << i << " Failed GradZX by: " << std::abs(actual_Hessian[2*dimension+0] - GMLS_GradZX) << std::endl;
+                }
+                if(std::abs(actual_Hessian[2*dimension+1] - GMLS_GradZY) > failure_tolerance) {
+                    all_passed = false;
+                    std::cout << i << " Failed GradZY by: " << std::abs(actual_Hessian[2*dimension+1] - GMLS_GradZY) << std::endl;
+                }
+                if(std::abs(actual_Hessian[2*dimension+2] - GMLS_GradZZ) > failure_tolerance) {
+                    all_passed = false;
+                    std::cout << i << " Failed GradZZ by: " << std::abs(actual_Hessian[2*dimension+2] - GMLS_GradZZ) << std::endl;
+                }
+            }
+        }
 
         // check curl
         if (order > 2) { // reconstructed solution not in basis unless order greater than 2 used

src/Compadre_GMLS_Targets.hpp

@@ -289,6 +289,38 @@ void GMLS::computeTargetFunctionals(const member_type& teamMember, scratch_vecto
                     }
                 });
                 additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
+            } else if (_operations(i) == TargetOperation::GradientOfScalarPointEvaluation) {
+                Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
+                    for (int e=0; e<num_evaluation_sites; ++e) { 
+                        int output_components = _basis_multiplier;
+                        for (int m2=0; m2<output_components; ++m2) { // output components 2
+                            int offset = getTargetOffsetIndexDevice(i, 0 /*in*/, m2 /*out*/, e);
+                            this->calcGradientPij(teamMember, t1.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, m2 /*partial_direction*/, _dimensions, _poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, e);
+                            for (int j=0; j<target_NP; ++j) {
+                                P_target_row(offset, j) = t1(j);
+                            }
+                        }
+                    }
+                });
+                additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
+            } else if (_operations(i) == TargetOperation::GradientOfVectorPointEvaluation) {
+                Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
+                    for (int e=0; e<num_evaluation_sites; ++e) { 
+                        for (int m0=0; m0<_sampling_multiplier; ++m0) { // input components
+                            int output_components = _basis_multiplier;
+                            for (int m1=0; m1<output_components; ++m1) { // output components 1
+                                for (int m2=0; m2<output_components; ++m2) { // output components 2
+                                    int offset = getTargetOffsetIndexDevice(i, m0 /*in*/, m1*output_components+m2 /*out*/, e);
+                                    this->calcGradientPij(teamMember, t1.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, m2 /*partial_direction*/, _dimensions, _poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, e);
+                                    for (int j=0; j<target_NP; ++j) {
+                                        P_target_row(offset, m1*target_NP + j) = t1(j);
+                                    }
+                                }
+                            }
+                        }
+                    }
+                });
+                additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
             } else if (_operations(i) == TargetOperation::DivergenceOfVectorPointEvaluation) {
                 if (_polynomial_sampling_functional == StaggeredEdgeIntegralSample) {
                     Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
@@ -478,6 +510,21 @@ void GMLS::computeTargetFunctionals(const member_type& teamMember, scratch_vecto
                     }
                 });
                 additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
+            } else if (_operations(i) == TargetOperation::GradientOfVectorPointEvaluation) {
+                Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
+                    for (int j=0; j<num_evaluation_sites; ++j) { 
+                        for (int m0=0; m0<_dimensions; ++m0) { // input components
+                            for (int m1=0; m1<_dimensions; ++m1) { // output components 1
+                                for (int m2=0; m2<_dimensions; ++m2) { // output componets 2
+                                    int offset = getTargetOffsetIndexDevice(i, m0 /*in*/, m1*_dimensions+m2 /*out*/, j);
+                                    auto row = Kokkos::subview(P_target_row, offset, Kokkos::ALL());
+                                    this->calcGradientPij(teamMember, row.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, m2 /*partial_direction*/, _dimensions, _poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, j);
+                                }
+                            }
+                        }
+                    }
+                });
+                additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
             } else if (_operations(i) == TargetOperation::PartialXOfScalarPointEvaluation) {
                 // copied from ScalarTaylorPolynomial
                 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {