Skip to content

Commit ba2cb04

Browse files
zatkins-devzatkins-dev
committed
mixed basis: add failing GPU test
1 parent c941d95 commit ba2cb04

3 files changed

Lines changed: 263 additions & 11 deletions

File tree

‎tests/t572-operator.c‎

Lines changed: 7 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -20,11 +20,11 @@ int main(int argc, char **argv) {
2020
CeedInt q_data_size = dim * (dim + 1) / 2;
2121
CeedInt num_elem = n_x * n_y;
2222
CeedInt num_dofs_u = (n_x * (p_u - 1) + 1) * (n_y * (p_u - 1) + 1), num_dofs_p = (n_x * (p_p - 1) + 1) * (n_y * (p_p - 1) + 1),
23-
num_qpts = num_elem * q * q;
24-
CeedInt ind_x[num_elem * p_u * p_u], ind_p[num_elem * p_p * p_p];
25-
CeedInt l_size = num_comp_u * num_dofs_u + num_comp_p * num_dofs_p;
26-
CeedScalar assembled_values[l_size * l_size];
27-
CeedScalar assembled_true[l_size * l_size];
23+
num_qpts = num_elem * q * q;
24+
CeedInt ind_x[num_elem * p_u * p_u], ind_p[num_elem * p_p * p_p];
25+
CeedInt l_size = num_comp_u * num_dofs_u + num_comp_p * num_dofs_p;
26+
CeedScalar assembled_values[l_size * l_size];
27+
CeedScalar assembled_true[l_size * l_size];
2828

2929
CeedInit(argv[1], &ceed);
3030

@@ -92,8 +92,8 @@ int main(int argc, char **argv) {
9292
CeedQFunctionAddInput(qf_diff, "qdata", q_data_size, CEED_EVAL_NONE);
9393
CeedQFunctionAddInput(qf_diff, "du", num_comp_u * dim, CEED_EVAL_GRAD);
9494
CeedQFunctionAddInput(qf_diff, "dp", num_comp_p * dim, CEED_EVAL_GRAD);
95-
CeedQFunctionAddOutput(qf_diff, "dv", num_comp_u * dim, CEED_EVAL_GRAD);
9695
CeedQFunctionAddOutput(qf_diff, "dq", num_comp_p * dim, CEED_EVAL_GRAD);
96+
CeedQFunctionAddOutput(qf_diff, "dv", num_comp_u * dim, CEED_EVAL_GRAD);
9797

9898
Context_t context = {num_comp_u, num_comp_p};
9999
CeedQFunctionContext ctx;
@@ -113,8 +113,8 @@ int main(int argc, char **argv) {
113113
CeedOperatorSetField(op_diff, "qdata", elem_restriction_q_data, CEED_BASIS_NONE, q_data);
114114
CeedOperatorSetField(op_diff, "du", elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE);
115115
CeedOperatorSetField(op_diff, "dp", elem_restriction_p, basis_p, CEED_VECTOR_ACTIVE);
116-
CeedOperatorSetField(op_diff, "dv", elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE);
117116
CeedOperatorSetField(op_diff, "dq", elem_restriction_p, basis_p, CEED_VECTOR_ACTIVE);
117+
CeedOperatorSetField(op_diff, "dv", elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE);
118118

119119
// Apply Setup Operator
120120
CeedOperatorApply(op_setup, x, q_data, CEED_REQUEST_IMMEDIATE);

‎tests/t572-operator.h‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -49,10 +49,10 @@ CEED_QFUNCTION(multi_basis_diff)(void *ctx, const CeedInt Q, const CeedScalar *c
4949
const CeedScalar *du = in[1];
5050
const CeedScalar *dp = in[2];
5151

52-
// out[0] is output to multiply against gradient v, shape [2, nc_u, Q]
53-
// out[1] is output to multiply against gradient q, shape [2, nc_p, Q]
54-
CeedScalar *dv = out[0];
55-
CeedScalar *dq = out[1];
52+
// out[0] is output to multiply against gradient q, shape [2, nc_p, Q]
53+
// out[1] is output to multiply against gradient v, shape [2, nc_u, Q]
54+
CeedScalar *dq = out[0];
55+
CeedScalar *dv = out[1];
5656

5757
const Context_t *context = (Context_t *)ctx;
5858

‎tests/t573-operator.c‎

Lines changed: 252 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,252 @@
1+
/// @file
2+
/// Test assembly of operator for operator with multiple active bases
3+
/// \test Test assembly of operator for operator with multiple active bases
4+
#include "t539-operator.h"
5+
6+
#include <ceed.h>
7+
#include <math.h>
8+
#include <stdio.h>
9+
10+
int main(int argc, char **argv) {
11+
Ceed ceed;
12+
CeedElemRestriction elem_restriction_x, elem_restriction_u_0, elem_restriction_u_1, elem_restr_q_data_mass, elem_restr_q_data_diff;
13+
CeedBasis basis_x, basis_u_0, basis_u_1;
14+
CeedQFunction qf_setup_mass, qf_setup_diff, qf_apply;
15+
CeedOperator op_setup_mass, op_setup_diff, op_apply;
16+
CeedVector q_data_mass, q_data_diff, x, assembled, u, v;
17+
CeedInt p_0 = 2, p_1 = 3, q = 4, dim = 2, num_comp_0 = 2, num_comp_1 = 1;
18+
CeedInt n_x = 1, n_y = 2, num_elem = n_x * n_y;
19+
CeedInt num_dofs_0 = (n_x * (p_0 - 1) + 1) * (n_y * (p_0 - 1) + 1), num_dofs_1 = (n_x * (p_1 - 1) + 1) * (n_y * (p_1 - 1) + 1);
20+
CeedInt num_qpts = num_elem * q * q;
21+
CeedInt ind_u_0[num_elem * p_0 * p_0], ind_u_1[num_elem * p_1 * p_1];
22+
CeedInt l_size = num_comp_0 * num_dofs_0 + num_comp_1 * num_dofs_1;
23+
CeedScalar assembled_values[l_size * l_size];
24+
CeedScalar assembled_true[l_size * l_size];
25+
26+
CeedInit(argv[1], &ceed);
27+
28+
// Vectors
29+
CeedVectorCreate(ceed, dim * num_dofs_0, &x);
30+
{
31+
CeedScalar x_array[dim * num_dofs_0];
32+
33+
for (CeedInt i = 0; i < n_x * (p_0 - 1) + 1; i++) {
34+
for (CeedInt j = 0; j < n_y * (p_0 - 1) + 1; j++) {
35+
x_array[i + j * (n_x * (p_0 - 1) + 1) + 0 * num_dofs_0] = (CeedScalar)i / ((p_0 - 1) * n_x);
36+
x_array[i + j * (n_x * (p_0 - 1) + 1) + 1 * num_dofs_0] = (CeedScalar)j / ((p_0 - 1) * n_y);
37+
}
38+
}
39+
CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, x_array);
40+
}
41+
CeedVectorCreate(ceed, num_comp_0 * num_dofs_0 + num_comp_1 * num_dofs_1, &u);
42+
CeedVectorCreate(ceed, num_comp_0 * num_dofs_0 + num_comp_1 * num_dofs_1, &v);
43+
CeedVectorCreate(ceed, num_qpts, &q_data_mass);
44+
CeedVectorCreate(ceed, num_qpts * dim * (dim + 1) / 2, &q_data_diff);
45+
46+
// Restrictions
47+
for (CeedInt i = 0; i < num_elem; i++) {
48+
CeedInt col, row, offset;
49+
col = i % n_x;
50+
row = i / n_x;
51+
offset = col * (p_0 - 1) + row * (n_x * (p_0 - 1) + 1) * (p_0 - 1);
52+
for (CeedInt j = 0; j < p_0; j++) {
53+
for (CeedInt k = 0; k < p_0; k++) ind_u_0[p_0 * (p_0 * i + k) + j] = offset + k * (n_x * (p_0 - 1) + 1) + j;
54+
}
55+
offset = col * (p_1 - 1) + row * (n_x * (p_1 - 1) + 1) * (p_1 - 1) + num_dofs_0 * num_comp_0;
56+
for (CeedInt j = 0; j < p_1; j++) {
57+
for (CeedInt k = 0; k < p_1; k++) ind_u_1[p_1 * (p_1 * i + k) + j] = offset + k * (n_x * (p_1 - 1) + 1) + j;
58+
}
59+
}
60+
CeedElemRestrictionCreate(ceed, num_elem, p_0 * p_0, dim, num_dofs_0, dim * num_dofs_0, CEED_MEM_HOST, CEED_USE_POINTER, ind_u_0,
61+
&elem_restriction_x);
62+
CeedElemRestrictionCreate(ceed, num_elem, p_0 * p_0, num_comp_0, num_dofs_0, num_comp_0 * num_dofs_0 + num_comp_1 * num_dofs_1, CEED_MEM_HOST,
63+
CEED_USE_POINTER, ind_u_0, &elem_restriction_u_0);
64+
CeedElemRestrictionCreate(ceed, num_elem, p_1 * p_1, num_comp_1, num_dofs_1, num_comp_0 * num_dofs_0 + num_comp_1 * num_dofs_1, CEED_MEM_HOST,
65+
CEED_USE_POINTER, ind_u_1, &elem_restriction_u_1);
66+
67+
CeedInt strides_q_data_mass[3] = {1, q * q, q * q};
68+
CeedElemRestrictionCreateStrided(ceed, num_elem, q * q, 1, num_qpts, strides_q_data_mass, &elem_restr_q_data_mass);
69+
70+
CeedInt strides_q_data_diff[3] = {1, q * q, dim * (dim + 1) / 2 * q * q};
71+
CeedElemRestrictionCreateStrided(ceed, num_elem, q * q, dim * (dim + 1) / 2, dim * (dim + 1) / 2 * num_qpts, strides_q_data_diff,
72+
&elem_restr_q_data_diff);
73+
74+
// Bases
75+
CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, p_0, q, CEED_GAUSS, &basis_x);
76+
CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_0, p_0, q, CEED_GAUSS, &basis_u_0);
77+
CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_1, p_1, q, CEED_GAUSS, &basis_u_1);
78+
79+
// QFunction - setup mass
80+
CeedQFunctionCreateInteriorByName(ceed, "Mass2DBuild", &qf_setup_mass);
81+
82+
// Operator - setup mass
83+
CeedOperatorCreate(ceed, qf_setup_mass, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_setup_mass);
84+
CeedOperatorSetField(op_setup_mass, "dx", elem_restriction_x, basis_x, CEED_VECTOR_ACTIVE);
85+
CeedOperatorSetField(op_setup_mass, "weights", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE);
86+
CeedOperatorSetField(op_setup_mass, "qdata", elem_restr_q_data_mass, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE);
87+
88+
// QFunction - setup diffusion
89+
CeedQFunctionCreateInteriorByName(ceed, "Poisson2DBuild", &qf_setup_diff);
90+
91+
// Operator - setup diffusion
92+
CeedOperatorCreate(ceed, qf_setup_diff, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_setup_diff);
93+
CeedOperatorSetField(op_setup_diff, "dx", elem_restriction_x, basis_x, CEED_VECTOR_ACTIVE);
94+
CeedOperatorSetField(op_setup_diff, "weights", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE);
95+
CeedOperatorSetField(op_setup_diff, "qdata", elem_restr_q_data_diff, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE);
96+
97+
// Apply Setup Operators
98+
CeedOperatorApply(op_setup_mass, x, q_data_mass, CEED_REQUEST_IMMEDIATE);
99+
CeedOperatorApply(op_setup_diff, x, q_data_diff, CEED_REQUEST_IMMEDIATE);
100+
101+
// QFunction - apply
102+
CeedQFunctionCreateInterior(ceed, 1, apply, apply_loc, &qf_apply);
103+
CeedQFunctionAddInput(qf_apply, "du_0", num_comp_0 * dim, CEED_EVAL_GRAD);
104+
CeedQFunctionAddInput(qf_apply, "mass qdata", 1, CEED_EVAL_NONE);
105+
CeedQFunctionAddInput(qf_apply, "diff qdata", dim * (dim + 1) / 2, CEED_EVAL_NONE);
106+
CeedQFunctionAddInput(qf_apply, "u_0", num_comp_0, CEED_EVAL_INTERP);
107+
CeedQFunctionAddInput(qf_apply, "u_1", num_comp_1, CEED_EVAL_INTERP);
108+
CeedQFunctionAddOutput(qf_apply, "v_0", num_comp_0, CEED_EVAL_INTERP);
109+
CeedQFunctionAddOutput(qf_apply, "v_1", num_comp_1, CEED_EVAL_INTERP);
110+
CeedQFunctionAddOutput(qf_apply, "dv_0", num_comp_0 * dim, CEED_EVAL_GRAD);
111+
112+
// Operator - apply
113+
CeedOperatorCreate(ceed, qf_apply, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_apply);
114+
CeedOperatorSetField(op_apply, "du_0", elem_restriction_u_0, basis_u_0, CEED_VECTOR_ACTIVE);
115+
CeedOperatorSetField(op_apply, "mass qdata", elem_restr_q_data_mass, CEED_BASIS_NONE, q_data_mass);
116+
CeedOperatorSetField(op_apply, "diff qdata", elem_restr_q_data_diff, CEED_BASIS_NONE, q_data_diff);
117+
CeedOperatorSetField(op_apply, "u_0", elem_restriction_u_0, basis_u_0, CEED_VECTOR_ACTIVE);
118+
CeedOperatorSetField(op_apply, "u_1", elem_restriction_u_1, basis_u_1, CEED_VECTOR_ACTIVE);
119+
CeedOperatorSetField(op_apply, "v_0", elem_restriction_u_0, basis_u_0, CEED_VECTOR_ACTIVE);
120+
CeedOperatorSetField(op_apply, "v_1", elem_restriction_u_1, basis_u_1, CEED_VECTOR_ACTIVE);
121+
CeedOperatorSetField(op_apply, "dv_0", elem_restriction_u_0, basis_u_0, CEED_VECTOR_ACTIVE);
122+
123+
// Fuly assemble operator
124+
CeedSize num_entries;
125+
CeedInt *rows;
126+
CeedInt *cols;
127+
128+
for (CeedInt k = 0; k < l_size * l_size; k++) {
129+
assembled_values[k] = 0.0;
130+
assembled_true[k] = 0.0;
131+
}
132+
CeedOperatorLinearAssembleSymbolic(op_apply, &num_entries, &rows, &cols);
133+
CeedVectorCreate(ceed, num_entries, &assembled);
134+
CeedOperatorLinearAssemble(op_apply, assembled);
135+
{
136+
const CeedScalar *assembled_array;
137+
138+
CeedVectorGetArrayRead(assembled, CEED_MEM_HOST, &assembled_array);
139+
for (CeedInt k = 0; k < num_entries; k++) assembled_values[rows[k] * l_size + cols[k]] += assembled_array[k];
140+
CeedVectorRestoreArrayRead(assembled, &assembled_array);
141+
}
142+
143+
// Manually assemble operator
144+
CeedInt old_index = -1;
145+
146+
CeedVectorSetValue(u, 0.0);
147+
for (CeedInt j = 0; j < l_size; j++) {
148+
CeedScalar *u_array;
149+
const CeedScalar *v_array;
150+
151+
// Set input
152+
CeedVectorGetArray(u, CEED_MEM_HOST, &u_array);
153+
u_array[j] = 1.0;
154+
if (j > 0) u_array[old_index] = 0.0;
155+
old_index = j;
156+
CeedVectorRestoreArray(u, &u_array);
157+
158+
// Compute effect of DoF ind
159+
CeedOperatorApply(op_apply, u, v, CEED_REQUEST_IMMEDIATE);
160+
161+
CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array);
162+
for (CeedInt i = 0; i < l_size; i++) assembled_true[i * l_size + j] = v_array[i];
163+
CeedVectorRestoreArrayRead(v, &v_array);
164+
}
165+
166+
// Check output
167+
for (CeedInt i = 0; i < l_size; i++) {
168+
for (CeedInt j = 0; j < l_size; j++) {
169+
const CeedScalar assembled_value = assembled_values[i * l_size + j];
170+
const CeedScalar assembled_true_value = assembled_true[i * l_size + j];
171+
const CeedScalar error = fabs(assembled_value - assembled_true_value) / (fmax(fabs(assembled_true_value), 1.0));
172+
173+
if (!(error < 1000. * CEED_EPSILON)) {
174+
// LCOV_EXCL_START
175+
const CeedInt node_out = (i < num_comp_0 * num_dofs_0) ? i % num_dofs_0 : (i - num_comp_0 * num_dofs_0) % num_dofs_0;
176+
const CeedInt comp_out = (i < num_comp_0 * num_dofs_0) ? i / num_dofs_0 : (i - num_comp_0 * num_dofs_0) / num_dofs_0;
177+
const CeedInt node_in = (j < num_comp_0 * num_dofs_0) ? j % num_dofs_0 : (j - num_comp_0 * num_dofs_0) % num_dofs_0;
178+
const CeedInt comp_in = (j < num_comp_0 * num_dofs_0) ? j / num_dofs_0 : (j - num_comp_0 * num_dofs_0) / num_dofs_0;
179+
180+
printf("[(%s, %" CeedInt_FMT ", %" CeedInt_FMT "), (%s, %" CeedInt_FMT ", %" CeedInt_FMT ")] Error in assembly: %g != %g (error: %g)\n",
181+
(i < num_comp_0 * num_dofs_0 ? "0" : "1"), node_out, comp_out, (j < num_comp_0 * num_dofs_0 ? "0" : "1"), node_in, comp_in,
182+
assembled_value, assembled_true_value, error);
183+
// LCOV_EXCL_STOP
184+
}
185+
}
186+
}
187+
// printf("\n");
188+
// printf("ASSEMBLED\n");
189+
// printf("\n");
190+
191+
// for (CeedInt i = 0; i < l_size; i++) {
192+
// if (i == num_comp_0 * num_dofs_0) {
193+
// for (CeedInt j = 0; j < l_size * 6 + 3; j++) printf("-");
194+
// printf("\n");
195+
// }
196+
// for (CeedInt j = 0; j < l_size; j++) {
197+
// // LCOV_EXCL_START
198+
// const CeedScalar assembled_value = assembled_values[i * l_size + j];
199+
200+
// if (j == num_comp_0 * num_dofs_0) printf(" | ");
201+
// if (assembled_value == 0) printf(" ");
202+
// else printf("%8.1e ", assembled_value);
203+
// // LCOV_EXCL_STOP
204+
// }
205+
// printf("\n");
206+
// }
207+
208+
// printf("\n");
209+
// printf("TRUE\n");
210+
// printf("\n");
211+
212+
// for (CeedInt i = 0; i < l_size; i++) {
213+
// if (i == num_comp_0 * num_dofs_0) {
214+
// for (CeedInt j = 0; j < l_size * 9 + 3; j++) printf("-");
215+
// printf("\n");
216+
// }
217+
// for (CeedInt j = 0; j < l_size; j++) {
218+
// // LCOV_EXCL_START
219+
// const CeedScalar assembled_value = assembled_true[i * l_size + j];
220+
221+
// if (j == num_comp_0 * num_dofs_0) printf(" | ");
222+
// if (assembled_value == 0) printf(" ");
223+
// else printf("%8.1e ", assembled_value);
224+
// // LCOV_EXCL_STOP
225+
// }
226+
// printf("\n");
227+
// }
228+
229+
// Cleanup
230+
CeedVectorDestroy(&x);
231+
CeedVectorDestroy(&assembled);
232+
CeedVectorDestroy(&u);
233+
CeedVectorDestroy(&v);
234+
CeedVectorDestroy(&q_data_mass);
235+
CeedVectorDestroy(&q_data_diff);
236+
CeedElemRestrictionDestroy(&elem_restriction_x);
237+
CeedElemRestrictionDestroy(&elem_restriction_u_0);
238+
CeedElemRestrictionDestroy(&elem_restriction_u_1);
239+
CeedElemRestrictionDestroy(&elem_restr_q_data_mass);
240+
CeedElemRestrictionDestroy(&elem_restr_q_data_diff);
241+
CeedBasisDestroy(&basis_x);
242+
CeedBasisDestroy(&basis_u_0);
243+
CeedBasisDestroy(&basis_u_1);
244+
CeedQFunctionDestroy(&qf_setup_mass);
245+
CeedQFunctionDestroy(&qf_setup_diff);
246+
CeedQFunctionDestroy(&qf_apply);
247+
CeedOperatorDestroy(&op_setup_mass);
248+
CeedOperatorDestroy(&op_setup_diff);
249+
CeedOperatorDestroy(&op_apply);
250+
CeedDestroy(&ceed);
251+
return 0;
252+
}

0 commit comments

Comments
 (0)