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_filter = diff_filter->num_comp_filter; 21 CeedInt num_comp_q, num_comp_qd, dim, num_qpts_1d, num_nodes_1d, num_comp_x; 22 CeedElemRestriction elem_restr_filter; 23 CeedBasis basis_filter; 24 25 PetscFunctionBeginUser; 26 PetscCall(DMGetDimension(user->dm, &dim)); 27 CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_x, &num_comp_x); 28 CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_q, &num_comp_q); 29 CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_qd_i, &num_comp_qd); 30 CeedBasisGetNumQuadraturePoints1D(ceed_data->basis_q, &num_qpts_1d); 31 CeedBasisGetNumNodes1D(ceed_data->basis_q, &num_nodes_1d); 32 33 PetscCall(GetRestrictionForDomain(ceed, dm_filter, 0, 0, 0, num_qpts_1d, 0, &elem_restr_filter, NULL, NULL)); 34 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_filter, num_nodes_1d, num_qpts_1d, CEED_GAUSS, &basis_filter); 35 36 { // -- Create RHS MatopApplyContext 37 CeedQFunction qf_rhs; 38 CeedOperator op_rhs; 39 switch (user->phys->state_var) { 40 case STATEVAR_PRIMITIVE: 41 CeedQFunctionCreateInterior(ceed, 1, DifferentialFilter_RHS_Prim, DifferentialFilter_RHS_Prim_loc, &qf_rhs); 42 break; 43 case STATEVAR_CONSERVATIVE: 44 CeedQFunctionCreateInterior(ceed, 1, DifferentialFilter_RHS_Conserv, DifferentialFilter_RHS_Conserv_loc, &qf_rhs); 45 break; 46 default: 47 SETERRQ(PetscObjectComm((PetscObject)user->dm), PETSC_ERR_SUP, "Differential filtering not available for chosen state variable"); 48 } 49 CeedQFunctionSetContext(qf_rhs, diff_filter_qfctx); 50 CeedQFunctionAddInput(qf_rhs, "q", num_comp_q, CEED_EVAL_INTERP); 51 CeedQFunctionAddInput(qf_rhs, "qdata", num_comp_qd, CEED_EVAL_NONE); 52 CeedQFunctionAddInput(qf_rhs, "x", num_comp_x, CEED_EVAL_INTERP); 53 CeedQFunctionAddOutput(qf_rhs, "v", num_comp_filter, CEED_EVAL_INTERP); 54 55 CeedOperatorCreate(ceed, qf_rhs, NULL, NULL, &op_rhs); 56 CeedOperatorSetField(op_rhs, "q", ceed_data->elem_restr_q, ceed_data->basis_q, CEED_VECTOR_ACTIVE); 57 CeedOperatorSetField(op_rhs, "qdata", ceed_data->elem_restr_qd_i, CEED_BASIS_COLLOCATED, ceed_data->q_data); 58 CeedOperatorSetField(op_rhs, "x", ceed_data->elem_restr_x, ceed_data->basis_x, ceed_data->x_coord); 59 CeedOperatorSetField(op_rhs, "v", elem_restr_filter, basis_filter, CEED_VECTOR_ACTIVE); 60 61 PetscCall(OperatorApplyContextCreate(user->dm, dm_filter, ceed, op_rhs, NULL, NULL, user->Q_loc, NULL, &diff_filter->op_rhs_ctx)); 62 63 CeedQFunctionDestroy(&qf_rhs); 64 CeedOperatorDestroy(&op_rhs); 65 } 66 67 { // Setup LHS Operator and KSP for the differential filtering solve 68 CeedQFunction qf_lhs; 69 CeedOperator op_lhs; 70 OperatorApplyContext mat_ctx; 71 Mat mat_lhs; 72 CeedInt num_comp_qd, dim, num_comp_grid_aniso; 73 CeedElemRestriction elem_restr_grid_aniso; 74 CeedVector grid_aniso_ceed; 75 76 PetscCall(DMGetDimension(user->dm, &dim)); 77 CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_qd_i, &num_comp_qd); 78 79 switch (num_comp_filter) { 80 case 1: 81 CeedQFunctionCreateInterior(ceed, 1, DifferentialFilter_LHS_1, DifferentialFilter_LHS_1_loc, &qf_lhs); 82 break; 83 case 5: 84 CeedQFunctionCreateInterior(ceed, 1, DifferentialFilter_LHS_5, DifferentialFilter_LHS_5_loc, &qf_lhs); 85 break; 86 case 11: 87 CeedQFunctionCreateInterior(ceed, 1, DifferentialFilter_LHS_11, DifferentialFilter_LHS_11_loc, &qf_lhs); 88 break; 89 default: 90 SETERRQ(PetscObjectComm((PetscObject)user->dm), PETSC_ERR_SUP, "Differential filtering not available for (%d) components", num_comp_filter); 91 } 92 93 // -- Get Grid anisotropy tensor 94 PetscCall(GridAnisotropyTensorCalculateCollocatedVector(ceed, user, ceed_data, &elem_restr_grid_aniso, &grid_aniso_ceed, &num_comp_grid_aniso)); 95 96 CeedQFunctionSetContext(qf_lhs, diff_filter_qfctx); 97 CeedQFunctionAddInput(qf_lhs, "q", num_comp_filter, CEED_EVAL_INTERP); 98 CeedQFunctionAddInput(qf_lhs, "Grad_q", num_comp_filter * dim, CEED_EVAL_GRAD); 99 CeedQFunctionAddInput(qf_lhs, "anisotropy tensor", num_comp_grid_aniso, CEED_EVAL_NONE); 100 CeedQFunctionAddInput(qf_lhs, "x", num_comp_x, CEED_EVAL_INTERP); 101 CeedQFunctionAddInput(qf_lhs, "qdata", num_comp_qd, CEED_EVAL_NONE); 102 CeedQFunctionAddOutput(qf_lhs, "v", num_comp_filter, CEED_EVAL_INTERP); 103 CeedQFunctionAddOutput(qf_lhs, "Grad_v", num_comp_filter * dim, CEED_EVAL_GRAD); 104 105 CeedOperatorCreate(ceed, qf_lhs, NULL, NULL, &op_lhs); 106 CeedOperatorSetField(op_lhs, "q", elem_restr_filter, basis_filter, CEED_VECTOR_ACTIVE); 107 CeedOperatorSetField(op_lhs, "Grad_q", elem_restr_filter, basis_filter, CEED_VECTOR_ACTIVE); 108 CeedOperatorSetField(op_lhs, "anisotropy tensor", elem_restr_grid_aniso, CEED_BASIS_COLLOCATED, grid_aniso_ceed); 109 CeedOperatorSetField(op_lhs, "x", ceed_data->elem_restr_x, ceed_data->basis_x, ceed_data->x_coord); 110 CeedOperatorSetField(op_lhs, "qdata", ceed_data->elem_restr_qd_i, CEED_BASIS_COLLOCATED, ceed_data->q_data); 111 CeedOperatorSetField(op_lhs, "v", elem_restr_filter, basis_filter, CEED_VECTOR_ACTIVE); 112 CeedOperatorSetField(op_lhs, "Grad_v", elem_restr_filter, basis_filter, CEED_VECTOR_ACTIVE); 113 114 PetscCall(OperatorApplyContextCreate(dm_filter, dm_filter, ceed, op_lhs, NULL, NULL, NULL, NULL, &mat_ctx)); 115 PetscCall(CreateMatShell_Ceed(mat_ctx, &mat_lhs)); 116 117 PetscCall(KSPCreate(PetscObjectComm((PetscObject)dm_filter), &diff_filter->ksp)); 118 PetscCall(KSPSetOptionsPrefix(diff_filter->ksp, "diff_filter_")); 119 { 120 PC pc; 121 PetscCall(KSPGetPC(diff_filter->ksp, &pc)); 122 PetscCall(PCSetType(pc, PCJACOBI)); 123 PetscCall(PCJacobiSetType(pc, PC_JACOBI_DIAGONAL)); 124 PetscCall(KSPSetType(diff_filter->ksp, KSPCG)); 125 PetscCall(KSPSetNormType(diff_filter->ksp, KSP_NORM_NATURAL)); 126 PetscCall(KSPSetTolerances(diff_filter->ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT)); 127 } 128 PetscCall(KSPSetOperators(diff_filter->ksp, mat_lhs, mat_lhs)); 129 PetscCall(KSPSetFromOptions(diff_filter->ksp)); 130 131 CeedQFunctionDestroy(&qf_lhs); 132 CeedOperatorDestroy(&op_lhs); 133 } 134 135 CeedElemRestrictionDestroy(&elem_restr_filter); 136 CeedBasisDestroy(&basis_filter); 137 PetscFunctionReturn(0); 138 } 139 140 // @brief Setup DM, operators, contexts, etc. for performing differential filtering 141 PetscErrorCode DifferentialFilterSetup(Ceed ceed, User user, CeedData ceed_data, ProblemData *problem) { 142 MPI_Comm comm = user->comm; 143 NewtonianIdealGasContext gas; 144 DifferentialFilterContext diff_filter_ctx; 145 CeedQFunctionContext diff_filter_qfctx; 146 147 PetscFunctionBeginUser; 148 PetscCall(PetscNew(&user->diff_filter)); 149 150 { // Create DM for filtered quantities 151 PetscFE fe; 152 PetscSection section; 153 PetscInt dim; 154 155 user->diff_filter->num_comp_filter = DIFF_FILTER_NUM_COMPONENTS; 156 157 PetscCall(DMClone(user->dm, &user->diff_filter->dm_filter)); 158 PetscCall(DMGetDimension(user->diff_filter->dm_filter, &dim)); 159 PetscCall(PetscObjectSetName((PetscObject)user->diff_filter->dm_filter, "Differential Filtering")); 160 161 PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, dim, user->diff_filter->num_comp_filter, PETSC_FALSE, user->app_ctx->degree, PETSC_DECIDE, &fe)); 162 PetscCall(PetscObjectSetName((PetscObject)fe, "Differential Filtering")); 163 PetscCall(DMAddField(user->diff_filter->dm_filter, NULL, (PetscObject)fe)); 164 PetscCall(DMCreateDS(user->diff_filter->dm_filter)); 165 PetscCall(DMPlexSetClosurePermutationTensor(user->diff_filter->dm_filter, PETSC_DETERMINE, NULL)); 166 167 PetscCall(DMGetLocalSection(user->diff_filter->dm_filter, §ion)); 168 PetscCall(PetscSectionSetFieldName(section, 0, "")); 169 PetscCall(PetscSectionSetComponentName(section, 0, DIFF_FILTER_PRESSURE, "FilteredPressure")); 170 PetscCall(PetscSectionSetComponentName(section, 0, DIFF_FILTER_VELOCITY_X, "FilteredVelocityX")); 171 PetscCall(PetscSectionSetComponentName(section, 0, DIFF_FILTER_VELOCITY_Y, "FilteredVelocityY")); 172 PetscCall(PetscSectionSetComponentName(section, 0, DIFF_FILTER_VELOCITY_Z, "FilteredVelocityZ")); 173 PetscCall(PetscSectionSetComponentName(section, 0, DIFF_FILTER_TEMPERATURE, "FilteredTemperature")); 174 PetscCall(PetscSectionSetComponentName(section, 0, DIFF_FILTER_VELOCITY_SQUARED_XX, "FilteredVelocitySquaredXX")); 175 PetscCall(PetscSectionSetComponentName(section, 0, DIFF_FILTER_VELOCITY_SQUARED_YY, "FilteredVelocitySquaredYY")); 176 PetscCall(PetscSectionSetComponentName(section, 0, DIFF_FILTER_VELOCITY_SQUARED_ZZ, "FilteredVelocitySquaredZZ")); 177 PetscCall(PetscSectionSetComponentName(section, 0, DIFF_FILTER_VELOCITY_SQUARED_YZ, "FilteredVelocitySquaredYZ")); 178 PetscCall(PetscSectionSetComponentName(section, 0, DIFF_FILTER_VELOCITY_SQUARED_XZ, "FilteredVelocitySquaredXZ")); 179 PetscCall(PetscSectionSetComponentName(section, 0, DIFF_FILTER_VELOCITY_SQUARED_XY, "FilteredVelocitySquaredXY")); 180 181 PetscCall(PetscFEDestroy(&fe)); 182 } 183 184 PetscCall(PetscNew(&diff_filter_ctx)); 185 diff_filter_ctx->grid_based_width = false; 186 for (int i = 0; i < 3; i++) diff_filter_ctx->width_scaling[i] = 1; 187 diff_filter_ctx->kernel_scaling = 0.1; 188 diff_filter_ctx->damping_function = DIFF_FILTER_DAMP_NONE; 189 diff_filter_ctx->friction_length = 0; 190 diff_filter_ctx->damping_constant = 25; 191 192 PetscOptionsBegin(comm, NULL, "Differential Filtering Options", NULL); 193 PetscInt narray = 3; 194 PetscCall(PetscOptionsBool("-diff_filter_grid_based_width", "Use filter width based on the grid size", NULL, diff_filter_ctx->grid_based_width, 195 (PetscBool *)&diff_filter_ctx->grid_based_width, NULL)); 196 PetscCall(PetscOptionsRealArray("-diff_filter_width_scaling", "Anisotropic scaling of filter width tensor", NULL, diff_filter_ctx->width_scaling, 197 &narray, NULL)); 198 PetscCall(PetscOptionsReal("-diff_filter_kernel_scaling", "Scaling to make differential kernel size \"equivalent\" to other filter kernels", NULL, 199 diff_filter_ctx->kernel_scaling, &diff_filter_ctx->kernel_scaling, NULL)); 200 PetscCall(PetscOptionsEnum("-diff_filter_wall_damping_function", "Damping function to use at the wall", NULL, DifferentialFilterDampingFunctions, 201 (PetscEnum)(diff_filter_ctx->damping_function), (PetscEnum *)&diff_filter_ctx->damping_function, NULL)); 202 PetscCall(PetscOptionsReal("-diff_filter_wall_damping_constant", "Contant for the wall-damping function", NULL, diff_filter_ctx->damping_constant, 203 &diff_filter_ctx->damping_constant, NULL)); 204 PetscCall(PetscOptionsReal("-diff_filter_friction_length", "Friction length associated with the flow, \\delta_\\nu. For wall-damping functions", 205 NULL, diff_filter_ctx->friction_length, &diff_filter_ctx->friction_length, NULL)); 206 PetscOptionsEnd(); 207 208 Units units = user->units; 209 for (int i = 0; i < 3; i++) diff_filter_ctx->width_scaling[i] *= units->meter; 210 diff_filter_ctx->kernel_scaling *= units->meter; 211 diff_filter_ctx->friction_length *= units->meter; 212 213 // -- Create QFContext 214 CeedQFunctionContextGetDataRead(problem->apply_vol_ifunction.qfunction_context, CEED_MEM_HOST, &gas); 215 diff_filter_ctx->gas = *gas; 216 CeedQFunctionContextRestoreDataRead(problem->apply_vol_ifunction.qfunction_context, &gas); 217 218 CeedQFunctionContextCreate(ceed, &diff_filter_qfctx); 219 CeedQFunctionContextSetData(diff_filter_qfctx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*diff_filter_ctx), diff_filter_ctx); 220 CeedQFunctionContextSetDataDestroy(diff_filter_qfctx, CEED_MEM_HOST, FreeContextPetsc); 221 222 // -- Setup Operators 223 PetscCall(DifferentialFilterCreateOperators(ceed, user, ceed_data, diff_filter_qfctx)); 224 225 CeedQFunctionContextDestroy(&diff_filter_qfctx); 226 PetscFunctionReturn(0); 227 } 228 229 // @brief Apply differential filter to the solution given by Q 230 PetscErrorCode DifferentialFilterApply(User user, const PetscReal solution_time, const Vec Q, Vec Filtered_Solution) { 231 DiffFilterData diff_filter = user->diff_filter; 232 233 PetscFunctionBeginUser; 234 PetscCall(UpdateBoundaryValues(user, diff_filter->op_rhs_ctx->X_loc, solution_time)); 235 ApplyCeedOperatorGlobalToGlobal(Q, Filtered_Solution, diff_filter->op_rhs_ctx); 236 PetscCall(VecViewFromOptions(Filtered_Solution, NULL, "-diff_filter_rhs_view")); 237 238 PetscCall(KSPSolve(diff_filter->ksp, Filtered_Solution, Filtered_Solution)); 239 240 PetscFunctionReturn(0); 241 } 242 243 // @brief TSMonitor for just applying differential filtering to the simulation 244 // This runs every time step and is primarily for testing purposes 245 PetscErrorCode TSMonitor_DifferentialFilter(TS ts, PetscInt steps, PetscReal solution_time, Vec Q, void *ctx) { 246 User user = (User)ctx; 247 DiffFilterData diff_filter = user->diff_filter; 248 Vec Filtered_Field; 249 250 PetscFunctionBeginUser; 251 PetscCall(DMGetGlobalVector(diff_filter->dm_filter, &Filtered_Field)); 252 253 PetscCall(DifferentialFilterApply(user, solution_time, Q, Filtered_Field)); 254 PetscCall(VecViewFromOptions(Filtered_Field, NULL, "-diff_filter_view")); 255 256 if (user->app_ctx->test_type == TESTTYPE_DIFF_FILTER) PetscCall(RegressionTests_NS(user->app_ctx, Filtered_Field)); 257 258 PetscCall(DMRestoreGlobalVector(diff_filter->dm_filter, &Filtered_Field)); 259 260 PetscFunctionReturn(0); 261 } 262 263 PetscErrorCode DifferentialFilterDataDestroy(DiffFilterData diff_filter) { 264 PetscFunctionBeginUser; 265 if (!diff_filter) PetscFunctionReturn(0); 266 267 OperatorApplyContextDestroy(diff_filter->op_rhs_ctx); 268 PetscCall(DMDestroy(&diff_filter->dm_filter)); 269 PetscCall(KSPDestroy(&diff_filter->ksp)); 270 271 PetscCall(PetscFree(diff_filter)); 272 273 PetscFunctionReturn(0); 274 } 275