1 // SPDX-FileCopyrightText: Copyright (c) 2017-2024, HONEE contributors. 2 // SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause 3 4 #include <navierstokes.h> 5 #include "petscerror.h" 6 7 /** 8 @brief Add `BCDefinition` to a `PetscSegBuffer` 9 10 @param[in] bc_def `BCDefinition` to add 11 @param[in,out] bc_defs_seg `PetscSegBuffer` to add to 12 **/ 13 static PetscErrorCode AddBCDefinitionToSegBuffer(BCDefinition bc_def, PetscSegBuffer bc_defs_seg) { 14 BCDefinition *bc_def_ptr; 15 16 PetscFunctionBeginUser; 17 if (bc_def == NULL) PetscFunctionReturn(PETSC_SUCCESS); 18 PetscCall(PetscSegBufferGet(bc_defs_seg, 1, &bc_def_ptr)); 19 *bc_def_ptr = bc_def; 20 PetscFunctionReturn(PETSC_SUCCESS); 21 } 22 23 /** 24 @brief Create and setup `BCDefinition`s from commandline options 25 26 @param[in] honee `Honee` 27 @param[in,out] problem `ProblemData` 28 @param[in] app_ctx `AppCtx` 29 **/ 30 PetscErrorCode BoundaryConditionSetUp(Honee honee, ProblemData problem, AppCtx app_ctx) { 31 PetscSegBuffer bc_defs_seg; 32 PetscBool flg; 33 BCDefinition bc_def; 34 35 PetscFunctionBeginUser; 36 PetscCall(PetscSegBufferCreate(sizeof(BCDefinition), 4, &bc_defs_seg)); 37 38 PetscOptionsBegin(honee->comm, NULL, "Boundary Condition Options", NULL); 39 40 PetscCall(PetscOptionsBCDefinition("-bc_wall", "Face IDs to apply wall BC", NULL, "wall", &bc_def, NULL)); 41 PetscCall(AddBCDefinitionToSegBuffer(bc_def, bc_defs_seg)); 42 if (bc_def) { 43 PetscInt num_essential_comps = 16, essential_comps[16]; 44 45 PetscCall(PetscOptionsIntArray("-wall_comps", "An array of constrained component numbers", NULL, essential_comps, &num_essential_comps, &flg)); 46 PetscCall(BCDefinitionSetEssential(bc_def, num_essential_comps, essential_comps)); 47 48 app_ctx->wall_forces.num_wall = bc_def->num_label_values; 49 PetscCall(PetscMalloc1(bc_def->num_label_values, &app_ctx->wall_forces.walls)); 50 PetscCall(PetscArraycpy(app_ctx->wall_forces.walls, bc_def->label_values, bc_def->num_label_values)); 51 } 52 53 { // Symmetry Boundary Conditions 54 const char *deprecated[3] = {"-bc_slip_x", "-bc_slip_y", "-bc_slip_z"}; 55 const char *flags[3] = {"-bc_symmetry_x", "-bc_symmetry_y", "-bc_symmetry_z"}; 56 57 for (PetscInt j = 0; j < 3; j++) { 58 PetscCall(PetscOptionsDeprecated(deprecated[j], flags[j], "libCEED 0.12.0", 59 "Use -bc_symmetry_[x,y,z] for direct equivalency, or -bc_slip for weak, Riemann-based, direction-invariant " 60 "slip/no-penatration boundary conditions")); 61 PetscCall(PetscOptionsBCDefinition(flags[j], "Face IDs to apply symmetry BC", NULL, "symmetry", &bc_def, NULL)); 62 if (!bc_def) { 63 PetscCall(PetscOptionsBCDefinition(deprecated[j], "Face IDs to apply symmetry BC", NULL, "symmetry", &bc_def, NULL)); 64 } 65 PetscCall(AddBCDefinitionToSegBuffer(bc_def, bc_defs_seg)); 66 if (bc_def) { 67 PetscInt essential_comps[1] = {j + 1}; 68 69 PetscCall(BCDefinitionSetEssential(bc_def, 1, essential_comps)); 70 } 71 } 72 } 73 74 PetscCall(PetscOptionsBCDefinition("-bc_inflow", "Face IDs to apply inflow BC", NULL, "inflow", &bc_def, NULL)); 75 PetscCall(AddBCDefinitionToSegBuffer(bc_def, bc_defs_seg)); 76 77 PetscCall(PetscOptionsBCDefinition("-bc_outflow", "Face IDs to apply outflow BC", NULL, "outflow", &bc_def, NULL)); 78 PetscCall(AddBCDefinitionToSegBuffer(bc_def, bc_defs_seg)); 79 80 PetscCall(PetscOptionsBCDefinition("-bc_freestream", "Face IDs to apply freestream BC", NULL, "freestream", &bc_def, NULL)); 81 PetscCall(AddBCDefinitionToSegBuffer(bc_def, bc_defs_seg)); 82 83 PetscCall(PetscOptionsBCDefinition("-bc_slip", "Face IDs to apply slip BC", NULL, "slip", &bc_def, NULL)); 84 PetscCall(AddBCDefinitionToSegBuffer(bc_def, bc_defs_seg)); 85 86 PetscOptionsEnd(); 87 88 PetscCall(PetscSegBufferGetSize(bc_defs_seg, &problem->num_bc_defs)); 89 PetscCall(PetscSegBufferExtractAlloc(bc_defs_seg, &problem->bc_defs)); 90 PetscCall(PetscSegBufferDestroy(&bc_defs_seg)); 91 92 //TODO: Verify that the BCDefinition don't have overlapping claims to boundary faces 93 94 PetscFunctionReturn(PETSC_SUCCESS); 95 } 96 97 /** 98 @brief Destroy `HoneeBCStruct` object 99 100 @param[in] ctx Pointer to `HoneeBCStruct` 101 **/ 102 PetscErrorCode HoneeBCDestroy(void **ctx) { 103 HoneeBCStruct honee_bc = *(HoneeBCStruct *)ctx; 104 Ceed ceed = honee_bc->honee->ceed; 105 106 PetscFunctionBeginUser; 107 if (honee_bc->qfctx) PetscCallCeed(ceed, CeedQFunctionContextDestroy(&honee_bc->qfctx)); 108 if (honee_bc->DestroyCtx) PetscCall((*(honee_bc->DestroyCtx))(&honee_bc->ctx)); 109 PetscCall(PetscFree(honee_bc)); 110 *ctx = NULL; 111 PetscFunctionReturn(PETSC_SUCCESS); 112 } 113 114 /** 115 @brief Create a QFunction matching the "standard" HONEE inputs for IFunctions 116 117 This assumes the `bc_def` context is a `HoneeBCStruct` and inputs/outputs of the IFunction QFunction are in the correct order. 118 119 @param[in] bc_def `BCDefinition` that hosts the IFunction 120 @param[in] qf_func_ptr `CeedQFunctionUser` for the IFunction 121 @param[in] qf_loc Absolute path to source of `CeedQFunctionUser` 122 @param[in] qfctx `CeedQFunctionContext` for the IFunction (also shared with the IJacobian, if applicable) 123 @param[out] qf_ifunc QFunction for the IFunction 124 **/ 125 PetscErrorCode HoneeBCCreateIFunctionQF(BCDefinition bc_def, CeedQFunctionUser qf_func_ptr, const char *qf_loc, CeedQFunctionContext qfctx, 126 CeedQFunction *qf_ifunc) { 127 Ceed ceed; 128 DM dm; 129 PetscInt dim, dim_sur, height = 1, num_comp_x, num_comp_q; 130 CeedInt q_data_size_sur, jac_data_size_sur; 131 HoneeBCStruct honee_bc; 132 133 PetscFunctionBeginUser; 134 PetscCall(BCDefinitionGetDM(bc_def, &dm)); 135 PetscCall(BCDefinitionGetContext(bc_def, &honee_bc)); 136 ceed = honee_bc->honee->ceed; 137 jac_data_size_sur = honee_bc->jac_data_size_sur; 138 139 PetscCall(DMGetDimension(dm, &dim)); 140 dim_sur = dim - height; 141 PetscCall(QDataBoundaryGetNumComponents(dm, &q_data_size_sur)); 142 PetscCall(DMGetCoordinateNumComps(dm, &num_comp_x)); 143 { 144 PetscSection section; 145 146 PetscCall(DMGetLocalSection(dm, §ion)); 147 PetscCall(PetscSectionGetFieldComponents(section, bc_def->dm_field, &num_comp_q)); 148 } 149 150 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, qf_func_ptr, qf_loc, qf_ifunc)); 151 PetscCallCeed(ceed, CeedQFunctionSetContext(*qf_ifunc, qfctx)); 152 PetscCallCeed(ceed, CeedQFunctionSetUserFlopsEstimate(*qf_ifunc, 0)); 153 PetscCallCeed(ceed, CeedQFunctionAddInput(*qf_ifunc, "q", num_comp_q, CEED_EVAL_INTERP)); 154 PetscCallCeed(ceed, CeedQFunctionAddInput(*qf_ifunc, "Grad_q", num_comp_q * dim_sur, CEED_EVAL_GRAD)); 155 PetscCallCeed(ceed, CeedQFunctionAddInput(*qf_ifunc, "surface qdata", q_data_size_sur, CEED_EVAL_NONE)); 156 PetscCallCeed(ceed, CeedQFunctionAddInput(*qf_ifunc, "x", num_comp_x, CEED_EVAL_INTERP)); 157 PetscCallCeed(ceed, CeedQFunctionAddOutput(*qf_ifunc, "v", num_comp_q, CEED_EVAL_INTERP)); 158 if (jac_data_size_sur) PetscCallCeed(ceed, CeedQFunctionAddOutput(*qf_ifunc, "surface jacobian data", jac_data_size_sur, CEED_EVAL_NONE)); 159 PetscFunctionReturn(PETSC_SUCCESS); 160 } 161 162 /** 163 @brief Create a QFunction matching the "standard" HONEE inputs for IJacobian 164 165 This assumes the `bc_def` context is a `HoneeBCStruct` and inputs/outputs of the IJacobian QFunction are in the correct order. 166 167 @param[in] bc_def `BCDefinition` that hosts the IJacobian 168 @param[in] qf_func_ptr `CeedQFunctionUser` for the IJacobian 169 @param[in] qf_loc Absolute path to source of `CeedQFunctionUser` 170 @param[in] qfctx `CeedQFunctionContext` for the IJacobian (also shared with the IFunction) 171 @param[out] qf_ijac QFunction for the IJacobian 172 **/ 173 PetscErrorCode HoneeBCCreateIJacobianQF(BCDefinition bc_def, CeedQFunctionUser qf_func_ptr, const char *qf_loc, CeedQFunctionContext qfctx, 174 CeedQFunction *qf_ijac) { 175 Ceed ceed; 176 DM dm; 177 PetscInt dim, dim_sur, height = 1, num_comp_x, num_comp_q; 178 CeedInt q_data_size_sur, jac_data_size_sur; 179 HoneeBCStruct honee_bc; 180 181 PetscFunctionBeginUser; 182 PetscCall(BCDefinitionGetDM(bc_def, &dm)); 183 PetscCall(BCDefinitionGetContext(bc_def, &honee_bc)); 184 ceed = honee_bc->honee->ceed; 185 jac_data_size_sur = honee_bc->jac_data_size_sur; 186 187 PetscCall(DMGetDimension(dm, &dim)); 188 dim_sur = dim - height; 189 PetscCall(QDataBoundaryGetNumComponents(dm, &q_data_size_sur)); 190 PetscCall(DMGetCoordinateNumComps(dm, &num_comp_x)); 191 { 192 PetscSection section; 193 194 PetscCall(DMGetLocalSection(dm, §ion)); 195 PetscCall(PetscSectionGetFieldComponents(section, bc_def->dm_field, &num_comp_q)); 196 } 197 198 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, qf_func_ptr, qf_loc, qf_ijac)); 199 PetscCallCeed(ceed, CeedQFunctionSetContext(*qf_ijac, qfctx)); 200 PetscCallCeed(ceed, CeedQFunctionSetUserFlopsEstimate(*qf_ijac, 0)); 201 PetscCallCeed(ceed, CeedQFunctionAddInput(*qf_ijac, "dq", num_comp_q, CEED_EVAL_INTERP)); 202 PetscCallCeed(ceed, CeedQFunctionAddInput(*qf_ijac, "Grad_dq", num_comp_q * dim_sur, CEED_EVAL_GRAD)); 203 PetscCallCeed(ceed, CeedQFunctionAddInput(*qf_ijac, "surface qdata", q_data_size_sur, CEED_EVAL_NONE)); 204 PetscCallCeed(ceed, CeedQFunctionAddInput(*qf_ijac, "x", num_comp_x, CEED_EVAL_INTERP)); 205 if (jac_data_size_sur) PetscCallCeed(ceed, CeedQFunctionAddInput(*qf_ijac, "surface jacobian data", jac_data_size_sur, CEED_EVAL_NONE)); 206 PetscCallCeed(ceed, CeedQFunctionAddOutput(*qf_ijac, "v", num_comp_q, CEED_EVAL_INTERP)); 207 PetscFunctionReturn(PETSC_SUCCESS); 208 } 209 210 /** 211 @brief Setups and adds IFunction operator to given composite operator 212 213 @param[in] bc_def `BCDefinition` for the boundary condition IFunction 214 @param[in] domain_label `DMLabel` for the face orientation label 215 @param[in] label_value Orientation value 216 @param[in] qf_ifunc QFunction for the IFunction 217 @param[in,out] op_ifunc Composite operator to be added to 218 @param[out] sub_op_ifunc Non-composite operator that is added in this function. To be passed to `HoneeBCAddIJacobianOp()` 219 **/ 220 PetscErrorCode HoneeBCAddIFunctionOp(BCDefinition bc_def, DMLabel domain_label, PetscInt label_value, CeedQFunction qf_ifunc, CeedOperator op_ifunc, 221 CeedOperator *sub_op_ifunc) { 222 Honee honee; 223 Ceed ceed; 224 DM dm, dm_coord; 225 HoneeBCStruct honee_bc; 226 PetscInt dim, height = 1, num_comp_x, num_comp_q; 227 CeedInt q_data_size, jac_data_size; 228 CeedVector q_data, jac_data; 229 CeedOperator sub_op_ifunc_; 230 CeedElemRestriction elem_restr_qd, elem_restr_q, elem_restr_x, elem_restr_jac; 231 CeedBasis basis_x, basis_q; 232 233 PetscFunctionBeginUser; 234 PetscCall(BCDefinitionGetDM(bc_def, &dm)); 235 PetscCall(BCDefinitionGetContext(bc_def, &honee_bc)); 236 honee = honee_bc->honee; 237 ceed = honee_bc->honee->ceed; 238 jac_data_size = honee_bc->jac_data_size_sur; 239 240 PetscCall(DMGetDimension(dm, &dim)); 241 PetscCall(QDataBoundaryGetNumComponents(dm, &q_data_size)); 242 PetscCall(DMGetCoordinateNumComps(dm, &num_comp_x)); 243 244 { 245 PetscSection section; 246 247 PetscCall(DMGetLocalSection(dm, §ion)); 248 PetscCall(PetscSectionGetFieldComponents(section, bc_def->dm_field, &num_comp_q)); 249 } 250 251 PetscCall(DMGetCoordinateDM(dm, &dm_coord)); 252 PetscCall(CreateBasisFromPlex(ceed, dm, domain_label, label_value, height, bc_def->dm_field, &basis_q)); 253 PetscCall(CreateBasisFromPlex(ceed, dm_coord, domain_label, label_value, height, bc_def->dm_field, &basis_x)); 254 PetscCall(DMPlexCeedElemRestrictionCreate(ceed, dm, domain_label, label_value, height, bc_def->dm_field, &elem_restr_q)); 255 PetscCall(DMPlexCeedElemRestrictionCoordinateCreate(ceed, dm, domain_label, label_value, height, &elem_restr_x)); 256 if (jac_data_size > 0) { 257 PetscCall(DMPlexCeedElemRestrictionQDataCreate(ceed, dm, domain_label, label_value, height, jac_data_size, &elem_restr_jac)); 258 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_jac, &jac_data, NULL)); 259 } 260 261 PetscCall(QDataBoundaryGet(ceed, dm, domain_label, label_value, elem_restr_x, basis_x, honee->x_coord, &elem_restr_qd, &q_data, &q_data_size)); 262 263 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_ifunc, NULL, NULL, &sub_op_ifunc_)); 264 PetscCallCeed(ceed, CeedOperatorSetField(sub_op_ifunc_, "q", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE)); 265 PetscCallCeed(ceed, CeedOperatorSetField(sub_op_ifunc_, "Grad_q", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE)); 266 PetscCallCeed(ceed, CeedOperatorSetField(sub_op_ifunc_, "surface qdata", elem_restr_qd, CEED_BASIS_NONE, q_data)); 267 PetscCallCeed(ceed, CeedOperatorSetField(sub_op_ifunc_, "x", elem_restr_x, basis_x, honee->x_coord)); 268 PetscCallCeed(ceed, CeedOperatorSetField(sub_op_ifunc_, "v", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE)); 269 if (jac_data_size > 0) PetscCallCeed(ceed, CeedOperatorSetField(sub_op_ifunc_, "surface jacobian data", elem_restr_jac, CEED_BASIS_NONE, jac_data)); 270 271 PetscCallCeed(ceed, CeedCompositeOperatorAddSub(op_ifunc, sub_op_ifunc_)); 272 *sub_op_ifunc = sub_op_ifunc_; 273 274 PetscCallCeed(ceed, CeedVectorDestroy(&q_data)); 275 PetscCallCeed(ceed, CeedBasisDestroy(&basis_q)); 276 PetscCallCeed(ceed, CeedBasisDestroy(&basis_x)); 277 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd)); 278 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_q)); 279 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_x)); 280 if (jac_data_size > 0) { 281 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_jac)); 282 PetscCallCeed(ceed, CeedVectorDestroy(&jac_data)); 283 } 284 PetscFunctionReturn(PETSC_SUCCESS); 285 } 286 287 /** 288 @brief Setups and adds IJacobian operator to given composite operator 289 290 The field data (element restriction, vector, etc) is read from `sub_op_ifunc` and used for the IJacobian. 291 292 @param[in] bc_def `BCDefinition` for the boundary condition IJacobian 293 @param[out] sub_op_ifunc IFunction operator corresponding to this IJacobian operator. 294 @param[in] domain_label `DMLabel` for the face orientation label 295 @param[in] label_value Orientation value 296 @param[in] qf_ijac QFunction for the IJacobian 297 @param[in,out] op_ijac Composite operator to be added to 298 **/ 299 PetscErrorCode HoneeBCAddIJacobianOp(BCDefinition bc_def, CeedOperator sub_op_ifunc, DMLabel domain_label, PetscInt label_value, 300 CeedQFunction qf_ijac, CeedOperator op_ijac) { 301 Honee honee; 302 Ceed ceed; 303 DM dm, dm_coord; 304 HoneeBCStruct honee_bc; 305 PetscInt dim, height = 1, num_comp_x, num_comp_q; 306 CeedInt q_data_size; 307 CeedVector q_data, jac_data; 308 CeedElemRestriction elem_restr_qd, elem_restr_q, elem_restr_x, elem_restr_jac; 309 CeedOperator sub_op_ijac; 310 CeedBasis basis_x, basis_q; 311 312 PetscFunctionBeginUser; 313 PetscCall(BCDefinitionGetContext(bc_def, &honee_bc)); 314 PetscBool use_jac_data = honee_bc->jac_data_size_sur > 0; 315 PetscCall(BCDefinitionGetDM(bc_def, &dm)); 316 honee = honee_bc->honee; 317 ceed = honee_bc->honee->ceed; 318 319 PetscCall(DMGetDimension(dm, &dim)); 320 PetscCall(QDataBoundaryGetNumComponents(dm, &q_data_size)); 321 PetscCall(DMGetCoordinateNumComps(dm, &num_comp_x)); 322 323 { 324 PetscSection section; 325 326 PetscCall(DMGetLocalSection(dm, §ion)); 327 PetscCall(PetscSectionGetFieldComponents(section, bc_def->dm_field, &num_comp_q)); 328 } 329 330 PetscCall(DMGetCoordinateDM(dm, &dm_coord)); 331 PetscCall(CreateBasisFromPlex(ceed, dm, domain_label, label_value, height, bc_def->dm_field, &basis_q)); 332 PetscCall(CreateBasisFromPlex(ceed, dm_coord, domain_label, label_value, height, bc_def->dm_field, &basis_x)); 333 334 { // Get restriction and basis from the RHS function 335 CeedOperatorField op_field; 336 337 PetscCallCeed(ceed, CeedOperatorGetFieldByName(sub_op_ifunc, "q", &op_field)); 338 PetscCallCeed(ceed, CeedOperatorFieldGetElemRestriction(op_field, &elem_restr_q)); 339 340 PetscCallCeed(ceed, CeedOperatorGetFieldByName(sub_op_ifunc, "x", &op_field)); 341 PetscCallCeed(ceed, CeedOperatorFieldGetElemRestriction(op_field, &elem_restr_x)); 342 343 if (use_jac_data) { 344 PetscCallCeed(ceed, CeedOperatorGetFieldByName(sub_op_ifunc, "surface jacobian data", &op_field)); 345 PetscCallCeed(ceed, CeedOperatorFieldGetData(op_field, NULL, &elem_restr_jac, NULL, &jac_data)); 346 } 347 } 348 349 PetscCall(QDataBoundaryGet(ceed, dm, domain_label, label_value, elem_restr_x, basis_x, honee->x_coord, &elem_restr_qd, &q_data, &q_data_size)); 350 351 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_ijac, NULL, NULL, &sub_op_ijac)); 352 PetscCallCeed(ceed, CeedOperatorSetField(sub_op_ijac, "dq", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE)); 353 PetscCallCeed(ceed, CeedOperatorSetField(sub_op_ijac, "Grad_dq", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE)); 354 PetscCallCeed(ceed, CeedOperatorSetField(sub_op_ijac, "surface qdata", elem_restr_qd, CEED_BASIS_NONE, q_data)); 355 PetscCallCeed(ceed, CeedOperatorSetField(sub_op_ijac, "x", elem_restr_x, basis_x, honee->x_coord)); 356 if (use_jac_data) PetscCallCeed(ceed, CeedOperatorSetField(sub_op_ijac, "surface jacobian data", elem_restr_jac, CEED_BASIS_NONE, jac_data)); 357 PetscCallCeed(ceed, CeedOperatorSetField(sub_op_ijac, "v", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE)); 358 359 PetscCallCeed(ceed, CeedCompositeOperatorAddSub(op_ijac, sub_op_ijac)); 360 361 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_q)); 362 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_x)); 363 if (use_jac_data) { 364 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_jac)); 365 PetscCallCeed(ceed, CeedVectorDestroy(&jac_data)); 366 } 367 PetscCallCeed(ceed, CeedVectorDestroy(&q_data)); 368 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd)); 369 PetscCallCeed(ceed, CeedBasisDestroy(&basis_q)); 370 PetscCallCeed(ceed, CeedBasisDestroy(&basis_x)); 371 PetscCallCeed(ceed, CeedOperatorDestroy(&sub_op_ijac)); 372 PetscFunctionReturn(PETSC_SUCCESS); 373 } 374