xref: /honee/src/differential_filter.c (revision 0d0968c28710ac8ae311048b492b03f8a4ab28eb) 
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 performing differential filtering
5 
6 #include "../qfunctions/differential_filter.h"
7 #include <ceed.h>
8 
9 #include <petscdmplex.h>
10 
11 #include <navierstokes.h>
12 
13 // @brief Create RHS and LHS operators for differential filtering
14 PetscErrorCode DifferentialFilterCreateOperators(Ceed ceed, Honee honee, CeedQFunctionContext diff_filter_qfctx) {
15   DiffFilterData      diff_filter = honee->diff_filter;
16   DM                  dm_filter   = diff_filter->dm_filter;
17   CeedInt             num_comp_q, q_data_size, num_comp_x;
18   PetscInt            dim;
19   CeedVector          q_data;
20   CeedElemRestriction elem_restr_qd;
21   DMLabel             domain_label = NULL;
22   PetscInt            label_value = 0, height = 0;
23 
24   PetscFunctionBeginUser;
25   PetscCall(DMGetDimension(honee->dm, &dim));
26   PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(honee->elem_restr_x, &num_comp_x));
27   PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(honee->elem_restr_q, &num_comp_q));
28 
29   PetscCall(QDataGet(ceed, dm_filter, domain_label, label_value, honee->elem_restr_x, honee->basis_x, honee->x_coord, &elem_restr_qd, &q_data,
30                      &q_data_size));
31 
32   {  // -- Create RHS MatopApplyContext
33     CeedQFunction qf_rhs;
34     CeedOperator  op_rhs;
35     switch (honee->phys->state_var) {
36       case STATEVAR_PRIMITIVE:
37         PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, DifferentialFilter_RHS_Prim, DifferentialFilter_RHS_Prim_loc, &qf_rhs));
38         break;
39       case STATEVAR_CONSERVATIVE:
40         PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, DifferentialFilter_RHS_Conserv, DifferentialFilter_RHS_Conserv_loc, &qf_rhs));
41         break;
42       case STATEVAR_ENTROPY:
43         PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, DifferentialFilter_RHS_Entropy, DifferentialFilter_RHS_Entropy_loc, &qf_rhs));
44         break;
45     }
46     if (diff_filter->do_mms_test) {
47       PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_rhs));
48       PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, DifferentialFilter_MMS_RHS, DifferentialFilter_MMS_RHS_loc, &qf_rhs));
49     }
50 
51     PetscCallCeed(ceed, CeedQFunctionSetContext(qf_rhs, diff_filter_qfctx));
52     PetscCallCeed(ceed, CeedQFunctionAddInput(qf_rhs, "q", num_comp_q, CEED_EVAL_INTERP));
53     PetscCallCeed(ceed, CeedQFunctionAddInput(qf_rhs, "qdata", q_data_size, CEED_EVAL_NONE));
54     PetscCallCeed(ceed, CeedQFunctionAddInput(qf_rhs, "x", num_comp_x, CEED_EVAL_INTERP));
55     for (PetscInt i = 0; i < diff_filter->num_filtered_fields; i++) {
56       char field_name[PETSC_MAX_PATH_LEN];
57 
58       PetscCall(PetscSNPrintf(field_name, PETSC_MAX_PATH_LEN, "v%" PetscInt_FMT, i));
59       PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_rhs, field_name, diff_filter->num_field_components[i], CEED_EVAL_INTERP));
60     }
61 
62     PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_rhs, NULL, NULL, &op_rhs));
63     PetscCallCeed(ceed, CeedOperatorSetField(op_rhs, "q", honee->elem_restr_q, honee->basis_q, CEED_VECTOR_ACTIVE));
64     PetscCallCeed(ceed, CeedOperatorSetField(op_rhs, "qdata", elem_restr_qd, CEED_BASIS_NONE, q_data));
65     PetscCallCeed(ceed, CeedOperatorSetField(op_rhs, "x", honee->elem_restr_x, honee->basis_x, honee->x_coord));
66     for (PetscInt dm_field = 0; dm_field < diff_filter->num_filtered_fields; dm_field++) {
67       char                field_name[PETSC_MAX_PATH_LEN];
68       CeedElemRestriction elem_restr_filter;
69       CeedBasis           basis_filter;
70 
71       PetscCall(DMPlexCeedElemRestrictionCreate(ceed, dm_filter, domain_label, label_value, height, dm_field, &elem_restr_filter));
72       PetscCall(CreateBasisFromPlex(ceed, dm_filter, domain_label, label_value, height, dm_field, &basis_filter));
73 
74       PetscCall(PetscSNPrintf(field_name, PETSC_MAX_PATH_LEN, "v%" PetscInt_FMT, dm_field));
75       PetscCallCeed(ceed, CeedOperatorSetField(op_rhs, field_name, elem_restr_filter, basis_filter, CEED_VECTOR_ACTIVE));
76 
77       PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_filter));
78       PetscCallCeed(ceed, CeedBasisDestroy(&basis_filter));
79     }
80 
81     PetscCall(OperatorApplyContextCreate(honee->dm, dm_filter, ceed, op_rhs, NULL, NULL, honee->Q_loc, NULL, &diff_filter->op_rhs_ctx));
82 
83     PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_rhs));
84     PetscCallCeed(ceed, CeedOperatorDestroy(&op_rhs));
85   }
86 
87   {  // Setup LHS Operator and KSP for the differential filtering solve
88     CeedOperator        op_lhs;
89     Mat                 mat_lhs;
90     PetscInt            dim, num_comp_grid_aniso;
91     CeedElemRestriction elem_restr_grid_aniso;
92     CeedVector          grid_aniso_ceed;
93 
94     PetscCall(DMGetDimension(honee->dm, &dim));
95     PetscCall(GridAnisotropyTensorCalculateCollocatedVector(ceed, honee, &elem_restr_grid_aniso, &grid_aniso_ceed, &num_comp_grid_aniso));
96 
97     PetscCallCeed(ceed, CeedCompositeOperatorCreate(ceed, &op_lhs));
98     for (PetscInt i = 0; i < diff_filter->num_filtered_fields; i++) {
99       CeedQFunction       qf_lhs;
100       PetscInt            num_comp_filter = diff_filter->num_field_components[i];
101       CeedOperator        op_lhs_sub;
102       CeedElemRestriction elem_restr_filter;
103       CeedBasis           basis_filter;
104 
105       switch (num_comp_filter) {
106         case 1:
107           PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, DifferentialFilter_LHS_1, DifferentialFilter_LHS_1_loc, &qf_lhs));
108           break;
109         case 5:
110           PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, DifferentialFilter_LHS_5, DifferentialFilter_LHS_5_loc, &qf_lhs));
111           break;
112         case 6:
113           PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, DifferentialFilter_LHS_6, DifferentialFilter_LHS_6_loc, &qf_lhs));
114           break;
115         case 11:
116           PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, DifferentialFilter_LHS_11, DifferentialFilter_LHS_11_loc, &qf_lhs));
117           break;
118         default:
119           SETERRQ(PetscObjectComm((PetscObject)honee->dm), PETSC_ERR_SUP, "Differential filtering not available for (%" PetscInt_FMT ") components",
120                   num_comp_filter);
121       }
122 
123       PetscCallCeed(ceed, CeedQFunctionSetContext(qf_lhs, diff_filter_qfctx));
124       PetscCallCeed(ceed, CeedQFunctionSetUserFlopsEstimate(qf_lhs, 0));
125       PetscCallCeed(ceed, CeedQFunctionAddInput(qf_lhs, "q", num_comp_filter, CEED_EVAL_INTERP));
126       PetscCallCeed(ceed, CeedQFunctionAddInput(qf_lhs, "Grad_q", num_comp_filter * dim, CEED_EVAL_GRAD));
127       PetscCallCeed(ceed, CeedQFunctionAddInput(qf_lhs, "anisotropy tensor", num_comp_grid_aniso, CEED_EVAL_NONE));
128       PetscCallCeed(ceed, CeedQFunctionAddInput(qf_lhs, "x", num_comp_x, CEED_EVAL_INTERP));
129       PetscCallCeed(ceed, CeedQFunctionAddInput(qf_lhs, "qdata", q_data_size, CEED_EVAL_NONE));
130       PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_lhs, "v", num_comp_filter, CEED_EVAL_INTERP));
131       PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_lhs, "Grad_v", num_comp_filter * dim, CEED_EVAL_GRAD));
132 
133       {
134         CeedOperatorField op_field;
135         char              field_name[PETSC_MAX_PATH_LEN];
136 
137         PetscCall(PetscSNPrintf(field_name, PETSC_MAX_PATH_LEN, "v%" PetscInt_FMT, i));
138         PetscCallCeed(ceed, CeedOperatorGetFieldByName(diff_filter->op_rhs_ctx->op, field_name, &op_field));
139         PetscCallCeed(ceed, CeedOperatorFieldGetData(op_field, NULL, &elem_restr_filter, &basis_filter, NULL));
140       }
141 
142       PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_lhs, NULL, NULL, &op_lhs_sub));
143       PetscCallCeed(ceed, CeedOperatorSetField(op_lhs_sub, "q", elem_restr_filter, basis_filter, CEED_VECTOR_ACTIVE));
144       PetscCallCeed(ceed, CeedOperatorSetField(op_lhs_sub, "Grad_q", elem_restr_filter, basis_filter, CEED_VECTOR_ACTIVE));
145       PetscCallCeed(ceed, CeedOperatorSetField(op_lhs_sub, "anisotropy tensor", elem_restr_grid_aniso, CEED_BASIS_NONE, grid_aniso_ceed));
146       PetscCallCeed(ceed, CeedOperatorSetField(op_lhs_sub, "x", honee->elem_restr_x, honee->basis_x, honee->x_coord));
147       PetscCallCeed(ceed, CeedOperatorSetField(op_lhs_sub, "qdata", elem_restr_qd, CEED_BASIS_NONE, q_data));
148       PetscCallCeed(ceed, CeedOperatorSetField(op_lhs_sub, "v", elem_restr_filter, basis_filter, CEED_VECTOR_ACTIVE));
149       PetscCallCeed(ceed, CeedOperatorSetField(op_lhs_sub, "Grad_v", elem_restr_filter, basis_filter, CEED_VECTOR_ACTIVE));
150 
151       PetscCallCeed(ceed, CeedCompositeOperatorAddSub(op_lhs, op_lhs_sub));
152       PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_filter));
153       PetscCallCeed(ceed, CeedBasisDestroy(&basis_filter));
154       PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_lhs));
155       PetscCallCeed(ceed, CeedOperatorDestroy(&op_lhs_sub));
156     }
157     PetscCallCeed(ceed, CeedVectorDestroy(&grid_aniso_ceed));
158     PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_grid_aniso));
159 
160     PetscCallCeed(ceed, CeedOperatorGetContextFieldLabel(op_lhs, "filter width scaling", &diff_filter->filter_width_scaling_label));
161     PetscCall(MatCreateCeed(dm_filter, dm_filter, op_lhs, NULL, &mat_lhs));
162 
163     PetscCall(KSPCreate(PetscObjectComm((PetscObject)dm_filter), &diff_filter->ksp));
164     PetscCall(KSPSetOptionsPrefix(diff_filter->ksp, "diff_filter_"));
165     {
166       PC pc;
167       PetscCall(KSPGetPC(diff_filter->ksp, &pc));
168       PetscCall(PCSetType(pc, PCJACOBI));
169       PetscCall(PCJacobiSetType(pc, PC_JACOBI_DIAGONAL));
170       PetscCall(KSPSetType(diff_filter->ksp, KSPCG));
171       PetscCall(KSPSetNormType(diff_filter->ksp, KSP_NORM_NATURAL));
172       PetscCall(KSPSetTolerances(diff_filter->ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT));
173     }
174     PetscCall(KSPSetFromOptions_WithMatCeed(diff_filter->ksp, mat_lhs));
175 
176     PetscCall(MatDestroy(&mat_lhs));
177     PetscCallCeed(ceed, CeedOperatorDestroy(&op_lhs));
178   }
179 
180   PetscCallCeed(ceed, CeedVectorDestroy(&q_data));
181   PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd));
182   PetscFunctionReturn(PETSC_SUCCESS);
183 }
184 
185 // @brief Setup DM, operators, contexts, etc. for performing differential filtering
186 PetscErrorCode DifferentialFilterSetup(Ceed ceed, Honee honee, ProblemData problem) {
187   MPI_Comm                  comm = honee->comm;
188   NewtonianIdealGasContext  gas;
189   DifferentialFilterContext diff_filter_ctx;
190   CeedQFunctionContext      diff_filter_qfctx;
191 
192   PetscFunctionBeginUser;
193   PetscCall(PetscNew(&honee->diff_filter));
194   DiffFilterData diff_filter = honee->diff_filter;
195   PetscCall(PetscOptionsGetBool(NULL, NULL, "-diff_filter_mms", &diff_filter->do_mms_test, NULL));
196 
197   {  // Create DM for filtered quantities
198     PetscSection section;
199 
200     PetscCall(DMClone(honee->dm, &diff_filter->dm_filter));
201     PetscCall(DMSetMatrixPreallocateSkip(diff_filter->dm_filter, PETSC_TRUE));
202     PetscCall(PetscObjectSetName((PetscObject)diff_filter->dm_filter, "Differential Filtering"));
203 
204     diff_filter->num_filtered_fields = diff_filter->do_mms_test ? 1 : 2;
205     PetscCall(PetscMalloc1(diff_filter->num_filtered_fields, &diff_filter->num_field_components));
206 
207     if (diff_filter->do_mms_test) {
208       PetscInt field_components;
209       diff_filter->num_field_components[0] = field_components = 1;
210       PetscCall(DMSetupByOrder_FEM(PETSC_TRUE, PETSC_TRUE, honee->app_ctx->degree, 1, honee->app_ctx->q_extra, diff_filter->num_filtered_fields,
211                                    &field_components, diff_filter->dm_filter));
212 
213       PetscCall(DMGetLocalSection(diff_filter->dm_filter, &section));
214       PetscCall(PetscSectionSetFieldName(section, 0, ""));
215       PetscCall(PetscSectionSetComponentName(section, 0, 0, "FilteredPhi"));
216     } else {
217       PetscInt field_components[2];
218       diff_filter->num_field_components[0] = field_components[0] = DIFF_FILTER_STATE_NUM;
219       diff_filter->num_field_components[1] = field_components[1] = DIFF_FILTER_VELOCITY_SQUARED_NUM;
220       PetscCall(DMSetupByOrder_FEM(PETSC_TRUE, PETSC_TRUE, honee->app_ctx->degree, 1, honee->app_ctx->q_extra, diff_filter->num_filtered_fields,
221                                    field_components, diff_filter->dm_filter));
222 
223       diff_filter->field_prim_state = 0;
224       diff_filter->field_velo_prod  = 1;
225       PetscCall(DMGetLocalSection(diff_filter->dm_filter, &section));
226       PetscCall(PetscSectionSetFieldName(section, diff_filter->field_prim_state, "Filtered Primitive State Variables"));
227       PetscCall(PetscSectionSetComponentName(section, 0, DIFF_FILTER_PRESSURE, "FilteredPressure"));
228       PetscCall(PetscSectionSetComponentName(section, 0, DIFF_FILTER_VELOCITY_X, "FilteredVelocityX"));
229       PetscCall(PetscSectionSetComponentName(section, 0, DIFF_FILTER_VELOCITY_Y, "FilteredVelocityY"));
230       PetscCall(PetscSectionSetComponentName(section, 0, DIFF_FILTER_VELOCITY_Z, "FilteredVelocityZ"));
231       PetscCall(PetscSectionSetComponentName(section, 0, DIFF_FILTER_TEMPERATURE, "FilteredTemperature"));
232       PetscCall(PetscSectionSetFieldName(section, diff_filter->field_velo_prod, "Filtered Velocity Products"));
233       PetscCall(PetscSectionSetComponentName(section, 1, DIFF_FILTER_VELOCITY_SQUARED_XX, "FilteredVelocitySquaredXX"));
234       PetscCall(PetscSectionSetComponentName(section, 1, DIFF_FILTER_VELOCITY_SQUARED_YY, "FilteredVelocitySquaredYY"));
235       PetscCall(PetscSectionSetComponentName(section, 1, DIFF_FILTER_VELOCITY_SQUARED_ZZ, "FilteredVelocitySquaredZZ"));
236       PetscCall(PetscSectionSetComponentName(section, 1, DIFF_FILTER_VELOCITY_SQUARED_YZ, "FilteredVelocitySquaredYZ"));
237       PetscCall(PetscSectionSetComponentName(section, 1, DIFF_FILTER_VELOCITY_SQUARED_XZ, "FilteredVelocitySquaredXZ"));
238       PetscCall(PetscSectionSetComponentName(section, 1, DIFF_FILTER_VELOCITY_SQUARED_XY, "FilteredVelocitySquaredXY"));
239     }
240   }
241 
242   PetscCall(PetscNew(&diff_filter_ctx));
243   diff_filter_ctx->grid_based_width = false;
244   for (int i = 0; i < 3; i++) diff_filter_ctx->width_scaling[i] = 1;
245   diff_filter_ctx->kernel_scaling   = 0.1;
246   diff_filter_ctx->damping_function = DIFF_FILTER_DAMP_NONE;
247   diff_filter_ctx->friction_length  = 0;
248   diff_filter_ctx->damping_constant = 25;
249 
250   PetscOptionsBegin(comm, NULL, "Differential Filtering Options", NULL);
251   PetscInt narray = 3;
252   PetscCall(PetscOptionsBool("-diff_filter_grid_based_width", "Use filter width based on the grid size", NULL, diff_filter_ctx->grid_based_width,
253                              (PetscBool *)&diff_filter_ctx->grid_based_width, NULL));
254   PetscCall(PetscOptionsRealArray("-diff_filter_width_scaling", "Anisotropic scaling of filter width tensor", NULL, diff_filter_ctx->width_scaling,
255                                   &narray, NULL));
256   PetscCall(PetscOptionsReal("-diff_filter_kernel_scaling", "Scaling to make differential kernel size \"equivalent\" to other filter kernels", NULL,
257                              diff_filter_ctx->kernel_scaling, &diff_filter_ctx->kernel_scaling, NULL));
258   PetscCall(PetscOptionsEnum("-diff_filter_wall_damping_function", "Damping function to use at the wall", NULL, DifferentialFilterDampingFunctions,
259                              (PetscEnum)(diff_filter_ctx->damping_function), (PetscEnum *)&diff_filter_ctx->damping_function, NULL));
260   PetscCall(PetscOptionsReal("-diff_filter_wall_damping_constant", "Contant for the wall-damping function", NULL, diff_filter_ctx->damping_constant,
261                              &diff_filter_ctx->damping_constant, NULL));
262   PetscCall(PetscOptionsReal("-diff_filter_friction_length", "Friction length associated with the flow, \\delta_\\nu. For wall-damping functions",
263                              NULL, diff_filter_ctx->friction_length, &diff_filter_ctx->friction_length, NULL));
264   PetscOptionsEnd();
265 
266   Units units = honee->units;
267   for (int i = 0; i < 3; i++) diff_filter_ctx->width_scaling[i] *= units->meter;
268   diff_filter_ctx->kernel_scaling *= units->meter;
269   diff_filter_ctx->friction_length *= units->meter;
270 
271   // -- Create QFContext
272   PetscCallCeed(ceed, CeedQFunctionContextGetDataRead(problem->apply_vol_ifunction.qfctx, CEED_MEM_HOST, &gas));
273   diff_filter_ctx->gas = *gas;
274   PetscCallCeed(ceed, CeedQFunctionContextRestoreDataRead(problem->apply_vol_ifunction.qfctx, &gas));
275 
276   PetscCallCeed(ceed, CeedQFunctionContextCreate(ceed, &diff_filter_qfctx));
277   PetscCallCeed(ceed, CeedQFunctionContextSetData(diff_filter_qfctx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*diff_filter_ctx), diff_filter_ctx));
278   PetscCallCeed(ceed, CeedQFunctionContextSetDataDestroy(diff_filter_qfctx, CEED_MEM_HOST, FreeContextPetsc));
279   PetscCallCeed(ceed, CeedQFunctionContextRegisterDouble(
280                           diff_filter_qfctx, "filter width scaling", offsetof(struct DifferentialFilterContext_, width_scaling),
281                           sizeof(diff_filter_ctx->width_scaling) / sizeof(diff_filter_ctx->width_scaling[0]), "Filter width scaling"));
282 
283   // -- Setup Operators
284   PetscCall(DifferentialFilterCreateOperators(ceed, honee, diff_filter_qfctx));
285 
286   PetscCallCeed(ceed, CeedQFunctionContextDestroy(&diff_filter_qfctx));
287   PetscFunctionReturn(PETSC_SUCCESS);
288 }
289 
290 // @brief Apply differential filter to the solution given by Q
291 PetscErrorCode DifferentialFilterApply(Honee honee, const PetscReal solution_time, const Vec Q, Vec Filtered_Solution) {
292   DiffFilterData   diff_filter = honee->diff_filter;
293   PetscObjectState X_loc_state;
294   Vec              RHS;
295 
296   PetscFunctionBeginUser;
297   PetscCall(PetscLogEventBegin(HONEE_DifferentialFilter, Q, Filtered_Solution, 0, 0));
298   PetscCall(DMGetNamedGlobalVector(diff_filter->dm_filter, "RHS", &RHS));
299   PetscCall(UpdateBoundaryValues(honee, diff_filter->op_rhs_ctx->X_loc, solution_time));
300   PetscCall(VecGetState(diff_filter->op_rhs_ctx->X_loc, &X_loc_state));
301   if (X_loc_state != diff_filter->X_loc_state) {
302     PetscCall(ApplyCeedOperatorGlobalToGlobal(Q, RHS, diff_filter->op_rhs_ctx));
303     PetscCall(VecGetState(diff_filter->op_rhs_ctx->X_loc, &X_loc_state));
304     diff_filter->X_loc_state = X_loc_state;
305   }
306   PetscCall(VecViewFromOptions(RHS, NULL, "-diff_filter_rhs_view"));
307 
308   PetscCall(KSPSolve(diff_filter->ksp, RHS, Filtered_Solution));
309   PetscCall(DMRestoreNamedGlobalVector(diff_filter->dm_filter, "RHS", &RHS));
310   PetscCall(PetscLogEventEnd(HONEE_DifferentialFilter, Q, Filtered_Solution, 0, 0));
311   PetscFunctionReturn(PETSC_SUCCESS);
312 }
313 
314 // @brief TSMonitor for just applying differential filtering to the simulation
315 // This runs every time step and is primarily for testing purposes
316 PetscErrorCode TSMonitor_DifferentialFilter(TS ts, PetscInt steps, PetscReal solution_time, Vec Q, void *ctx) {
317   Honee          honee       = (Honee)ctx;
318   DiffFilterData diff_filter = honee->diff_filter;
319   Vec            Filtered_Field;
320 
321   PetscFunctionBeginUser;
322   PetscCall(DMGetGlobalVector(diff_filter->dm_filter, &Filtered_Field));
323 
324   PetscCall(DifferentialFilterApply(honee, solution_time, Q, Filtered_Field));
325   PetscCall(VecViewFromOptions(Filtered_Field, NULL, "-diff_filter_view"));
326   if (honee->app_ctx->test_type == TESTTYPE_DIFF_FILTER) PetscCall(RegressionTest(honee->app_ctx, Filtered_Field));
327 
328   PetscCall(DMRestoreGlobalVector(diff_filter->dm_filter, &Filtered_Field));
329   PetscFunctionReturn(PETSC_SUCCESS);
330 }
331 
332 PetscErrorCode DifferentialFilterDataDestroy(DiffFilterData diff_filter) {
333   PetscFunctionBeginUser;
334   if (!diff_filter) PetscFunctionReturn(PETSC_SUCCESS);
335 
336   PetscCall(OperatorApplyContextDestroy(diff_filter->op_rhs_ctx));
337   PetscCall(DMDestroy(&diff_filter->dm_filter));
338   PetscCall(KSPDestroy(&diff_filter->ksp));
339 
340   PetscCall(PetscFree(diff_filter->num_field_components));
341   PetscCall(PetscFree(diff_filter));
342   PetscFunctionReturn(PETSC_SUCCESS);
343 }
344 
345 PetscErrorCode DifferentialFilterMmsICSetup(ProblemData problem) {
346   PetscFunctionBeginUser;
347   problem->ics.qf_func_ptr = DifferentialFilter_MMS_IC;
348   problem->ics.qf_loc      = DifferentialFilter_MMS_IC_loc;
349   PetscFunctionReturn(PETSC_SUCCESS);
350 }
351