xref: /honee/src/diff_flux_projection.c (revision 36038bbcf8b5827c2ae723a5a7c52ea5c4047506)
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 Setup direct projection of divergence of diffusive flux
91 
92   @param[in] ceed      `Ceed` context
93   @param[in] user      `User` context
94   @param[in] ceed_data `CeedData` context
95   @param[in] problem   `ProblemData` context
96 **/
97 static PetscErrorCode DivDiffFluxProjectionSetup_Direct(Ceed ceed, User user, CeedData ceed_data, ProblemData problem) {
98   DivDiffFluxProjectionData diff_flux_proj = user->diff_flux_proj;
99   NodalProjectionData       projection     = diff_flux_proj->projection;
100   CeedOperator              op_rhs;
101   CeedBasis                 basis_diff_flux;
102   CeedElemRestriction       elem_restr_diff_flux_volume, elem_restr_qd;
103   CeedVector                q_data;
104   CeedInt                   num_comp_q, q_data_size;
105   PetscInt                  dim, label_value = 0;
106   DMLabel                   domain_label = NULL;
107   MPI_Comm                  comm         = PetscObjectComm((PetscObject)projection->dm);
108 
109   PetscFunctionBeginUser;
110   // -- Get Pre-requisite things
111   PetscCall(DMGetDimension(projection->dm, &dim));
112   PetscCallCeed(ceed, CeedBasisGetNumComponents(ceed_data->basis_q, &num_comp_q));
113 
114   {  // Get elem_restr_diff_flux and basis_diff_flux
115     CeedOperator     *sub_ops;
116     CeedOperatorField op_field;
117     PetscInt          sub_op_index = 0;  // will be 0 for the volume op
118 
119     PetscCallCeed(ceed, CeedCompositeOperatorGetSubList(user->op_ifunction, &sub_ops));
120     PetscCallCeed(ceed, CeedOperatorGetFieldByName(sub_ops[sub_op_index], "div F_diff", &op_field));
121     PetscCallCeed(ceed, CeedOperatorFieldGetElemRestriction(op_field, &elem_restr_diff_flux_volume));
122     PetscCallCeed(ceed, CeedOperatorFieldGetBasis(op_field, &basis_diff_flux));
123   }
124   PetscCall(QDataGet(ceed, projection->dm, domain_label, label_value, ceed_data->elem_restr_x, ceed_data->basis_x, ceed_data->x_coord, &elem_restr_qd,
125                      &q_data, &q_data_size));
126 
127   {  // Create Projection RHS OperatorApplyContext
128     CeedOperator op_rhs;
129 
130     PetscCheck(diff_flux_proj->CreateRHSOperator_Direct, comm, PETSC_ERR_ARG_WRONGSTATE,
131                "Must define CreateRHSOperator_Direct to use indirect div_diff_flux projection");
132     PetscCall(diff_flux_proj->CreateRHSOperator_Direct(user, ceed_data, diff_flux_proj, &op_rhs));
133     PetscCall(DMCreateLocalVector(projection->dm, &diff_flux_proj->DivDiffFlux_loc));
134     diff_flux_proj->ceed_vec_has_array = PETSC_FALSE;
135     PetscCall(OperatorApplyContextCreate(user->dm, projection->dm, ceed, op_rhs, NULL, NULL, NULL, diff_flux_proj->DivDiffFlux_loc,
136                                          &projection->l2_rhs_ctx));
137     PetscCallCeed(ceed, CeedOperatorDestroy(&op_rhs));
138   }
139 
140   {  // -- Build Mass operator
141     CeedQFunction qf_mass;
142     CeedOperator  op_mass;
143 
144     PetscCall(CreateMassQFunction(ceed, projection->num_comp, q_data_size, &qf_mass));
145     PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_mass, NULL, NULL, &op_mass));
146     PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "u", elem_restr_diff_flux_volume, basis_diff_flux, CEED_VECTOR_ACTIVE));
147     PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "qdata", elem_restr_qd, CEED_BASIS_NONE, q_data));
148     PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "v", elem_restr_diff_flux_volume, basis_diff_flux, CEED_VECTOR_ACTIVE));
149 
150     {  // -- Setup KSP for L^2 projection
151       Mat mat_mass;
152 
153       PetscCall(MatCreateCeed(projection->dm, projection->dm, op_mass, NULL, &mat_mass));
154 
155       PetscCall(KSPCreate(comm, &projection->ksp));
156       PetscCall(KSPSetOptionsPrefix(projection->ksp, "div_diff_flux_projection_"));
157       {  // lumped by default
158         PC pc;
159         PetscCall(KSPGetPC(projection->ksp, &pc));
160         PetscCall(PCSetType(pc, PCJACOBI));
161         PetscCall(PCJacobiSetType(pc, PC_JACOBI_ROWSUM));
162         PetscCall(KSPSetType(projection->ksp, KSPPREONLY));
163       }
164       PetscCall(KSPSetFromOptions_WithMatCeed(projection->ksp, mat_mass));
165       PetscCall(MatDestroy(&mat_mass));
166     }
167     PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_mass));
168     PetscCallCeed(ceed, CeedOperatorDestroy(&op_mass));
169   }
170   PetscCallCeed(ceed, CeedVectorDestroy(&q_data));
171   PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd));
172   PetscCallCeed(ceed, CeedOperatorDestroy(&op_rhs));
173   PetscFunctionReturn(PETSC_SUCCESS);
174 }
175 
176 /**
177   @brief Setup indirect projection of divergence of diffusive flux
178 
179   @param[in]     ceed      `Ceed` context
180   @param[in,out] user      `User` context
181   @param[in]     ceed_data `CeedData` context
182   @param[in]     problem   `ProblemData` context
183 **/
184 static PetscErrorCode DivDiffFluxProjectionSetup_Indirect(Ceed ceed, User user, CeedData ceed_data, ProblemData problem) {
185   DivDiffFluxProjectionData diff_flux_proj = user->diff_flux_proj;
186   NodalProjectionData       projection     = diff_flux_proj->projection;
187   CeedBasis                 basis_diff_flux;
188   CeedElemRestriction       elem_restr_diff_flux, elem_restr_qd;
189   CeedVector                q_data;
190   CeedInt                   num_comp_q, q_data_size;
191   PetscInt                  dim;
192   PetscInt                  label_value = 0, height = 0, dm_field = 0;
193   DMLabel                   domain_label = NULL;
194   MPI_Comm                  comm         = PetscObjectComm((PetscObject)projection->dm);
195 
196   PetscFunctionBeginUser;
197   PetscCall(DMGetDimension(projection->dm, &dim));
198   PetscCallCeed(ceed, CeedBasisGetNumComponents(ceed_data->basis_q, &num_comp_q));
199 
200   PetscCall(DMPlexCeedElemRestrictionCreate(ceed, projection->dm, domain_label, label_value, height, dm_field, &elem_restr_diff_flux));
201   PetscCall(CreateBasisFromPlex(ceed, projection->dm, domain_label, label_value, height, dm_field, &basis_diff_flux));
202   PetscCall(QDataGet(ceed, projection->dm, domain_label, label_value, ceed_data->elem_restr_x, ceed_data->basis_x, ceed_data->x_coord, &elem_restr_qd,
203                      &q_data, &q_data_size));
204 
205   {
206     CeedOperator op_rhs;
207 
208     PetscCheck(diff_flux_proj->CreateRHSOperator_Indirect, comm, PETSC_ERR_ARG_WRONGSTATE,
209                "Must define CreateRHSOperator_Indirect to use indirect div_diff_flux projection");
210     PetscCall(diff_flux_proj->CreateRHSOperator_Indirect(user, ceed_data, diff_flux_proj, &op_rhs));
211     PetscCall(OperatorApplyContextCreate(user->dm, projection->dm, ceed, op_rhs, NULL, NULL, NULL, NULL, &projection->l2_rhs_ctx));
212     PetscCallCeed(ceed, CeedOperatorDestroy(&op_rhs));
213   }
214 
215   {  // -- Build Mass operator
216     CeedQFunction qf_mass;
217     CeedOperator  op_mass;
218 
219     PetscCall(CreateMassQFunction(ceed, projection->num_comp, q_data_size, &qf_mass));
220     PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_mass, NULL, NULL, &op_mass));
221     PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "u", elem_restr_diff_flux, basis_diff_flux, CEED_VECTOR_ACTIVE));
222     PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "qdata", elem_restr_qd, CEED_BASIS_NONE, q_data));
223     PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "v", elem_restr_diff_flux, basis_diff_flux, CEED_VECTOR_ACTIVE));
224 
225     {  // -- Setup KSP for L^2 projection
226       Mat mat_mass;
227 
228       PetscCall(MatCreateCeed(projection->dm, projection->dm, op_mass, NULL, &mat_mass));
229 
230       PetscCall(KSPCreate(comm, &projection->ksp));
231       PetscCall(KSPSetOptionsPrefix(projection->ksp, "div_diff_flux_projection_"));
232       {  // lumped by default
233         PC pc;
234         PetscCall(KSPGetPC(projection->ksp, &pc));
235         PetscCall(PCSetType(pc, PCJACOBI));
236         PetscCall(PCJacobiSetType(pc, PC_JACOBI_ROWSUM));
237         PetscCall(KSPSetType(projection->ksp, KSPPREONLY));
238       }
239       PetscCall(KSPSetFromOptions_WithMatCeed(projection->ksp, mat_mass));
240       PetscCall(MatDestroy(&mat_mass));
241     }
242     PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_mass));
243     PetscCallCeed(ceed, CeedOperatorDestroy(&op_mass));
244   }
245 
246   {  // Create OperatorApplyContext to calculate divergence at quadrature points
247     CeedQFunction       qf_calc_divergence;
248     CeedOperator        op_calc_divergence;
249     CeedElemRestriction elem_restr_div_diff_flux;
250 
251     {  // Get elem_restr_div_diff_flux
252       CeedOperator     *sub_ops;
253       CeedOperatorField op_field;
254       PetscInt          sub_op_index = 0;  // will be 0 for the volume op
255 
256       PetscCallCeed(ceed, CeedCompositeOperatorGetSubList(user->op_ifunction, &sub_ops));
257       PetscCallCeed(ceed, CeedOperatorGetFieldByName(sub_ops[sub_op_index], "div F_diff", &op_field));
258       PetscCallCeed(ceed, CeedOperatorFieldGetElemRestriction(op_field, &elem_restr_div_diff_flux));
259     }
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, CeedQFunctionDestroy(&qf_calc_divergence));
277     PetscCallCeed(ceed, CeedOperatorDestroy(&op_calc_divergence));
278   }
279   PetscCallCeed(ceed, CeedBasisDestroy(&basis_diff_flux));
280   PetscCallCeed(ceed, CeedVectorDestroy(&q_data));
281   PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd));
282   PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_diff_flux));
283   PetscFunctionReturn(PETSC_SUCCESS);
284 }
285 
286 /**
287   @brief Setup projection of divergence of diffusive flux
288 
289   @param[in]     ceed      `Ceed` context
290   @param[in,out] user      `User` context
291   @param[in]     ceed_data `CeedData` context
292   @param[in]     problem   `ProblemData` context
293 **/
294 PetscErrorCode DivDiffFluxProjectionSetup(Ceed ceed, User user, CeedData ceed_data, ProblemData problem) {
295   PetscFunctionBeginUser;
296   switch (user->app_ctx->divFdiffproj_method) {
297     case DIV_DIFF_FLUX_PROJ_DIRECT:
298       PetscCall(DivDiffFluxProjectionSetup_Direct(ceed, user, ceed_data, problem));
299       break;
300     case DIV_DIFF_FLUX_PROJ_INDIRECT:
301       PetscCall(DivDiffFluxProjectionSetup_Indirect(ceed, user, ceed_data, problem));
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   PetscCall(NodalProjectionDataDestroy(diff_flux_proj->projection));
377   PetscCall(OperatorApplyContextDestroy(diff_flux_proj->calc_div_diff_flux));
378   if (diff_flux_proj->ceed_vec_has_array) {
379     PetscCall(VecReadCeedToPetsc(diff_flux_proj->div_diff_flux_ceed, diff_flux_proj->DivDiffFlux_memtype, diff_flux_proj->DivDiffFlux_loc));
380     diff_flux_proj->ceed_vec_has_array = PETSC_FALSE;
381   }
382   PetscCall(CeedVectorDestroy(&diff_flux_proj->div_diff_flux_ceed));
383   PetscCall(VecDestroy(&diff_flux_proj->DivDiffFlux_loc));
384   PetscCall(PetscFree(diff_flux_proj));
385   PetscFunctionReturn(PETSC_SUCCESS);
386 }
387