1 /// @file
2 /// Test full assembly of Poisson operator (see t534)
3 /// \test Test full assembly of Poisson operator
4 #include <ceed.h>
5 #include <math.h>
6 #include <stdio.h>
7 #include <stdlib.h>
8
9 #include "t534-operator.h"
10
main(int argc,char ** argv)11 int main(int argc, char **argv) {
12 Ceed ceed;
13 CeedElemRestriction elem_restriction_x, elem_restriction_u, elem_restriction_q_data;
14 CeedBasis basis_x, basis_u;
15 CeedQFunction qf_setup, qf_diff;
16 CeedOperator op_setup, op_diff;
17 CeedVector q_data, x, u, v;
18 CeedInt p = 3, q = 4, dim = 2;
19 CeedInt n_x = 3, n_y = 2;
20 CeedInt num_elem = n_x * n_y;
21 CeedInt num_dofs = (n_x * 2 + 1) * (n_y * 2 + 1), num_qpts = num_elem * q * q;
22 CeedInt ind_x[num_elem * p * p];
23 CeedScalar assembled_values[num_dofs * num_dofs];
24 CeedScalar assembled_true[num_dofs * num_dofs];
25
26 CeedInit(argv[1], &ceed);
27
28 // Vectors
29 CeedVectorCreate(ceed, dim * num_dofs, &x);
30 {
31 CeedScalar x_array[dim * num_dofs];
32
33 for (CeedInt i = 0; i < n_x * 2 + 1; i++) {
34 for (CeedInt j = 0; j < n_y * 2 + 1; j++) {
35 x_array[i + j * (n_x * 2 + 1) + 0 * num_dofs] = (CeedScalar)i / (2 * n_x);
36 x_array[i + j * (n_x * 2 + 1) + 1 * num_dofs] = (CeedScalar)j / (2 * n_y);
37 }
38 }
39 CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, x_array);
40 }
41 CeedVectorCreate(ceed, num_dofs, &u);
42 CeedVectorCreate(ceed, num_dofs, &v);
43 CeedVectorCreate(ceed, num_qpts * dim * (dim + 1) / 2, &q_data);
44
45 // Restrictions
46 for (CeedInt i = 0; i < num_elem; i++) {
47 CeedInt col, row, offset;
48 col = i % n_x;
49 row = i / n_x;
50 offset = col * (p - 1) + row * (n_x * 2 + 1) * (p - 1);
51 for (CeedInt j = 0; j < p; j++) {
52 for (CeedInt k = 0; k < p; k++) ind_x[p * (p * i + k) + j] = offset + k * (n_x * 2 + 1) + j;
53 }
54 }
55 CeedElemRestrictionCreate(ceed, num_elem, p * p, dim, num_dofs, dim * num_dofs, CEED_MEM_HOST, CEED_USE_POINTER, ind_x, &elem_restriction_x);
56 CeedElemRestrictionCreate(ceed, num_elem, p * p, 1, 1, num_dofs, CEED_MEM_HOST, CEED_USE_POINTER, ind_x, &elem_restriction_u);
57
58 CeedInt strides_q_data[3] = {1, q * q, q * q * dim * (dim + 1) / 2};
59 CeedElemRestrictionCreateStrided(ceed, num_elem, q * q, dim * (dim + 1) / 2, dim * (dim + 1) / 2 * num_qpts, strides_q_data,
60 &elem_restriction_q_data);
61
62 // Bases
63 CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, p, q, CEED_GAUSS, &basis_x);
64 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, p, q, CEED_GAUSS, &basis_u);
65
66 // QFunction - setup
67 CeedQFunctionCreateInterior(ceed, 1, setup, setup_loc, &qf_setup);
68 CeedQFunctionAddInput(qf_setup, "dx", dim * dim, CEED_EVAL_GRAD);
69 CeedQFunctionAddInput(qf_setup, "weight", 1, CEED_EVAL_WEIGHT);
70 CeedQFunctionAddOutput(qf_setup, "q data", dim * (dim + 1) / 2, CEED_EVAL_NONE);
71
72 // Operator - setup
73 CeedOperatorCreate(ceed, qf_setup, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_setup);
74 CeedOperatorSetField(op_setup, "dx", elem_restriction_x, basis_x, CEED_VECTOR_ACTIVE);
75 CeedOperatorSetField(op_setup, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE);
76 CeedOperatorSetField(op_setup, "q data", elem_restriction_q_data, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE);
77
78 // Apply Setup Operator
79 CeedOperatorApply(op_setup, x, q_data, CEED_REQUEST_IMMEDIATE);
80
81 // QFunction - apply
82 CeedQFunctionCreateInterior(ceed, 1, diff, diff_loc, &qf_diff);
83 CeedQFunctionAddInput(qf_diff, "du", dim, CEED_EVAL_GRAD);
84 CeedQFunctionAddInput(qf_diff, "q data", dim * (dim + 1) / 2, CEED_EVAL_NONE);
85 CeedQFunctionAddOutput(qf_diff, "dv", dim, CEED_EVAL_GRAD);
86
87 // Operator - apply
88 CeedOperatorCreate(ceed, qf_diff, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_diff);
89 CeedOperatorSetField(op_diff, "du", elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE);
90 CeedOperatorSetField(op_diff, "q data", elem_restriction_q_data, CEED_BASIS_NONE, q_data);
91 CeedOperatorSetField(op_diff, "dv", elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE);
92
93 // Fully assemble operator
94 CeedSize num_entries;
95 CeedInt *rows;
96 CeedInt *cols;
97 CeedVector assembled;
98
99 for (CeedInt k = 0; k < num_dofs * num_dofs; ++k) {
100 assembled_values[k] = 0.0;
101 assembled_true[k] = 0.0;
102 }
103 CeedOperatorLinearAssembleSymbolic(op_diff, &num_entries, &rows, &cols);
104 CeedVectorCreate(ceed, num_entries, &assembled);
105 CeedOperatorLinearAssemble(op_diff, assembled);
106 {
107 const CeedScalar *assembled_array;
108
109 CeedVectorGetArrayRead(assembled, CEED_MEM_HOST, &assembled_array);
110 for (CeedInt k = 0; k < num_entries; ++k) assembled_values[rows[k] * num_dofs + cols[k]] += assembled_array[k];
111 CeedVectorRestoreArrayRead(assembled, &assembled_array);
112 }
113
114 // Manually assemble operator
115 CeedVectorSetValue(u, 0.0);
116 for (CeedInt j = 0; j < num_dofs; j++) {
117 CeedScalar *u_array;
118 const CeedScalar *v_array;
119
120 // Set input
121 CeedVectorGetArray(u, CEED_MEM_HOST, &u_array);
122 u_array[j] = 1.0;
123 if (j) u_array[j - 1] = 0.0;
124 CeedVectorRestoreArray(u, &u_array);
125
126 // Compute entries for column j
127 CeedOperatorApply(op_diff, u, v, CEED_REQUEST_IMMEDIATE);
128
129 CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array);
130 for (CeedInt i = 0; i < num_dofs; i++) assembled_true[i * num_dofs + j] = v_array[i];
131 CeedVectorRestoreArrayRead(v, &v_array);
132 }
133
134 // Check output
135 for (CeedInt i = 0; i < num_dofs; i++) {
136 for (CeedInt j = 0; j < num_dofs; j++) {
137 if (fabs(assembled_values[i * num_dofs + j] - assembled_true[i * num_dofs + j]) > 100. * CEED_EPSILON) {
138 // LCOV_EXCL_START
139 printf("[%" CeedInt_FMT ", %" CeedInt_FMT "] Error in assembly: %f != %f\n", i, j, assembled_values[i * num_dofs + j],
140 assembled_true[i * num_dofs + j]);
141 // LCOV_EXCL_STOP
142 }
143 }
144 }
145
146 // Cleanup
147 free(rows);
148 free(cols);
149 CeedVectorDestroy(&x);
150 CeedVectorDestroy(&q_data);
151 CeedVectorDestroy(&u);
152 CeedVectorDestroy(&v);
153 CeedVectorDestroy(&assembled);
154 CeedElemRestrictionDestroy(&elem_restriction_u);
155 CeedElemRestrictionDestroy(&elem_restriction_x);
156 CeedElemRestrictionDestroy(&elem_restriction_q_data);
157 CeedBasisDestroy(&basis_u);
158 CeedBasisDestroy(&basis_x);
159 CeedQFunctionDestroy(&qf_setup);
160 CeedQFunctionDestroy(&qf_diff);
161 CeedOperatorDestroy(&op_setup);
162 CeedOperatorDestroy(&op_diff);
163 CeedDestroy(&ceed);
164 return 0;
165 }
166