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