Fix order of tensor contractions for CeedBasisApplyAtPoints (CEED#1324

) * Fix order of tensor contractions such that x is the inner dimension * Fix typo in docstring (arbirtary -> arbitrary)
argonne-lcf · Sep 5, 2023 · 53ef286 · 53ef286
1 parent c6c10e0
commit 53ef286
Show file tree

Hide file tree

Showing 8 changed files with 20 additions and 23 deletions.
diff --git a/interface/ceed-basis.c b/interface/ceed-basis.c
@@ -1609,12 +1609,11 @@ int CeedBasisApplyAtPoints(CeedBasis basis, CeedInt num_points, CeedTransposeMod
           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--) {
+            for (CeedInt d = 0; d < dim; 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]));
+                                               d == 0 ? chebyshev_coeffs : tmp[d % 2], d == (dim - 1) ? &v_array[p * num_comp] : tmp[(d + 1) % 2]));
               pre /= Q_1d;
               post *= 1;
             }
@@ -1626,18 +1625,17 @@ int CeedBasisApplyAtPoints(CeedBasis basis, CeedInt num_points, CeedTransposeMod
 
           // ---- 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--) {
+            for (CeedInt pass = 0; pass < dim; pass++) {
               CeedInt pre = num_comp * CeedIntPow(Q_1d, dim - 1), post = 1;
 
-              for (CeedInt d = dim - 1; d >= 0; d--) {
+              for (CeedInt d = 0; d < dim; 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]));
+                                                 d == 0 ? chebyshev_coeffs : tmp[d % 2],
+                                                 d == (dim - 1) ? &v_array[p * num_comp * dim + pass] : tmp[(d + 1) % 2]));
                 pre /= Q_1d;
                 post *= 1;
               }
@@ -1671,12 +1669,11 @@ 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--) {
+          for (CeedInt d = 0; d < dim; 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, 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]));
+            CeedCall(CeedTensorContractApply(basis->contract, pre, 1, post, Q_1d, chebyshev_x, t_mode, p > 0 && d == (dim - 1),
+                                             d == 0 ? &u_array[p * num_comp] : tmp[d % 2], d == (dim - 1) ? chebyshev_coeffs : tmp[(d + 1) % 2]));
             pre /= 1;
             post *= Q_1d;
           }

diff --git a/tests/t350-basis.c b/tests/t350-basis.c
@@ -1,5 +1,5 @@
 /// @file
-/// Test polynomial interpolation to arbirtary points in 1D
+/// Test polynomial interpolation to arbitrary points in 1D
 /// \test Test polynomial interpolation to arbitrary points in 1D
 #include <ceed.h>
 #include <math.h>

diff --git a/tests/t351-basis.c b/tests/t351-basis.c
@@ -1,5 +1,5 @@
 /// @file
-/// Test polynomial interpolation to arbirtary points in multiple dimensions
+/// Test polynomial interpolation to arbitrary points in multiple dimensions
 /// \test Test polynomial interpolation to arbitrary points in multiple dimensions
 #include <ceed.h>
 #include <math.h>
@@ -36,7 +36,7 @@ int main(int argc, char **argv) {
       CeedScalar x_array[x_dim * dim];
 
       for (CeedInt d = 0; d < dim; d++) {
-        for (CeedInt i = 0; i < x_dim; i++) x_array[d * x_dim + i] = (i % CeedIntPow(2, dim - d)) / CeedIntPow(2, dim - d - 1) ? 1 : -1;
+        for (CeedInt i = 0; i < x_dim; i++) x_array[d * x_dim + i] = (i % CeedIntPow(2, d + 1)) / CeedIntPow(2, d) ? 1 : -1;
       }
       CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, x_array);
     }

diff --git a/tests/t352-basis.c b/tests/t352-basis.c
@@ -1,5 +1,5 @@
 /// @file
-/// Test polynomial interpolation to arbirtary points with multiple components in multiple dimensions
+/// Test polynomial interpolation to arbitrary points with multiple components in multiple dimensions
 /// \test Test polynomial interpolation to arbitrary points with multiple components in multiple dimensions
 #include <ceed.h>
 #include <math.h>
@@ -36,7 +36,7 @@ int main(int argc, char **argv) {
       CeedScalar x_array[x_dim * dim];
 
       for (CeedInt d = 0; d < dim; d++) {
-        for (CeedInt i = 0; i < x_dim; i++) x_array[d * x_dim + i] = (i % CeedIntPow(2, dim - d)) / CeedIntPow(2, dim - d - 1) ? 1 : -1;
+        for (CeedInt i = 0; i < x_dim; i++) x_array[d * x_dim + i] = (i % CeedIntPow(2, d + 1)) / CeedIntPow(2, d) ? 1 : -1;
       }
       CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, x_array);
     }

diff --git a/tests/t353-basis.c b/tests/t353-basis.c
@@ -1,5 +1,5 @@
 /// @file
-/// Test polynomial interpolation transpose from arbirtary points in 1D
+/// Test polynomial interpolation transpose from arbitrary points in 1D
 /// \test Test polynomial interpolation transpose from arbitrary points in 1D
 #include <ceed.h>
 #include <math.h>

diff --git a/tests/t354-basis.c b/tests/t354-basis.c
@@ -1,5 +1,5 @@
 /// @file
-/// Test polynomial interpolation to arbirtary points in multiple dimensions
+/// Test polynomial interpolation to arbitrary points in multiple dimensions
 /// \test Test polynomial interpolation to arbitrary points in multiple dimensions
 #include <ceed.h>
 #include <math.h>
@@ -40,7 +40,7 @@ int main(int argc, char **argv) {
       CeedScalar x_array[x_dim * dim];
 
       for (CeedInt d = 0; d < dim; d++) {
-        for (CeedInt i = 0; i < x_dim; i++) x_array[d * x_dim + i] = (i % CeedIntPow(2, dim - d)) / CeedIntPow(2, dim - d - 1) ? 1 : -1;
+        for (CeedInt i = 0; i < x_dim; i++) x_array[d * x_dim + i] = (i % CeedIntPow(2, d + 1)) / CeedIntPow(2, d) ? 1 : -1;
       }
       CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, x_array);
     }

diff --git a/tests/t355-basis.c b/tests/t355-basis.c
@@ -1,5 +1,5 @@
 /// @file
-/// Test polynomial gradient to arbirtary points in 1D
+/// Test polynomial gradient to arbitrary points in 1D
 /// \test Test polynomial gradient to arbitrary points in 1D
 #include <ceed.h>
 #include <math.h>

diff --git a/tests/t356-basis.c b/tests/t356-basis.c
@@ -1,5 +1,5 @@
 /// @file
-/// Test polynomial gradient to arbirtary points in multiple dimensions
+/// Test polynomial gradient to arbitrary points in multiple dimensions
 /// \test Test polynomial graient to arbitrary points in multiple dimensions
 #include <ceed.h>
 #include <math.h>
@@ -46,7 +46,7 @@ int main(int argc, char **argv) {
       CeedScalar x_array[x_dim * dim];
 
       for (CeedInt d = 0; d < dim; d++) {
-        for (CeedInt i = 0; i < x_dim; i++) x_array[d * x_dim + i] = (i % CeedIntPow(2, dim - d)) / CeedIntPow(2, dim - d - 1) ? 1 : -1;
+        for (CeedInt i = 0; i < x_dim; i++) x_array[d * x_dim + i] = (i % CeedIntPow(2, d + 1)) / CeedIntPow(2, d) ? 1 : -1;
       }
       CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, x_array);
     }