1 /// @file
2 /// Test creation, action, and destruction for mass matrix operator with multigrid level, non-tensor basis
3 /// \test Test creation, action, and destruction for mass matrix operator with multigrid level, non-tensor basis
4 #include <ceed.h>
5 #include <math.h>
6 #include <stdio.h>
7 #include <stdlib.h>
8
main(int argc,char ** argv)9 int main(int argc, char **argv) {
10 Ceed ceed;
11 CeedElemRestriction elem_restriction_x, elem_restriction_q_data, elem_restriction_u_coarse, elem_restriction_u_fine;
12 CeedBasis basis_x, basis_u;
13 CeedQFunction qf_setup, qf_mass;
14 CeedOperator op_setup, op_mass_coarse, op_mass_fine, op_prolong, op_restrict;
15 CeedVector q_data, x, u_coarse, u_fine, v_coarse, v_fine, p_mult_fine;
16 CeedInt num_elem = 15, p_coarse = 3, p_fine = 5, q = 8;
17 CeedInt num_dofs_x = num_elem + 1, num_dofs_u_coarse = num_elem * (p_coarse - 1) + 1, num_dofs_u_fine = num_elem * (p_fine - 1) + 1;
18 CeedInt ind_u_coarse[num_elem * p_coarse], ind_u_fine[num_elem * p_fine], ind_x[num_elem * 2];
19
20 CeedInit(argv[1], &ceed);
21
22 CeedVectorCreate(ceed, num_dofs_x, &x);
23 {
24 CeedScalar x_array[num_dofs_x];
25
26 for (CeedInt i = 0; i < num_dofs_x; i++) x_array[i] = (CeedScalar)i / (num_dofs_x - 1);
27 CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, x_array);
28 }
29 CeedVectorCreate(ceed, num_dofs_u_fine, &p_mult_fine);
30 CeedVectorCreate(ceed, num_dofs_u_coarse, &u_coarse);
31 CeedVectorCreate(ceed, num_dofs_u_coarse, &v_coarse);
32 CeedVectorCreate(ceed, num_dofs_u_fine, &u_fine);
33 CeedVectorCreate(ceed, num_dofs_u_fine, &v_fine);
34 CeedVectorCreate(ceed, num_elem * q, &q_data);
35
36 // Restrictions
37 for (CeedInt i = 0; i < num_elem; i++) {
38 ind_x[2 * i + 0] = i;
39 ind_x[2 * i + 1] = i + 1;
40 }
41 CeedElemRestrictionCreate(ceed, num_elem, 2, 1, 1, num_dofs_x, CEED_MEM_HOST, CEED_USE_POINTER, ind_x, &elem_restriction_x);
42
43 for (CeedInt i = 0; i < num_elem; i++) {
44 for (CeedInt j = 0; j < p_coarse; j++) {
45 ind_u_coarse[p_coarse * i + j] = i * (p_coarse - 1) + j;
46 }
47 }
48 CeedElemRestrictionCreate(ceed, num_elem, p_coarse, 1, 1, num_dofs_u_coarse, CEED_MEM_HOST, CEED_USE_POINTER, ind_u_coarse,
49 &elem_restriction_u_coarse);
50
51 for (CeedInt i = 0; i < num_elem; i++) {
52 for (CeedInt j = 0; j < p_fine; j++) {
53 ind_u_fine[p_fine * i + j] = i * (p_fine - 1) + j;
54 }
55 }
56 CeedElemRestrictionCreate(ceed, num_elem, p_fine, 1, 1, num_dofs_u_fine, CEED_MEM_HOST, CEED_USE_POINTER, ind_u_fine, &elem_restriction_u_fine);
57
58 CeedInt strides_q_data[3] = {1, q, q};
59 CeedElemRestrictionCreateStrided(ceed, num_elem, q, 1, q * num_elem, strides_q_data, &elem_restriction_q_data);
60
61 // Bases
62 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, 2, q, CEED_GAUSS, &basis_x);
63 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, p_fine, q, CEED_GAUSS, &basis_u);
64
65 // QFunctions
66 CeedQFunctionCreateInteriorByName(ceed, "Mass1DBuild", &qf_setup);
67 CeedQFunctionCreateInteriorByName(ceed, "MassApply", &qf_mass);
68
69 // Operators
70 CeedOperatorCreate(ceed, qf_setup, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_setup);
71
72 CeedOperatorSetField(op_setup, "weights", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE);
73 CeedOperatorSetField(op_setup, "dx", elem_restriction_x, basis_x, CEED_VECTOR_ACTIVE);
74 CeedOperatorSetField(op_setup, "qdata", elem_restriction_q_data, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE);
75
76 CeedOperatorCreate(ceed, qf_mass, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_mass_fine);
77 CeedOperatorSetField(op_mass_fine, "qdata", elem_restriction_q_data, CEED_BASIS_NONE, q_data);
78 CeedOperatorSetField(op_mass_fine, "u", elem_restriction_u_fine, basis_u, CEED_VECTOR_ACTIVE);
79 CeedOperatorSetField(op_mass_fine, "v", elem_restriction_u_fine, basis_u, CEED_VECTOR_ACTIVE);
80
81 CeedOperatorApply(op_setup, x, q_data, CEED_REQUEST_IMMEDIATE);
82
83 // Create multigrid level
84 CeedVectorSetValue(p_mult_fine, 1.0);
85 CeedBasis basis_u_coarse, basis_coarse_to_fine;
86 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, p_coarse, q, CEED_GAUSS, &basis_u_coarse);
87 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, p_coarse, p_fine, CEED_GAUSS_LOBATTO, &basis_coarse_to_fine);
88 const CeedScalar *interp_coarse_to_fine;
89 CeedBasisGetInterp1D(basis_coarse_to_fine, &interp_coarse_to_fine);
90 CeedOperatorMultigridLevelCreateH1(op_mass_fine, p_mult_fine, elem_restriction_u_coarse, basis_u_coarse, interp_coarse_to_fine, &op_mass_coarse,
91 &op_prolong, &op_restrict);
92
93 // Coarse problem
94 CeedVectorSetValue(u_coarse, 1.0);
95 CeedOperatorApply(op_mass_coarse, u_coarse, v_coarse, CEED_REQUEST_IMMEDIATE);
96
97 // Check output
98 {
99 const CeedScalar *v_array;
100 CeedScalar sum = 0.;
101
102 CeedVectorGetArrayRead(v_coarse, CEED_MEM_HOST, &v_array);
103 for (CeedInt i = 0; i < num_dofs_u_coarse; i++) {
104 sum += v_array[i];
105 }
106 CeedVectorRestoreArrayRead(v_coarse, &v_array);
107 if (fabs(sum - 1.) > 1000. * CEED_EPSILON) printf("Computed Area Coarse Grid: %f != True Area: 1.0\n", sum);
108 }
109
110 // Prolong coarse u
111 CeedOperatorApply(op_prolong, u_coarse, u_fine, CEED_REQUEST_IMMEDIATE);
112
113 // Fine problem
114 CeedOperatorApply(op_mass_fine, u_fine, v_fine, CEED_REQUEST_IMMEDIATE);
115
116 // Check output
117 {
118 const CeedScalar *v_array;
119 CeedScalar sum = 0.;
120
121 CeedVectorGetArrayRead(v_fine, CEED_MEM_HOST, &v_array);
122 for (CeedInt i = 0; i < num_dofs_u_fine; i++) {
123 sum += v_array[i];
124 }
125 CeedVectorRestoreArrayRead(v_fine, &v_array);
126 if (fabs(sum - 1.) > 1000. * CEED_EPSILON) printf("Computed Area Fine Grid: %f != True Area: 1.0\n", sum);
127 }
128
129 // Restrict state to coarse grid
130 CeedOperatorApply(op_restrict, v_fine, v_coarse, CEED_REQUEST_IMMEDIATE);
131
132 // Check output
133 {
134 const CeedScalar *v_array;
135 CeedScalar sum = 0.;
136
137 CeedVectorGetArrayRead(v_coarse, CEED_MEM_HOST, &v_array);
138 for (CeedInt i = 0; i < num_dofs_u_coarse; i++) {
139 sum += v_array[i];
140 }
141 CeedVectorRestoreArrayRead(v_coarse, &v_array);
142 if (fabs(sum - 1.) > 1000. * CEED_EPSILON) printf("Computed Area Coarse Grid: %f != True Area: 1.0\n", sum);
143 }
144
145 // Cleanup
146 CeedVectorDestroy(&x);
147 CeedVectorDestroy(&u_coarse);
148 CeedVectorDestroy(&u_fine);
149 CeedVectorDestroy(&v_coarse);
150 CeedVectorDestroy(&v_fine);
151 CeedVectorDestroy(&p_mult_fine);
152 CeedVectorDestroy(&q_data);
153 CeedElemRestrictionDestroy(&elem_restriction_u_coarse);
154 CeedElemRestrictionDestroy(&elem_restriction_u_fine);
155 CeedElemRestrictionDestroy(&elem_restriction_x);
156 CeedElemRestrictionDestroy(&elem_restriction_q_data);
157 CeedBasisDestroy(&basis_u);
158 CeedBasisDestroy(&basis_u_coarse);
159 CeedBasisDestroy(&basis_coarse_to_fine);
160 CeedBasisDestroy(&basis_x);
161 CeedQFunctionDestroy(&qf_setup);
162 CeedQFunctionDestroy(&qf_mass);
163 CeedOperatorDestroy(&op_setup);
164 CeedOperatorDestroy(&op_mass_coarse);
165 CeedOperatorDestroy(&op_mass_fine);
166 CeedOperatorDestroy(&op_prolong);
167 CeedOperatorDestroy(&op_restrict);
168 CeedDestroy(&ceed);
169 return 0;
170 }
171