diff --git a/doc/sphinx/source/releasenotes.md b/doc/sphinx/source/releasenotes.md index ab88370d1a..ad6a97e509 100644 --- a/doc/sphinx/source/releasenotes.md +++ b/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. diff --git a/interface/ceed-basis.c b/interface/ceed-basis.c index 2bf3483a18..26891a688b 100644 --- a/interface/ceed-basis.c +++ b/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,29 +1592,53 @@ 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)); @@ -1576,6 +1646,7 @@ int CeedBasisApplyAtPoints(CeedBasis basis, CeedInt num_points, CeedTransposeMod 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; diff --git a/tests/t355-basis.c b/tests/t355-basis.c new file mode 100644 index 0000000000..1d1fe86835 --- /dev/null +++ b/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 +#include +#include + +#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; +} diff --git a/tests/t356-basis.c b/tests/t356-basis.c new file mode 100644 index 0000000000..b04584901f --- /dev/null +++ b/tests/t356-basis.c @@ -0,0 +1,115 @@ +/// @file +/// Test polynomial gradient to arbirtary points in multiple dimensions +/// \test Test polynomial graient to arbitrary points in multiple dimensions +#include +#include +#include + +static CeedScalar Eval(CeedInt dim, const CeedScalar x[]) { + CeedScalar result = 1, center = 0.1; + for (CeedInt d = 0; d < dim; d++) { + result *= tanh(x[d] - center); + center += 0.1; + } + return result; +} + +static CeedScalar EvalGrad(CeedInt dim, const CeedScalar x[], CeedInt direction) { + CeedScalar result = 1, center = 0.1; + for (CeedInt d = 0; d < dim; d++) { + if (d == direction) result *= 1.0 / cosh(x[d] - center) / cosh(x[d] - center); + else result *= tanh(x[d] - center); + center += 0.1; + } + return result; +} + +int main(int argc, char **argv) { + Ceed ceed; + + CeedInit(argv[1], &ceed); + + for (CeedInt dim = 1; dim <= 3; dim++) { + CeedVector x, x_nodes, x_points, u, v; + CeedBasis basis_x, basis_u; + const CeedInt p = 9, q = 9, num_points = 4, x_dim = CeedIntPow(2, dim), p_dim = CeedIntPow(p, dim); + + CeedVectorCreate(ceed, x_dim * dim, &x); + CeedVectorCreate(ceed, p_dim * dim, &x_nodes); + CeedVectorCreate(ceed, num_points * dim, &x_points); + CeedVectorCreate(ceed, p_dim, &u); + CeedVectorCreate(ceed, num_points * dim, &v); + + // Get nodal coordinates + CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, 2, p, CEED_GAUSS_LOBATTO, &basis_x); + { + 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; + } + 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_dim]; + + CeedVectorGetArrayRead(x_nodes, CEED_MEM_HOST, &x_array); + for (CeedInt i = 0; i < p_dim; i++) { + CeedScalar coord[dim]; + + for (CeedInt d = 0; d < dim; d++) coord[d] = x_array[d * p_dim + i]; + u_array[i] = Eval(dim, coord); + } + CeedVectorRestoreArrayRead(x_nodes, &x_array); + CeedVectorSetArray(u, CEED_MEM_HOST, CEED_COPY_VALUES, (CeedScalar *)&u_array); + } + + // Interpolate to arbitrary points + CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, p, q, CEED_GAUSS, &basis_u); + { + CeedScalar x_array[12] = {-0.33, -0.65, 0.16, 0.99, -0.65, 0.16, 0.99, -0.33, 0.16, 0.99, -0.33, -0.65}; + + 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 coord[dim]; + + for (CeedInt d = 0; d < dim; d++) coord[d] = x_array[d + i * dim]; + for (CeedInt d = 0; d < dim; d++) { + const CeedScalar dfx = EvalGrad(dim, coord, d); + if (fabs(v_array[i * dim + d] - dfx) > 1E-3) { + // LCOV_EXCL_START + printf("[%" CeedInt_FMT "] %f != %f = df(%f", dim, v_array[i * dim + d], dfx, coord[0]); + for (CeedInt d = 1; d < dim; d++) printf(", %f", coord[d]); + puts(")"); + // LCOV_EXCL_STOP + } + } + } + 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; +}