Skip to content

Commit 1f5047a

Browse files
zatkins-devzatkins-dev
committed
tests: Add multi-basis assembly tests
1 parent 942790e commit 1f5047a

4 files changed

Lines changed: 537 additions & 0 deletions

File tree

‎tests/t571-operator.c‎

Lines changed: 206 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,206 @@
1+
/// @file
2+
/// Test full assembly of multi-basis asymmetric mass-like operator
3+
/// \test Test full assembly of multi-basis asymmetric mass-like operator
4+
#include "t571-operator.h"
5+
6+
#include <ceed.h>
7+
#include <math.h>
8+
#include <stdio.h>
9+
#include <stdlib.h>
10+
11+
int main(int argc, char **argv) {
12+
Ceed ceed;
13+
CeedElemRestriction elem_restriction_x, elem_restriction_u, elem_restriction_p, elem_restriction_q_data;
14+
CeedBasis basis_x, basis_u, basis_p;
15+
CeedQFunction qf_setup, qf_mass;
16+
CeedOperator op_setup, op_mass;
17+
CeedVector q_data, x, u, v;
18+
CeedInt p_u = 3, p_p = 2, q = 3, dim = 2, num_comp_u = 3, num_comp_p = 2;
19+
CeedInt n_x = 3, n_y = 2;
20+
CeedInt num_elem = n_x * n_y;
21+
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),
22+
num_qpts = num_elem * q * q;
23+
CeedInt ind_x[num_elem * p_u * p_u], ind_p[num_elem * p_p * p_p];
24+
CeedInt l_size = num_comp_u * num_dofs_u + num_comp_p * num_dofs_p;
25+
CeedScalar assembled_values[l_size * l_size];
26+
CeedScalar assembled_true[l_size * l_size];
27+
28+
CeedInit(argv[1], &ceed);
29+
30+
// Vectors
31+
CeedVectorCreate(ceed, dim * num_dofs_u, &x);
32+
{
33+
CeedScalar x_array[dim * num_dofs_u];
34+
35+
for (CeedInt i = 0; i < n_x * (p_u - 1) + 1; i++) {
36+
for (CeedInt j = 0; j < n_y * (p_u - 1) + 1; j++) {
37+
x_array[i + j * (n_x * 2 + 1) + 0 * num_dofs_u] = (CeedScalar)i / (n_x * (p_u - 1));
38+
x_array[i + j * (n_x * 2 + 1) + 1 * num_dofs_u] = (CeedScalar)j / (n_y * (p_u - 1));
39+
}
40+
}
41+
CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, x_array);
42+
}
43+
CeedVectorCreate(ceed, l_size, &u);
44+
CeedVectorCreate(ceed, l_size, &v);
45+
CeedVectorCreate(ceed, num_qpts, &q_data);
46+
47+
// Restrictions
48+
for (CeedInt i = 0; i < num_elem; i++) {
49+
CeedInt col, row, offset;
50+
51+
col = i % n_x;
52+
row = i / n_x;
53+
offset = col * (p_u - 1) + row * (n_x * (p_u - 1) + 1) * (p_u - 1);
54+
for (CeedInt j = 0; j < p_u; j++) {
55+
for (CeedInt k = 0; k < p_u; k++) ind_x[p_u * (p_u * i + k) + j] = offset + k * p_u + j;
56+
}
57+
}
58+
const CeedInt offset_p = num_comp_u * num_dofs_u;
59+
60+
for (CeedInt k = 0; k < num_elem; k++) {
61+
CeedInt col, row, offset;
62+
63+
col = k % n_x;
64+
row = k / n_x;
65+
offset = col * (p_p - 1) + row * (n_x * (p_p - 1) + 1) * (p_p - 1);
66+
// Data for node i, component j, element k can be found in the L-vector at index offsets[i + k*elem_size] + j*comp_stride.
67+
for (CeedInt j = 0; j < p_p; j++) {
68+
for (CeedInt i = 0; i < p_p; i++) ind_p[k * p_p * p_p + (p_p * i + j)] = offset + i * p_p + j + offset_p;
69+
}
70+
}
71+
CeedElemRestrictionCreate(ceed, num_elem, p_u * p_u, dim, num_dofs_u, dim * num_dofs_u, CEED_MEM_HOST, CEED_USE_POINTER, ind_x,
72+
&elem_restriction_x);
73+
CeedElemRestrictionCreate(ceed, num_elem, p_u * p_u, num_comp_u, num_dofs_u, l_size, CEED_MEM_HOST, CEED_USE_POINTER, ind_x, &elem_restriction_u);
74+
CeedElemRestrictionCreate(ceed, num_elem, p_p * p_p, num_comp_p, num_dofs_p, l_size, CEED_MEM_HOST, CEED_USE_POINTER, ind_p, &elem_restriction_p);
75+
76+
CeedInt strides_q_data[3] = {1, q * q * num_elem, q * q};
77+
CeedElemRestrictionCreateStrided(ceed, num_elem, q * q, 1, num_qpts * 1, strides_q_data, &elem_restriction_q_data);
78+
79+
// Bases
80+
CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, p_u, q, CEED_GAUSS, &basis_x);
81+
CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_u, p_u, q, CEED_GAUSS, &basis_u);
82+
CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_p, p_p, q, CEED_GAUSS, &basis_p);
83+
84+
// QFunctions
85+
CeedQFunctionCreateInterior(ceed, 1, setup, setup_loc, &qf_setup);
86+
CeedQFunctionAddInput(qf_setup, "weight", 1, CEED_EVAL_WEIGHT);
87+
CeedQFunctionAddInput(qf_setup, "dx", dim * dim, CEED_EVAL_GRAD);
88+
CeedQFunctionAddOutput(qf_setup, "qdata", 1, CEED_EVAL_NONE);
89+
90+
CeedQFunctionCreateInterior(ceed, 1, multi_basis, multi_basis_loc, &qf_mass);
91+
CeedQFunctionAddInput(qf_mass, "qdata", 1, CEED_EVAL_NONE);
92+
CeedQFunctionAddInput(qf_mass, "du", num_comp_u, CEED_EVAL_INTERP);
93+
CeedQFunctionAddInput(qf_mass, "dp", num_comp_p, CEED_EVAL_INTERP);
94+
CeedQFunctionAddOutput(qf_mass, "dv", num_comp_u, CEED_EVAL_INTERP);
95+
CeedQFunctionAddOutput(qf_mass, "dq", num_comp_p, CEED_EVAL_INTERP);
96+
97+
Context_t context = {num_comp_u, num_comp_p};
98+
CeedQFunctionContext ctx;
99+
100+
CeedQFunctionContextCreate(ceed, &ctx);
101+
CeedQFunctionContextSetData(ctx, CEED_MEM_HOST, CEED_COPY_VALUES, sizeof(Context_t), &context);
102+
CeedQFunctionSetContext(qf_mass, ctx);
103+
CeedQFunctionContextDestroy(&ctx);
104+
105+
// Operators
106+
CeedOperatorCreate(ceed, qf_setup, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_setup);
107+
CeedOperatorSetField(op_setup, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE);
108+
CeedOperatorSetField(op_setup, "dx", elem_restriction_x, basis_x, CEED_VECTOR_ACTIVE);
109+
CeedOperatorSetField(op_setup, "qdata", elem_restriction_q_data, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE);
110+
111+
CeedOperatorCreate(ceed, qf_mass, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_mass);
112+
CeedOperatorSetField(op_mass, "qdata", elem_restriction_q_data, CEED_BASIS_NONE, q_data);
113+
CeedOperatorSetField(op_mass, "du", elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE);
114+
CeedOperatorSetField(op_mass, "dp", elem_restriction_p, basis_p, CEED_VECTOR_ACTIVE);
115+
CeedOperatorSetField(op_mass, "dv", elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE);
116+
CeedOperatorSetField(op_mass, "dq", elem_restriction_p, basis_p, CEED_VECTOR_ACTIVE);
117+
118+
// Apply Setup Operator
119+
CeedOperatorApply(op_setup, x, q_data, CEED_REQUEST_IMMEDIATE);
120+
121+
// Fuly assemble operator
122+
CeedSize num_entries;
123+
CeedInt *rows;
124+
CeedInt *cols;
125+
CeedVector assembled;
126+
127+
for (CeedInt k = 0; k < l_size * l_size; k++) {
128+
assembled_values[k] = 0.0;
129+
assembled_true[k] = 0.0;
130+
}
131+
CeedOperatorLinearAssembleSymbolic(op_mass, &num_entries, &rows, &cols);
132+
CeedVectorCreate(ceed, num_entries, &assembled);
133+
CeedOperatorLinearAssemble(op_mass, assembled);
134+
{
135+
const CeedScalar *assembled_array;
136+
137+
CeedVectorGetArrayRead(assembled, CEED_MEM_HOST, &assembled_array);
138+
for (CeedInt k = 0; k < num_entries; k++) assembled_values[rows[k] * l_size + cols[k]] += assembled_array[k];
139+
CeedVectorRestoreArrayRead(assembled, &assembled_array);
140+
}
141+
142+
// Manually assemble operator
143+
CeedInt old_index = -1;
144+
145+
CeedVectorSetValue(u, 0.0);
146+
for (CeedInt j = 0; j < l_size; j++) {
147+
CeedScalar *u_array;
148+
const CeedScalar *v_array;
149+
150+
// Set input
151+
CeedVectorGetArray(u, CEED_MEM_HOST, &u_array);
152+
u_array[j] = 1.0;
153+
if (j > 0) u_array[old_index] = 0.0;
154+
old_index = j;
155+
CeedVectorRestoreArray(u, &u_array);
156+
157+
// Compute effect of DoF ind
158+
CeedOperatorApply(op_mass, u, v, CEED_REQUEST_IMMEDIATE);
159+
160+
CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array);
161+
for (CeedInt i = 0; i < l_size; i++) assembled_true[i * l_size + j] = v_array[i];
162+
CeedVectorRestoreArrayRead(v, &v_array);
163+
}
164+
165+
// Check output
166+
for (CeedInt i = 0; i < l_size; i++) {
167+
for (CeedInt j = 0; j < l_size; j++) {
168+
const CeedScalar assembled_value = assembled_values[i * l_size + j];
169+
const CeedScalar assembled_true_value = assembled_true[i * l_size + j];
170+
if (!(fabs(assembled_value - assembled_true_value) / (fmax(fabs(assembled_true_value), 1.0)) < 100. * CEED_EPSILON)) {
171+
// LCOV_EXCL_START
172+
const CeedInt node_out = (i < num_comp_u * num_dofs_u) ? i % num_dofs_u : (i - num_comp_u * num_dofs_u) % num_dofs_p;
173+
const CeedInt comp_out = (i < num_comp_u * num_dofs_u) ? i / num_dofs_u : (i - num_comp_u * num_dofs_u) / num_dofs_p;
174+
const CeedInt node_in = (j < num_comp_u * num_dofs_u) ? j % num_dofs_u : (j - num_comp_u * num_dofs_u) % num_dofs_p;
175+
const CeedInt comp_in = (j < num_comp_u * num_dofs_u) ? j / num_dofs_u : (j - num_comp_u * num_dofs_u) / num_dofs_p;
176+
177+
printf("[(%s, %" CeedInt_FMT ", %" CeedInt_FMT "), (%s, %" CeedInt_FMT ", %" CeedInt_FMT ")] Error in assembly: %g != %g (error: %g)\n",
178+
(i < num_comp_u * num_dofs_u ? "u" : "p"), node_out, comp_out, (j < num_comp_u * num_dofs_u ? "u" : "p"), node_in, comp_in,
179+
assembled_value, assembled_true_value, fabs(assembled_value - assembled_true_value) / (fmax(fabs(assembled_true_value), 1.0)));
180+
// LCOV_EXCL_STOP
181+
}
182+
}
183+
}
184+
185+
// Cleanup
186+
free(rows);
187+
free(cols);
188+
CeedVectorDestroy(&x);
189+
CeedVectorDestroy(&q_data);
190+
CeedVectorDestroy(&u);
191+
CeedVectorDestroy(&v);
192+
CeedVectorDestroy(&assembled);
193+
CeedElemRestrictionDestroy(&elem_restriction_u);
194+
CeedElemRestrictionDestroy(&elem_restriction_p);
195+
CeedElemRestrictionDestroy(&elem_restriction_x);
196+
CeedElemRestrictionDestroy(&elem_restriction_q_data);
197+
CeedBasisDestroy(&basis_u);
198+
CeedBasisDestroy(&basis_p);
199+
CeedBasisDestroy(&basis_x);
200+
CeedQFunctionDestroy(&qf_setup);
201+
CeedQFunctionDestroy(&qf_mass);
202+
CeedOperatorDestroy(&op_setup);
203+
CeedOperatorDestroy(&op_mass);
204+
CeedDestroy(&ceed);
205+
return 0;
206+
}

‎tests/t571-operator.h‎

Lines changed: 44 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,44 @@
1+
// Copyright (c) 2017-2026, Lawrence Livermore National Security, LLC and other CEED contributors.
2+
// All Rights Reserved. See the top-level LICENSE and NOTICE files for details.
3+
//
4+
// SPDX-License-Identifier: BSD-2-Clause
5+
//
6+
// This file is part of CEED: http://github.com/ceed
7+
8+
#include <ceed/types.h>
9+
10+
typedef struct {
11+
CeedInt num_comp_u;
12+
CeedInt num_comp_p;
13+
} Context_t;
14+
15+
CEED_QFUNCTION(setup)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
16+
const CeedScalar *weight = in[0], *J = in[1];
17+
CeedScalar *rho = out[0];
18+
for (CeedInt i = 0; i < Q; i++) {
19+
rho[i] = weight[i] * (J[i + Q * 0] * J[i + Q * 3] - J[i + Q * 1] * J[i + Q * 2]);
20+
}
21+
return 0;
22+
}
23+
24+
CEED_QFUNCTION(multi_basis)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
25+
const CeedScalar(*q_data) = (const CeedScalar(*))in[0];
26+
const CeedScalar(*u)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1];
27+
const CeedScalar(*p)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2];
28+
CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
29+
CeedScalar(*q)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[1];
30+
31+
const Context_t *context = (Context_t *)ctx;
32+
33+
CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
34+
CeedScalar mass_u = 0;
35+
for (CeedInt j = 0; j < context->num_comp_u; j++) {
36+
v[j][i] = q_data[i] * u[j][i];
37+
mass_u += v[j][i];
38+
}
39+
for (CeedInt j = 0; j < context->num_comp_p; j++) {
40+
q[j][i] = q_data[i] * p[j][i] + (j + 1) * mass_u;
41+
}
42+
}
43+
return 0;
44+
}

0 commit comments

Comments
 (0)