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