From 53ef2869b49ed7dd91aa8ced0c36f3ed2fcf2517 Mon Sep 17 00:00:00 2001 From: Zach Atkins Date: Tue, 5 Sep 2023 12:07:35 -0600 Subject: [PATCH] Fix order of tensor contractions for `CeedBasisApplyAtPoints` (#1324) * Fix order of tensor contractions such that x is the inner dimension * Fix typo in docstring (arbirtary -> arbitrary) --- interface/ceed-basis.c | 21 +++++++++------------ tests/t350-basis.c | 2 +- tests/t351-basis.c | 4 ++-- tests/t352-basis.c | 4 ++-- tests/t353-basis.c | 2 +- tests/t354-basis.c | 4 ++-- tests/t355-basis.c | 2 +- tests/t356-basis.c | 4 ++-- 8 files changed, 20 insertions(+), 23 deletions(-) diff --git a/interface/ceed-basis.c b/interface/ceed-basis.c index 015806d99b..d88eddb90d 100644 --- 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 index 4e8b9f5388..54979bb9a0 100644 --- 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 #include diff --git a/tests/t351-basis.c b/tests/t351-basis.c index 955cd4aec2..27fa6a3c8c 100644 --- 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 #include @@ -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 index 3a00356c67..b84c3f510a 100644 --- 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 #include @@ -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 index e44b773acc..22f80ddcdd 100644 --- 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 #include diff --git a/tests/t354-basis.c b/tests/t354-basis.c index fe4c0aa82b..e0b730158b 100644 --- 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 #include @@ -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 index 1d1fe86835..7fd7906dcb 100644 --- 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 #include diff --git a/tests/t356-basis.c b/tests/t356-basis.c index b04584901f..8f7fc0eb49 100644 --- 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 #include @@ -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); }