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