From 9aa6d808f7e75d10548254ac7f7619382bb1b20a Mon Sep 17 00:00:00 2001 From: Jeremy L Thompson Date: Tue, 6 Aug 2024 13:21:29 -0600 Subject: [PATCH] basis - add CeedBasisApplyAdd + CPU impl --- backends/ref/ceed-ref-basis.c | 47 ++++--- include/ceed-impl.h | 1 + include/ceed/ceed.h | 1 + .../jit-source/cuda/cuda-ref-basis-tensor.h | 12 +- interface/ceed-basis.c | 62 +++++++++- interface/ceed.c | 1 + tests/README.md | 3 +- tests/t360-basis.c | 56 +++++++++ tests/t361-basis.c | 116 ++++++++++++++++++ tests/t362-basis.c | 59 +++++++++ tests/t363-basis.c | 54 ++++++++ 11 files changed, 382 insertions(+), 30 deletions(-) create mode 100644 tests/t360-basis.c create mode 100644 tests/t361-basis.c create mode 100644 tests/t362-basis.c create mode 100644 tests/t363-basis.c diff --git a/backends/ref/ceed-ref-basis.c b/backends/ref/ceed-ref-basis.c index b82e8bb278..dd395eae1e 100644 --- a/backends/ref/ceed-ref-basis.c +++ b/backends/ref/ceed-ref-basis.c @@ -16,11 +16,11 @@ //------------------------------------------------------------------------------ // Basis Apply //------------------------------------------------------------------------------ -static int CeedBasisApply_Ref(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector U, CeedVector V) { +static int CeedBasisApplyCore_Ref(CeedBasis basis, bool apply_add, CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector U, + CeedVector V) { Ceed ceed; - bool is_tensor_basis; + bool is_tensor_basis, add = apply_add || (t_mode == CEED_TRANSPOSE); CeedInt dim, num_comp, q_comp, num_nodes, num_qpts; - const CeedInt add = (t_mode == CEED_TRANSPOSE); const CeedScalar *u; CeedScalar *v; CeedTensorContract contract; @@ -36,14 +36,9 @@ static int CeedBasisApply_Ref(CeedBasis basis, CeedInt num_elem, CeedTransposeMo CeedCallBackend(CeedBasisGetTensorContract(basis, &contract)); if (U != CEED_VECTOR_NONE) CeedCallBackend(CeedVectorGetArrayRead(U, CEED_MEM_HOST, &u)); else CeedCheck(eval_mode == CEED_EVAL_WEIGHT, ceed, CEED_ERROR_BACKEND, "An input vector is required for this CeedEvalMode"); - CeedCallBackend(CeedVectorGetArrayWrite(V, CEED_MEM_HOST, &v)); - // Clear v if operating in transpose - if (t_mode == CEED_TRANSPOSE) { - const CeedInt v_size = num_elem * num_comp * num_nodes; - - for (CeedInt i = 0; i < v_size; i++) v[i] = (CeedScalar)0.0; - } + if (t_mode == CEED_TRANSPOSE && !apply_add) CeedCallBackend(CeedVectorSetValue(V, 0.0)); + CeedCallBackend(CeedVectorGetArrayWrite(V, CEED_MEM_HOST, &v)); CeedCallBackend(CeedBasisIsTensor(basis, &is_tensor_basis)); if (is_tensor_basis) { @@ -101,12 +96,13 @@ static int CeedBasisApply_Ref(CeedBasis basis, CeedInt num_elem, CeedTransposeMo // or Grad to quadrature points (Transpose) for (CeedInt d = 0; d < dim; d++) { CeedCallBackend(CeedTensorContractApply(contract, pre, P, post, Q, (t_mode == CEED_NOTRANSPOSE ? interp_1d : impl->collo_grad_1d), t_mode, - add && (d > 0), - (t_mode == CEED_NOTRANSPOSE ? (d == 0 ? u : tmp[d % 2]) : u + d * num_qpts * num_comp * num_elem), + (t_mode == CEED_TRANSPOSE) && (d > 0), + (t_mode == CEED_NOTRANSPOSE ? (d == 0 ? u : tmp[d % 2]) : &u[d * num_qpts * num_comp * num_elem]), (t_mode == CEED_NOTRANSPOSE ? (d == dim - 1 ? interp : tmp[(d + 1) % 2]) : interp))); pre /= P; post *= Q; } + // Grad to quadrature points (NoTranspose) // or Interpolate to nodes (Transpose) P = Q_1d, Q = Q_1d; @@ -117,9 +113,10 @@ static int CeedBasisApply_Ref(CeedBasis basis, CeedInt num_elem, CeedTransposeMo pre = num_comp * CeedIntPow(P, dim - 1), post = num_elem; for (CeedInt d = 0; d < dim; d++) { CeedCallBackend(CeedTensorContractApply( - contract, pre, P, post, Q, (t_mode == CEED_NOTRANSPOSE ? impl->collo_grad_1d : interp_1d), t_mode, add && (d == dim - 1), + contract, pre, P, post, Q, (t_mode == CEED_NOTRANSPOSE ? impl->collo_grad_1d : interp_1d), t_mode, + (t_mode == CEED_NOTRANSPOSE && apply_add) || (t_mode == CEED_TRANSPOSE && (d == dim - 1)), (t_mode == CEED_NOTRANSPOSE ? interp : (d == 0 ? interp : tmp[d % 2])), - (t_mode == CEED_NOTRANSPOSE ? v + d * num_qpts * num_comp * num_elem : (d == dim - 1 ? v : tmp[(d + 1) % 2])))); + (t_mode == CEED_NOTRANSPOSE ? &v[d * num_qpts * num_comp * num_elem] : (d == dim - 1 ? v : tmp[(d + 1) % 2])))); pre /= P; post *= Q; } @@ -133,8 +130,8 @@ static int CeedBasisApply_Ref(CeedBasis basis, CeedInt num_elem, CeedTransposeMo for (CeedInt d = 0; d < dim; d++) { CeedCallBackend(CeedTensorContractApply(contract, pre, P, post, Q, grad_1d, t_mode, add && (d > 0), - t_mode == CEED_NOTRANSPOSE ? u : u + d * num_comp * num_qpts * num_elem, - t_mode == CEED_TRANSPOSE ? v : v + d * num_comp * num_qpts * num_elem)); + t_mode == CEED_NOTRANSPOSE ? u : &u[d * num_comp * num_qpts * num_elem], + t_mode == CEED_TRANSPOSE ? v : &v[d * num_comp * num_qpts * num_elem])); pre /= P; post *= Q; } @@ -156,8 +153,8 @@ static int CeedBasisApply_Ref(CeedBasis basis, CeedInt num_elem, CeedTransposeMo for (CeedInt d = 0; d < dim; d++) { CeedCallBackend(CeedTensorContractApply( contract, pre, P, post, Q, (p == d) ? grad_1d : interp_1d, t_mode, add && (d == dim - 1), - (d == 0 ? (t_mode == CEED_NOTRANSPOSE ? u : u + p * num_comp * num_qpts * num_elem) : tmp[d % 2]), - (d == dim - 1 ? (t_mode == CEED_TRANSPOSE ? v : v + p * num_comp * num_qpts * num_elem) : tmp[(d + 1) % 2]))); + (d == 0 ? (t_mode == CEED_NOTRANSPOSE ? u : &u[p * num_comp * num_qpts * num_elem]) : tmp[d % 2]), + (d == dim - 1 ? (t_mode == CEED_TRANSPOSE ? v : &v[p * num_comp * num_qpts * num_elem]) : tmp[(d + 1) % 2]))); pre /= P; post *= Q; } @@ -249,6 +246,16 @@ static int CeedBasisApply_Ref(CeedBasis basis, CeedInt num_elem, CeedTransposeMo return CEED_ERROR_SUCCESS; } +static int CeedBasisApply_Ref(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector U, CeedVector V) { + CeedCallBackend(CeedBasisApplyCore_Ref(basis, false, num_elem, t_mode, eval_mode, U, V)); + return CEED_ERROR_SUCCESS; +} + +static int CeedBasisApplyAdd_Ref(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector U, CeedVector V) { + CeedCallBackend(CeedBasisApplyCore_Ref(basis, true, num_elem, t_mode, eval_mode, U, V)); + return CEED_ERROR_SUCCESS; +} + //------------------------------------------------------------------------------ // Basis Destroy Tensor //------------------------------------------------------------------------------ @@ -297,6 +304,7 @@ int CeedBasisCreateTensorH1_Ref(CeedInt dim, CeedInt P_1d, CeedInt Q_1d, const C CeedCallBackend(CeedBasisSetTensorContract(basis, contract)); CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Apply", CeedBasisApply_Ref)); + CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "ApplyAdd", CeedBasisApplyAdd_Ref)); CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Destroy", CeedBasisDestroyTensor_Ref)); return CEED_ERROR_SUCCESS; } @@ -316,6 +324,7 @@ int CeedBasisCreateH1_Ref(CeedElemTopology topo, CeedInt dim, CeedInt num_nodes, CeedCallBackend(CeedBasisSetTensorContract(basis, contract)); CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Apply", CeedBasisApply_Ref)); + CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "ApplyAdd", CeedBasisApplyAdd_Ref)); return CEED_ERROR_SUCCESS; } @@ -334,6 +343,7 @@ int CeedBasisCreateHdiv_Ref(CeedElemTopology topo, CeedInt dim, CeedInt num_node CeedCallBackend(CeedBasisSetTensorContract(basis, contract)); CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Apply", CeedBasisApply_Ref)); + CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "ApplyAdd", CeedBasisApplyAdd_Ref)); return CEED_ERROR_SUCCESS; } @@ -352,6 +362,7 @@ int CeedBasisCreateHcurl_Ref(CeedElemTopology topo, CeedInt dim, CeedInt num_nod CeedCallBackend(CeedBasisSetTensorContract(basis, contract)); CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Apply", CeedBasisApply_Ref)); + CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "ApplyAdd", CeedBasisApplyAdd_Ref)); return CEED_ERROR_SUCCESS; } diff --git a/include/ceed-impl.h b/include/ceed-impl.h index de206e9e59..889d8a5d31 100644 --- a/include/ceed-impl.h +++ b/include/ceed-impl.h @@ -186,6 +186,7 @@ struct CeedElemRestriction_private { struct CeedBasis_private { Ceed ceed; int (*Apply)(CeedBasis, CeedInt, CeedTransposeMode, CeedEvalMode, CeedVector, CeedVector); + int (*ApplyAdd)(CeedBasis, CeedInt, CeedTransposeMode, CeedEvalMode, CeedVector, CeedVector); int (*ApplyAtPoints)(CeedBasis, CeedInt, const CeedInt *, CeedTransposeMode, CeedEvalMode, CeedVector, CeedVector, CeedVector); int (*Destroy)(CeedBasis); int ref_count; diff --git a/include/ceed/ceed.h b/include/ceed/ceed.h index f55aa9c6f3..fb0d37a35d 100644 --- a/include/ceed/ceed.h +++ b/include/ceed/ceed.h @@ -305,6 +305,7 @@ CEED_EXTERN int CeedBasisCreateProjection(CeedBasis basis_from, CeedBasis basis_ CEED_EXTERN int CeedBasisReferenceCopy(CeedBasis basis, CeedBasis *basis_copy); CEED_EXTERN int CeedBasisView(CeedBasis basis, FILE *stream); CEED_EXTERN int CeedBasisApply(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector u, CeedVector v); +CEED_EXTERN int CeedBasisApplyAdd(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector u, CeedVector v); CEED_EXTERN int CeedBasisApplyAtPoints(CeedBasis basis, CeedInt num_elem, const CeedInt *num_points, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector x_ref, CeedVector u, CeedVector v); CEED_EXTERN int CeedBasisGetCeed(CeedBasis basis, Ceed *ceed); diff --git a/include/ceed/jit-source/cuda/cuda-ref-basis-tensor.h b/include/ceed/jit-source/cuda/cuda-ref-basis-tensor.h index 468ec978c0..4709d86760 100644 --- a/include/ceed/jit-source/cuda/cuda-ref-basis-tensor.h +++ b/include/ceed/jit-source/cuda/cuda-ref-basis-tensor.h @@ -57,13 +57,13 @@ extern "C" __global__ void Interp(const CeedInt num_elem, const CeedInt is_trans // Contract along middle index for (CeedInt k = i; k < writeLen; k += blockDim.x) { - const CeedInt c = k % post; - const CeedInt j = (k / post) % Q; - const CeedInt a = k / (post * Q); - CeedScalar vk = 0; + const CeedInt c = k % post; + const CeedInt j = (k / post) % Q; + const CeedInt a = k / (post * Q); + CeedScalar v_k = 0; - for (CeedInt b = 0; b < P; b++) vk += s_interp_1d[j * stride_0 + b * stride_1] * in[(a * P + b) * post + c]; - out[k] = vk; + for (CeedInt b = 0; b < P; b++) v_k += s_interp_1d[j * stride_0 + b * stride_1] * in[(a * P + b) * post + c]; + out[k] = v_k; } post *= Q; } diff --git a/interface/ceed-basis.c b/interface/ceed-basis.c index bc16186459..52db44d6ef 100644 --- a/interface/ceed-basis.c +++ b/interface/ceed-basis.c @@ -1486,7 +1486,7 @@ int CeedBasisView(CeedBasis basis, FILE *stream) { } /** - @brief Apply basis evaluation from nodes to quadrature points or vice versa + @brief Check input vector dimensions for CeedBasisApply[Add] @param[in] basis `CeedBasis` to evaluate @param[in] num_elem The number of elements to apply the basis evaluation to; @@ -1504,9 +1504,9 @@ int CeedBasisView(CeedBasis basis, FILE *stream) { @return An error code: 0 - success, otherwise - failure - @ref User + @ref Developer **/ -int CeedBasisApply(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector u, CeedVector v) { +static int CeedBasisApplyCheckDims(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector u, CeedVector v) { CeedInt dim, num_comp, q_comp, num_nodes, num_qpts; CeedSize u_length = 0, v_length; Ceed ceed; @@ -1520,8 +1520,6 @@ int CeedBasisApply(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode, CeedCall(CeedVectorGetLength(v, &v_length)); if (u) CeedCall(CeedVectorGetLength(u, &u_length)); - CeedCheck(basis->Apply, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not support CeedBasisApply"); - // Check compatibility of topological and geometrical dimensions CeedCheck((t_mode == CEED_TRANSPOSE && v_length % num_nodes == 0 && u_length % num_qpts == 0) || (t_mode == CEED_NOTRANSPOSE && u_length % num_nodes == 0 && v_length % num_qpts == 0), @@ -1544,11 +1542,65 @@ int CeedBasisApply(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode, break; } CeedCheck(has_good_dims, ceed, CEED_ERROR_DIMENSION, "Input/output vectors too short for basis and evaluation mode"); + return CEED_ERROR_SUCCESS; +} +/** + @brief Apply basis evaluation from nodes to quadrature points or vice versa + + @param[in] basis `CeedBasis` to evaluate + @param[in] num_elem The number of elements to apply the basis evaluation to; + the backend will specify the ordering in @ref CeedElemRestrictionCreate() + @param[in] t_mode @ref CEED_NOTRANSPOSE to evaluate from nodes to quadrature points; + @ref CEED_TRANSPOSE to apply the transpose, mapping from quadrature points to nodes + @param[in] eval_mode @ref CEED_EVAL_NONE to use values directly, + @ref CEED_EVAL_INTERP to use interpolated values, + @ref CEED_EVAL_GRAD to use gradients, + @ref CEED_EVAL_DIV to use divergence, + @ref CEED_EVAL_CURL to use curl, + @ref CEED_EVAL_WEIGHT to use quadrature weights + @param[in] u Input `CeedVector` + @param[out] v Output `CeedVector` + + @return An error code: 0 - success, otherwise - failure + + @ref User +**/ +int CeedBasisApply(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector u, CeedVector v) { + CeedCall(CeedBasisApplyCheckDims(basis, num_elem, t_mode, eval_mode, u, v)); + CeedCheck(basis->Apply, CeedBasisReturnCeed(basis), CEED_ERROR_UNSUPPORTED, "Backend does not support CeedBasisApply"); CeedCall(basis->Apply(basis, num_elem, t_mode, eval_mode, u, v)); return CEED_ERROR_SUCCESS; } +/** + @brief Apply basis evaluation from quadrature points to nodes and sum into target vector + + @param[in] basis `CeedBasis` to evaluate + @param[in] num_elem The number of elements to apply the basis evaluation to; + the backend will specify the ordering in @ref CeedElemRestrictionCreate() + @param[in] t_mode @ref CEED_TRANSPOSE to apply the transpose, mapping from quadrature points to nodes + @param[in] eval_mode @ref CEED_EVAL_NONE to use values directly, + @ref CEED_EVAL_INTERP to use interpolated values, + @ref CEED_EVAL_GRAD to use gradients, + @ref CEED_EVAL_DIV to use divergence, + @ref CEED_EVAL_CURL to use curl, + @ref CEED_EVAL_WEIGHT to use quadrature weights + @param[in] u Input `CeedVector` + @param[out] v Output `CeedVector` to sum into + + @return An error code: 0 - success, otherwise - failure + + @ref User +**/ +int CeedBasisApplyAdd(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector u, CeedVector v) { + CeedCheck(t_mode == CEED_TRANSPOSE, CeedBasisReturnCeed(basis), CEED_ERROR_UNSUPPORTED, "CeedBasisApplyAdd only supports CEED_TRANSPOSE"); + CeedCall(CeedBasisApplyCheckDims(basis, num_elem, t_mode, eval_mode, u, v)); + CeedCheck(basis->ApplyAdd, CeedBasisReturnCeed(basis), CEED_ERROR_UNSUPPORTED, "Backend does not implement CeedBasisApplyAdd"); + CeedCall(basis->ApplyAdd(basis, num_elem, t_mode, eval_mode, u, v)); + return CEED_ERROR_SUCCESS; +} + /** @brief Apply basis evaluation from nodes to arbitrary points diff --git a/interface/ceed.c b/interface/ceed.c index 3c79c59b33..3a9be0c5dc 100644 --- a/interface/ceed.c +++ b/interface/ceed.c @@ -952,6 +952,7 @@ int CeedInit(const char *resource, Ceed *ceed) { CEED_FTABLE_ENTRY(CeedElemRestriction, GetAtPointsElementOffset), CEED_FTABLE_ENTRY(CeedElemRestriction, Destroy), CEED_FTABLE_ENTRY(CeedBasis, Apply), + CEED_FTABLE_ENTRY(CeedBasis, ApplyAdd), CEED_FTABLE_ENTRY(CeedBasis, ApplyAtPoints), CEED_FTABLE_ENTRY(CeedBasis, Destroy), CEED_FTABLE_ENTRY(CeedTensorContract, Apply), diff --git a/tests/README.md b/tests/README.md index fd6e426420..031ff5a030 100644 --- a/tests/README.md +++ b/tests/README.md @@ -15,7 +15,8 @@ The tests are organized by API object, and some tests are further organized, as     2. CeedBasis simplex basis tests\     3. CeedBasis non-tensor H(div) basis tests\     4. CeedBasis non-tensor H(curl) basis tests\ -    5. CeedBasis evaluation at arbitrary points tests +    5. CeedBasis evaluation at arbitrary points tests\ + 6. CeedBasis ApplyAdd tests 4. CeedQFunction Tests\     0. CeedQFunction user code tests\     1. CeedQFunction gallery code tests diff --git a/tests/t360-basis.c b/tests/t360-basis.c new file mode 100644 index 0000000000..987550cc63 --- /dev/null +++ b/tests/t360-basis.c @@ -0,0 +1,56 @@ +/// @file +/// Test interpolation ApplyAdd in multiple dimensions +/// \test Test interpolation ApplyAdd in multiple dimensions +#include +#include +#include + +int main(int argc, char **argv) { + Ceed ceed; + + CeedInit(argv[1], &ceed); + + for (CeedInt dim = 1; dim <= 3; dim++) { + CeedVector u, u_q, v, v_q, w_q; + CeedBasis basis; + CeedInt p = 4, q = 5, p_dim = CeedIntPow(p, dim), q_dim = CeedIntPow(q, dim); + + CeedVectorCreate(ceed, p_dim, &u); + CeedVectorCreate(ceed, p_dim, &v); + CeedVectorSetValue(u, 1.0); + CeedVectorSetValue(v, 0.0); + CeedVectorCreate(ceed, q_dim, &u_q); + CeedVectorCreate(ceed, q_dim, &v_q); + CeedVectorCreate(ceed, q_dim, &w_q); + + CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, p, q, CEED_GAUSS, &basis); + + // Compute area + CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, u, u_q); + CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, CEED_VECTOR_NONE, w_q); + CeedVectorPointwiseMult(v_q, u_q, w_q); + CeedBasisApply(basis, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, v_q, v); + // Double area computed + CeedBasisApplyAdd(basis, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, v_q, v); + + // Check area + { + const CeedScalar *v_array; + CeedScalar area = 0.0; + + CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array); + for (CeedInt i = 0; i < p_dim; i++) area += v_array[i]; + if (fabs(area - 2.0 * CeedIntPow(2, dim)) > 1E-6) printf("Incorrect area computed %f != %f\n", area, 2.0 * CeedIntPow(2, dim)); + CeedVectorRestoreArrayRead(v, &v_array); + } + + CeedVectorDestroy(&u); + CeedVectorDestroy(&v); + CeedVectorDestroy(&u_q); + CeedVectorDestroy(&v_q); + CeedVectorDestroy(&w_q); + CeedBasisDestroy(&basis); + } + CeedDestroy(&ceed); + return 0; +} diff --git a/tests/t361-basis.c b/tests/t361-basis.c new file mode 100644 index 0000000000..6671a39ae5 --- /dev/null +++ b/tests/t361-basis.c @@ -0,0 +1,116 @@ +/// @file +/// Test grad ApplyAdd in multiple dimensions +/// \test Test grad ApplyAdd in multiple dimensions +#include +#include +#include + +static CeedScalar Eval(CeedInt dim, const CeedScalar x[]) { + CeedScalar result = tanh(x[0] + 0.1); + if (dim > 1) result += atan(x[1] + 0.2); + if (dim > 2) result += exp(-(x[2] + 0.3) * (x[2] + 0.3)); + return result; +} + +static CeedScalar GetTolerance(CeedScalarType scalar_type, int dim) { + CeedScalar tol; + if (scalar_type == CEED_SCALAR_FP32) { + if (dim == 3) tol = 0.05; + else tol = 1.e-3; + } else { + tol = 1.e-10; + } + return 2.0 * tol; +} + +int main(int argc, char **argv) { + Ceed ceed; + + CeedInit(argv[1], &ceed); + + for (CeedInt dim = 1; dim <= 3; dim++) { + CeedVector x, x_q, u, u_q, ones, v; + CeedBasis basis_x_lobatto, basis_u_gauss; + CeedInt p = 8, q = 10, p_dim = CeedIntPow(p, dim), q_dim = CeedIntPow(q, dim), x_dim = CeedIntPow(2, dim); + CeedScalar sum_1 = 0, sum_2 = 0; + + CeedVectorCreate(ceed, x_dim * dim, &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, d + 1)) / CeedIntPow(2, d) ? 1 : -1; + } + CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, x_array); + } + CeedVectorCreate(ceed, p_dim * dim, &x_q); + CeedVectorSetValue(x_q, 0); + CeedVectorCreate(ceed, p_dim, &u); + CeedVectorCreate(ceed, q_dim * dim, &u_q); + CeedVectorSetValue(u_q, 0); + CeedVectorCreate(ceed, q_dim * dim, &ones); + CeedVectorSetValue(ones, 1); + CeedVectorCreate(ceed, p_dim, &v); + CeedVectorSetValue(v, 0); + + // Get function values at quadrature points + CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, 2, p, CEED_GAUSS_LOBATTO, &basis_x_lobatto); + CeedBasisApply(basis_x_lobatto, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x, x_q); + + { + const CeedScalar *x_q_array; + CeedScalar u_array[p_dim]; + + CeedVectorGetArrayRead(x_q, CEED_MEM_HOST, &x_q_array); + for (CeedInt i = 0; i < p_dim; i++) { + CeedScalar coord[dim]; + + for (CeedInt d = 0; d < dim; d++) coord[d] = x_q_array[d * p_dim + i]; + u_array[i] = Eval(dim, coord); + } + CeedVectorRestoreArrayRead(x_q, &x_q_array); + CeedVectorSetArray(u, CEED_MEM_HOST, CEED_COPY_VALUES, u_array); + } + + // Calculate G u at quadrature points, G' * 1 at dofs + CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, p, q, CEED_GAUSS, &basis_u_gauss); + CeedBasisApply(basis_u_gauss, 1, CEED_NOTRANSPOSE, CEED_EVAL_GRAD, u, u_q); + CeedVectorScale(u_q, 2.0); + CeedBasisApply(basis_u_gauss, 1, CEED_TRANSPOSE, CEED_EVAL_GRAD, ones, v); + CeedBasisApplyAdd(basis_u_gauss, 1, CEED_TRANSPOSE, CEED_EVAL_GRAD, ones, v); + + // Check if 1' * G * u = u' * (G' * 1) + { + const CeedScalar *v_array, *u_array, *u_q_array; + + CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array); + CeedVectorGetArrayRead(u, CEED_MEM_HOST, &u_array); + CeedVectorGetArrayRead(u_q, CEED_MEM_HOST, &u_q_array); + for (CeedInt i = 0; i < p_dim; i++) sum_1 += v_array[i] * u_array[i]; + for (CeedInt i = 0; i < dim * q_dim; i++) sum_2 += u_q_array[i]; + CeedVectorRestoreArrayRead(v, &v_array); + CeedVectorRestoreArrayRead(u, &u_array); + CeedVectorRestoreArrayRead(u_q, &u_q_array); + } + { + CeedScalarType scalar_type; + + CeedGetScalarType(&scalar_type); + + CeedScalar tol = GetTolerance(scalar_type, dim); + + if (fabs(sum_1 - sum_2) > tol) printf("[%" CeedInt_FMT "] %0.12f != %0.12f\n", dim, sum_1, sum_2); + } + + CeedVectorDestroy(&x); + CeedVectorDestroy(&x_q); + CeedVectorDestroy(&u); + CeedVectorDestroy(&u_q); + CeedVectorDestroy(&ones); + CeedVectorDestroy(&v); + CeedBasisDestroy(&basis_x_lobatto); + CeedBasisDestroy(&basis_u_gauss); + } + CeedDestroy(&ceed); + return 0; +} diff --git a/tests/t362-basis.c b/tests/t362-basis.c new file mode 100644 index 0000000000..bff1937d66 --- /dev/null +++ b/tests/t362-basis.c @@ -0,0 +1,59 @@ +/// @file +/// Test integration ApplyAdd with a 2D Simplex non-tensor H^1 basis +/// \test Test integration ApplyAdd with a 2D Simplex non-tensor H^1 basis +#include +#include +#include + +#include "t320-basis.h" + +// main test +int main(int argc, char **argv) { + Ceed ceed; + CeedVector u, v, u_q, v_q, w_q; + const CeedInt p = 6, q = 4, dim = 2; + CeedBasis basis; + CeedScalar q_ref[dim * q], q_weight[q]; + CeedScalar interp[p * q], grad[dim * p * q]; + + CeedInit(argv[1], &ceed); + + CeedVectorCreate(ceed, p, &u); + CeedVectorCreate(ceed, p, &v); + CeedVectorSetValue(u, 1.0); + CeedVectorSetValue(v, 0.0); + CeedVectorCreate(ceed, q, &u_q); + CeedVectorCreate(ceed, q, &v_q); + CeedVectorCreate(ceed, q, &w_q); + + Build2DSimplex(q_ref, q_weight, interp, grad); + CeedBasisCreateH1(ceed, CEED_TOPOLOGY_TRIANGLE, 1, p, q, interp, grad, q_ref, q_weight, &basis); + + // Compute area + CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, u, u_q); + CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, CEED_VECTOR_NONE, w_q); + CeedVectorPointwiseMult(v_q, u_q, w_q); + CeedBasisApply(basis, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, v_q, v); + // Double area computed + CeedBasisApplyAdd(basis, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, v_q, v); + + // Check area + { + const CeedScalar *v_array; + CeedScalar area = 0.0; + + CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array); + for (CeedInt i = 0; i < p; i++) area += v_array[i]; + if (fabs(area - 1.0) > 1E-6) printf("Incorrect area computed %f != %f\n", area, 1.0); + CeedVectorRestoreArrayRead(v, &v_array); + } + + CeedVectorDestroy(&u); + CeedVectorDestroy(&v); + CeedVectorDestroy(&u_q); + CeedVectorDestroy(&v_q); + CeedVectorDestroy(&w_q); + CeedBasisDestroy(&basis); + CeedDestroy(&ceed); + return 0; +} diff --git a/tests/t363-basis.c b/tests/t363-basis.c new file mode 100644 index 0000000000..6c19f34027 --- /dev/null +++ b/tests/t363-basis.c @@ -0,0 +1,54 @@ +/// @file +/// Test grad transpose ApplyAdd with a 2D Simplex non-tensor H^1 basis +/// \test Test grad transpose ApplyAdd with a 2D Simplex non-tensor H^1 basis +#include +#include +#include + +#include "t320-basis.h" + +int main(int argc, char **argv) { + Ceed ceed; + CeedVector u, v; + const CeedInt p = 6, q = 4, dim = 2; + CeedBasis basis; + CeedScalar q_ref[dim * q], q_weight[q]; + CeedScalar interp[p * q], grad[dim * p * q]; + CeedScalar column_sum[p]; + + CeedInit(argv[1], &ceed); + + CeedVectorCreate(ceed, q * dim, &u); + CeedVectorSetValue(u, 1); + CeedVectorCreate(ceed, p, &v); + CeedVectorSetValue(v, 0); + + Build2DSimplex(q_ref, q_weight, interp, grad); + CeedBasisCreateH1(ceed, CEED_TOPOLOGY_TRIANGLE, 1, p, q, interp, grad, q_ref, q_weight, &basis); + + CeedBasisApply(basis, 1, CEED_TRANSPOSE, CEED_EVAL_GRAD, u, v); + CeedBasisApplyAdd(basis, 1, CEED_TRANSPOSE, CEED_EVAL_GRAD, u, v); + + // Check values at quadrature points + for (int i = 0; i < p; i++) { + column_sum[i] = 0; + for (int j = 0; j < q * dim; j++) { + column_sum[i] += grad[i + j * p]; + } + } + { + const CeedScalar *v_array; + + CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array); + for (int i = 0; i < p; i++) { + if (fabs(column_sum[i] - v_array[i] / 2.0) > 100. * CEED_EPSILON) printf("[%" CeedInt_FMT "] %f != %f\n", i, v_array[i] / 2.0, column_sum[i]); + } + CeedVectorRestoreArrayRead(v, &v_array); + } + + CeedVectorDestroy(&u); + CeedVectorDestroy(&v); + CeedBasisDestroy(&basis); + CeedDestroy(&ceed); + return 0; +}