Merge pull request #1262 from CEED/jeremy/grad-at-points

Add CEED_EVAL_GRAD support for CeedBasisApplyAtPoints

Author: jeremylt

doc/sphinx/source/releasenotes.md

@@ -17,6 +17,7 @@ For example, `CeedOperatorContextGetFieldLabel` was renamed to `CeedOperatorGetC
 - Update `/cpu/self/memcheck/*` backends to help verify `CeedVector` array access assumptions and `CeedQFunction` user output assumptions.
 - Update {c:func}`CeedOperatorLinearAssembleDiagonal` to provide default implementation that supports `CeedOperator` with multiple active bases.
 - Added Sycl backends `/gpu/sycl/ref` and `/gpu/sycl/shared`
+- Added {c:func}`CeedBasisApplyAtPoints` for evalution of values and derivaties at arbitrary points inside elements
 
 ### Examples
 
@@ -60,7 +61,7 @@ This issue will be fixed in a future OCCA release.
 - Support for primitive variables for more accurate boundary layers and all-speed flow.
 - Added $YZ\beta$ shock capturing scheme and Shock Tube example.
 - Added Channel example, with comparison to analytic solutions.
-- Added Flat Plate with boundary layer mesh and compressible Blasius inflow condition based on Chebyshev collocation solution of the Blasius equations. 
+- Added Flat Plate with boundary layer mesh and compressible Blasius inflow condition based on Chebyshev collocation solution of the Blasius equations.
 - Added strong and weak synthetic turbulence generation (STG) inflow boundary conditions.
 - Added "freestream" boundary conditions based on HLLC Riemann solver.
 - Automated stabilization coefficients for different basis degree.

interface/ceed-basis.c

@@ -34,6 +34,53 @@ const CeedBasis CEED_BASIS_COLLOCATED = &ceed_basis_collocated;
 /// @addtogroup CeedBasisDeveloper
 /// @{
 
+/**
+  @brief Compute Chebyshev polynomial values at a point
+
+  @param[in]  x           Coordinate to evaluate Chebyshev polynomials at
+  @param[in]  n           Number of Chebyshev polynomials to evaluate, n >= 2
+  @param[out] chebyshev_x Array of Chebyshev polynomial values
+
+  @return An error code: 0 - success, otherwise - failure
+
+  @ref Developer
+**/
+static int CeedChebyshevPolynomialsAtPoint(CeedScalar x, CeedInt n, CeedScalar *chebyshev_x) {
+  chebyshev_x[0] = 1.0;
+  chebyshev_x[1] = 2 * x;
+  for (CeedInt i = 2; i < n; i++) chebyshev_x[i] = 2 * x * chebyshev_x[i - 1] - chebyshev_x[i - 2];
+
+  return CEED_ERROR_SUCCESS;
+}
+
+/**
+  @brief Compute values of the derivative of Chebyshev polynomials at a point
+
+  @param[in]  x           Coordinate to evaluate derivative of Chebyshev polynomials at
+  @param[in]  n           Number of Chebyshev polynomials to evaluate, n >= 2
+  @param[out] chebyshev_x Array of Chebyshev polynomial derivative values
+
+  @return An error code: 0 - success, otherwise - failure
+
+  @ref Developer
+**/
+static int CeedChebyshevDerivativeAtPoint(CeedScalar x, CeedInt n, CeedScalar *chebyshev_dx) {
+  CeedScalar chebyshev_x[3];
+
+  chebyshev_x[1]  = 1.0;
+  chebyshev_x[2]  = 2 * x;
+  chebyshev_dx[0] = 0.0;
+  chebyshev_dx[1] = 2.0;
+  for (CeedInt i = 2; i < n; i++) {
+    chebyshev_x[0]  = chebyshev_x[1];
+    chebyshev_x[1]  = chebyshev_x[2];
+    chebyshev_x[2]  = 2 * x * chebyshev_x[1] - chebyshev_x[0];
+    chebyshev_dx[i] = 2 * x * chebyshev_dx[i - 1] + 2 * chebyshev_x[1] - chebyshev_dx[i - 2];
+  }
+
+  return CEED_ERROR_SUCCESS;
+}
+
 /**
   @brief Compute Householder reflection
 
@@ -1438,6 +1485,9 @@ int CeedBasisApplyAtPoints(CeedBasis basis, CeedInt num_points, CeedTransposeMod
                    (t_mode == CEED_NOTRANSPOSE && (v_length >= num_points * num_q_comp || u_length >= num_nodes * num_comp)));
       break;
     case CEED_EVAL_GRAD:
+      good_dims = ((t_mode == CEED_TRANSPOSE && (u_length >= num_points * num_q_comp * dim || v_length >= num_nodes * num_comp)) ||
+                   (t_mode == CEED_NOTRANSPOSE && (v_length >= num_points * num_q_comp * dim || u_length >= num_nodes * num_comp)));
+      break;
     case CEED_EVAL_NONE:
     case CEED_EVAL_WEIGHT:
     case CEED_EVAL_DIV:
@@ -1456,6 +1506,8 @@ int CeedBasisApplyAtPoints(CeedBasis basis, CeedInt num_points, CeedTransposeMod
 
   // Default implementation
   CeedCheck(basis->is_tensor_basis, basis->ceed, CEED_ERROR_UNSUPPORTED, "Evaluation at arbitrary points only supported for tensor product bases");
+  CeedCheck(eval_mode == CEED_EVAL_INTERP || t_mode == CEED_NOTRANSPOSE, basis->ceed, CEED_ERROR_UNSUPPORTED, "%s evaluation only supported for %s",
+            CeedEvalModes[eval_mode], CeedTransposeModes[CEED_NOTRANSPOSE]);
   if (!basis->basis_chebyshev) {
     // Build matrix mapping from quadrature point values to Chebyshev coefficients
     CeedScalar       *tau, *C, *I, *chebyshev_coeffs_1d;
@@ -1466,13 +1518,7 @@ int CeedBasisApplyAtPoints(CeedBasis basis, CeedInt num_points, CeedTransposeMod
     CeedCheck(P_1d > 0 && Q_1d > 0, basis->ceed, CEED_ERROR_INCOMPATIBLE, "Basis dimensions are malformed");
     CeedCall(CeedCalloc(Q_1d * Q_1d, &C));
     CeedCall(CeedBasisGetQRef(basis, &q_ref_1d));
-    for (CeedInt i = 0; i < Q_1d; i++) {
-      const CeedScalar x = q_ref_1d[i];
-
-      C[i * Q_1d + 0] = 1.0;
-      C[i * Q_1d + 1] = 2 * x;
-      for (CeedInt j = 2; j < Q_1d; j++) C[i * Q_1d + j] = 2 * x * C[i * Q_1d + j - 1] - C[i * Q_1d + j - 2];
-    }
+    for (CeedInt i = 0; i < Q_1d; i++) CeedCall(CeedChebyshevPolynomialsAtPoint(q_ref_1d[i], Q_1d, &C[i * Q_1d]));
 
     // Inverse of coefficient matrix
     CeedCall(CeedCalloc(Q_1d * Q_1d, &chebyshev_coeffs_1d));
@@ -1546,36 +1592,61 @@ int CeedBasisApplyAtPoints(CeedBasis basis, CeedInt num_points, CeedTransposeMod
       CeedCall(CeedVectorGetArrayRead(basis->vec_chebyshev, CEED_MEM_HOST, &chebyshev_coeffs));
       CeedCall(CeedVectorGetArrayRead(x_ref, CEED_MEM_HOST, &x_array_read));
       CeedCall(CeedVectorGetArrayWrite(v, CEED_MEM_HOST, &v_array));
-      {
-        CeedScalar tmp[2][num_comp * CeedIntPow(Q_1d, dim)], chebyshev_x[Q_1d];
-
-        // ---- Values at point
-        for (CeedInt p = 0; p < num_points; p++) {
-          CeedInt pre = num_comp * CeedIntPow(Q_1d, dim - 1), post = 1;
-
-          for (CeedInt d = dim - 1; d >= 0; d--) {
-            // ------ Compute Chebyshev polynomial values
-            {
-              const CeedScalar x = x_array_read[p * dim + d];
-
-              chebyshev_x[0] = 1.0;
-              chebyshev_x[1] = 2 * x;
-              for (CeedInt j = 2; j < Q_1d; j++) chebyshev_x[j] = 2 * x * chebyshev_x[j - 1] - chebyshev_x[j - 2];
+      switch (eval_mode) {
+        case CEED_EVAL_INTERP: {
+          CeedScalar tmp[2][num_comp * CeedIntPow(Q_1d, dim)], chebyshev_x[Q_1d];
+
+          // ---- Values at point
+          for (CeedInt p = 0; p < num_points; p++) {
+            CeedInt pre = num_comp * CeedIntPow(Q_1d, dim - 1), post = 1;
+
+            // Note: stepping "backwards" through the tensor contractions to agree with the ordering of the Chebyshev coefficients
+            for (CeedInt d = dim - 1; d >= 0; d--) {
+              // ------ Tensor contract with current Chebyshev polynomial values
+              CeedCall(CeedChebyshevPolynomialsAtPoint(x_array_read[p * dim + d], Q_1d, chebyshev_x));
+              CeedCall(CeedTensorContractApply(basis->contract, pre, Q_1d, post, 1, chebyshev_x, t_mode, false,
+                                               d == (dim - 1) ? chebyshev_coeffs : tmp[d % 2], d == 0 ? &v_array[p * num_comp] : tmp[(d + 1) % 2]));
+              pre /= Q_1d;
+              post *= 1;
             }
-            // ------ Tensor contract
-            CeedCall(CeedTensorContractApply(basis->contract, pre, Q_1d, post, 1, chebyshev_x, t_mode, false,
-                                             d == (dim - 1) ? chebyshev_coeffs : tmp[d % 2], d == 0 ? &v_array[p * num_comp] : tmp[(d + 1) % 2]));
-            pre /= Q_1d;
-            post *= 1;
           }
+          break;
         }
+        case CEED_EVAL_GRAD: {
+          CeedScalar tmp[2][num_comp * CeedIntPow(Q_1d, dim)], chebyshev_x[Q_1d];
+
+          // ---- Values at point
+          for (CeedInt p = 0; p < num_points; p++) {
+            // Note: stepping "backwards" through the tensor contractions to agree with the ordering of the Chebyshev coefficients
+            // Dim**2 contractions, apply grad when pass == dim
+            for (CeedInt pass = dim - 1; pass >= 0; pass--) {
+              CeedInt pre = num_comp * CeedIntPow(Q_1d, dim - 1), post = 1;
+
+              for (CeedInt d = dim - 1; d >= 0; d--) {
+                // ------ Tensor contract with current Chebyshev polynomial values
+                if (pass == d) CeedCall(CeedChebyshevDerivativeAtPoint(x_array_read[p * dim + d], Q_1d, chebyshev_x));
+                else CeedCall(CeedChebyshevPolynomialsAtPoint(x_array_read[p * dim + d], Q_1d, chebyshev_x));
+                CeedCall(CeedTensorContractApply(basis->contract, pre, Q_1d, post, 1, chebyshev_x, t_mode, false,
+                                                 d == (dim - 1) ? chebyshev_coeffs : tmp[d % 2],
+                                                 d == 0 ? &v_array[p * num_comp * dim + pass] : tmp[(d + 1) % 2]));
+                pre /= Q_1d;
+                post *= 1;
+              }
+            }
+          }
+          break;
+        }
+        default:
+          // Nothing to do, this won't occur
+          break;
       }
       CeedCall(CeedVectorRestoreArrayRead(basis->vec_chebyshev, &chebyshev_coeffs));
       CeedCall(CeedVectorRestoreArrayRead(x_ref, &x_array_read));
       CeedCall(CeedVectorRestoreArray(v, &v_array));
       break;
     }
     case CEED_TRANSPOSE: {
+      // Note: No switch on e_mode here because only CEED_EVAL_INTERP is supported at this time
       // Arbitrary points to nodes
       CeedScalar       *chebyshev_coeffs;
       const CeedScalar *u_array, *x_array_read;
@@ -1591,16 +1662,10 @@ int CeedBasisApplyAtPoints(CeedBasis basis, CeedInt num_points, CeedTransposeMod
         for (CeedInt p = 0; p < num_points; p++) {
           CeedInt pre = num_comp * 1, post = 1;
 
+          // Note: stepping "backwards" through the tensor contractions to agree with the ordering of the Chebyshev coefficients
           for (CeedInt d = dim - 1; d >= 0; d--) {
-            // ------ Compute Chebyshev polynomial values
-            {
-              const CeedScalar x = x_array_read[p * dim + d];
-
-              chebyshev_x[0] = 1.0;
-              chebyshev_x[1] = 2 * x;
-              for (CeedInt j = 2; j < Q_1d; j++) chebyshev_x[j] = 2 * x * chebyshev_x[j - 1] - chebyshev_x[j - 2];
-            }
-            // ------ Tensor contract
+            // ------ Tensor contract with current Chebyshev polynomial values
+            CeedCall(CeedChebyshevPolynomialsAtPoint(x_array_read[p * dim + d], Q_1d, chebyshev_x));
             CeedCall(CeedTensorContractApply(basis->contract, pre, 1, post, Q_1d, chebyshev_x, t_mode, p > 0 && d == 0,
                                              d == (dim - 1) ? &u_array[p * num_comp] : tmp[d % 2], d == 0 ? chebyshev_coeffs : tmp[(d + 1) % 2]));
             pre /= 1;

tests/t355-basis.c

@@ -0,0 +1,89 @@
+/// @file
+/// Test polynomial gradient to arbirtary points in 1D
+/// \test Test polynomial gradient to arbitrary points in 1D
+#include <ceed.h>
+#include <math.h>
+#include <stdio.h>
+
+#define ALEN(a) (sizeof(a) / sizeof((a)[0]))
+
+static CeedScalar Eval(CeedScalar x, CeedInt n, const CeedScalar *c) {
+  CeedScalar y = c[n - 1];
+  for (CeedInt i = n - 2; i >= 0; i--) y = y * x + c[i];
+  return y;
+}
+
+static CeedScalar EvalGrad(CeedScalar x, CeedInt n, const CeedScalar *c) {
+  CeedScalar y = (n - 1) * c[n - 1];
+  for (CeedInt i = n - 2; i >= 1; i--) y = y * x + i * c[i];
+  return y;
+}
+
+int main(int argc, char **argv) {
+  Ceed             ceed;
+  CeedVector       x, x_nodes, x_points, u, v;
+  CeedBasis        basis_x, basis_u;
+  const CeedInt    p = 5, q = 5, num_points = 4;
+  const CeedScalar c[4] = {1, 2, 3, 4};  // f = 1 + 2x + 3x^2 + ..., df = 2 + 6x + 12x^2 + ...
+
+  CeedInit(argv[1], &ceed);
+
+  CeedVectorCreate(ceed, 2, &x);
+  CeedVectorCreate(ceed, p, &x_nodes);
+  CeedVectorCreate(ceed, num_points, &x_points);
+  CeedVectorCreate(ceed, p, &u);
+  CeedVectorCreate(ceed, num_points, &v);
+
+  // Get nodal coordinates
+  CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, 2, p, CEED_GAUSS_LOBATTO, &basis_x);
+  {
+    CeedScalar x_array[2];
+
+    for (CeedInt i = 0; i < 2; i++) x_array[i] = CeedIntPow(-1, i + 1);
+    CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, x_array);
+  }
+  CeedBasisApply(basis_x, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x, x_nodes);
+
+  // Set values of u at nodes
+  {
+    const CeedScalar *x_array;
+    CeedScalar        u_array[p];
+
+    CeedVectorGetArrayRead(x_nodes, CEED_MEM_HOST, &x_array);
+    for (CeedInt i = 0; i < p; i++) u_array[i] = Eval(x_array[i], ALEN(c), c);
+    CeedVectorRestoreArrayRead(x_nodes, &x_array);
+    CeedVectorSetArray(u, CEED_MEM_HOST, CEED_COPY_VALUES, (CeedScalar *)&u_array);
+  }
+
+  // Interpolate to arbitrary points
+  CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, p, q, CEED_GAUSS, &basis_u);
+  {
+    CeedScalar x_array[4] = {-0.33, -0.65, 0.16, 0.99};
+
+    CeedVectorSetArray(x_points, CEED_MEM_HOST, CEED_COPY_VALUES, x_array);
+  }
+  CeedBasisApplyAtPoints(basis_u, num_points, CEED_NOTRANSPOSE, CEED_EVAL_GRAD, x_points, u, v);
+
+  {
+    const CeedScalar *x_array, *v_array;
+
+    CeedVectorGetArrayRead(x_points, CEED_MEM_HOST, &x_array);
+    CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array);
+    for (CeedInt i = 0; i < num_points; i++) {
+      CeedScalar dfx = EvalGrad(x_array[i], ALEN(c), c);
+      if (fabs(v_array[i] - dfx) > 500. * CEED_EPSILON) printf("%f != %f = df(%f)\n", v_array[i], dfx, x_array[i]);
+    }
+    CeedVectorRestoreArrayRead(x_points, &x_array);
+    CeedVectorRestoreArrayRead(v, &v_array);
+  }
+
+  CeedVectorDestroy(&x);
+  CeedVectorDestroy(&x_nodes);
+  CeedVectorDestroy(&x_points);
+  CeedVectorDestroy(&u);
+  CeedVectorDestroy(&v);
+  CeedBasisDestroy(&basis_x);
+  CeedBasisDestroy(&basis_u);
+  CeedDestroy(&ceed);
+  return 0;
+}