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