xref: /honee/src/diff_flux_projection.c (revision 0880fbb64d64d53bf2220d9ebe5cb405355858f6)
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] ceed      `Ceed` context
114   @param[in] user      `User` context
115   @param[in] ceed_data `CeedData` context
116   @param[in] problem   `ProblemData` context
117 **/
118 static PetscErrorCode DivDiffFluxProjectionSetup_Direct(Ceed ceed, User user, CeedData ceed_data, ProblemData problem) {
119   DivDiffFluxProjectionData diff_flux_proj = user->diff_flux_proj;
120   NodalProjectionData       projection     = diff_flux_proj->projection;
121   MPI_Comm                  comm           = PetscObjectComm((PetscObject)projection->dm);
122 
123   PetscFunctionBeginUser;
124   {  // Create Projection RHS OperatorApplyContext
125     CeedOperator op_rhs;
126 
127     PetscCheck(diff_flux_proj->CreateRHSOperator_Direct, comm, PETSC_ERR_ARG_WRONGSTATE,
128                "Must define CreateRHSOperator_Direct to use indirect div_diff_flux projection");
129     PetscCall(diff_flux_proj->CreateRHSOperator_Direct(user, ceed_data, diff_flux_proj, &op_rhs));
130     PetscCall(DMCreateLocalVector(projection->dm, &diff_flux_proj->DivDiffFlux_loc));
131     diff_flux_proj->ceed_vec_has_array = PETSC_FALSE;
132     PetscCall(OperatorApplyContextCreate(user->dm, projection->dm, ceed, op_rhs, NULL, NULL, NULL, diff_flux_proj->DivDiffFlux_loc,
133                                          &projection->l2_rhs_ctx));
134     PetscCallCeed(ceed, CeedOperatorDestroy(&op_rhs));
135   }
136 
137   {  // -- Build Mass operator
138     CeedQFunction       qf_mass;
139     CeedOperator        op_mass;
140     CeedBasis           basis_diff_flux             = NULL;
141     CeedElemRestriction elem_restr_diff_flux_volume = NULL, elem_restr_qd;
142     CeedVector          q_data;
143     CeedInt             q_data_size;
144     PetscInt            label_value  = 0;
145     DMLabel             domain_label = NULL;
146 
147     PetscCall(DivDiffFluxProjectionGetOperatorFieldData(diff_flux_proj, &elem_restr_diff_flux_volume, &basis_diff_flux, NULL, NULL));
148     PetscCall(QDataGet(ceed, projection->dm, domain_label, label_value, ceed_data->elem_restr_x, ceed_data->basis_x, ceed_data->x_coord,
149                        &elem_restr_qd, &q_data, &q_data_size));
150 
151     PetscCall(CreateMassQFunction(ceed, projection->num_comp, q_data_size, &qf_mass));
152     PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_mass, NULL, NULL, &op_mass));
153     PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "u", elem_restr_diff_flux_volume, basis_diff_flux, CEED_VECTOR_ACTIVE));
154     PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "qdata", elem_restr_qd, CEED_BASIS_NONE, q_data));
155     PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "v", elem_restr_diff_flux_volume, basis_diff_flux, CEED_VECTOR_ACTIVE));
156 
157     {  // -- Setup KSP for L^2 projection
158       Mat mat_mass;
159 
160       PetscCall(MatCreateCeed(projection->dm, projection->dm, op_mass, NULL, &mat_mass));
161 
162       PetscCall(KSPCreate(comm, &projection->ksp));
163       PetscCall(KSPSetOptionsPrefix(projection->ksp, "div_diff_flux_projection_"));
164       {  // lumped by default
165         PC pc;
166         PetscCall(KSPGetPC(projection->ksp, &pc));
167         PetscCall(PCSetType(pc, PCJACOBI));
168         PetscCall(PCJacobiSetType(pc, PC_JACOBI_ROWSUM));
169         PetscCall(KSPSetType(projection->ksp, KSPPREONLY));
170       }
171       PetscCall(KSPSetFromOptions_WithMatCeed(projection->ksp, mat_mass));
172       PetscCall(MatDestroy(&mat_mass));
173     }
174 
175     PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_diff_flux_volume));
176     PetscCallCeed(ceed, CeedBasisDestroy(&basis_diff_flux));
177     PetscCallCeed(ceed, CeedVectorDestroy(&q_data));
178     PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd));
179     PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_mass));
180     PetscCallCeed(ceed, CeedOperatorDestroy(&op_mass));
181   }
182   PetscFunctionReturn(PETSC_SUCCESS);
183 }
184 
185 /**
186   @brief Setup indirect projection of divergence of diffusive flux
187 
188   @param[in]     ceed      `Ceed` context
189   @param[in,out] user      `User` context
190   @param[in]     ceed_data `CeedData` context
191   @param[in]     problem   `ProblemData` context
192 **/
193 static PetscErrorCode DivDiffFluxProjectionSetup_Indirect(Ceed ceed, User user, CeedData ceed_data, ProblemData problem) {
194   DivDiffFluxProjectionData diff_flux_proj = user->diff_flux_proj;
195   NodalProjectionData       projection     = diff_flux_proj->projection;
196   CeedBasis                 basis_diff_flux;
197   CeedElemRestriction       elem_restr_diff_flux, elem_restr_qd;
198   CeedVector                q_data;
199   CeedInt                   q_data_size;
200   MPI_Comm                  comm = PetscObjectComm((PetscObject)projection->dm);
201 
202   PetscFunctionBeginUser;
203   {
204     PetscInt label_value = 0, height = 0, dm_field = 0;
205     DMLabel  domain_label = NULL;
206 
207     PetscCall(DMPlexCeedElemRestrictionCreate(ceed, projection->dm, domain_label, label_value, height, dm_field, &elem_restr_diff_flux));
208     PetscCall(CreateBasisFromPlex(ceed, projection->dm, domain_label, label_value, height, dm_field, &basis_diff_flux));
209     PetscCall(QDataGet(ceed, projection->dm, domain_label, label_value, ceed_data->elem_restr_x, ceed_data->basis_x, ceed_data->x_coord,
210                        &elem_restr_qd, &q_data, &q_data_size));
211   }
212 
213   {
214     CeedOperator op_rhs;
215 
216     PetscCheck(diff_flux_proj->CreateRHSOperator_Indirect, comm, PETSC_ERR_ARG_WRONGSTATE,
217                "Must define CreateRHSOperator_Indirect to use indirect div_diff_flux projection");
218     PetscCall(diff_flux_proj->CreateRHSOperator_Indirect(user, ceed_data, diff_flux_proj, &op_rhs));
219     PetscCall(OperatorApplyContextCreate(user->dm, projection->dm, ceed, op_rhs, NULL, NULL, NULL, NULL, &projection->l2_rhs_ctx));
220     PetscCallCeed(ceed, CeedOperatorDestroy(&op_rhs));
221   }
222 
223   {  // -- Build Mass operator
224     CeedQFunction qf_mass;
225     CeedOperator  op_mass;
226 
227     PetscCall(CreateMassQFunction(ceed, projection->num_comp, q_data_size, &qf_mass));
228     PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_mass, NULL, NULL, &op_mass));
229     PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "u", elem_restr_diff_flux, basis_diff_flux, CEED_VECTOR_ACTIVE));
230     PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "qdata", elem_restr_qd, CEED_BASIS_NONE, q_data));
231     PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "v", elem_restr_diff_flux, basis_diff_flux, CEED_VECTOR_ACTIVE));
232 
233     {  // -- Setup KSP for L^2 projection
234       Mat mat_mass;
235 
236       PetscCall(MatCreateCeed(projection->dm, projection->dm, op_mass, NULL, &mat_mass));
237 
238       PetscCall(KSPCreate(comm, &projection->ksp));
239       PetscCall(KSPSetOptionsPrefix(projection->ksp, "div_diff_flux_projection_"));
240       {  // lumped by default
241         PC pc;
242         PetscCall(KSPGetPC(projection->ksp, &pc));
243         PetscCall(PCSetType(pc, PCJACOBI));
244         PetscCall(PCJacobiSetType(pc, PC_JACOBI_ROWSUM));
245         PetscCall(KSPSetType(projection->ksp, KSPPREONLY));
246       }
247       PetscCall(KSPSetFromOptions_WithMatCeed(projection->ksp, mat_mass));
248       PetscCall(MatDestroy(&mat_mass));
249     }
250     PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_mass));
251     PetscCallCeed(ceed, CeedOperatorDestroy(&op_mass));
252   }
253 
254   {  // Create OperatorApplyContext to calculate divergence at quadrature points
255     CeedQFunction       qf_calc_divergence;
256     CeedOperator        op_calc_divergence;
257     CeedElemRestriction elem_restr_div_diff_flux = NULL;
258     PetscInt            dim;
259 
260     PetscCall(DMGetDimension(projection->dm, &dim));
261     PetscCall(DivDiffFluxProjectionGetOperatorFieldData(diff_flux_proj, &elem_restr_div_diff_flux, NULL, NULL, NULL));
262 
263     PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, ComputeDivDiffusiveFlux3D_4, ComputeDivDiffusiveFlux3D_4_loc, &qf_calc_divergence));
264 
265     PetscCallCeed(ceed, CeedQFunctionAddInput(qf_calc_divergence, "Grad F_diff", projection->num_comp * dim, CEED_EVAL_GRAD));
266     PetscCallCeed(ceed, CeedQFunctionAddInput(qf_calc_divergence, "qdata", q_data_size, CEED_EVAL_NONE));
267     PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_calc_divergence, "Div F_diff", 4, CEED_EVAL_NONE));
268 
269     PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_calc_divergence, NULL, NULL, &op_calc_divergence));
270     PetscCallCeed(ceed, CeedOperatorSetField(op_calc_divergence, "Grad F_diff", elem_restr_diff_flux, basis_diff_flux, CEED_VECTOR_ACTIVE));
271     PetscCallCeed(ceed, CeedOperatorSetField(op_calc_divergence, "qdata", elem_restr_qd, CEED_BASIS_NONE, q_data));
272     PetscCallCeed(
273         ceed, CeedOperatorSetField(op_calc_divergence, "Div F_diff", elem_restr_div_diff_flux, CEED_BASIS_NONE, diff_flux_proj->div_diff_flux_ceed));
274 
275     PetscCall(
276         OperatorApplyContextCreate(projection->dm, NULL, ceed, op_calc_divergence, NULL, NULL, NULL, NULL, &diff_flux_proj->calc_div_diff_flux));
277 
278     PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_div_diff_flux));
279     PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_calc_divergence));
280     PetscCallCeed(ceed, CeedOperatorDestroy(&op_calc_divergence));
281   }
282   PetscCallCeed(ceed, CeedBasisDestroy(&basis_diff_flux));
283   PetscCallCeed(ceed, CeedVectorDestroy(&q_data));
284   PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd));
285   PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_diff_flux));
286   PetscFunctionReturn(PETSC_SUCCESS);
287 }
288 
289 /**
290   @brief Setup projection of divergence of diffusive flux
291 
292   @param[in]     ceed      `Ceed` context
293   @param[in,out] user      `User` context
294   @param[in]     ceed_data `CeedData` context
295   @param[in]     problem   `ProblemData` context
296 **/
297 PetscErrorCode DivDiffFluxProjectionSetup(Ceed ceed, User user, CeedData ceed_data, ProblemData problem) {
298   PetscFunctionBeginUser;
299   switch (user->app_ctx->divFdiffproj_method) {
300     case DIV_DIFF_FLUX_PROJ_DIRECT:
301       PetscCall(DivDiffFluxProjectionSetup_Direct(ceed, user, ceed_data, problem));
302       break;
303     case DIV_DIFF_FLUX_PROJ_INDIRECT:
304       PetscCall(DivDiffFluxProjectionSetup_Indirect(ceed, user, ceed_data, problem));
305       break;
306     case DIV_DIFF_FLUX_PROJ_NONE:
307       SETERRQ(PetscObjectComm((PetscObject)user->dm), PETSC_ERR_ARG_WRONG, "Should not reach here with div_diff_flux_projection_method %s",
308               DivDiffFluxProjectionMethods[user->app_ctx->divFdiffproj_method]);
309       break;
310   }
311   PetscFunctionReturn(PETSC_SUCCESS);
312 }
313 
314 /**
315   @brief Project the divergence of diffusive flux
316 
317   This implicitly sets the `CeedVector` input (`div_diff_flux_ceed`) to the divergence of diffusive flux.
318 
319   @param[in]  diff_flux_proj `NodalProjectionData` for the projection
320   @param[in]  Q_loc          Localized solution vector
321 **/
322 PetscErrorCode DivDiffFluxProjectionApply(DivDiffFluxProjectionData diff_flux_proj, Vec Q_loc) {
323   NodalProjectionData projection = diff_flux_proj->projection;
324 
325   PetscFunctionBeginUser;
326   PetscCall(PetscLogEventBegin(FLUIDS_DivDiffFluxProjection, Q_loc, 0, 0, 0));
327   switch (diff_flux_proj->method) {
328     case DIV_DIFF_FLUX_PROJ_DIRECT: {
329       Vec DivDiffFlux;
330 
331       PetscCall(DMGetGlobalVector(projection->dm, &DivDiffFlux));
332       if (diff_flux_proj->ceed_vec_has_array) {
333         PetscCall(VecReadCeedToPetsc(diff_flux_proj->div_diff_flux_ceed, diff_flux_proj->DivDiffFlux_memtype, diff_flux_proj->DivDiffFlux_loc));
334         diff_flux_proj->ceed_vec_has_array = PETSC_FALSE;
335       }
336       PetscCall(ApplyCeedOperatorLocalToGlobal(Q_loc, DivDiffFlux, projection->l2_rhs_ctx));
337       PetscCall(VecViewFromOptions(DivDiffFlux, NULL, "-div_diff_flux_projection_rhs_view"));
338 
339       PetscCall(KSPSolve(projection->ksp, DivDiffFlux, DivDiffFlux));
340       PetscCall(VecViewFromOptions(DivDiffFlux, NULL, "-div_diff_flux_projection_view"));
341 
342       PetscCall(DMGlobalToLocal(projection->dm, DivDiffFlux, INSERT_VALUES, diff_flux_proj->DivDiffFlux_loc));
343       PetscCall(VecReadPetscToCeed(diff_flux_proj->DivDiffFlux_loc, &diff_flux_proj->DivDiffFlux_memtype, diff_flux_proj->div_diff_flux_ceed));
344       diff_flux_proj->ceed_vec_has_array = PETSC_TRUE;
345 
346       PetscCall(DMRestoreGlobalVector(projection->dm, &DivDiffFlux));
347       break;
348     }
349     case DIV_DIFF_FLUX_PROJ_INDIRECT: {
350       Vec DiffFlux;
351 
352       PetscCall(DMGetGlobalVector(projection->dm, &DiffFlux));
353       PetscCall(ApplyCeedOperatorLocalToGlobal(Q_loc, DiffFlux, projection->l2_rhs_ctx));
354       PetscCall(VecViewFromOptions(DiffFlux, NULL, "-div_diff_flux_projection_rhs_view"));
355 
356       PetscCall(KSPSolve(projection->ksp, DiffFlux, DiffFlux));
357       PetscCall(VecViewFromOptions(DiffFlux, NULL, "-div_diff_flux_projection_view"));
358 
359       PetscCall(ApplyCeedOperatorGlobalToLocal(DiffFlux, NULL, diff_flux_proj->calc_div_diff_flux));
360       PetscCall(DMRestoreGlobalVector(projection->dm, &DiffFlux));
361     } break;
362     case DIV_DIFF_FLUX_PROJ_NONE:
363       SETERRQ(PetscObjectComm((PetscObject)projection->dm), PETSC_ERR_ARG_WRONG, "Should not reach here with div_diff_flux_projection_method %s",
364               DivDiffFluxProjectionMethods[diff_flux_proj->method]);
365       break;
366   }
367   PetscCall(PetscLogEventEnd(FLUIDS_DivDiffFluxProjection, Q_loc, 0, 0, 0));
368   PetscFunctionReturn(PETSC_SUCCESS);
369 }
370 
371 /**
372   @brief Destroy `DivDiffFluxProjectionData` object
373 
374   @param[in,out] diff_flux_proj Object to destroy
375 **/
376 PetscErrorCode DivDiffFluxProjectionDataDestroy(DivDiffFluxProjectionData diff_flux_proj) {
377   PetscFunctionBeginUser;
378   if (diff_flux_proj == NULL) PetscFunctionReturn(PETSC_SUCCESS);
379   Ceed ceed = CeedVectorReturnCeed(diff_flux_proj->div_diff_flux_ceed);
380 
381   PetscCall(NodalProjectionDataDestroy(diff_flux_proj->projection));
382   PetscCall(OperatorApplyContextDestroy(diff_flux_proj->calc_div_diff_flux));
383   if (diff_flux_proj->ceed_vec_has_array) {
384     PetscCall(VecReadCeedToPetsc(diff_flux_proj->div_diff_flux_ceed, diff_flux_proj->DivDiffFlux_memtype, diff_flux_proj->DivDiffFlux_loc));
385     diff_flux_proj->ceed_vec_has_array = PETSC_FALSE;
386   }
387   PetscCallCeed(ceed, CeedVectorDestroy(&diff_flux_proj->div_diff_flux_ceed));
388   PetscCallCeed(ceed, CeedElemRestrictionDestroy(&diff_flux_proj->elem_restr_div_diff_flux));
389   PetscCallCeed(ceed, CeedBasisDestroy(&diff_flux_proj->basis_div_diff_flux));
390   PetscCall(VecDestroy(&diff_flux_proj->DivDiffFlux_loc));
391   PetscCall(PetscFree(diff_flux_proj));
392   PetscFunctionReturn(PETSC_SUCCESS);
393 }
394