-
Notifications
You must be signed in to change notification settings - Fork 45
/
t533-operator.c
143 lines (122 loc) · 5.4 KB
/
t533-operator.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
/// @file
/// Test assembly of mass matrix operator diagonal
/// \test Test assembly of mass matrix operator diagonal
#include <ceed.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "t510-operator.h"
int main(int argc, char **argv) {
Ceed ceed;
CeedElemRestriction elem_restriction_x, elem_restriction_u, elem_restriction_q_data;
CeedBasis basis_x, basis_u;
CeedQFunction qf_setup, qf_mass;
CeedOperator op_setup, op_mass;
CeedVector q_data, x, assembled, u, v;
CeedInt num_elem = 6, p = 3, q = 4, dim = 2;
CeedInt nx = 3, ny = 2;
CeedInt num_dofs = (nx * 2 + 1) * (ny * 2 + 1), num_qpts = num_elem * q * q;
CeedInt ind_x[num_elem * p * p];
CeedScalar assembled_true[num_dofs];
CeedInit(argv[1], &ceed);
// Vectors
CeedVectorCreate(ceed, dim * num_dofs, &x);
{
CeedScalar x_array[dim * num_dofs];
for (CeedInt i = 0; i < nx * 2 + 1; i++)
for (CeedInt j = 0; j < ny * 2 + 1; j++) {
x_array[i + j * (nx * 2 + 1) + 0 * num_dofs] = (CeedScalar)i / (2 * nx);
x_array[i + j * (nx * 2 + 1) + 1 * num_dofs] = (CeedScalar)j / (2 * ny);
}
CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, x_array);
}
CeedVectorCreate(ceed, num_dofs, &u);
CeedVectorCreate(ceed, num_dofs, &v);
CeedVectorCreate(ceed, num_qpts, &q_data);
// Restrictions
for (CeedInt i = 0; i < num_elem; i++) {
CeedInt col, row, offset;
col = i % nx;
row = i / nx;
offset = col * (p - 1) + row * (nx * 2 + 1) * (p - 1);
for (CeedInt j = 0; j < p; j++)
for (CeedInt k = 0; k < p; k++) ind_x[p * (p * i + k) + j] = offset + k * (nx * 2 + 1) + j;
}
CeedElemRestrictionCreate(ceed, num_elem, p * p, dim, num_dofs, dim * num_dofs, CEED_MEM_HOST, CEED_USE_POINTER, ind_x, &elem_restriction_x);
CeedElemRestrictionCreate(ceed, num_elem, p * p, 1, 1, num_dofs, CEED_MEM_HOST, CEED_USE_POINTER, ind_x, &elem_restriction_u);
CeedInt strides_q_data[3] = {1, q * q, q * q};
CeedElemRestrictionCreateStrided(ceed, num_elem, q * q, 1, num_qpts, strides_q_data, &elem_restriction_q_data);
// Bases
CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, p, q, CEED_GAUSS, &basis_x);
CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, p, q, CEED_GAUSS, &basis_u);
// QFunctions
CeedQFunctionCreateInterior(ceed, 1, setup, setup_loc, &qf_setup);
CeedQFunctionAddInput(qf_setup, "weight", 1, CEED_EVAL_WEIGHT);
CeedQFunctionAddInput(qf_setup, "dx", dim * dim, CEED_EVAL_GRAD);
CeedQFunctionAddOutput(qf_setup, "rho", 1, CEED_EVAL_NONE);
CeedQFunctionCreateInterior(ceed, 1, mass, mass_loc, &qf_mass);
CeedQFunctionAddInput(qf_mass, "rho", 1, CEED_EVAL_NONE);
CeedQFunctionAddInput(qf_mass, "u", 1, CEED_EVAL_INTERP);
CeedQFunctionAddOutput(qf_mass, "v", 1, CEED_EVAL_INTERP);
// Operators
CeedOperatorCreate(ceed, qf_setup, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_setup);
CeedOperatorSetField(op_setup, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE);
CeedOperatorSetField(op_setup, "dx", elem_restriction_x, basis_x, CEED_VECTOR_ACTIVE);
CeedOperatorSetField(op_setup, "rho", elem_restriction_q_data, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE);
CeedOperatorCreate(ceed, qf_mass, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_mass);
CeedOperatorSetField(op_mass, "rho", elem_restriction_q_data, CEED_BASIS_NONE, q_data);
CeedOperatorSetField(op_mass, "u", elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE);
CeedOperatorSetField(op_mass, "v", elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE);
// Apply Setup Operator
CeedOperatorApply(op_setup, x, q_data, CEED_REQUEST_IMMEDIATE);
// Assemble diagonal
CeedVectorCreate(ceed, num_dofs, &assembled);
CeedOperatorLinearAssembleDiagonal(op_mass, assembled, CEED_REQUEST_IMMEDIATE);
// Manually assemble diagonal
CeedVectorSetValue(u, 0.0);
for (CeedInt i = 0; i < num_dofs; i++) {
CeedScalar *u_array;
const CeedScalar *v_array;
// Set input
CeedVectorGetArray(u, CEED_MEM_HOST, &u_array);
u_array[i] = 1.0;
if (i) u_array[i - 1] = 0.0;
CeedVectorRestoreArray(u, &u_array);
// Compute diag entry for DoF i
CeedOperatorApply(op_mass, u, v, CEED_REQUEST_IMMEDIATE);
// Retrieve entry
CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array);
assembled_true[i] = v_array[i];
CeedVectorRestoreArrayRead(v, &v_array);
}
// Check output
{
const CeedScalar *assembled_array;
CeedVectorGetArrayRead(assembled, CEED_MEM_HOST, &assembled_array);
for (CeedInt i = 0; i < num_dofs; i++) {
if (fabs(assembled_array[i] - assembled_true[i]) > 100. * CEED_EPSILON) {
// LCOV_EXCL_START
printf("[%" CeedInt_FMT "] Error in assembly: %f != %f\n", i, assembled_array[i], assembled_true[i]);
// LCOV_EXCL_STOP
}
}
CeedVectorRestoreArrayRead(assembled, &assembled_array);
}
// Cleanup
CeedVectorDestroy(&x);
CeedVectorDestroy(&assembled);
CeedVectorDestroy(&q_data);
CeedVectorDestroy(&u);
CeedVectorDestroy(&v);
CeedElemRestrictionDestroy(&elem_restriction_u);
CeedElemRestrictionDestroy(&elem_restriction_x);
CeedElemRestrictionDestroy(&elem_restriction_q_data);
CeedBasisDestroy(&basis_u);
CeedBasisDestroy(&basis_x);
CeedQFunctionDestroy(&qf_setup);
CeedQFunctionDestroy(&qf_mass);
CeedOperatorDestroy(&op_setup);
CeedOperatorDestroy(&op_mass);
CeedDestroy(&ceed);
return 0;
}