1 #include "../include/matops.h"
2
3 #include "../include/petscutils.h"
4
5 // -----------------------------------------------------------------------------
6 // Setup apply operator context data
7 // -----------------------------------------------------------------------------
SetupApplyOperatorCtx(MPI_Comm comm,DM dm,Ceed ceed,CeedData ceed_data,Vec X_loc,OperatorApplyContext op_apply_ctx)8 PetscErrorCode SetupApplyOperatorCtx(MPI_Comm comm, DM dm, Ceed ceed, CeedData ceed_data, Vec X_loc, OperatorApplyContext op_apply_ctx) {
9 PetscFunctionBeginUser;
10 op_apply_ctx->comm = comm;
11 op_apply_ctx->dm = dm;
12 op_apply_ctx->X_loc = X_loc;
13 PetscCall(VecDuplicate(X_loc, &op_apply_ctx->Y_loc));
14 op_apply_ctx->x_ceed = ceed_data->x_ceed;
15 op_apply_ctx->y_ceed = ceed_data->y_ceed;
16 op_apply_ctx->op = ceed_data->op_apply;
17 op_apply_ctx->ceed = ceed;
18 PetscFunctionReturn(PETSC_SUCCESS);
19 }
20
21 // -----------------------------------------------------------------------------
22 // Setup error operator context data
23 // -----------------------------------------------------------------------------
SetupErrorOperatorCtx(MPI_Comm comm,DM dm,Ceed ceed,CeedData ceed_data,Vec X_loc,CeedOperator op_error,OperatorApplyContext op_error_ctx)24 PetscErrorCode SetupErrorOperatorCtx(MPI_Comm comm, DM dm, Ceed ceed, CeedData ceed_data, Vec X_loc, CeedOperator op_error,
25 OperatorApplyContext op_error_ctx) {
26 PetscFunctionBeginUser;
27 op_error_ctx->comm = comm;
28 op_error_ctx->dm = dm;
29 op_error_ctx->X_loc = X_loc;
30 PetscCall(VecDuplicate(X_loc, &op_error_ctx->Y_loc));
31 op_error_ctx->x_ceed = ceed_data->x_ceed;
32 op_error_ctx->y_ceed = ceed_data->y_ceed;
33 op_error_ctx->op = op_error;
34 op_error_ctx->ceed = ceed;
35 PetscFunctionReturn(PETSC_SUCCESS);
36 }
37
38 // -----------------------------------------------------------------------------
39 // This function returns the computed diagonal of the operator
40 // -----------------------------------------------------------------------------
MatGetDiag(Mat A,Vec D)41 PetscErrorCode MatGetDiag(Mat A, Vec D) {
42 OperatorApplyContext op_apply_ctx;
43
44 PetscFunctionBeginUser;
45 PetscCall(MatShellGetContext(A, &op_apply_ctx));
46
47 // Compute Diagonal via libCEED
48 PetscMemType mem_type;
49
50 // Place PETSc vector in libCEED vector
51 PetscCall(VecP2C(op_apply_ctx->Y_loc, &mem_type, op_apply_ctx->y_ceed));
52
53 // Compute Diagonal
54 CeedOperatorLinearAssembleDiagonal(op_apply_ctx->op, op_apply_ctx->y_ceed, CEED_REQUEST_IMMEDIATE);
55
56 // Local-to-Global
57 PetscCall(VecC2P(op_apply_ctx->y_ceed, mem_type, op_apply_ctx->Y_loc));
58 PetscCall(VecZeroEntries(D));
59 PetscCall(DMLocalToGlobal(op_apply_ctx->dm, op_apply_ctx->Y_loc, ADD_VALUES, D));
60 PetscFunctionReturn(PETSC_SUCCESS);
61 }
62
63 // -----------------------------------------------------------------------------
64 // This function uses libCEED to compute the action of the Laplacian with Dirichlet boundary conditions
65 // -----------------------------------------------------------------------------
ApplyLocal_Ceed(Vec X,Vec Y,OperatorApplyContext op_apply_ctx)66 PetscErrorCode ApplyLocal_Ceed(Vec X, Vec Y, OperatorApplyContext op_apply_ctx) {
67 PetscMemType x_mem_type, y_mem_type;
68
69 PetscFunctionBeginUser;
70 // Global-to-local
71 PetscCall(DMGlobalToLocal(op_apply_ctx->dm, X, INSERT_VALUES, op_apply_ctx->X_loc));
72
73 // Setup libCEED vectors
74 PetscCall(VecReadP2C(op_apply_ctx->X_loc, &x_mem_type, op_apply_ctx->x_ceed));
75 PetscCall(VecP2C(op_apply_ctx->Y_loc, &y_mem_type, op_apply_ctx->y_ceed));
76
77 // Apply libCEED operator
78 CeedOperatorApply(op_apply_ctx->op, op_apply_ctx->x_ceed, op_apply_ctx->y_ceed, CEED_REQUEST_IMMEDIATE);
79
80 // Restore PETSc vectors
81 PetscCall(VecReadC2P(op_apply_ctx->x_ceed, x_mem_type, op_apply_ctx->X_loc));
82 PetscCall(VecC2P(op_apply_ctx->y_ceed, y_mem_type, op_apply_ctx->Y_loc));
83
84 // Local-to-global
85 PetscCall(VecZeroEntries(Y));
86 PetscCall(DMLocalToGlobal(op_apply_ctx->dm, op_apply_ctx->Y_loc, ADD_VALUES, Y));
87 PetscFunctionReturn(PETSC_SUCCESS);
88 }
89
90 // -----------------------------------------------------------------------------
91 // This function wraps the libCEED operator for a MatShell
92 // -----------------------------------------------------------------------------
MatMult_Ceed(Mat A,Vec X,Vec Y)93 PetscErrorCode MatMult_Ceed(Mat A, Vec X, Vec Y) {
94 OperatorApplyContext op_apply_ctx;
95
96 PetscFunctionBeginUser;
97 PetscCall(MatShellGetContext(A, &op_apply_ctx));
98
99 // libCEED for local action of residual evaluator
100 PetscCall(ApplyLocal_Ceed(X, Y, op_apply_ctx));
101 PetscFunctionReturn(PETSC_SUCCESS);
102 }
103
104 // -----------------------------------------------------------------------------
105 // This function uses libCEED to compute the action of the prolongation operator
106 // -----------------------------------------------------------------------------
MatMult_Prolong(Mat A,Vec X,Vec Y)107 PetscErrorCode MatMult_Prolong(Mat A, Vec X, Vec Y) {
108 ProlongRestrContext pr_restr_ctx;
109 PetscMemType c_mem_type, f_mem_type;
110
111 PetscFunctionBeginUser;
112 PetscCall(MatShellGetContext(A, &pr_restr_ctx));
113
114 // Global-to-local
115 PetscCall(VecZeroEntries(pr_restr_ctx->loc_vec_c));
116 PetscCall(DMGlobalToLocal(pr_restr_ctx->dmc, X, INSERT_VALUES, pr_restr_ctx->loc_vec_c));
117
118 // Setup libCEED vectors
119 PetscCall(VecReadP2C(pr_restr_ctx->loc_vec_c, &c_mem_type, pr_restr_ctx->ceed_vec_c));
120 PetscCall(VecP2C(pr_restr_ctx->loc_vec_f, &f_mem_type, pr_restr_ctx->ceed_vec_f));
121
122 // Apply libCEED operator
123 CeedOperatorApply(pr_restr_ctx->op_prolong, pr_restr_ctx->ceed_vec_c, pr_restr_ctx->ceed_vec_f, CEED_REQUEST_IMMEDIATE);
124
125 // Restore PETSc vectors
126 PetscCall(VecReadC2P(pr_restr_ctx->ceed_vec_c, c_mem_type, pr_restr_ctx->loc_vec_c));
127 PetscCall(VecC2P(pr_restr_ctx->ceed_vec_f, f_mem_type, pr_restr_ctx->loc_vec_f));
128
129 // Multiplicity
130 PetscCall(VecPointwiseMult(pr_restr_ctx->loc_vec_f, pr_restr_ctx->loc_vec_f, pr_restr_ctx->mult_vec));
131
132 // Local-to-global
133 PetscCall(VecZeroEntries(Y));
134 PetscCall(DMLocalToGlobal(pr_restr_ctx->dmf, pr_restr_ctx->loc_vec_f, ADD_VALUES, Y));
135 PetscFunctionReturn(PETSC_SUCCESS);
136 }
137
138 // -----------------------------------------------------------------------------
139 // This function uses libCEED to compute the action of the restriction operator
140 // -----------------------------------------------------------------------------
MatMult_Restrict(Mat A,Vec X,Vec Y)141 PetscErrorCode MatMult_Restrict(Mat A, Vec X, Vec Y) {
142 ProlongRestrContext pr_restr_ctx;
143 PetscMemType c_mem_type, f_mem_type;
144
145 PetscFunctionBeginUser;
146 PetscCall(MatShellGetContext(A, &pr_restr_ctx));
147
148 // Global-to-local
149 PetscCall(VecZeroEntries(pr_restr_ctx->loc_vec_f));
150 PetscCall(DMGlobalToLocal(pr_restr_ctx->dmf, X, INSERT_VALUES, pr_restr_ctx->loc_vec_f));
151
152 // Multiplicity
153 PetscCall(VecPointwiseMult(pr_restr_ctx->loc_vec_f, pr_restr_ctx->loc_vec_f, pr_restr_ctx->mult_vec));
154
155 // Setup libCEED vectors
156 PetscCall(VecReadP2C(pr_restr_ctx->loc_vec_f, &f_mem_type, pr_restr_ctx->ceed_vec_f));
157 PetscCall(VecP2C(pr_restr_ctx->loc_vec_c, &c_mem_type, pr_restr_ctx->ceed_vec_c));
158
159 // Apply CEED operator
160 CeedOperatorApply(pr_restr_ctx->op_restrict, pr_restr_ctx->ceed_vec_f, pr_restr_ctx->ceed_vec_c, CEED_REQUEST_IMMEDIATE);
161
162 // Restore PETSc vectors
163 PetscCall(VecReadC2P(pr_restr_ctx->ceed_vec_f, f_mem_type, pr_restr_ctx->loc_vec_f));
164 PetscCall(VecC2P(pr_restr_ctx->ceed_vec_c, c_mem_type, pr_restr_ctx->loc_vec_c));
165
166 // Local-to-global
167 PetscCall(VecZeroEntries(Y));
168 PetscCall(DMLocalToGlobal(pr_restr_ctx->dmc, pr_restr_ctx->loc_vec_c, ADD_VALUES, Y));
169 PetscFunctionReturn(PETSC_SUCCESS);
170 }
171
172 // -----------------------------------------------------------------------------
173 // This function calculates the error in the final solution
174 // -----------------------------------------------------------------------------
ComputeL2Error(Vec X,PetscScalar * l2_error,OperatorApplyContext op_error_ctx)175 PetscErrorCode ComputeL2Error(Vec X, PetscScalar *l2_error, OperatorApplyContext op_error_ctx) {
176 Vec E;
177
178 PetscFunctionBeginUser;
179 PetscCall(VecDuplicate(X, &E));
180 PetscCall(ApplyLocal_Ceed(X, E, op_error_ctx));
181 PetscCall(VecViewFromOptions(E, NULL, "-error_view"));
182 PetscScalar error_sq = 1.0;
183 PetscCall(VecSum(E, &error_sq));
184 *l2_error = sqrt(error_sq);
185 PetscCall(VecDestroy(&E));
186 PetscFunctionReturn(PETSC_SUCCESS);
187 }
188
189 // -----------------------------------------------------------------------------
190