1 /// @file
2 /// Test assembly of Poisson operator QFunction
3 /// \test Test assembly of Poisson operator QFunction
4 #include "t531-operator.h"
5
6 #include <ceed.h>
7 #include <math.h>
8 #include <stdio.h>
9 #include <stdlib.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, elem_restriction_assembled = NULL;
14 CeedBasis basis_x, basis_u;
15 CeedQFunction qf_setup, qf_diff, qf_diff_assembled;
16 CeedOperator op_setup, op_diff, op_diff_assembled;
17 CeedVector q_data, x, assembled = NULL, u, v, v_assembled;
18 CeedInt num_elem = 6, p = 3, q = 4, dim = 2;
19 CeedInt nx = 3, ny = 2;
20 CeedInt num_dofs = (nx * 2 + 1) * (ny * 2 + 1), num_qpts = num_elem * q * q;
21 CeedInt ind_x[num_elem * p * p];
22
23 CeedInit(argv[1], &ceed);
24
25 // Vectors
26 CeedVectorCreate(ceed, dim * num_dofs, &x);
27 {
28 CeedScalar x_array[dim * num_dofs];
29
30 for (CeedInt i = 0; i < nx * 2 + 1; i++) {
31 for (CeedInt j = 0; j < ny * 2 + 1; j++) {
32 x_array[i + j * (nx * 2 + 1) + 0 * num_dofs] = (CeedScalar)i / (2 * nx) + 0.5 * j;
33 x_array[i + j * (nx * 2 + 1) + 1 * num_dofs] = (CeedScalar)j / (2 * ny) + 0.5 * i;
34 }
35 }
36 CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, x_array);
37 }
38 CeedVectorCreate(ceed, num_dofs, &u);
39 {
40 CeedScalar *u_array;
41
42 CeedVectorGetArrayWrite(u, CEED_MEM_HOST, &u_array);
43 for (CeedInt i = 0; i < nx * 2 + 1; i++) {
44 for (CeedInt j = 0; j < ny * 2 + 1; j++) {
45 u_array[i + j * (nx * 2 + 1)] = i * nx + j * ny;
46 }
47 }
48 CeedVectorRestoreArray(u, &u_array);
49 }
50 CeedVectorCreate(ceed, num_dofs, &v);
51 CeedVectorCreate(ceed, num_dofs, &v_assembled);
52 CeedVectorCreate(ceed, num_qpts * dim * (dim + 1) / 2, &q_data);
53
54 // Restrictions
55 for (CeedInt i = 0; i < num_elem; i++) {
56 CeedInt col, row, offset;
57 col = i % nx;
58 row = i / nx;
59 offset = col * (p - 1) + row * (nx * 2 + 1) * (p - 1);
60 for (CeedInt j = 0; j < p; j++) {
61 for (CeedInt k = 0; k < p; k++) ind_x[p * (p * i + k) + j] = offset + k * (nx * 2 + 1) + j;
62 }
63 }
64 CeedElemRestrictionCreate(ceed, num_elem, p * p, dim, num_dofs, dim * num_dofs, CEED_MEM_HOST, CEED_USE_POINTER, ind_x, &elem_restriction_x);
65 CeedElemRestrictionCreate(ceed, num_elem, p * p, 1, 1, num_dofs, CEED_MEM_HOST, CEED_USE_POINTER, ind_x, &elem_restriction_u);
66
67 CeedInt strides_q_data[3] = {1, q * q, q * q * dim * (dim + 1) / 2};
68 CeedElemRestrictionCreateStrided(ceed, num_elem, q * q, dim * (dim + 1) / 2, dim * (dim + 1) / 2 * num_qpts, strides_q_data,
69 &elem_restriction_q_data);
70
71 // Bases
72 CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, p, q, CEED_GAUSS, &basis_x);
73 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, p, q, CEED_GAUSS, &basis_u);
74
75 // QFunction - setup
76 CeedQFunctionCreateInterior(ceed, 1, setup, setup_loc, &qf_setup);
77 CeedQFunctionAddInput(qf_setup, "dx", dim * dim, CEED_EVAL_GRAD);
78 CeedQFunctionAddInput(qf_setup, "weight", 1, CEED_EVAL_WEIGHT);
79 CeedQFunctionAddOutput(qf_setup, "q data", dim * (dim + 1) / 2, CEED_EVAL_NONE);
80
81 // Operator - setup
82 CeedOperatorCreate(ceed, qf_setup, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_setup);
83 CeedOperatorSetField(op_setup, "dx", elem_restriction_x, basis_x, CEED_VECTOR_ACTIVE);
84 CeedOperatorSetField(op_setup, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE);
85 CeedOperatorSetField(op_setup, "q data", elem_restriction_q_data, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE);
86
87 // Apply Setup Operator
88 CeedOperatorApply(op_setup, x, q_data, CEED_REQUEST_IMMEDIATE);
89
90 // QFunction - apply
91 CeedQFunctionCreateInterior(ceed, 1, diff, diff_loc, &qf_diff);
92 CeedQFunctionAddInput(qf_diff, "du", dim, CEED_EVAL_GRAD);
93 CeedQFunctionAddInput(qf_diff, "q data", dim * (dim + 1) / 2, CEED_EVAL_NONE);
94 CeedQFunctionAddOutput(qf_diff, "dv", dim, CEED_EVAL_GRAD);
95
96 // Operator - apply
97 CeedOperatorCreate(ceed, qf_diff, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_diff);
98 CeedOperatorSetField(op_diff, "du", elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE);
99 CeedOperatorSetField(op_diff, "q data", elem_restriction_q_data, CEED_BASIS_NONE, q_data);
100 CeedOperatorSetField(op_diff, "dv", elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE);
101
102 // Apply original Poisson Operator
103 CeedOperatorApply(op_diff, u, v, CEED_REQUEST_IMMEDIATE);
104
105 // Assemble QFunction
106 CeedOperatorSetQFunctionAssemblyReuse(op_diff, true);
107 CeedOperatorLinearAssembleQFunctionBuildOrUpdate(op_diff, &assembled, &elem_restriction_assembled, CEED_REQUEST_IMMEDIATE);
108 // Second call will be no-op since SetQFunctionUpdated was not called
109 CeedOperatorSetQFunctionAssemblyDataUpdateNeeded(op_diff, false);
110 CeedOperatorLinearAssembleQFunctionBuildOrUpdate(op_diff, &assembled, &elem_restriction_assembled, CEED_REQUEST_IMMEDIATE);
111
112 // QFunction - apply assembled
113 CeedQFunctionCreateInterior(ceed, 1, diff_lin, diff_lin_loc, &qf_diff_assembled);
114 CeedQFunctionAddInput(qf_diff_assembled, "du", dim, CEED_EVAL_GRAD);
115 CeedQFunctionAddInput(qf_diff_assembled, "q data", dim * dim, CEED_EVAL_NONE);
116 CeedQFunctionAddOutput(qf_diff_assembled, "dv", dim, CEED_EVAL_GRAD);
117
118 // Operator - apply assembled
119 CeedOperatorCreate(ceed, qf_diff_assembled, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_diff_assembled);
120 CeedOperatorSetField(op_diff_assembled, "du", elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE);
121 CeedOperatorSetField(op_diff_assembled, "q data", elem_restriction_assembled, CEED_BASIS_NONE, assembled);
122 CeedOperatorSetField(op_diff_assembled, "dv", elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE);
123
124 // Apply new Poisson Operator
125 CeedOperatorApply(op_diff_assembled, u, v_assembled, CEED_REQUEST_IMMEDIATE);
126
127 // Check output
128 {
129 const CeedScalar *v_array, *v_assembled_array;
130
131 CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array);
132 CeedVectorGetArrayRead(v_assembled, CEED_MEM_HOST, &v_assembled_array);
133 for (CeedInt i = 0; i < num_dofs; i++) {
134 if (fabs(v_array[i] - v_assembled_array[i]) > 100. * CEED_EPSILON) {
135 // LCOV_EXCL_START
136 printf("Error: Linearized operator computed v[%d] = %f != %f\n", i, v_assembled_array[i], v_array[i]);
137 // LCOV_EXCL_STOP
138 }
139 }
140 CeedVectorRestoreArrayRead(v, &v_array);
141 CeedVectorRestoreArrayRead(v_assembled, &v_assembled_array);
142 }
143
144 // Cleanup
145 CeedVectorDestroy(&x);
146 CeedVectorDestroy(&assembled);
147 CeedVectorDestroy(&q_data);
148 CeedVectorDestroy(&u);
149 CeedVectorDestroy(&v);
150 CeedVectorDestroy(&v_assembled);
151 CeedElemRestrictionDestroy(&elem_restriction_u);
152 CeedElemRestrictionDestroy(&elem_restriction_x);
153 CeedElemRestrictionDestroy(&elem_restriction_q_data);
154 CeedElemRestrictionDestroy(&elem_restriction_assembled);
155 CeedBasisDestroy(&basis_u);
156 CeedBasisDestroy(&basis_x);
157 CeedQFunctionDestroy(&qf_setup);
158 CeedQFunctionDestroy(&qf_diff);
159 CeedQFunctionDestroy(&qf_diff_assembled);
160 CeedOperatorDestroy(&op_setup);
161 CeedOperatorDestroy(&op_diff);
162 CeedOperatorDestroy(&op_diff_assembled);
163 CeedDestroy(&ceed);
164 return 0;
165 }
166