1 /// @file
2 /// Test gradient transpose ApplyAdd in multiple dimensions at arbitrary points
3 /// \test Test gradient transpose ApplyAdd in multiple dimensions at arbitrary points
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.005;
19 else tol = 1.e-4;
20 } else {
21 tol = 1.e-11;
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_nodes, x_points, u, u_points, v, ones;
33 CeedBasis basis_x, basis_u;
34 const CeedInt p = 9, q = 9, num_points = 4, x_dim = CeedIntPow(2, dim), p_dim = CeedIntPow(p, dim);
35 CeedScalar sum_1 = 0, sum_2 = 0;
36
37 CeedVectorCreate(ceed, x_dim * dim, &x);
38 CeedVectorCreate(ceed, p_dim * dim, &x_nodes);
39 CeedVectorCreate(ceed, num_points * dim, &x_points);
40 CeedVectorCreate(ceed, p_dim, &u);
41 CeedVectorCreate(ceed, num_points * dim, &u_points);
42 CeedVectorCreate(ceed, p_dim, &v);
43 CeedVectorCreate(ceed, num_points * dim, &ones);
44
45 CeedVectorSetValue(ones, 1);
46 CeedVectorSetValue(v, 0);
47
48 // Get nodal coordinates
49 CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, 2, p, CEED_GAUSS_LOBATTO, &basis_x);
50 {
51 CeedScalar x_array[x_dim * dim];
52
53 for (CeedInt d = 0; d < dim; d++) {
54 for (CeedInt i = 0; i < x_dim; i++) x_array[d * x_dim + i] = (i % CeedIntPow(2, d + 1)) / CeedIntPow(2, d) ? 1 : -1;
55 }
56 CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, x_array);
57 }
58 CeedBasisApply(basis_x, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x, x_nodes);
59
60 // Set values of u at nodes
61 {
62 const CeedScalar *x_array;
63 CeedScalar u_array[p_dim];
64
65 CeedVectorGetArrayRead(x_nodes, CEED_MEM_HOST, &x_array);
66 for (CeedInt i = 0; i < p_dim; i++) {
67 CeedScalar coord[dim];
68
69 for (CeedInt d = 0; d < dim; d++) coord[d] = x_array[d * p_dim + i];
70 u_array[i] = Eval(dim, coord);
71 }
72 CeedVectorRestoreArrayRead(x_nodes, &x_array);
73 CeedVectorSetArray(u, CEED_MEM_HOST, CEED_COPY_VALUES, (CeedScalar *)&u_array);
74 }
75
76 // Interpolate to arbitrary points
77 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, p, q, CEED_GAUSS, &basis_u);
78 {
79 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};
80
81 CeedVectorSetArray(x_points, CEED_MEM_HOST, CEED_COPY_VALUES, x_array);
82 }
83
84 // Calculate G u at arbitrary points, G' * 1 at dofs
85 CeedBasisApplyAtPoints(basis_u, 1, &num_points, CEED_NOTRANSPOSE, CEED_EVAL_GRAD, x_points, u, u_points);
86 CeedBasisApplyAtPoints(basis_u, 1, &num_points, CEED_TRANSPOSE, CEED_EVAL_GRAD, x_points, ones, v);
87 // Double it
88 CeedBasisApplyAddAtPoints(basis_u, 1, &num_points, CEED_TRANSPOSE, CEED_EVAL_GRAD, x_points, ones, v);
89 {
90 const CeedScalar *u_array, *v_array, *u_points_array;
91
92 CeedVectorGetArrayRead(u, CEED_MEM_HOST, &u_array);
93 CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array);
94 CeedVectorGetArrayRead(u_points, CEED_MEM_HOST, &u_points_array);
95 for (CeedInt i = 0; i < p_dim; i++) sum_1 += v_array[i] * u_array[i];
96 for (CeedInt i = 0; i < num_points * dim; i++) sum_2 += u_points_array[i];
97 CeedVectorRestoreArrayRead(u, &u_array);
98 CeedVectorRestoreArrayRead(v, &v_array);
99 CeedVectorRestoreArrayRead(u_points, &u_points_array);
100 }
101 {
102 CeedScalarType scalar_type;
103
104 CeedGetScalarType(&scalar_type);
105
106 CeedScalar tol = GetTolerance(scalar_type, dim);
107
108 if (fabs(sum_1 - 2.0 * sum_2) > tol) printf("[%" CeedInt_FMT "] %f != %f\n", dim, sum_1, 2.0 * sum_2);
109 }
110
111 CeedVectorDestroy(&x);
112 CeedVectorDestroy(&x_nodes);
113 CeedVectorDestroy(&x_points);
114 CeedVectorDestroy(&u);
115 CeedVectorDestroy(&u_points);
116 CeedVectorDestroy(&ones);
117 CeedVectorDestroy(&v);
118 CeedBasisDestroy(&basis_x);
119 CeedBasisDestroy(&basis_u);
120 }
121 CeedDestroy(&ceed);
122 return 0;
123 }
124