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
Eval(CeedInt dim,const CeedScalar x[])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
GetTolerance(CeedScalarType scalar_type,int dim)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
main(int argc,char ** argv)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
68 for (CeedInt d = 0; d < dim; d++) coord[d] = x_q_array[d * p_dim + i];
69 u_array[i] = Eval(dim, coord);
70 }
71 CeedVectorRestoreArrayRead(x_q, &x_q_array);
72 CeedVectorSetArray(u, CEED_MEM_HOST, CEED_COPY_VALUES, u_array);
73 }
74
75 // Calculate G u at quadrature points, G' * 1 at dofs
76 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, p, q, CEED_GAUSS, &basis_u_gauss);
77 CeedBasisApply(basis_u_gauss, 1, CEED_NOTRANSPOSE, CEED_EVAL_GRAD, u, u_q);
78 CeedBasisApply(basis_u_gauss, 1, CEED_TRANSPOSE, CEED_EVAL_GRAD, ones, v);
79
80 // Check if 1' * G * u = u' * (G' * 1)
81 {
82 const CeedScalar *v_array, *u_array, *u_q_array;
83
84 CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array);
85 CeedVectorGetArrayRead(u, CEED_MEM_HOST, &u_array);
86 CeedVectorGetArrayRead(u_q, CEED_MEM_HOST, &u_q_array);
87 for (CeedInt i = 0; i < p_dim; i++) sum_1 += v_array[i] * u_array[i];
88 for (CeedInt i = 0; i < dim * q_dim; i++) sum_2 += u_q_array[i];
89 CeedVectorRestoreArrayRead(v, &v_array);
90 CeedVectorRestoreArrayRead(u, &u_array);
91 CeedVectorRestoreArrayRead(u_q, &u_q_array);
92 }
93 {
94 CeedScalarType scalar_type;
95
96 CeedGetScalarType(&scalar_type);
97
98 CeedScalar tol = GetTolerance(scalar_type, dim);
99
100 if (fabs(sum_1 - sum_2) > tol) printf("[%" CeedInt_FMT "] %0.12f != %0.12f\n", dim, sum_1, sum_2);
101 }
102
103 CeedVectorDestroy(&x);
104 CeedVectorDestroy(&x_q);
105 CeedVectorDestroy(&u);
106 CeedVectorDestroy(&u_q);
107 CeedVectorDestroy(&ones);
108 CeedVectorDestroy(&v);
109 CeedBasisDestroy(&basis_x_lobatto);
110 CeedBasisDestroy(&basis_u_gauss);
111 }
112 CeedDestroy(&ceed);
113 return 0;
114 }
115