1 /// @file 2 /// Test grad in multiple dimensions 3 /// \test Test grad in multiple dimensions 4 #include <ceed.h> 5 #include <math.h> 6 7 static CeedScalar Eval(CeedInt dim, const CeedScalar x[]) { 8 CeedScalar result = tanh(x[0] + 0.1); 9 if (dim > 1) result += atan(x[1] + 0.2); 10 if (dim > 2) result += exp(-(x[2] + 0.3) * (x[2] + 0.3)); 11 return result; 12 } 13 14 static CeedScalar GetTolerance(CeedScalarType scalar_type, int dim) { 15 CeedScalar tol; 16 if (scalar_type == CEED_SCALAR_FP32) { 17 if (dim == 3) tol = 0.05; 18 else tol = 1.e-3; 19 } else { 20 tol = 1.e-11; 21 } 22 return tol; 23 } 24 25 int main(int argc, char **argv) { 26 Ceed ceed; 27 28 CeedInit(argv[1], &ceed); 29 30 for (CeedInt dim = 1; dim <= 3; dim++) { 31 CeedVector X, X_q, U, U_q, ones, grad_T_ones; 32 CeedBasis basis_x_lobatto, basis_u_gauss; 33 CeedInt P = 8, Q = 10, P_dim = CeedIntPow(P, dim), Q_dim = CeedIntPow(Q, dim), X_dim = CeedIntPow(2, dim); 34 CeedScalar x[X_dim * dim], u[P_dim]; 35 const CeedScalar *x_q, *u_q, *grad_t_ones_array; 36 CeedScalar sum_1 = 0, sum_2 = 0; 37 38 for (CeedInt d = 0; d < dim; d++) { 39 for (CeedInt i = 0; i < X_dim; i++) x[d * X_dim + i] = (i % CeedIntPow(2, dim - d)) / CeedIntPow(2, dim - d - 1) ? 1 : -1; 40 } 41 42 CeedVectorCreate(ceed, X_dim * dim, &X); 43 CeedVectorSetArray(X, CEED_MEM_HOST, CEED_USE_POINTER, x); 44 CeedVectorCreate(ceed, P_dim * dim, &X_q); 45 CeedVectorSetValue(X_q, 0); 46 CeedVectorCreate(ceed, P_dim, &U); 47 CeedVectorCreate(ceed, Q_dim * dim, &U_q); 48 CeedVectorSetValue(U_q, 0); 49 CeedVectorCreate(ceed, Q_dim * dim, &ones); 50 CeedVectorSetValue(ones, 1); 51 CeedVectorCreate(ceed, P_dim, &grad_T_ones); 52 CeedVectorSetValue(grad_T_ones, 0); 53 54 // Get function values at quadrature points 55 CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, 2, P, CEED_GAUSS_LOBATTO, &basis_x_lobatto); 56 CeedBasisApply(basis_x_lobatto, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, X, X_q); 57 58 CeedVectorGetArrayRead(X_q, CEED_MEM_HOST, &x_q); 59 for (CeedInt i = 0; i < P_dim; i++) { 60 CeedScalar xx[dim]; 61 for (CeedInt d = 0; d < dim; d++) xx[d] = x_q[d * P_dim + i]; 62 u[i] = Eval(dim, xx); 63 } 64 CeedVectorRestoreArrayRead(X_q, &x_q); 65 CeedVectorSetArray(U, CEED_MEM_HOST, CEED_USE_POINTER, u); 66 67 // Calculate G u at quadrature points, G' * 1 at dofs 68 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, P, Q, CEED_GAUSS, &basis_u_gauss); 69 CeedBasisApply(basis_u_gauss, 1, CEED_NOTRANSPOSE, CEED_EVAL_GRAD, U, U_q); 70 CeedBasisApply(basis_u_gauss, 1, CEED_TRANSPOSE, CEED_EVAL_GRAD, ones, grad_T_ones); 71 72 // Check if 1' * G * u = u' * (G' * 1) 73 CeedVectorGetArrayRead(grad_T_ones, CEED_MEM_HOST, &grad_t_ones_array); 74 CeedVectorGetArrayRead(U_q, CEED_MEM_HOST, &u_q); 75 for (CeedInt i = 0; i < P_dim; i++) sum_1 += grad_t_ones_array[i] * u[i]; 76 for (CeedInt i = 0; i < dim * Q_dim; i++) sum_2 += u_q[i]; 77 CeedVectorRestoreArrayRead(grad_T_ones, &grad_t_ones_array); 78 CeedVectorRestoreArrayRead(U_q, &u_q); 79 CeedScalar tol = GetTolerance(CEED_SCALAR_TYPE, dim); 80 if (fabs(sum_1 - sum_2) > tol) printf("[%" CeedInt_FMT "] %f != %f\n", dim, sum_1, sum_2); 81 82 CeedVectorDestroy(&X); 83 CeedVectorDestroy(&X_q); 84 CeedVectorDestroy(&U); 85 CeedVectorDestroy(&U_q); 86 CeedVectorDestroy(&ones); 87 CeedVectorDestroy(&grad_T_ones); 88 CeedBasisDestroy(&basis_x_lobatto); 89 CeedBasisDestroy(&basis_u_gauss); 90 } 91 CeedDestroy(&ceed); 92 return 0; 93 } 94