1 // SPDX-FileCopyrightText: Copyright (c) 2017-2024, HONEE contributors.
2 // SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause
3 /// @file
4 /// Functions for setting up and projecting the divergence of the diffusive flux
5
6 #include "../qfunctions/diff_flux_projection.h"
7
8 #include <petscdmplex.h>
9
10 #include <navierstokes.h>
11
12 const char *const DivDiffFluxProjectionMethods[] = {"NONE", "DIRECT", "INDIRECT", "DivDiffFluxProjectionMethod", "DIV_DIFF_FLUX_PROJ_", NULL};
13
14 /**
15 @brief Create `DivDiffFluxProjectionData` for solution DM in `honee`
16
17 @param[in] honee `Honee` context
18 @param[in] divFdiffproj_method Method used to perform the divergence of diffusive flux method (can be `DIV_DIFF_FLUX_PROJ_NONE`)
19 @param[in] num_diff_flux_comps Number of components that makes up the diffusive flux (e.g. 1 for scalar advection-diffusion)
20 @param[out] diff_flux_proj The `DivDiffFluxProjectionData` object created, or set to `NULL` if `divFdiffproj_method = DIV_DIFF_FLUX_PROJ_NONE`
21 **/
DivDiffFluxProjectionCreate(Honee honee,DivDiffFluxProjectionMethod divFdiffproj_method,PetscInt num_diff_flux_comps,DivDiffFluxProjectionData * diff_flux_proj)22 PetscErrorCode DivDiffFluxProjectionCreate(Honee honee, DivDiffFluxProjectionMethod divFdiffproj_method, PetscInt num_diff_flux_comps,
23 DivDiffFluxProjectionData *diff_flux_proj) {
24 PetscInt height = 0, dm_field = 0, dim, degree = honee->app_ctx->degree, q_extra = honee->app_ctx->q_extra;
25 DivDiffFluxProjectionData diff_flux_proj_;
26 NodalProjectionData projection;
27
28 PetscFunctionBeginUser;
29 if (divFdiffproj_method == DIV_DIFF_FLUX_PROJ_NONE) {
30 *diff_flux_proj = NULL;
31 PetscFunctionReturn(PETSC_SUCCESS);
32 }
33 PetscCall(PetscNew(&projection));
34 PetscCall(PetscNew(&diff_flux_proj_));
35 *diff_flux_proj_ = (struct DivDiffFluxProjectionData_){
36 .method = divFdiffproj_method,
37 .num_diff_flux_comps = num_diff_flux_comps,
38 .projection = projection,
39 };
40
41 PetscCall(DMClone(honee->dm, &projection->dm));
42 PetscCall(DMSetMatrixPreallocateSkip(projection->dm, PETSC_TRUE));
43 PetscCall(DMGetDimension(projection->dm, &dim));
44 switch (diff_flux_proj_->method) {
45 case DIV_DIFF_FLUX_PROJ_DIRECT: {
46 projection->num_comp = diff_flux_proj_->num_diff_flux_comps;
47 PetscCall(PetscObjectSetName((PetscObject)projection->dm, "DivDiffFluxProj"));
48 PetscCall(DMSetupByOrder_FEM(PETSC_TRUE, PETSC_TRUE, degree, 1, q_extra, 1, &projection->num_comp, projection->dm));
49
50 PetscCall(DMPlexCeedElemRestrictionCreate(honee->ceed, projection->dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, height, dm_field,
51 &diff_flux_proj_->elem_restr_div_diff_flux));
52 PetscCallCeed(honee->ceed,
53 CeedElemRestrictionCreateVector(diff_flux_proj_->elem_restr_div_diff_flux, &diff_flux_proj_->div_diff_flux_ceed, NULL));
54 PetscCall(DMPlexCeedBasisCreate(honee->ceed, projection->dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, height, dm_field,
55 &diff_flux_proj_->basis_div_diff_flux));
56 diff_flux_proj_->eval_mode_div_diff_flux = CEED_EVAL_INTERP;
57
58 { // Create face labels on projection->dm for boundary integrals
59 DMLabel face_sets_label;
60 PetscInt num_face_set_values, *face_set_values;
61
62 PetscCall(DMGetLabel(honee->dm, "Face Sets", &face_sets_label));
63 PetscCall(DMLabelCreateGlobalValueArray(honee->dm, face_sets_label, &num_face_set_values, &face_set_values));
64 for (PetscInt f = 0; f < num_face_set_values; f++) {
65 DMLabel face_orientation_label;
66 char *face_orientation_label_name;
67
68 PetscCall(DMPlexCreateFaceLabel(honee->dm, face_set_values[f], &face_orientation_label_name));
69 PetscCall(DMGetLabel(honee->dm, face_orientation_label_name, &face_orientation_label));
70 PetscCall(DMAddLabel(projection->dm, face_orientation_label));
71 PetscCall(PetscFree(face_orientation_label_name));
72 }
73 PetscCall(PetscFree(face_set_values));
74 }
75 } break;
76 case DIV_DIFF_FLUX_PROJ_INDIRECT: {
77 projection->num_comp = diff_flux_proj_->num_diff_flux_comps * dim;
78 PetscCall(PetscObjectSetName((PetscObject)projection->dm, "DiffFluxProj"));
79 PetscCall(DMSetupByOrder_FEM(PETSC_TRUE, PETSC_TRUE, degree, 1, q_extra, 1, &projection->num_comp, projection->dm));
80
81 PetscCall(DMPlexCeedElemRestrictionQDataCreate(honee->ceed, projection->dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, height,
82 diff_flux_proj_->num_diff_flux_comps, &diff_flux_proj_->elem_restr_div_diff_flux));
83 PetscCallCeed(honee->ceed,
84 CeedElemRestrictionCreateVector(diff_flux_proj_->elem_restr_div_diff_flux, &diff_flux_proj_->div_diff_flux_ceed, NULL));
85 diff_flux_proj_->basis_div_diff_flux = CEED_BASIS_NONE;
86 diff_flux_proj_->eval_mode_div_diff_flux = CEED_EVAL_NONE;
87 } break;
88 case DIV_DIFF_FLUX_PROJ_NONE:
89 SETERRQ(honee->comm, PETSC_ERR_ARG_WRONG, "Should not reach here with div_diff_flux_projection_method %s",
90 DivDiffFluxProjectionMethods[divFdiffproj_method]);
91 break;
92 }
93 *diff_flux_proj = diff_flux_proj_;
94 PetscFunctionReturn(PETSC_SUCCESS);
95 };
96
97 /**
98 @brief Return the objects required for the Divergence of Diffusive flux to be read by a `CeedOperator`
99
100 @param[in] diff_flux_proj Projection object
101 @param[out] elem_restr Element restriction for the divergence of diffusive flux, or `NULL`
102 @param[out] basis Basis for the divergence of diffusive flux, or `NULL`
103 @param[out] vector Vector for the divergence of diffusive flux, or `NULL`
104 @param[out] eval_mode Eval mode for the divergence of diffusive flux, or `NULL`
105 **/
DivDiffFluxProjectionGetOperatorFieldData(DivDiffFluxProjectionData diff_flux_proj,CeedElemRestriction * elem_restr,CeedBasis * basis,CeedVector * vector,CeedEvalMode * eval_mode)106 PetscErrorCode DivDiffFluxProjectionGetOperatorFieldData(DivDiffFluxProjectionData diff_flux_proj, CeedElemRestriction *elem_restr, CeedBasis *basis,
107 CeedVector *vector, CeedEvalMode *eval_mode) {
108 Ceed ceed = CeedVectorReturnCeed(diff_flux_proj->div_diff_flux_ceed);
109
110 PetscFunctionBeginUser;
111 if (elem_restr) PetscCallCeed(ceed, CeedElemRestrictionReferenceCopy(diff_flux_proj->elem_restr_div_diff_flux, elem_restr));
112 if (basis) PetscCallCeed(ceed, CeedBasisReferenceCopy(diff_flux_proj->basis_div_diff_flux, basis));
113 if (vector) PetscCallCeed(ceed, CeedVectorReferenceCopy(diff_flux_proj->div_diff_flux_ceed, vector));
114 if (eval_mode) *eval_mode = diff_flux_proj->eval_mode_div_diff_flux;
115 PetscFunctionReturn(PETSC_SUCCESS);
116 }
117
118 /**
119 @brief Setup direct projection of divergence of diffusive flux
120
121 @param[in] honee `Honee` context
122 @param[in,out] diff_flux_proj Flux projection object to setup
123 **/
DivDiffFluxProjectionSetup_Direct(Honee honee,DivDiffFluxProjectionData diff_flux_proj)124 static PetscErrorCode DivDiffFluxProjectionSetup_Direct(Honee honee, DivDiffFluxProjectionData diff_flux_proj) {
125 Ceed ceed = honee->ceed;
126 NodalProjectionData projection = diff_flux_proj->projection;
127 MPI_Comm comm = PetscObjectComm((PetscObject)projection->dm);
128
129 PetscFunctionBeginUser;
130 { // Create Projection RHS OperatorApplyContext
131 CeedOperator op_rhs;
132
133 PetscCheck(diff_flux_proj->CreateRHSOperator_Direct, comm, PETSC_ERR_ARG_WRONGSTATE,
134 "Must define CreateRHSOperator_Direct to use indirect div_diff_flux projection");
135 PetscCall(diff_flux_proj->CreateRHSOperator_Direct(honee, diff_flux_proj, &op_rhs));
136 PetscCall(DMCreateLocalVector(projection->dm, &diff_flux_proj->DivDiffFlux_loc));
137 diff_flux_proj->ceed_vec_has_array = PETSC_FALSE;
138 PetscCall(OperatorApplyContextCreate(honee->dm, projection->dm, ceed, op_rhs, NULL, NULL, NULL, diff_flux_proj->DivDiffFlux_loc,
139 &projection->l2_rhs_ctx));
140 PetscCallCeed(ceed, CeedOperatorDestroy(&op_rhs));
141 }
142
143 { // -- Build Mass operator
144 CeedQFunction qf_mass;
145 CeedOperator op_mass;
146 CeedBasis basis_div_diff_flux = NULL;
147 CeedElemRestriction elem_restr_div_diff_flux_volume = NULL, elem_restr_qd;
148 CeedVector q_data;
149 CeedInt q_data_size;
150
151 PetscCall(DivDiffFluxProjectionGetOperatorFieldData(diff_flux_proj, &elem_restr_div_diff_flux_volume, &basis_div_diff_flux, NULL, NULL));
152 PetscCall(QDataGet(ceed, projection->dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, &elem_restr_qd, &q_data, &q_data_size));
153
154 PetscCall(HoneeMassQFunctionCreate(ceed, projection->num_comp, q_data_size, &qf_mass));
155 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_mass, NULL, NULL, &op_mass));
156 PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "u", elem_restr_div_diff_flux_volume, basis_div_diff_flux, CEED_VECTOR_ACTIVE));
157 PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "qdata", elem_restr_qd, CEED_BASIS_NONE, q_data));
158 PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "v", elem_restr_div_diff_flux_volume, basis_div_diff_flux, CEED_VECTOR_ACTIVE));
159
160 { // -- Setup KSP for L^2 projection
161 Mat mat_mass;
162
163 PetscCall(MatCreateCeed(projection->dm, projection->dm, op_mass, NULL, &mat_mass));
164
165 PetscCall(KSPCreate(comm, &projection->ksp));
166 PetscCall(KSPSetOptionsPrefix(projection->ksp, "div_diff_flux_projection_"));
167 { // lumped by default
168 PC pc;
169 PetscCall(KSPGetPC(projection->ksp, &pc));
170 PetscCall(PCSetType(pc, PCJACOBI));
171 PetscCall(PCJacobiSetType(pc, PC_JACOBI_ROWSUM));
172 PetscCall(KSPSetType(projection->ksp, KSPPREONLY));
173 }
174 PetscCall(KSPSetFromOptions_WithMatCeed(projection->ksp, mat_mass));
175 PetscCall(MatDestroy(&mat_mass));
176 }
177
178 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_div_diff_flux_volume));
179 PetscCallCeed(ceed, CeedBasisDestroy(&basis_div_diff_flux));
180 PetscCallCeed(ceed, CeedVectorDestroy(&q_data));
181 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd));
182 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_mass));
183 PetscCallCeed(ceed, CeedOperatorDestroy(&op_mass));
184 }
185 PetscFunctionReturn(PETSC_SUCCESS);
186 }
187
188 /**
189 @brief Setup indirect projection of divergence of diffusive flux
190
191 @param[in] honee `Honee` context
192 @param[in,out] diff_flux_proj Flux projection object to setup
193 **/
DivDiffFluxProjectionSetup_Indirect(Honee honee,DivDiffFluxProjectionData diff_flux_proj)194 static PetscErrorCode DivDiffFluxProjectionSetup_Indirect(Honee honee, DivDiffFluxProjectionData diff_flux_proj) {
195 Ceed ceed = honee->ceed;
196 NodalProjectionData projection = diff_flux_proj->projection;
197 CeedBasis basis_diff_flux;
198 CeedElemRestriction elem_restr_diff_flux, elem_restr_qd;
199 CeedVector q_data;
200 CeedInt q_data_size;
201 MPI_Comm comm = PetscObjectComm((PetscObject)projection->dm);
202
203 PetscFunctionBeginUser;
204 {
205 PetscInt height = 0, dm_field = 0;
206
207 PetscCall(DMPlexCeedElemRestrictionCreate(ceed, projection->dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, height, dm_field, &elem_restr_diff_flux));
208 PetscCall(DMPlexCeedBasisCreate(ceed, projection->dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, height, dm_field, &basis_diff_flux));
209 PetscCall(QDataGet(ceed, projection->dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, &elem_restr_qd, &q_data, &q_data_size));
210 }
211
212 {
213 CeedOperator op_rhs;
214
215 PetscCheck(diff_flux_proj->CreateRHSOperator_Indirect, comm, PETSC_ERR_ARG_WRONGSTATE,
216 "Must define CreateRHSOperator_Indirect to use indirect div_diff_flux projection");
217 PetscCall(diff_flux_proj->CreateRHSOperator_Indirect(honee, diff_flux_proj, &op_rhs));
218 PetscCall(OperatorApplyContextCreate(honee->dm, projection->dm, ceed, op_rhs, NULL, NULL, NULL, NULL, &projection->l2_rhs_ctx));
219 PetscCallCeed(ceed, CeedOperatorDestroy(&op_rhs));
220 }
221
222 { // -- Build Mass operator
223 CeedQFunction qf_mass;
224 CeedOperator op_mass;
225
226 PetscCall(HoneeMassQFunctionCreate(ceed, projection->num_comp, q_data_size, &qf_mass));
227 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_mass, NULL, NULL, &op_mass));
228 PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "u", elem_restr_diff_flux, basis_diff_flux, CEED_VECTOR_ACTIVE));
229 PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "qdata", elem_restr_qd, CEED_BASIS_NONE, q_data));
230 PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "v", elem_restr_diff_flux, basis_diff_flux, CEED_VECTOR_ACTIVE));
231
232 { // -- Setup KSP for L^2 projection
233 Mat mat_mass;
234
235 PetscCall(MatCreateCeed(projection->dm, projection->dm, op_mass, NULL, &mat_mass));
236
237 PetscCall(KSPCreate(comm, &projection->ksp));
238 PetscCall(KSPSetOptionsPrefix(projection->ksp, "div_diff_flux_projection_"));
239 { // lumped by default
240 PC pc;
241 PetscCall(KSPGetPC(projection->ksp, &pc));
242 PetscCall(PCSetType(pc, PCJACOBI));
243 PetscCall(PCJacobiSetType(pc, PC_JACOBI_ROWSUM));
244 PetscCall(KSPSetType(projection->ksp, KSPPREONLY));
245 }
246 PetscCall(KSPSetFromOptions_WithMatCeed(projection->ksp, mat_mass));
247 PetscCall(MatDestroy(&mat_mass));
248 }
249 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_mass));
250 PetscCallCeed(ceed, CeedOperatorDestroy(&op_mass));
251 }
252
253 { // Create OperatorApplyContext to calculate divergence at quadrature points
254 CeedQFunction qf_calc_divergence = NULL;
255 CeedOperator op_calc_divergence;
256 CeedElemRestriction elem_restr_div_diff_flux = NULL;
257 PetscInt dim;
258
259 PetscCall(DMGetDimension(projection->dm, &dim));
260 PetscCall(DivDiffFluxProjectionGetOperatorFieldData(diff_flux_proj, &elem_restr_div_diff_flux, NULL, NULL, NULL));
261
262 switch (dim) {
263 case 2:
264 switch (diff_flux_proj->num_diff_flux_comps) {
265 case 1:
266 PetscCallCeed(ceed,
267 CeedQFunctionCreateInterior(ceed, 1, ComputeDivDiffusiveFlux2D_1, ComputeDivDiffusiveFlux2D_1_loc, &qf_calc_divergence));
268 break;
269 }
270 break;
271 case 3:
272 switch (diff_flux_proj->num_diff_flux_comps) {
273 case 1:
274 PetscCallCeed(ceed,
275 CeedQFunctionCreateInterior(ceed, 1, ComputeDivDiffusiveFlux3D_1, ComputeDivDiffusiveFlux3D_1_loc, &qf_calc_divergence));
276 break;
277 case 4:
278 PetscCallCeed(ceed,
279 CeedQFunctionCreateInterior(ceed, 1, ComputeDivDiffusiveFlux3D_4, ComputeDivDiffusiveFlux3D_4_loc, &qf_calc_divergence));
280 break;
281 }
282 break;
283 }
284 PetscCheck(qf_calc_divergence, comm, PETSC_ERR_SUP,
285 "QFunction for calculating divergence of diffusive flux does not exist for"
286 " %" PetscInt_FMT " dimensional grid and %" PetscInt_FMT
287 " number of components.\nA new qfunction can be easily added; see source code for pattern.",
288 dim, diff_flux_proj->num_diff_flux_comps);
289
290 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_calc_divergence, "Grad F_diff", projection->num_comp * dim, CEED_EVAL_GRAD));
291 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_calc_divergence, "qdata", q_data_size, CEED_EVAL_NONE));
292 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_calc_divergence, "Div F_diff", diff_flux_proj->num_diff_flux_comps, CEED_EVAL_NONE));
293
294 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_calc_divergence, NULL, NULL, &op_calc_divergence));
295 PetscCallCeed(ceed, CeedOperatorSetField(op_calc_divergence, "Grad F_diff", elem_restr_diff_flux, basis_diff_flux, CEED_VECTOR_ACTIVE));
296 PetscCallCeed(ceed, CeedOperatorSetField(op_calc_divergence, "qdata", elem_restr_qd, CEED_BASIS_NONE, q_data));
297 PetscCallCeed(ceed, CeedOperatorSetField(op_calc_divergence, "Div F_diff", elem_restr_div_diff_flux, CEED_BASIS_NONE,
298 diff_flux_proj->div_diff_flux_ceed));
299
300 PetscCall(OperatorApplyContextCreate(projection->dm, NULL, ceed, op_calc_divergence, NULL, CEED_VECTOR_NONE, NULL, NULL,
301 &diff_flux_proj->calc_div_diff_flux));
302
303 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_div_diff_flux));
304 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_calc_divergence));
305 PetscCallCeed(ceed, CeedOperatorDestroy(&op_calc_divergence));
306 }
307 PetscCallCeed(ceed, CeedBasisDestroy(&basis_diff_flux));
308 PetscCallCeed(ceed, CeedVectorDestroy(&q_data));
309 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd));
310 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_diff_flux));
311 PetscFunctionReturn(PETSC_SUCCESS);
312 }
313
314 /**
315 @brief Setup projection of divergence of diffusive flux
316
317 @param[in] honee `Honee` context
318 @param[in,out] diff_flux_proj Flux projection object to setup
319 **/
DivDiffFluxProjectionSetup(Honee honee,DivDiffFluxProjectionData diff_flux_proj)320 PetscErrorCode DivDiffFluxProjectionSetup(Honee honee, DivDiffFluxProjectionData diff_flux_proj) {
321 PetscFunctionBeginUser;
322 switch (honee->app_ctx->divFdiffproj_method) {
323 case DIV_DIFF_FLUX_PROJ_DIRECT:
324 PetscCall(DivDiffFluxProjectionSetup_Direct(honee, diff_flux_proj));
325 break;
326 case DIV_DIFF_FLUX_PROJ_INDIRECT:
327 PetscCall(DivDiffFluxProjectionSetup_Indirect(honee, diff_flux_proj));
328 break;
329 case DIV_DIFF_FLUX_PROJ_NONE:
330 SETERRQ(PetscObjectComm((PetscObject)honee->dm), PETSC_ERR_ARG_WRONG, "Should not reach here with div_diff_flux_projection_method %s",
331 DivDiffFluxProjectionMethods[honee->app_ctx->divFdiffproj_method]);
332 break;
333 }
334 PetscFunctionReturn(PETSC_SUCCESS);
335 }
336
337 /**
338 @brief Project the divergence of diffusive flux
339
340 This implicitly sets the `CeedVector` input (`div_diff_flux_ceed`) to the divergence of diffusive flux.
341
342 @param[in] diff_flux_proj `NodalProjectionData` for the projection
343 @param[in] Q_loc Localized solution vector
344 **/
DivDiffFluxProjectionApply(DivDiffFluxProjectionData diff_flux_proj,Vec Q_loc)345 PetscErrorCode DivDiffFluxProjectionApply(DivDiffFluxProjectionData diff_flux_proj, Vec Q_loc) {
346 NodalProjectionData projection = diff_flux_proj->projection;
347
348 PetscFunctionBeginUser;
349 PetscCall(PetscLogEventBegin(HONEE_DivDiffFluxProjection, Q_loc, 0, 0, 0));
350 switch (diff_flux_proj->method) {
351 case DIV_DIFF_FLUX_PROJ_DIRECT: {
352 Vec DivDiffFlux, RHS;
353
354 PetscCall(DMGetGlobalVector(projection->dm, &DivDiffFlux));
355 PetscCall(DMGetGlobalVector(projection->dm, &RHS));
356 if (diff_flux_proj->ceed_vec_has_array) {
357 PetscCall(VecReadCeedToPetsc(diff_flux_proj->div_diff_flux_ceed, diff_flux_proj->DivDiffFlux_memtype, diff_flux_proj->DivDiffFlux_loc));
358 diff_flux_proj->ceed_vec_has_array = PETSC_FALSE;
359 }
360 PetscCall(ApplyCeedOperatorLocalToGlobal(Q_loc, RHS, projection->l2_rhs_ctx));
361 PetscCall(VecViewFromOptions(DivDiffFlux, NULL, "-div_diff_flux_projection_rhs_view"));
362
363 {
364 // Run PCApply manually if using ksp_type preonly -pc_type jacobi
365 // This is to avoid an AllReduce call in KSPSolve_Preonly, which causes significant slowdowns for lumped mass matrix solves.
366 // See https://gitlab.com/petsc/petsc/-/merge_requests/8048 for more details and a possible fix
367 PC pc;
368 PetscBool ispreonly, isjacobi;
369 PetscCall(KSPGetPC(projection->ksp, &pc));
370 PetscCall(PetscObjectTypeCompare((PetscObject)projection->ksp, KSPPREONLY, &ispreonly));
371 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCJACOBI, &isjacobi));
372 if (ispreonly && isjacobi) PetscCall(PCApply(pc, RHS, DivDiffFlux));
373 else PetscCall(KSPSolve(projection->ksp, RHS, DivDiffFlux));
374 }
375 PetscCall(VecViewFromOptions(DivDiffFlux, NULL, "-div_diff_flux_projection_view"));
376
377 PetscCall(DMGlobalToLocal(projection->dm, DivDiffFlux, INSERT_VALUES, diff_flux_proj->DivDiffFlux_loc));
378 PetscCall(VecReadPetscToCeed(diff_flux_proj->DivDiffFlux_loc, &diff_flux_proj->DivDiffFlux_memtype, diff_flux_proj->div_diff_flux_ceed));
379 diff_flux_proj->ceed_vec_has_array = PETSC_TRUE;
380
381 PetscCall(DMRestoreGlobalVector(projection->dm, &RHS));
382 PetscCall(DMRestoreGlobalVector(projection->dm, &DivDiffFlux));
383 break;
384 }
385 case DIV_DIFF_FLUX_PROJ_INDIRECT: {
386 Vec DiffFlux, RHS;
387
388 PetscCall(DMGetGlobalVector(projection->dm, &DiffFlux));
389 PetscCall(DMGetGlobalVector(projection->dm, &RHS));
390 PetscCall(ApplyCeedOperatorLocalToGlobal(Q_loc, RHS, projection->l2_rhs_ctx));
391 PetscCall(VecViewFromOptions(DiffFlux, NULL, "-div_diff_flux_projection_rhs_view"));
392
393 {
394 // Run PCApply manually if using -ksp_type preonly -pc_type jacobi
395 // This is to avoid an AllReduce call in KSPSolve_Preonly, which causes significant slowdowns for lumped mass matrix solves.
396 // See https://gitlab.com/petsc/petsc/-/merge_requests/8048 for more details and a possible fix
397 PC pc;
398 PetscBool ispreonly, isjacobi;
399 PetscCall(KSPGetPC(projection->ksp, &pc));
400 PetscCall(PetscObjectTypeCompare((PetscObject)projection->ksp, KSPPREONLY, &ispreonly));
401 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCJACOBI, &isjacobi));
402 if (ispreonly && isjacobi) PetscCall(PCApply(pc, RHS, DiffFlux));
403 else PetscCall(KSPSolve(projection->ksp, RHS, DiffFlux));
404 }
405 PetscCall(VecViewFromOptions(DiffFlux, NULL, "-div_diff_flux_projection_view"));
406
407 PetscCall(ApplyCeedOperatorGlobalToLocal(DiffFlux, NULL, diff_flux_proj->calc_div_diff_flux));
408 PetscCall(DMRestoreGlobalVector(projection->dm, &RHS));
409 PetscCall(DMRestoreGlobalVector(projection->dm, &DiffFlux));
410 } break;
411 case DIV_DIFF_FLUX_PROJ_NONE:
412 SETERRQ(PetscObjectComm((PetscObject)projection->dm), PETSC_ERR_ARG_WRONG, "Should not reach here with div_diff_flux_projection_method %s",
413 DivDiffFluxProjectionMethods[diff_flux_proj->method]);
414 break;
415 }
416 PetscCall(PetscLogEventEnd(HONEE_DivDiffFluxProjection, Q_loc, 0, 0, 0));
417 PetscFunctionReturn(PETSC_SUCCESS);
418 }
419
420 /**
421 @brief Destroy `DivDiffFluxProjectionData` object
422
423 @param[in,out] diff_flux_proj Object to destroy
424 **/
DivDiffFluxProjectionDataDestroy(DivDiffFluxProjectionData diff_flux_proj)425 PetscErrorCode DivDiffFluxProjectionDataDestroy(DivDiffFluxProjectionData diff_flux_proj) {
426 PetscFunctionBeginUser;
427 if (diff_flux_proj == NULL) PetscFunctionReturn(PETSC_SUCCESS);
428 Ceed ceed = CeedVectorReturnCeed(diff_flux_proj->div_diff_flux_ceed);
429
430 PetscCall(NodalProjectionDataDestroy(&diff_flux_proj->projection));
431 PetscCall(OperatorApplyContextDestroy(diff_flux_proj->calc_div_diff_flux));
432 if (diff_flux_proj->ceed_vec_has_array) {
433 PetscCall(VecReadCeedToPetsc(diff_flux_proj->div_diff_flux_ceed, diff_flux_proj->DivDiffFlux_memtype, diff_flux_proj->DivDiffFlux_loc));
434 diff_flux_proj->ceed_vec_has_array = PETSC_FALSE;
435 }
436 PetscCallCeed(ceed, CeedVectorDestroy(&diff_flux_proj->div_diff_flux_ceed));
437 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&diff_flux_proj->elem_restr_div_diff_flux));
438 PetscCallCeed(ceed, CeedBasisDestroy(&diff_flux_proj->basis_div_diff_flux));
439 PetscCall(VecDestroy(&diff_flux_proj->DivDiffFlux_loc));
440 PetscCall(PetscFree(diff_flux_proj));
441 PetscFunctionReturn(PETSC_SUCCESS);
442 }
443