1 // Copyright (c) 2017-2024, 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_NONE, 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 dm_field = 0; dm_field < diff_filter->num_filtered_fields; dm_field++) { 62 char field_name[PETSC_MAX_PATH_LEN]; 63 CeedElemRestriction elem_restr_filter; 64 CeedBasis basis_filter; 65 DMLabel domain_label = NULL; 66 PetscInt label_value = 0, height = 0; 67 PetscCall(DMPlexCeedElemRestrictionCreate(ceed, dm_filter, domain_label, label_value, height, dm_field, &elem_restr_filter)); 68 PetscCall(CreateBasisFromPlex(ceed, dm_filter, domain_label, label_value, height, dm_field, &basis_filter)); 69 70 PetscCall(PetscSNPrintf(field_name, PETSC_MAX_PATH_LEN, "v%" PetscInt_FMT, dm_field)); 71 PetscCallCeed(ceed, CeedOperatorSetField(op_rhs, field_name, elem_restr_filter, basis_filter, CEED_VECTOR_ACTIVE)); 72 } 73 74 PetscCall(OperatorApplyContextCreate(user->dm, dm_filter, ceed, op_rhs, NULL, NULL, user->Q_loc, NULL, &diff_filter->op_rhs_ctx)); 75 76 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_rhs)); 77 PetscCallCeed(ceed, CeedOperatorDestroy(&op_rhs)); 78 } 79 80 { // Setup LHS Operator and KSP for the differential filtering solve 81 CeedOperator op_lhs; 82 Mat mat_lhs; 83 CeedInt num_comp_qd; 84 PetscInt dim, num_comp_grid_aniso; 85 CeedElemRestriction elem_restr_grid_aniso; 86 CeedVector grid_aniso_ceed; 87 88 PetscCall(DMGetDimension(user->dm, &dim)); 89 PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_qd_i, &num_comp_qd)); 90 91 // -- Get Grid anisotropy tensor 92 PetscCall(GridAnisotropyTensorCalculateCollocatedVector(ceed, user, ceed_data, &elem_restr_grid_aniso, &grid_aniso_ceed, &num_comp_grid_aniso)); 93 94 PetscCallCeed(ceed, CeedCompositeOperatorCreate(ceed, &op_lhs)); 95 for (PetscInt i = 0; i < diff_filter->num_filtered_fields; i++) { 96 CeedQFunction qf_lhs; 97 PetscInt num_comp_filter = diff_filter->num_field_components[i]; 98 CeedOperator op_lhs_sub; 99 CeedElemRestriction elem_restr_filter; 100 CeedBasis basis_filter; 101 102 switch (num_comp_filter) { 103 case 1: 104 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, DifferentialFilter_LHS_1, DifferentialFilter_LHS_1_loc, &qf_lhs)); 105 break; 106 case 5: 107 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, DifferentialFilter_LHS_5, DifferentialFilter_LHS_5_loc, &qf_lhs)); 108 break; 109 case 6: 110 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, DifferentialFilter_LHS_6, DifferentialFilter_LHS_6_loc, &qf_lhs)); 111 break; 112 case 11: 113 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, DifferentialFilter_LHS_11, DifferentialFilter_LHS_11_loc, &qf_lhs)); 114 break; 115 default: 116 SETERRQ(PetscObjectComm((PetscObject)user->dm), PETSC_ERR_SUP, "Differential filtering not available for (%" PetscInt_FMT ") components", 117 num_comp_filter); 118 } 119 120 PetscCallCeed(ceed, CeedQFunctionSetContext(qf_lhs, diff_filter_qfctx)); 121 PetscCallCeed(ceed, CeedQFunctionSetUserFlopsEstimate(qf_lhs, 0)); 122 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_lhs, "q", num_comp_filter, CEED_EVAL_INTERP)); 123 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_lhs, "Grad_q", num_comp_filter * dim, CEED_EVAL_GRAD)); 124 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_lhs, "anisotropy tensor", num_comp_grid_aniso, CEED_EVAL_NONE)); 125 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_lhs, "x", num_comp_x, CEED_EVAL_INTERP)); 126 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_lhs, "qdata", num_comp_qd, CEED_EVAL_NONE)); 127 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_lhs, "v", num_comp_filter, CEED_EVAL_INTERP)); 128 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_lhs, "Grad_v", num_comp_filter * dim, CEED_EVAL_GRAD)); 129 130 { 131 CeedOperatorField op_field; 132 char field_name[PETSC_MAX_PATH_LEN]; 133 PetscCall(PetscSNPrintf(field_name, PETSC_MAX_PATH_LEN, "v%" PetscInt_FMT, i)); 134 PetscCallCeed(ceed, CeedOperatorGetFieldByName(diff_filter->op_rhs_ctx->op, field_name, &op_field)); 135 PetscCallCeed(ceed, CeedOperatorFieldGetElemRestriction(op_field, &elem_restr_filter)); 136 PetscCallCeed(ceed, CeedOperatorFieldGetBasis(op_field, &basis_filter)); 137 } 138 139 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_lhs, NULL, NULL, &op_lhs_sub)); 140 PetscCallCeed(ceed, CeedOperatorSetField(op_lhs_sub, "q", elem_restr_filter, basis_filter, CEED_VECTOR_ACTIVE)); 141 PetscCallCeed(ceed, CeedOperatorSetField(op_lhs_sub, "Grad_q", elem_restr_filter, basis_filter, CEED_VECTOR_ACTIVE)); 142 PetscCallCeed(ceed, CeedOperatorSetField(op_lhs_sub, "anisotropy tensor", elem_restr_grid_aniso, CEED_BASIS_NONE, grid_aniso_ceed)); 143 PetscCallCeed(ceed, CeedOperatorSetField(op_lhs_sub, "x", ceed_data->elem_restr_x, ceed_data->basis_x, ceed_data->x_coord)); 144 PetscCallCeed(ceed, CeedOperatorSetField(op_lhs_sub, "qdata", ceed_data->elem_restr_qd_i, CEED_BASIS_NONE, ceed_data->q_data)); 145 PetscCallCeed(ceed, CeedOperatorSetField(op_lhs_sub, "v", elem_restr_filter, basis_filter, CEED_VECTOR_ACTIVE)); 146 PetscCallCeed(ceed, CeedOperatorSetField(op_lhs_sub, "Grad_v", elem_restr_filter, basis_filter, CEED_VECTOR_ACTIVE)); 147 148 PetscCallCeed(ceed, CeedCompositeOperatorAddSub(op_lhs, op_lhs_sub)); 149 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_lhs)); 150 PetscCallCeed(ceed, CeedOperatorDestroy(&op_lhs_sub)); 151 } 152 PetscCall(MatCeedCreate(dm_filter, dm_filter, op_lhs, NULL, &mat_lhs)); 153 154 PetscCall(KSPCreate(PetscObjectComm((PetscObject)dm_filter), &diff_filter->ksp)); 155 PetscCall(KSPSetOptionsPrefix(diff_filter->ksp, "diff_filter_")); 156 { 157 PC pc; 158 PetscCall(KSPGetPC(diff_filter->ksp, &pc)); 159 PetscCall(PCSetType(pc, PCJACOBI)); 160 PetscCall(PCJacobiSetType(pc, PC_JACOBI_DIAGONAL)); 161 PetscCall(KSPSetType(diff_filter->ksp, KSPCG)); 162 PetscCall(KSPSetNormType(diff_filter->ksp, KSP_NORM_NATURAL)); 163 PetscCall(KSPSetTolerances(diff_filter->ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT)); 164 } 165 PetscCall(KSPSetFromOptions_WithMatCeed(diff_filter->ksp, mat_lhs)); 166 167 PetscCall(MatDestroy(&mat_lhs)); 168 PetscCallCeed(ceed, CeedOperatorDestroy(&op_lhs)); 169 } 170 PetscFunctionReturn(PETSC_SUCCESS); 171 } 172 173 // @brief Setup DM, operators, contexts, etc. for performing differential filtering 174 PetscErrorCode DifferentialFilterSetup(Ceed ceed, User user, CeedData ceed_data, ProblemData *problem) { 175 MPI_Comm comm = user->comm; 176 NewtonianIdealGasContext gas; 177 DifferentialFilterContext diff_filter_ctx; 178 CeedQFunctionContext diff_filter_qfctx; 179 180 PetscFunctionBeginUser; 181 PetscCall(PetscNew(&user->diff_filter)); 182 DiffFilterData diff_filter = user->diff_filter; 183 PetscCall(PetscOptionsGetBool(NULL, NULL, "-diff_filter_mms", &diff_filter->do_mms_test, NULL)); 184 185 { // Create DM for filtered quantities 186 PetscSection section; 187 188 PetscCall(DMClone(user->dm, &diff_filter->dm_filter)); 189 PetscCall(PetscObjectSetName((PetscObject)diff_filter->dm_filter, "Differential Filtering")); 190 191 diff_filter->num_filtered_fields = diff_filter->do_mms_test ? 1 : 2; 192 PetscCall(PetscMalloc1(diff_filter->num_filtered_fields, &diff_filter->num_field_components)); 193 194 if (diff_filter->do_mms_test) { 195 PetscInt field_components; 196 diff_filter->num_field_components[0] = field_components = 1; 197 PetscCall(DMSetupByOrder_FEM(PETSC_TRUE, PETSC_TRUE, user->app_ctx->degree, 1, user->app_ctx->q_extra, diff_filter->num_filtered_fields, 198 &field_components, diff_filter->dm_filter)); 199 200 PetscCall(DMGetLocalSection(diff_filter->dm_filter, §ion)); 201 PetscCall(PetscSectionSetFieldName(section, 0, "")); 202 PetscCall(PetscSectionSetComponentName(section, 0, 0, "FilteredPhi")); 203 } else { 204 PetscInt field_components[2]; 205 diff_filter->num_field_components[0] = field_components[0] = DIFF_FILTER_STATE_NUM; 206 diff_filter->num_field_components[1] = field_components[1] = DIFF_FILTER_VELOCITY_SQUARED_NUM; 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 diff_filter->field_prim_state = 0; 211 diff_filter->field_velo_prod = 1; 212 PetscCall(DMGetLocalSection(diff_filter->dm_filter, §ion)); 213 PetscCall(PetscSectionSetFieldName(section, diff_filter->field_prim_state, "Filtered Primitive State Variables")); 214 PetscCall(PetscSectionSetComponentName(section, 0, DIFF_FILTER_PRESSURE, "FilteredPressure")); 215 PetscCall(PetscSectionSetComponentName(section, 0, DIFF_FILTER_VELOCITY_X, "FilteredVelocityX")); 216 PetscCall(PetscSectionSetComponentName(section, 0, DIFF_FILTER_VELOCITY_Y, "FilteredVelocityY")); 217 PetscCall(PetscSectionSetComponentName(section, 0, DIFF_FILTER_VELOCITY_Z, "FilteredVelocityZ")); 218 PetscCall(PetscSectionSetComponentName(section, 0, DIFF_FILTER_TEMPERATURE, "FilteredTemperature")); 219 PetscCall(PetscSectionSetFieldName(section, diff_filter->field_velo_prod, "Filtered Velocity Products")); 220 PetscCall(PetscSectionSetComponentName(section, 1, DIFF_FILTER_VELOCITY_SQUARED_XX, "FilteredVelocitySquaredXX")); 221 PetscCall(PetscSectionSetComponentName(section, 1, DIFF_FILTER_VELOCITY_SQUARED_YY, "FilteredVelocitySquaredYY")); 222 PetscCall(PetscSectionSetComponentName(section, 1, DIFF_FILTER_VELOCITY_SQUARED_ZZ, "FilteredVelocitySquaredZZ")); 223 PetscCall(PetscSectionSetComponentName(section, 1, DIFF_FILTER_VELOCITY_SQUARED_YZ, "FilteredVelocitySquaredYZ")); 224 PetscCall(PetscSectionSetComponentName(section, 1, DIFF_FILTER_VELOCITY_SQUARED_XZ, "FilteredVelocitySquaredXZ")); 225 PetscCall(PetscSectionSetComponentName(section, 1, DIFF_FILTER_VELOCITY_SQUARED_XY, "FilteredVelocitySquaredXY")); 226 } 227 } 228 229 PetscCall(PetscNew(&diff_filter_ctx)); 230 diff_filter_ctx->grid_based_width = false; 231 for (int i = 0; i < 3; i++) diff_filter_ctx->width_scaling[i] = 1; 232 diff_filter_ctx->kernel_scaling = 0.1; 233 diff_filter_ctx->damping_function = DIFF_FILTER_DAMP_NONE; 234 diff_filter_ctx->friction_length = 0; 235 diff_filter_ctx->damping_constant = 25; 236 237 PetscOptionsBegin(comm, NULL, "Differential Filtering Options", NULL); 238 PetscInt narray = 3; 239 PetscCall(PetscOptionsBool("-diff_filter_grid_based_width", "Use filter width based on the grid size", NULL, diff_filter_ctx->grid_based_width, 240 (PetscBool *)&diff_filter_ctx->grid_based_width, NULL)); 241 PetscCall(PetscOptionsRealArray("-diff_filter_width_scaling", "Anisotropic scaling of filter width tensor", NULL, diff_filter_ctx->width_scaling, 242 &narray, NULL)); 243 PetscCall(PetscOptionsReal("-diff_filter_kernel_scaling", "Scaling to make differential kernel size \"equivalent\" to other filter kernels", NULL, 244 diff_filter_ctx->kernel_scaling, &diff_filter_ctx->kernel_scaling, NULL)); 245 PetscCall(PetscOptionsEnum("-diff_filter_wall_damping_function", "Damping function to use at the wall", NULL, DifferentialFilterDampingFunctions, 246 (PetscEnum)(diff_filter_ctx->damping_function), (PetscEnum *)&diff_filter_ctx->damping_function, NULL)); 247 PetscCall(PetscOptionsReal("-diff_filter_wall_damping_constant", "Contant for the wall-damping function", NULL, diff_filter_ctx->damping_constant, 248 &diff_filter_ctx->damping_constant, NULL)); 249 PetscCall(PetscOptionsReal("-diff_filter_friction_length", "Friction length associated with the flow, \\delta_\\nu. For wall-damping functions", 250 NULL, diff_filter_ctx->friction_length, &diff_filter_ctx->friction_length, NULL)); 251 PetscOptionsEnd(); 252 253 Units units = user->units; 254 for (int i = 0; i < 3; i++) diff_filter_ctx->width_scaling[i] *= units->meter; 255 diff_filter_ctx->kernel_scaling *= units->meter; 256 diff_filter_ctx->friction_length *= units->meter; 257 258 // -- Create QFContext 259 PetscCallCeed(ceed, CeedQFunctionContextGetDataRead(problem->apply_vol_ifunction.qfunction_context, CEED_MEM_HOST, &gas)); 260 diff_filter_ctx->gas = *gas; 261 PetscCallCeed(ceed, CeedQFunctionContextRestoreDataRead(problem->apply_vol_ifunction.qfunction_context, &gas)); 262 263 PetscCallCeed(ceed, CeedQFunctionContextCreate(ceed, &diff_filter_qfctx)); 264 PetscCallCeed(ceed, CeedQFunctionContextSetData(diff_filter_qfctx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*diff_filter_ctx), diff_filter_ctx)); 265 PetscCallCeed(ceed, CeedQFunctionContextSetDataDestroy(diff_filter_qfctx, CEED_MEM_HOST, FreeContextPetsc)); 266 267 // -- Setup Operators 268 PetscCall(DifferentialFilterCreateOperators(ceed, user, ceed_data, diff_filter_qfctx)); 269 270 PetscCallCeed(ceed, CeedQFunctionContextDestroy(&diff_filter_qfctx)); 271 PetscFunctionReturn(PETSC_SUCCESS); 272 } 273 274 // @brief Apply differential filter to the solution given by Q 275 PetscErrorCode DifferentialFilterApply(User user, const PetscReal solution_time, const Vec Q, Vec Filtered_Solution) { 276 DiffFilterData diff_filter = user->diff_filter; 277 278 PetscFunctionBeginUser; 279 PetscCall(PetscLogEventBegin(FLUIDS_DifferentialFilter, Q, Filtered_Solution, 0, 0)); 280 PetscCall(UpdateBoundaryValues(user, diff_filter->op_rhs_ctx->X_loc, solution_time)); 281 ApplyCeedOperatorGlobalToGlobal(Q, Filtered_Solution, diff_filter->op_rhs_ctx); 282 PetscCall(VecViewFromOptions(Filtered_Solution, NULL, "-diff_filter_rhs_view")); 283 284 PetscCall(KSPSolve(diff_filter->ksp, Filtered_Solution, Filtered_Solution)); 285 PetscCall(PetscLogEventEnd(FLUIDS_DifferentialFilter, Q, Filtered_Solution, 0, 0)); 286 PetscFunctionReturn(PETSC_SUCCESS); 287 } 288 289 // @brief TSMonitor for just applying differential filtering to the simulation 290 // This runs every time step and is primarily for testing purposes 291 PetscErrorCode TSMonitor_DifferentialFilter(TS ts, PetscInt steps, PetscReal solution_time, Vec Q, void *ctx) { 292 User user = (User)ctx; 293 DiffFilterData diff_filter = user->diff_filter; 294 Vec Filtered_Field; 295 296 PetscFunctionBeginUser; 297 PetscCall(DMGetGlobalVector(diff_filter->dm_filter, &Filtered_Field)); 298 299 PetscCall(DifferentialFilterApply(user, solution_time, Q, Filtered_Field)); 300 PetscCall(VecViewFromOptions(Filtered_Field, NULL, "-diff_filter_view")); 301 if (user->app_ctx->test_type == TESTTYPE_DIFF_FILTER) PetscCall(RegressionTest(user->app_ctx, Filtered_Field)); 302 303 PetscCall(DMRestoreGlobalVector(diff_filter->dm_filter, &Filtered_Field)); 304 PetscFunctionReturn(PETSC_SUCCESS); 305 } 306 307 PetscErrorCode DifferentialFilterDataDestroy(DiffFilterData diff_filter) { 308 PetscFunctionBeginUser; 309 if (!diff_filter) PetscFunctionReturn(PETSC_SUCCESS); 310 311 PetscCall(OperatorApplyContextDestroy(diff_filter->op_rhs_ctx)); 312 PetscCall(DMDestroy(&diff_filter->dm_filter)); 313 PetscCall(KSPDestroy(&diff_filter->ksp)); 314 315 PetscCall(PetscFree(diff_filter->num_field_components)); 316 PetscCall(PetscFree(diff_filter)); 317 PetscFunctionReturn(PETSC_SUCCESS); 318 } 319 320 PetscErrorCode DifferentialFilterMmsICSetup(ProblemData *problem) { 321 PetscFunctionBeginUser; 322 problem->ics.qfunction = DifferentialFilter_MMS_IC; 323 problem->ics.qfunction_loc = DifferentialFilter_MMS_IC_loc; 324 PetscFunctionReturn(PETSC_SUCCESS); 325 } 326