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,21 +15,22 @@
 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];
 }
 
 template <int Q_1D>
 inline __device__ void ChebyshevDerivativeAtPoint(const CeedScalar x, CeedScalar *chebyshev_dx) {
-  CeedScalar chebyshev_x[3];
+  CeedScalar chebyshev_x[2];
 
-  chebyshev_x[1]  = 1.0;
-  chebyshev_x[2]  = 2 * x;
+  chebyshev_x[0]  = 1.0;
+  chebyshev_x[1]  = 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++) {
-    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];
+    // dT_i/dx = i * dU_{i-1}/dx
+    chebyshev_dx[i]    = i * chebyshev_x[(i + 1) % 2];
+    chebyshev_x[i % 2] = 2 * x * chebyshev_x[(i + 1) % 2] - chebyshev_x[i % 2];
   }
 }
 

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

@@ -15,21 +15,22 @@
 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];
 }
 
 template <int Q_1D>
 inline __device__ void ChebyshevDerivativeAtPoint(const CeedScalar x, CeedScalar *chebyshev_dx) {
-  CeedScalar chebyshev_x[3];
+  CeedScalar chebyshev_x[2];
 
-  chebyshev_x[1]  = 1.0;
-  chebyshev_x[2]  = 2 * x;
+  chebyshev_x[0]  = 1.0;
+  chebyshev_x[1]  = 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++) {
-    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];
+    // dT_i/dx = i * dU_{i-1}/dx
+    chebyshev_dx[i]    = i * chebyshev_x[(i + 1) % 2];
+    chebyshev_x[i % 2] = 2 * x * chebyshev_x[(i + 1) % 2] - chebyshev_x[i % 2];
   }
 }
 

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

@@ -15,21 +15,22 @@
 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];
 }
 
 template <int Q_1D>
 inline __device__ void ChebyshevDerivativeAtPoint(const CeedScalar x, CeedScalar *chebyshev_dx) {
-  CeedScalar chebyshev_x[3];
+  CeedScalar chebyshev_x[2];
 
-  chebyshev_x[1]  = 1.0;
-  chebyshev_x[2]  = 2 * x;
+  chebyshev_x[0]  = 1.0;
+  chebyshev_x[1]  = 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++) {
-    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];
+    // dT_i/dx = i * dU_{i-1}/dx
+    chebyshev_dx[i]    = i * chebyshev_x[(i + 1) % 2];
+    chebyshev_x[i % 2] = 2 * x * chebyshev_x[(i + 1) % 2] - chebyshev_x[i % 2];
   }
 }
 

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

@@ -15,21 +15,22 @@
 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];
 }
 
 template <int Q_1D>
 inline __device__ void ChebyshevDerivativeAtPoint(const CeedScalar x, CeedScalar *chebyshev_dx) {
-  CeedScalar chebyshev_x[3];
+  CeedScalar chebyshev_x[2];
 
-  chebyshev_x[1]  = 1.0;
-  chebyshev_x[2]  = 2 * x;
+  chebyshev_x[0]  = 1.0;
+  chebyshev_x[1]  = 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++) {
-    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];
+    // dT_i/dx = i * dU_{i-1}/dx
+    chebyshev_dx[i]    = i * chebyshev_x[(i + 1) % 2];
+    chebyshev_x[i % 2] = 2 * x * chebyshev_x[(i + 1) % 2] - chebyshev_x[i % 2];
   }
 }
 

interface/ceed-basis.c

@@ -46,8 +46,10 @@ const CeedBasis CEED_BASIS_NONE = &ceed_basis_none;
   @ref Developer
 **/
 static int CeedChebyshevPolynomialsAtPoint(CeedScalar x, CeedInt n, CeedScalar *chebyshev_x) {
+  // Chebyshev polynomials of the first kind, 3 term recurrence: T_0(x) = 1, T_1(x) = x, T_n = 2 x T_{n-1} - T_{n-2}
+  // For more info, see e.g. https://en.wikipedia.org/wiki/Chebyshev_polynomials#Recurrence_definition
   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;
 }
@@ -64,17 +66,20 @@ static int CeedChebyshevPolynomialsAtPoint(CeedScalar x, CeedInt n, CeedScalar *
   @ref Developer
 **/
 static int CeedChebyshevDerivativeAtPoint(CeedScalar x, CeedInt n, CeedScalar *chebyshev_dx) {
+  // Note, these are Chebyshev polynomials of the 2nd kind, used for derivatives
+  // See e.g. https://en.wikipedia.org/wiki/Chebyshev_polynomials#Differentiation_and_integration for the recurrence
   CeedScalar chebyshev_x[3];
 
   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 +612,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;