at-points: change at-points operators to use Chebyshev polynomials of the first kind (normalized)

Author: zatkins-dev

CHANGELOG.md

@@ -30,7 +30,7 @@ Specifically, directories set with `CeedAddJitSourceRoot(ceed, "foo/bar")` will
 - Added support to code generation backends `/gpu/cuda/gen` and `/gpu/hip/gen` for operators with both tensor and non-tensor bases.
 - Add `CeedGetGitVersion()` to access the Git commit and dirty state of the repository at build time.
 - Add `CeedGetBuildConfiguration()` to access compilers, flags, and related information about the build environment.
-- Add support for full `CeedOperator` assembly for operators with multiple active fields with different bases for CPU backends and `/gpu/cuda/ref` and `/gpu/hip/gen` backends. 
+- Add support for full `CeedOperator` assembly for operators with multiple active fields with different bases for CPU backends and `/gpu/cuda/ref` and `/gpu/hip/gen` backends.
 
 ### Examples
 

include/ceed/jit-source/cuda/cuda-ref-basis-tensor-at-points.h

@@ -15,7 +15,7 @@
 template <int Q_1D>
 inline __device__ void ChebyshevPolynomialsAtPoint(const CeedScalar x, CeedScalar *chebyshev_x) {
   chebyshev_x[0] = 1.0;
-  chebyshev_x[1] = 2 * x;
+  chebyshev_x[1] = x;
   for (CeedInt i = 2; i < Q_1D; i++) chebyshev_x[i] = 2 * x * chebyshev_x[i - 1] - chebyshev_x[i - 2];
 }
 
@@ -26,10 +26,11 @@ inline __device__ void ChebyshevDerivativeAtPoint(const CeedScalar x, CeedScalar
   chebyshev_x[1]  = 1.0;
   chebyshev_x[2]  = 2 * x;
   chebyshev_dx[0] = 0.0;
-  chebyshev_dx[1] = 2.0;
+  chebyshev_dx[1] = 1.0;
   for (CeedInt i = 2; i < Q_1D; i++) {
+    // dT_i/dx = i * dU_{i-1}/dx
+    chebyshev_dx[i]          = i * chebyshev_x[i % 3];
     chebyshev_x[(i + 1) % 3] = 2 * x * chebyshev_x[(i + 0) % 3] - chebyshev_x[(i + 2) % 3];
-    chebyshev_dx[i]          = 2 * x * chebyshev_dx[i - 1] + 2 * chebyshev_x[(i + 0) % 3] - chebyshev_dx[i - 2];
   }
 }
 

include/ceed/jit-source/cuda/cuda-shared-basis-tensor-at-points-templates.h

@@ -15,7 +15,7 @@
 template <int Q_1D>
 inline __device__ void ChebyshevPolynomialsAtPoint(const CeedScalar x, CeedScalar *chebyshev_x) {
   chebyshev_x[0] = 1.0;
-  chebyshev_x[1] = 2 * x;
+  chebyshev_x[1] = x;
   for (CeedInt i = 2; i < Q_1D; i++) chebyshev_x[i] = 2 * x * chebyshev_x[i - 1] - chebyshev_x[i - 2];
 }
 
@@ -26,10 +26,11 @@ inline __device__ void ChebyshevDerivativeAtPoint(const CeedScalar x, CeedScalar
   chebyshev_x[1]  = 1.0;
   chebyshev_x[2]  = 2 * x;
   chebyshev_dx[0] = 0.0;
-  chebyshev_dx[1] = 2.0;
+  chebyshev_dx[1] = 1.0;
   for (CeedInt i = 2; i < Q_1D; i++) {
+    // dT_i/dx = i * dU_{i-1}/dx
+    chebyshev_dx[i]          = i * chebyshev_x[i % 3];
     chebyshev_x[(i + 1) % 3] = 2 * x * chebyshev_x[(i + 0) % 3] - chebyshev_x[(i + 2) % 3];
-    chebyshev_dx[i]          = 2 * x * chebyshev_dx[i - 1] + 2 * chebyshev_x[(i + 0) % 3] - chebyshev_dx[i - 2];
   }
 }
 

interface/ceed-basis.c

@@ -47,7 +47,7 @@ const CeedBasis CEED_BASIS_NONE = &ceed_basis_none;
 **/
 static int CeedChebyshevPolynomialsAtPoint(CeedScalar x, CeedInt n, CeedScalar *chebyshev_x) {
   chebyshev_x[0] = 1.0;
-  chebyshev_x[1] = 2 * x;
+  chebyshev_x[1] = 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;
 }
@@ -69,12 +69,13 @@ static int CeedChebyshevDerivativeAtPoint(CeedScalar x, CeedInt n, CeedScalar *c
   chebyshev_x[1]  = 1.0;
   chebyshev_x[2]  = 2 * x;
   chebyshev_dx[0] = 0.0;
-  chebyshev_dx[1] = 2.0;
+  chebyshev_dx[1] = 1.0;
   for (CeedInt i = 2; i < n; i++) {
+    // dT_i/dx = i * dU_{i-1}/dx
     chebyshev_x[0]  = chebyshev_x[1];
     chebyshev_x[1]  = chebyshev_x[2];
+    chebyshev_dx[i] = i * chebyshev_x[1];
     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;
 }
@@ -607,7 +608,6 @@ static int CeedBasisApplyAtPoints_Core(CeedBasis basis, bool apply_add, CeedInt
       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;