xref: /honee/src/diff_flux_projection.c (revision 8561fee26e1c8104caa6d917d7ecea939d705f29)
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 /**
13   @brief Create `DivDiffFluxProjectionData` for solution DM in `user`
14 
15   @param[in]  user                `User` context
16   @param[in]  num_diff_flux_comps Number of components that makes up the diffusive flux (e.g. 1 for scalar advection-diffusion)
17   @param[out] pdiff_flux_proj     The `DivDiffFluxProjectionData` object created, or `NULL` if not created
18 **/
19 PetscErrorCode DivDiffFluxProjectionCreate(User user, PetscInt num_diff_flux_comps, DivDiffFluxProjectionData *pdiff_flux_proj) {
20   PetscInt                  label_value = 0, height = 0, dm_field = 0, dim, degree = user->app_ctx->degree, q_extra = user->app_ctx->q_extra;
21   DMLabel                   domain_label = NULL;
22   DivDiffFluxProjectionData diff_flux_proj;
23   NodalProjectionData       projection;
24 
25   PetscFunctionBeginUser;
26   if (user->app_ctx->divFdiffproj_method == DIV_DIFF_FLUX_PROJ_NONE) {
27     *pdiff_flux_proj = NULL;
28     PetscFunctionReturn(PETSC_SUCCESS);
29   }
30   PetscCall(PetscNew(pdiff_flux_proj));
31   diff_flux_proj = *pdiff_flux_proj;
32   PetscCall(PetscNew(&user->diff_flux_proj->projection));
33   projection                          = user->diff_flux_proj->projection;
34   diff_flux_proj->method              = user->app_ctx->divFdiffproj_method;
35   diff_flux_proj->num_diff_flux_comps = num_diff_flux_comps;
36 
37   PetscCall(DMClone(user->dm, &projection->dm));
38   PetscCall(DMSetMatrixPreallocateSkip(projection->dm, PETSC_TRUE));
39   PetscCall(DMGetDimension(projection->dm, &dim));
40   switch (diff_flux_proj->method) {
41     case DIV_DIFF_FLUX_PROJ_DIRECT: {
42       projection->num_comp = diff_flux_proj->num_diff_flux_comps;
43       PetscCall(PetscObjectSetName((PetscObject)projection->dm, "DivDiffFluxProj"));
44       PetscCall(DMSetupByOrder_FEM(PETSC_TRUE, PETSC_TRUE, degree, 1, q_extra, 1, &projection->num_comp, projection->dm));
45 
46       PetscCall(DMPlexCeedElemRestrictionCreate(user->ceed, projection->dm, domain_label, label_value, height, dm_field,
47                                                 &diff_flux_proj->elem_restr_div_diff_flux));
48       PetscCallCeed(user->ceed, CeedElemRestrictionCreateVector(diff_flux_proj->elem_restr_div_diff_flux, &diff_flux_proj->div_diff_flux_ceed, NULL));
49       PetscCall(CreateBasisFromPlex(user->ceed, projection->dm, domain_label, label_value, height, dm_field, &diff_flux_proj->basis_div_diff_flux));
50       diff_flux_proj->eval_mode_div_diff_flux = CEED_EVAL_INTERP;
51 
52       {  // Create face labels on projection->dm for boundary integrals
53         DMLabel  face_sets_label;
54         PetscInt num_face_set_values, *face_set_values;
55 
56         PetscCall(DMGetLabel(user->dm, "Face Sets", &face_sets_label));
57         PetscCall(DMLabelCreateGlobalValueArray(user->dm, face_sets_label, &num_face_set_values, &face_set_values));
58         for (PetscInt f = 0; f < num_face_set_values; f++) {
59           DMLabel face_orientation_label;
60           char   *face_orientation_label_name;
61 
62           PetscCall(DMPlexCreateFaceLabel(user->dm, face_set_values[f], &face_orientation_label_name));
63           PetscCall(DMGetLabel(user->dm, face_orientation_label_name, &face_orientation_label));
64           PetscCall(DMAddLabel(projection->dm, face_orientation_label));
65           PetscCall(PetscFree(face_orientation_label_name));
66         }
67         PetscCall(PetscFree(face_set_values));
68       }
69     } break;
70     case DIV_DIFF_FLUX_PROJ_INDIRECT: {
71       projection->num_comp = diff_flux_proj->num_diff_flux_comps * dim;
72       PetscCall(PetscObjectSetName((PetscObject)projection->dm, "DiffFluxProj"));
73       PetscCall(DMSetupByOrder_FEM(PETSC_TRUE, PETSC_TRUE, degree, 1, q_extra, 1, &projection->num_comp, projection->dm));
74 
75       PetscCall(DMPlexCeedElemRestrictionQDataCreate(user->ceed, projection->dm, domain_label, label_value, height,
76                                                      diff_flux_proj->num_diff_flux_comps, &diff_flux_proj->elem_restr_div_diff_flux));
77       PetscCallCeed(user->ceed, CeedElemRestrictionCreateVector(diff_flux_proj->elem_restr_div_diff_flux, &diff_flux_proj->div_diff_flux_ceed, NULL));
78       diff_flux_proj->basis_div_diff_flux     = CEED_BASIS_NONE;
79       diff_flux_proj->eval_mode_div_diff_flux = CEED_EVAL_NONE;
80     } break;
81     case DIV_DIFF_FLUX_PROJ_NONE:
82       SETERRQ(PetscObjectComm((PetscObject)user->dm), PETSC_ERR_ARG_WRONG, "Should not reach here with div_diff_flux_projection_method %s",
83               DivDiffFluxProjectionMethods[user->app_ctx->divFdiffproj_method]);
84       break;
85   }
86   PetscFunctionReturn(PETSC_SUCCESS);
87 };
88 
89 /**
90   @brief Return the objects required for the Divergence of Diffusive flux to be read by a `CeedOperator`
91 
92   @param[in]  diff_flux_proj Projection object
93   @param[out] elem_restr     Element restriction for the divergence of diffusive flux, or `NULL`
94   @param[out] basis          Basis for the divergence of diffusive flux, or `NULL`
95   @param[out] vector         Vector for the divergence of diffusive flux, or `NULL`
96   @param[out] eval_mode      Eval mode for the divergence of diffusive flux, or `NULL`
97 **/
98 PetscErrorCode DivDiffFluxProjectionGetOperatorFieldData(DivDiffFluxProjectionData diff_flux_proj, CeedElemRestriction *elem_restr, CeedBasis *basis,
99                                                          CeedVector *vector, CeedEvalMode *eval_mode) {
100   Ceed ceed = CeedVectorReturnCeed(diff_flux_proj->div_diff_flux_ceed);
101 
102   PetscFunctionBeginUser;
103   if (elem_restr) PetscCallCeed(ceed, CeedElemRestrictionReferenceCopy(diff_flux_proj->elem_restr_div_diff_flux, elem_restr));
104   if (basis) PetscCallCeed(ceed, CeedBasisReferenceCopy(diff_flux_proj->basis_div_diff_flux, basis));
105   if (vector) PetscCallCeed(ceed, CeedVectorReferenceCopy(diff_flux_proj->div_diff_flux_ceed, vector));
106   if (eval_mode) *eval_mode = diff_flux_proj->eval_mode_div_diff_flux;
107   PetscFunctionReturn(PETSC_SUCCESS);
108 }
109 
110 /**
111   @brief Setup direct projection of divergence of diffusive flux
112 
113   @param[in]     user           `User` context
114   @param[in]     ceed_data      `CeedData` context
115   @param[in,out] diff_flux_proj Flux projection object to setup
116 **/
117 static PetscErrorCode DivDiffFluxProjectionSetup_Direct(User user, CeedData ceed_data, DivDiffFluxProjectionData diff_flux_proj) {
118   Ceed                ceed       = user->ceed;
119   NodalProjectionData projection = diff_flux_proj->projection;
120   MPI_Comm            comm       = PetscObjectComm((PetscObject)projection->dm);
121 
122   PetscFunctionBeginUser;
123   {  // Create Projection RHS OperatorApplyContext
124     CeedOperator op_rhs;
125 
126     PetscCheck(diff_flux_proj->CreateRHSOperator_Direct, comm, PETSC_ERR_ARG_WRONGSTATE,
127                "Must define CreateRHSOperator_Direct to use indirect div_diff_flux projection");
128     PetscCall(diff_flux_proj->CreateRHSOperator_Direct(user, ceed_data, diff_flux_proj, &op_rhs));
129     PetscCall(DMCreateLocalVector(projection->dm, &diff_flux_proj->DivDiffFlux_loc));
130     diff_flux_proj->ceed_vec_has_array = PETSC_FALSE;
131     PetscCall(OperatorApplyContextCreate(user->dm, projection->dm, ceed, op_rhs, NULL, NULL, NULL, diff_flux_proj->DivDiffFlux_loc,
132                                          &projection->l2_rhs_ctx));
133     PetscCallCeed(ceed, CeedOperatorDestroy(&op_rhs));
134   }
135 
136   {  // -- Build Mass operator
137     CeedQFunction       qf_mass;
138     CeedOperator        op_mass;
139     CeedBasis           basis_div_diff_flux             = NULL;
140     CeedElemRestriction elem_restr_div_diff_flux_volume = NULL, elem_restr_qd;
141     CeedVector          q_data;
142     CeedInt             q_data_size;
143     PetscInt            label_value  = 0;
144     DMLabel             domain_label = NULL;
145 
146     PetscCall(DivDiffFluxProjectionGetOperatorFieldData(diff_flux_proj, &elem_restr_div_diff_flux_volume, &basis_div_diff_flux, NULL, NULL));
147     PetscCall(QDataGet(ceed, projection->dm, domain_label, label_value, ceed_data->elem_restr_x, ceed_data->basis_x, ceed_data->x_coord,
148                        &elem_restr_qd, &q_data, &q_data_size));
149 
150     PetscCall(CreateMassQFunction(ceed, projection->num_comp, q_data_size, &qf_mass));
151     PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_mass, NULL, NULL, &op_mass));
152     PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "u", elem_restr_div_diff_flux_volume, basis_div_diff_flux, CEED_VECTOR_ACTIVE));
153     PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "qdata", elem_restr_qd, CEED_BASIS_NONE, q_data));
154     PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "v", elem_restr_div_diff_flux_volume, basis_div_diff_flux, CEED_VECTOR_ACTIVE));
155 
156     {  // -- Setup KSP for L^2 projection
157       Mat mat_mass;
158 
159       PetscCall(MatCreateCeed(projection->dm, projection->dm, op_mass, NULL, &mat_mass));
160 
161       PetscCall(KSPCreate(comm, &projection->ksp));
162       PetscCall(KSPSetOptionsPrefix(projection->ksp, "div_diff_flux_projection_"));
163       {  // lumped by default
164         PC pc;
165         PetscCall(KSPGetPC(projection->ksp, &pc));
166         PetscCall(PCSetType(pc, PCJACOBI));
167         PetscCall(PCJacobiSetType(pc, PC_JACOBI_ROWSUM));
168         PetscCall(KSPSetType(projection->ksp, KSPPREONLY));
169       }
170       PetscCall(KSPSetFromOptions_WithMatCeed(projection->ksp, mat_mass));
171       PetscCall(MatDestroy(&mat_mass));
172     }
173 
174     PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_div_diff_flux_volume));
175     PetscCallCeed(ceed, CeedBasisDestroy(&basis_div_diff_flux));
176     PetscCallCeed(ceed, CeedVectorDestroy(&q_data));
177     PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd));
178     PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_mass));
179     PetscCallCeed(ceed, CeedOperatorDestroy(&op_mass));
180   }
181   PetscFunctionReturn(PETSC_SUCCESS);
182 }
183 
184 /**
185   @brief Setup indirect projection of divergence of diffusive flux
186 
187   @param[in]     user           `User` context
188   @param[in]     ceed_data      `CeedData` context
189   @param[in,out] diff_flux_proj Flux projection object to setup
190 **/
191 static PetscErrorCode DivDiffFluxProjectionSetup_Indirect(User user, CeedData ceed_data, DivDiffFluxProjectionData diff_flux_proj) {
192   Ceed                ceed       = user->ceed;
193   NodalProjectionData projection = diff_flux_proj->projection;
194   CeedBasis           basis_diff_flux;
195   CeedElemRestriction elem_restr_diff_flux, elem_restr_qd;
196   CeedVector          q_data;
197   CeedInt             q_data_size;
198   MPI_Comm            comm = PetscObjectComm((PetscObject)projection->dm);
199 
200   PetscFunctionBeginUser;
201   {
202     PetscInt label_value = 0, height = 0, dm_field = 0;
203     DMLabel  domain_label = NULL;
204 
205     PetscCall(DMPlexCeedElemRestrictionCreate(ceed, projection->dm, domain_label, label_value, height, dm_field, &elem_restr_diff_flux));
206     PetscCall(CreateBasisFromPlex(ceed, projection->dm, domain_label, label_value, height, dm_field, &basis_diff_flux));
207     PetscCall(QDataGet(ceed, projection->dm, domain_label, label_value, ceed_data->elem_restr_x, ceed_data->basis_x, ceed_data->x_coord,
208                        &elem_restr_qd, &q_data, &q_data_size));
209   }
210 
211   {
212     CeedOperator op_rhs;
213 
214     PetscCheck(diff_flux_proj->CreateRHSOperator_Indirect, comm, PETSC_ERR_ARG_WRONGSTATE,
215                "Must define CreateRHSOperator_Indirect to use indirect div_diff_flux projection");
216     PetscCall(diff_flux_proj->CreateRHSOperator_Indirect(user, ceed_data, diff_flux_proj, &op_rhs));
217     PetscCall(OperatorApplyContextCreate(user->dm, projection->dm, ceed, op_rhs, NULL, NULL, NULL, NULL, &projection->l2_rhs_ctx));
218     PetscCallCeed(ceed, CeedOperatorDestroy(&op_rhs));
219   }
220 
221   {  // -- Build Mass operator
222     CeedQFunction qf_mass;
223     CeedOperator  op_mass;
224 
225     PetscCall(CreateMassQFunction(ceed, projection->num_comp, q_data_size, &qf_mass));
226     PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_mass, NULL, NULL, &op_mass));
227     PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "u", elem_restr_diff_flux, basis_diff_flux, CEED_VECTOR_ACTIVE));
228     PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "qdata", elem_restr_qd, CEED_BASIS_NONE, q_data));
229     PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "v", elem_restr_diff_flux, basis_diff_flux, CEED_VECTOR_ACTIVE));
230 
231     {  // -- Setup KSP for L^2 projection
232       Mat mat_mass;
233 
234       PetscCall(MatCreateCeed(projection->dm, projection->dm, op_mass, NULL, &mat_mass));
235 
236       PetscCall(KSPCreate(comm, &projection->ksp));
237       PetscCall(KSPSetOptionsPrefix(projection->ksp, "div_diff_flux_projection_"));
238       {  // lumped by default
239         PC pc;
240         PetscCall(KSPGetPC(projection->ksp, &pc));
241         PetscCall(PCSetType(pc, PCJACOBI));
242         PetscCall(PCJacobiSetType(pc, PC_JACOBI_ROWSUM));
243         PetscCall(KSPSetType(projection->ksp, KSPPREONLY));
244       }
245       PetscCall(KSPSetFromOptions_WithMatCeed(projection->ksp, mat_mass));
246       PetscCall(MatDestroy(&mat_mass));
247     }
248     PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_mass));
249     PetscCallCeed(ceed, CeedOperatorDestroy(&op_mass));
250   }
251 
252   {  // Create OperatorApplyContext to calculate divergence at quadrature points
253     CeedQFunction       qf_calc_divergence;
254     CeedOperator        op_calc_divergence;
255     CeedElemRestriction elem_restr_div_diff_flux = NULL;
256     PetscInt            dim;
257 
258     PetscCall(DMGetDimension(projection->dm, &dim));
259     PetscCall(DivDiffFluxProjectionGetOperatorFieldData(diff_flux_proj, &elem_restr_div_diff_flux, NULL, NULL, NULL));
260 
261     PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, ComputeDivDiffusiveFlux3D_4, ComputeDivDiffusiveFlux3D_4_loc, &qf_calc_divergence));
262 
263     PetscCallCeed(ceed, CeedQFunctionAddInput(qf_calc_divergence, "Grad F_diff", projection->num_comp * dim, CEED_EVAL_GRAD));
264     PetscCallCeed(ceed, CeedQFunctionAddInput(qf_calc_divergence, "qdata", q_data_size, CEED_EVAL_NONE));
265     PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_calc_divergence, "Div F_diff", 4, CEED_EVAL_NONE));
266 
267     PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_calc_divergence, NULL, NULL, &op_calc_divergence));
268     PetscCallCeed(ceed, CeedOperatorSetField(op_calc_divergence, "Grad F_diff", elem_restr_diff_flux, basis_diff_flux, CEED_VECTOR_ACTIVE));
269     PetscCallCeed(ceed, CeedOperatorSetField(op_calc_divergence, "qdata", elem_restr_qd, CEED_BASIS_NONE, q_data));
270     PetscCallCeed(
271         ceed, CeedOperatorSetField(op_calc_divergence, "Div F_diff", elem_restr_div_diff_flux, CEED_BASIS_NONE, diff_flux_proj->div_diff_flux_ceed));
272 
273     PetscCall(
274         OperatorApplyContextCreate(projection->dm, NULL, ceed, op_calc_divergence, NULL, NULL, NULL, NULL, &diff_flux_proj->calc_div_diff_flux));
275 
276     PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_div_diff_flux));
277     PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_calc_divergence));
278     PetscCallCeed(ceed, CeedOperatorDestroy(&op_calc_divergence));
279   }
280   PetscCallCeed(ceed, CeedBasisDestroy(&basis_diff_flux));
281   PetscCallCeed(ceed, CeedVectorDestroy(&q_data));
282   PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd));
283   PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_diff_flux));
284   PetscFunctionReturn(PETSC_SUCCESS);
285 }
286 
287 /**
288   @brief Setup projection of divergence of diffusive flux
289 
290   @param[in]     user           `User` context
291   @param[in]     ceed_data      `CeedData` context
292   @param[in,out] diff_flux_proj Flux projection object to setup
293 **/
294 PetscErrorCode DivDiffFluxProjectionSetup(User user, CeedData ceed_data, DivDiffFluxProjectionData diff_flux_proj) {
295   PetscFunctionBeginUser;
296   switch (user->app_ctx->divFdiffproj_method) {
297     case DIV_DIFF_FLUX_PROJ_DIRECT:
298       PetscCall(DivDiffFluxProjectionSetup_Direct(user, ceed_data, diff_flux_proj));
299       break;
300     case DIV_DIFF_FLUX_PROJ_INDIRECT:
301       PetscCall(DivDiffFluxProjectionSetup_Indirect(user, ceed_data, diff_flux_proj));
302       break;
303     case DIV_DIFF_FLUX_PROJ_NONE:
304       SETERRQ(PetscObjectComm((PetscObject)user->dm), PETSC_ERR_ARG_WRONG, "Should not reach here with div_diff_flux_projection_method %s",
305               DivDiffFluxProjectionMethods[user->app_ctx->divFdiffproj_method]);
306       break;
307   }
308   PetscFunctionReturn(PETSC_SUCCESS);
309 }
310 
311 /**
312   @brief Project the divergence of diffusive flux
313 
314   This implicitly sets the `CeedVector` input (`div_diff_flux_ceed`) to the divergence of diffusive flux.
315 
316   @param[in]  diff_flux_proj `NodalProjectionData` for the projection
317   @param[in]  Q_loc          Localized solution vector
318 **/
319 PetscErrorCode DivDiffFluxProjectionApply(DivDiffFluxProjectionData diff_flux_proj, Vec Q_loc) {
320   NodalProjectionData projection = diff_flux_proj->projection;
321 
322   PetscFunctionBeginUser;
323   PetscCall(PetscLogEventBegin(FLUIDS_DivDiffFluxProjection, Q_loc, 0, 0, 0));
324   switch (diff_flux_proj->method) {
325     case DIV_DIFF_FLUX_PROJ_DIRECT: {
326       Vec DivDiffFlux;
327 
328       PetscCall(DMGetGlobalVector(projection->dm, &DivDiffFlux));
329       if (diff_flux_proj->ceed_vec_has_array) {
330         PetscCall(VecReadCeedToPetsc(diff_flux_proj->div_diff_flux_ceed, diff_flux_proj->DivDiffFlux_memtype, diff_flux_proj->DivDiffFlux_loc));
331         diff_flux_proj->ceed_vec_has_array = PETSC_FALSE;
332       }
333       PetscCall(ApplyCeedOperatorLocalToGlobal(Q_loc, DivDiffFlux, projection->l2_rhs_ctx));
334       PetscCall(VecViewFromOptions(DivDiffFlux, NULL, "-div_diff_flux_projection_rhs_view"));
335 
336       PetscCall(KSPSolve(projection->ksp, DivDiffFlux, DivDiffFlux));
337       PetscCall(VecViewFromOptions(DivDiffFlux, NULL, "-div_diff_flux_projection_view"));
338 
339       PetscCall(DMGlobalToLocal(projection->dm, DivDiffFlux, INSERT_VALUES, diff_flux_proj->DivDiffFlux_loc));
340       PetscCall(VecReadPetscToCeed(diff_flux_proj->DivDiffFlux_loc, &diff_flux_proj->DivDiffFlux_memtype, diff_flux_proj->div_diff_flux_ceed));
341       diff_flux_proj->ceed_vec_has_array = PETSC_TRUE;
342 
343       PetscCall(DMRestoreGlobalVector(projection->dm, &DivDiffFlux));
344       break;
345     }
346     case DIV_DIFF_FLUX_PROJ_INDIRECT: {
347       Vec DiffFlux;
348 
349       PetscCall(DMGetGlobalVector(projection->dm, &DiffFlux));
350       PetscCall(ApplyCeedOperatorLocalToGlobal(Q_loc, DiffFlux, projection->l2_rhs_ctx));
351       PetscCall(VecViewFromOptions(DiffFlux, NULL, "-div_diff_flux_projection_rhs_view"));
352 
353       PetscCall(KSPSolve(projection->ksp, DiffFlux, DiffFlux));
354       PetscCall(VecViewFromOptions(DiffFlux, NULL, "-div_diff_flux_projection_view"));
355 
356       PetscCall(ApplyCeedOperatorGlobalToLocal(DiffFlux, NULL, diff_flux_proj->calc_div_diff_flux));
357       PetscCall(DMRestoreGlobalVector(projection->dm, &DiffFlux));
358     } break;
359     case DIV_DIFF_FLUX_PROJ_NONE:
360       SETERRQ(PetscObjectComm((PetscObject)projection->dm), PETSC_ERR_ARG_WRONG, "Should not reach here with div_diff_flux_projection_method %s",
361               DivDiffFluxProjectionMethods[diff_flux_proj->method]);
362       break;
363   }
364   PetscCall(PetscLogEventEnd(FLUIDS_DivDiffFluxProjection, Q_loc, 0, 0, 0));
365   PetscFunctionReturn(PETSC_SUCCESS);
366 }
367 
368 /**
369   @brief Destroy `DivDiffFluxProjectionData` object
370 
371   @param[in,out] diff_flux_proj Object to destroy
372 **/
373 PetscErrorCode DivDiffFluxProjectionDataDestroy(DivDiffFluxProjectionData diff_flux_proj) {
374   PetscFunctionBeginUser;
375   if (diff_flux_proj == NULL) PetscFunctionReturn(PETSC_SUCCESS);
376   Ceed ceed = CeedVectorReturnCeed(diff_flux_proj->div_diff_flux_ceed);
377 
378   PetscCall(NodalProjectionDataDestroy(diff_flux_proj->projection));
379   PetscCall(OperatorApplyContextDestroy(diff_flux_proj->calc_div_diff_flux));
380   if (diff_flux_proj->ceed_vec_has_array) {
381     PetscCall(VecReadCeedToPetsc(diff_flux_proj->div_diff_flux_ceed, diff_flux_proj->DivDiffFlux_memtype, diff_flux_proj->DivDiffFlux_loc));
382     diff_flux_proj->ceed_vec_has_array = PETSC_FALSE;
383   }
384   PetscCallCeed(ceed, CeedVectorDestroy(&diff_flux_proj->div_diff_flux_ceed));
385   PetscCallCeed(ceed, CeedElemRestrictionDestroy(&diff_flux_proj->elem_restr_div_diff_flux));
386   PetscCallCeed(ceed, CeedBasisDestroy(&diff_flux_proj->basis_div_diff_flux));
387   PetscCall(VecDestroy(&diff_flux_proj->DivDiffFlux_loc));
388   PetscCall(PetscFree(diff_flux_proj));
389   PetscFunctionReturn(PETSC_SUCCESS);
390 }
391