xref: /libCEED/examples/fluids/src/differential_filter.c (revision d4cc18453651bd0f94c1a2e078b2646a92dafdcc) 
1 // Copyright (c) 2017-2026, 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 #include <ceed.h>
12 
13 #include <petscdmplex.h>
14 
15 #include "../navierstokes.h"
16 
17 // @brief Create RHS and LHS operators for differential filtering
DifferentialFilterCreateOperators(Ceed ceed,User user,CeedData ceed_data,CeedQFunctionContext diff_filter_qfctx)18 PetscErrorCode DifferentialFilterCreateOperators(Ceed ceed, User user, CeedData ceed_data, CeedQFunctionContext diff_filter_qfctx) {
19   DiffFilterData diff_filter = user->diff_filter;
20   DM             dm_filter   = diff_filter->dm_filter;
21   CeedInt        num_comp_q, num_comp_qd, num_comp_x;
22   PetscInt       dim;
23 
24   PetscFunctionBeginUser;
25   PetscCall(DMGetDimension(user->dm, &dim));
26   PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_x, &num_comp_x));
27   PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_q, &num_comp_q));
28   PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_qd_i, &num_comp_qd));
29 
30   {  // -- Create RHS MatopApplyContext
31     CeedQFunction qf_rhs;
32     CeedOperator  op_rhs;
33     switch (user->phys->state_var) {
34       case STATEVAR_PRIMITIVE:
35         PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, DifferentialFilter_RHS_Prim, DifferentialFilter_RHS_Prim_loc, &qf_rhs));
36         break;
37       case STATEVAR_CONSERVATIVE:
38         PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, DifferentialFilter_RHS_Conserv, DifferentialFilter_RHS_Conserv_loc, &qf_rhs));
39         break;
40       case STATEVAR_ENTROPY:
41         PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, DifferentialFilter_RHS_Entropy, DifferentialFilter_RHS_Entropy_loc, &qf_rhs));
42         break;
43     }
44     if (diff_filter->do_mms_test) {
45       PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_rhs));
46       PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, DifferentialFilter_MMS_RHS, DifferentialFilter_MMS_RHS_loc, &qf_rhs));
47     }
48 
49     PetscCallCeed(ceed, CeedQFunctionSetContext(qf_rhs, diff_filter_qfctx));
50     PetscCallCeed(ceed, CeedQFunctionAddInput(qf_rhs, "q", num_comp_q, CEED_EVAL_INTERP));
51     PetscCallCeed(ceed, CeedQFunctionAddInput(qf_rhs, "qdata", num_comp_qd, CEED_EVAL_NONE));
52     PetscCallCeed(ceed, CeedQFunctionAddInput(qf_rhs, "x", num_comp_x, CEED_EVAL_INTERP));
53     for (PetscInt i = 0; i < diff_filter->num_filtered_fields; i++) {
54       char field_name[PETSC_MAX_PATH_LEN];
55       PetscCall(PetscSNPrintf(field_name, PETSC_MAX_PATH_LEN, "v%" PetscInt_FMT, i));
56       PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_rhs, field_name, diff_filter->num_field_components[i], CEED_EVAL_INTERP));
57     }
58 
59     PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_rhs, NULL, NULL, &op_rhs));
60     PetscCallCeed(ceed, CeedOperatorSetField(op_rhs, "q", ceed_data->elem_restr_q, ceed_data->basis_q, CEED_VECTOR_ACTIVE));
61     PetscCallCeed(ceed, CeedOperatorSetField(op_rhs, "qdata", ceed_data->elem_restr_qd_i, CEED_BASIS_NONE, ceed_data->q_data));
62     PetscCallCeed(ceed, CeedOperatorSetField(op_rhs, "x", ceed_data->elem_restr_x, ceed_data->basis_x, ceed_data->x_coord));
63     for (PetscInt dm_field = 0; dm_field < diff_filter->num_filtered_fields; dm_field++) {
64       char                field_name[PETSC_MAX_PATH_LEN];
65       CeedElemRestriction elem_restr_filter;
66       CeedBasis           basis_filter;
67       DMLabel             domain_label = NULL;
68       PetscInt            label_value = 0, height = 0;
69       PetscCall(DMPlexCeedElemRestrictionCreate(ceed, dm_filter, domain_label, label_value, height, dm_field, &elem_restr_filter));
70       PetscCall(CreateBasisFromPlex(ceed, dm_filter, domain_label, label_value, height, dm_field, &basis_filter));
71 
72       PetscCall(PetscSNPrintf(field_name, PETSC_MAX_PATH_LEN, "v%" PetscInt_FMT, dm_field));
73       PetscCallCeed(ceed, CeedOperatorSetField(op_rhs, field_name, elem_restr_filter, basis_filter, CEED_VECTOR_ACTIVE));
74 
75       PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_filter));
76       PetscCallCeed(ceed, CeedBasisDestroy(&basis_filter));
77     }
78 
79     PetscCall(OperatorApplyContextCreate(user->dm, dm_filter, ceed, op_rhs, NULL, NULL, user->Q_loc, NULL, &diff_filter->op_rhs_ctx));
80 
81     PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_rhs));
82     PetscCallCeed(ceed, CeedOperatorDestroy(&op_rhs));
83   }
84 
85   {  // Setup LHS Operator and KSP for the differential filtering solve
86     CeedOperator        op_lhs;
87     Mat                 mat_lhs;
88     CeedInt             num_comp_qd;
89     PetscInt            dim, num_comp_grid_aniso;
90     CeedElemRestriction elem_restr_grid_aniso;
91     CeedVector          grid_aniso_ceed;
92 
93     PetscCall(DMGetDimension(user->dm, &dim));
94     PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_qd_i, &num_comp_qd));
95 
96     // -- Get Grid anisotropy tensor
97     PetscCall(GridAnisotropyTensorCalculateCollocatedVector(ceed, user, ceed_data, &elem_restr_grid_aniso, &grid_aniso_ceed, &num_comp_grid_aniso));
98 
99     PetscCallCeed(ceed, CeedOperatorCreateComposite(ceed, &op_lhs));
100     for (PetscInt i = 0; i < diff_filter->num_filtered_fields; i++) {
101       CeedQFunction       qf_lhs;
102       PetscInt            num_comp_filter = diff_filter->num_field_components[i];
103       CeedOperator        op_lhs_sub;
104       CeedElemRestriction elem_restr_filter;
105       CeedBasis           basis_filter;
106 
107       switch (num_comp_filter) {
108         case 1:
109           PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, DifferentialFilter_LHS_1, DifferentialFilter_LHS_1_loc, &qf_lhs));
110           break;
111         case 5:
112           PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, DifferentialFilter_LHS_5, DifferentialFilter_LHS_5_loc, &qf_lhs));
113           break;
114         case 6:
115           PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, DifferentialFilter_LHS_6, DifferentialFilter_LHS_6_loc, &qf_lhs));
116           break;
117         case 11:
118           PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, DifferentialFilter_LHS_11, DifferentialFilter_LHS_11_loc, &qf_lhs));
119           break;
120         default:
121           SETERRQ(PetscObjectComm((PetscObject)user->dm), PETSC_ERR_SUP, "Differential filtering not available for (%" PetscInt_FMT ") components",
122                   num_comp_filter);
123       }
124 
125       PetscCallCeed(ceed, CeedQFunctionSetContext(qf_lhs, diff_filter_qfctx));
126       PetscCallCeed(ceed, CeedQFunctionSetUserFlopsEstimate(qf_lhs, 0));
127       PetscCallCeed(ceed, CeedQFunctionAddInput(qf_lhs, "q", num_comp_filter, CEED_EVAL_INTERP));
128       PetscCallCeed(ceed, CeedQFunctionAddInput(qf_lhs, "Grad_q", num_comp_filter * dim, CEED_EVAL_GRAD));
129       PetscCallCeed(ceed, CeedQFunctionAddInput(qf_lhs, "anisotropy tensor", num_comp_grid_aniso, CEED_EVAL_NONE));
130       PetscCallCeed(ceed, CeedQFunctionAddInput(qf_lhs, "x", num_comp_x, CEED_EVAL_INTERP));
131       PetscCallCeed(ceed, CeedQFunctionAddInput(qf_lhs, "qdata", num_comp_qd, CEED_EVAL_NONE));
132       PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_lhs, "v", num_comp_filter, CEED_EVAL_INTERP));
133       PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_lhs, "Grad_v", num_comp_filter * dim, CEED_EVAL_GRAD));
134 
135       {
136         CeedOperatorField op_field;
137         char              field_name[PETSC_MAX_PATH_LEN];
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", ceed_data->elem_restr_x, ceed_data->basis_x, ceed_data->x_coord));
148       PetscCallCeed(ceed, CeedOperatorSetField(op_lhs_sub, "qdata", ceed_data->elem_restr_qd_i, CEED_BASIS_NONE, ceed_data->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   PetscFunctionReturn(PETSC_SUCCESS);
181 }
182 
183 // @brief Setup DM, operators, contexts, etc. for performing differential filtering
DifferentialFilterSetup(Ceed ceed,User user,CeedData ceed_data,ProblemData problem)184 PetscErrorCode DifferentialFilterSetup(Ceed ceed, User user, CeedData ceed_data, ProblemData problem) {
185   MPI_Comm                  comm = user->comm;
186   NewtonianIdealGasContext  gas;
187   DifferentialFilterContext diff_filter_ctx;
188   CeedQFunctionContext      diff_filter_qfctx;
189 
190   PetscFunctionBeginUser;
191   PetscCall(PetscNew(&user->diff_filter));
192   DiffFilterData diff_filter = user->diff_filter;
193   PetscCall(PetscOptionsGetBool(NULL, NULL, "-diff_filter_mms", &diff_filter->do_mms_test, NULL));
194 
195   {  // Create DM for filtered quantities
196     PetscSection section;
197 
198     PetscCall(DMClone(user->dm, &diff_filter->dm_filter));
199     PetscCall(PetscObjectSetName((PetscObject)diff_filter->dm_filter, "Differential Filtering"));
200 
201     diff_filter->num_filtered_fields = diff_filter->do_mms_test ? 1 : 2;
202     PetscCall(PetscMalloc1(diff_filter->num_filtered_fields, &diff_filter->num_field_components));
203 
204     if (diff_filter->do_mms_test) {
205       PetscInt field_components;
206       diff_filter->num_field_components[0] = field_components = 1;
207       PetscCall(DMSetupByOrder_FEM(PETSC_TRUE, PETSC_TRUE, user->app_ctx->degree, 1, user->app_ctx->q_extra, diff_filter->num_filtered_fields,
208                                    &field_components, diff_filter->dm_filter));
209 
210       PetscCall(DMGetLocalSection(diff_filter->dm_filter, &section));
211       PetscCall(PetscSectionSetFieldName(section, 0, ""));
212       PetscCall(PetscSectionSetComponentName(section, 0, 0, "FilteredPhi"));
213     } else {
214       PetscInt field_components[2];
215       diff_filter->num_field_components[0] = field_components[0] = DIFF_FILTER_STATE_NUM;
216       diff_filter->num_field_components[1] = field_components[1] = DIFF_FILTER_VELOCITY_SQUARED_NUM;
217       PetscCall(DMSetupByOrder_FEM(PETSC_TRUE, PETSC_TRUE, user->app_ctx->degree, 1, user->app_ctx->q_extra, diff_filter->num_filtered_fields,
218                                    field_components, diff_filter->dm_filter));
219 
220       diff_filter->field_prim_state = 0;
221       diff_filter->field_velo_prod  = 1;
222       PetscCall(DMGetLocalSection(diff_filter->dm_filter, &section));
223       PetscCall(PetscSectionSetFieldName(section, diff_filter->field_prim_state, "Filtered Primitive State Variables"));
224       PetscCall(PetscSectionSetComponentName(section, 0, DIFF_FILTER_PRESSURE, "FilteredPressure"));
225       PetscCall(PetscSectionSetComponentName(section, 0, DIFF_FILTER_VELOCITY_X, "FilteredVelocityX"));
226       PetscCall(PetscSectionSetComponentName(section, 0, DIFF_FILTER_VELOCITY_Y, "FilteredVelocityY"));
227       PetscCall(PetscSectionSetComponentName(section, 0, DIFF_FILTER_VELOCITY_Z, "FilteredVelocityZ"));
228       PetscCall(PetscSectionSetComponentName(section, 0, DIFF_FILTER_TEMPERATURE, "FilteredTemperature"));
229       PetscCall(PetscSectionSetFieldName(section, diff_filter->field_velo_prod, "Filtered Velocity Products"));
230       PetscCall(PetscSectionSetComponentName(section, 1, DIFF_FILTER_VELOCITY_SQUARED_XX, "FilteredVelocitySquaredXX"));
231       PetscCall(PetscSectionSetComponentName(section, 1, DIFF_FILTER_VELOCITY_SQUARED_YY, "FilteredVelocitySquaredYY"));
232       PetscCall(PetscSectionSetComponentName(section, 1, DIFF_FILTER_VELOCITY_SQUARED_ZZ, "FilteredVelocitySquaredZZ"));
233       PetscCall(PetscSectionSetComponentName(section, 1, DIFF_FILTER_VELOCITY_SQUARED_YZ, "FilteredVelocitySquaredYZ"));
234       PetscCall(PetscSectionSetComponentName(section, 1, DIFF_FILTER_VELOCITY_SQUARED_XZ, "FilteredVelocitySquaredXZ"));
235       PetscCall(PetscSectionSetComponentName(section, 1, DIFF_FILTER_VELOCITY_SQUARED_XY, "FilteredVelocitySquaredXY"));
236     }
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   PetscCallCeed(ceed, CeedQFunctionContextRegisterDouble(diff_filter_qfctx, "filter width scaling",
277                                                          offsetof(struct DifferentialFilterContext_, width_scaling),
278                                                          sizeof(diff_filter_ctx->width_scaling) / sizeof(diff_filter_ctx->width_scaling[0]),
279                                                          "Filter width scaling"));
280 
281   // -- Setup Operators
282   PetscCall(DifferentialFilterCreateOperators(ceed, user, ceed_data, diff_filter_qfctx));
283 
284   PetscCallCeed(ceed, CeedQFunctionContextDestroy(&diff_filter_qfctx));
285   PetscFunctionReturn(PETSC_SUCCESS);
286 }
287 
288 // @brief Apply differential filter to the solution given by Q
DifferentialFilterApply(User user,const PetscReal solution_time,const Vec Q,Vec Filtered_Solution)289 PetscErrorCode DifferentialFilterApply(User user, const PetscReal solution_time, const Vec Q, Vec Filtered_Solution) {
290   DiffFilterData   diff_filter = user->diff_filter;
291   PetscObjectState X_loc_state;
292   Vec              RHS;
293 
294   PetscFunctionBeginUser;
295   PetscCall(PetscLogEventBegin(FLUIDS_DifferentialFilter, Q, Filtered_Solution, 0, 0));
296   PetscCall(DMGetNamedGlobalVector(diff_filter->dm_filter, "RHS", &RHS));
297   PetscCall(UpdateBoundaryValues(user, diff_filter->op_rhs_ctx->X_loc, solution_time));
298   PetscCall(VecGetState(diff_filter->op_rhs_ctx->X_loc, &X_loc_state));
299   if (X_loc_state != diff_filter->X_loc_state) {
300     PetscCall(ApplyCeedOperatorGlobalToGlobal(Q, RHS, diff_filter->op_rhs_ctx));
301     PetscCall(VecGetState(diff_filter->op_rhs_ctx->X_loc, &X_loc_state));
302     diff_filter->X_loc_state = X_loc_state;
303   }
304   PetscCall(VecViewFromOptions(RHS, NULL, "-diff_filter_rhs_view"));
305 
306   PetscCall(KSPSolve(diff_filter->ksp, RHS, Filtered_Solution));
307   PetscCall(DMRestoreNamedGlobalVector(diff_filter->dm_filter, "RHS", &RHS));
308   PetscCall(PetscLogEventEnd(FLUIDS_DifferentialFilter, Q, Filtered_Solution, 0, 0));
309   PetscFunctionReturn(PETSC_SUCCESS);
310 }
311 
312 // @brief TSMonitor for just applying differential filtering to the simulation
313 // This runs every time step and is primarily for testing purposes
TSMonitor_DifferentialFilter(TS ts,PetscInt steps,PetscReal solution_time,Vec Q,void * ctx)314 PetscErrorCode TSMonitor_DifferentialFilter(TS ts, PetscInt steps, PetscReal solution_time, Vec Q, void *ctx) {
315   User           user        = (User)ctx;
316   DiffFilterData diff_filter = user->diff_filter;
317   Vec            Filtered_Field;
318 
319   PetscFunctionBeginUser;
320   PetscCall(DMGetGlobalVector(diff_filter->dm_filter, &Filtered_Field));
321 
322   PetscCall(DifferentialFilterApply(user, solution_time, Q, Filtered_Field));
323   PetscCall(VecViewFromOptions(Filtered_Field, NULL, "-diff_filter_view"));
324   if (user->app_ctx->test_type == TESTTYPE_DIFF_FILTER) PetscCall(RegressionTest(user->app_ctx, Filtered_Field));
325 
326   PetscCall(DMRestoreGlobalVector(diff_filter->dm_filter, &Filtered_Field));
327   PetscFunctionReturn(PETSC_SUCCESS);
328 }
329 
DifferentialFilterDataDestroy(DiffFilterData diff_filter)330 PetscErrorCode DifferentialFilterDataDestroy(DiffFilterData diff_filter) {
331   PetscFunctionBeginUser;
332   if (!diff_filter) PetscFunctionReturn(PETSC_SUCCESS);
333 
334   PetscCall(OperatorApplyContextDestroy(diff_filter->op_rhs_ctx));
335   PetscCall(DMDestroy(&diff_filter->dm_filter));
336   PetscCall(KSPDestroy(&diff_filter->ksp));
337 
338   PetscCall(PetscFree(diff_filter->num_field_components));
339   PetscCall(PetscFree(diff_filter));
340   PetscFunctionReturn(PETSC_SUCCESS);
341 }
342 
DifferentialFilterMmsICSetup(ProblemData problem)343 PetscErrorCode DifferentialFilterMmsICSetup(ProblemData problem) {
344   PetscFunctionBeginUser;
345   problem->ics.qfunction     = DifferentialFilter_MMS_IC;
346   problem->ics.qfunction_loc = DifferentialFilter_MMS_IC_loc;
347   PetscFunctionReturn(PETSC_SUCCESS);
348 }
349