1 // SPDX-FileCopyrightText: Copyright (c) 2017-2024, HONEE contributors. 2 // SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause 3 4 /// @file 5 /// Utility functions for setting up ADVECTION 6 7 #include "../qfunctions/advection.h" 8 9 #include <ceed.h> 10 #include <petscdm.h> 11 12 #include <navierstokes.h> 13 14 const char *const AdvDifWindTypes[] = {"ROTATION", "TRANSLATION", "BOUNDARY_LAYER", "WindType", "ADVDIF_WIND_", NULL}; 15 const char *const AdvDifICTypes[] = {"SPHERE", "CYLINDER", "COSINE_HILL", "SKEW", "WAVE", "BOUNDARY_LAYER", "AdvectionICType", "ADVDIF_IC_", NULL}; 16 const char *const AdvDifWaveTypes[] = {"SINE", "SQUARE", "AdvDifWaveType", "ADVDIF_WAVE_", NULL}; 17 const char *const AdvDifBubbleContinuityTypes[] = {"SMOOTH", "BACK_SHARP", "THICK", "COSINE", "BubbleContinuityType", "ADVDIF_BUBBLE_CONTINUITY_", 18 NULL}; 19 const char *const StabilizationTauTypes[] = {"CTAU", "ADVDIFF_SHAKIB", "StabilizationTauType", "STAB_TAU_", NULL}; 20 21 // @brief Create CeedOperator for stabilized mass KSP for explicit timestepping 22 // 23 // Only used for SUPG stabilization 24 PetscErrorCode CreateKSPMassOperator_AdvectionStabilized(Honee honee, CeedOperator *op_mass) { 25 Ceed ceed = honee->ceed; 26 CeedInt num_comp_q, q_data_size; 27 CeedQFunction qf_mass = NULL; 28 CeedElemRestriction elem_restr_q, elem_restr_qd; 29 CeedBasis basis_q; 30 CeedVector q_data; 31 CeedQFunctionContext qfctx = NULL; 32 PetscInt dim; 33 34 PetscFunctionBeginUser; 35 PetscCall(DMGetDimension(honee->dm, &dim)); 36 { // Get restriction and basis from the RHS function 37 CeedOperator *sub_ops; 38 CeedOperatorField op_field; 39 PetscInt sub_op_index = 0; // will be 0 for the volume op 40 41 PetscCallCeed(ceed, CeedCompositeOperatorGetSubList(honee->op_rhs_ctx->op, &sub_ops)); 42 PetscCallCeed(ceed, CeedOperatorGetFieldByName(sub_ops[sub_op_index], "q", &op_field)); 43 PetscCallCeed(ceed, CeedOperatorFieldGetData(op_field, NULL, &elem_restr_q, &basis_q, NULL)); 44 PetscCallCeed(ceed, CeedOperatorGetFieldByName(sub_ops[sub_op_index], "qdata", &op_field)); 45 PetscCallCeed(ceed, CeedOperatorFieldGetData(op_field, NULL, &elem_restr_qd, NULL, &q_data)); 46 47 PetscCallCeed(ceed, CeedOperatorGetContext(sub_ops[sub_op_index], &qfctx)); 48 } 49 50 PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(elem_restr_q, &num_comp_q)); 51 PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(elem_restr_qd, &q_data_size)); 52 53 switch (dim) { 54 case 2: 55 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, MassFunction_Advection2D, MassFunction_Advection2D_loc, &qf_mass)); 56 break; 57 case 3: 58 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, MassFunction_Advection, MassFunction_Advection_loc, &qf_mass)); 59 break; 60 } 61 62 PetscCallCeed(ceed, CeedQFunctionSetContext(qf_mass, qfctx)); 63 PetscCallCeed(ceed, CeedQFunctionSetUserFlopsEstimate(qf_mass, 0)); 64 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_mass, "q_dot", 5, CEED_EVAL_INTERP)); 65 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_mass, "q", 5, CEED_EVAL_INTERP)); 66 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_mass, "qdata", q_data_size, CEED_EVAL_NONE)); 67 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_mass, "v", 5, CEED_EVAL_INTERP)); 68 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_mass, "Grad_v", 5 * dim, CEED_EVAL_GRAD)); 69 70 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_mass, NULL, NULL, op_mass)); 71 PetscCallCeed(ceed, CeedOperatorSetName(*op_mass, "RHS Mass Operator, Advection-Diffusion Stabilized")); 72 PetscCallCeed(ceed, CeedOperatorSetField(*op_mass, "q_dot", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE)); 73 PetscCallCeed(ceed, CeedOperatorSetField(*op_mass, "q", elem_restr_q, basis_q, honee->q_ceed)); 74 PetscCallCeed(ceed, CeedOperatorSetField(*op_mass, "qdata", elem_restr_qd, CEED_BASIS_NONE, q_data)); 75 PetscCallCeed(ceed, CeedOperatorSetField(*op_mass, "v", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE)); 76 PetscCallCeed(ceed, CeedOperatorSetField(*op_mass, "Grad_v", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE)); 77 78 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_q)); 79 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd)); 80 PetscCallCeed(ceed, CeedVectorDestroy(&q_data)); 81 PetscCallCeed(ceed, CeedBasisDestroy(&basis_q)); 82 PetscCallCeed(ceed, CeedQFunctionContextDestroy(&qfctx)); 83 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_mass)); 84 PetscFunctionReturn(PETSC_SUCCESS); 85 } 86 87 /** 88 @brief Create RHS CeedOperator for direct projection of divergence of diffusive flux 89 90 @param[in] honee `Honee` context 91 @param[in] diff_flux_proj `DivDiffFluxProjectionData` object 92 @param[out] op_rhs Operator to calculate the RHS of the L^2 projection 93 **/ 94 static PetscErrorCode DivDiffFluxProjectionCreateRHS_Direct_AdvDif(Honee honee, DivDiffFluxProjectionData diff_flux_proj, CeedOperator *op_rhs) { 95 Ceed ceed = honee->ceed; 96 NodalProjectionData projection = diff_flux_proj->projection; 97 CeedInt num_comp_q; 98 PetscInt dim, label_value = 0; 99 DMLabel domain_label = NULL; 100 CeedQFunctionContext advection_qfctx = NULL; 101 102 PetscFunctionBeginUser; 103 // -- Get Pre-requisite things 104 PetscCall(DMGetDimension(projection->dm, &dim)); 105 PetscCallCeed(ceed, CeedBasisGetNumComponents(honee->basis_q, &num_comp_q)); 106 107 { // Get advection-diffusion QF context 108 CeedOperator *sub_ops; 109 PetscInt sub_op_index = 0; // will be 0 for the volume op 110 111 if (honee->op_ifunction) PetscCallCeed(ceed, CeedCompositeOperatorGetSubList(honee->op_ifunction, &sub_ops)); 112 else PetscCallCeed(ceed, CeedCompositeOperatorGetSubList(honee->op_rhs_ctx->op, &sub_ops)); 113 PetscCallCeed(ceed, CeedOperatorGetContext(sub_ops[sub_op_index], &advection_qfctx)); 114 } 115 PetscCallCeed(ceed, CeedCompositeOperatorCreate(ceed, op_rhs)); 116 { // Add the volume integral CeedOperator 117 CeedQFunction qf_rhs_volume = NULL; 118 CeedOperator op_rhs_volume; 119 CeedVector q_data; 120 CeedElemRestriction elem_restr_qd, elem_restr_diff_flux_volume = NULL; 121 CeedBasis basis_diff_flux = NULL; 122 CeedInt q_data_size; 123 124 PetscCall(DivDiffFluxProjectionGetOperatorFieldData(diff_flux_proj, &elem_restr_diff_flux_volume, &basis_diff_flux, NULL, NULL)); 125 PetscCall(QDataGet(ceed, projection->dm, domain_label, label_value, honee->elem_restr_x, honee->basis_x, honee->x_coord, &elem_restr_qd, &q_data, 126 &q_data_size)); 127 switch (dim) { 128 case 2: 129 PetscCallCeed( 130 ceed, CeedQFunctionCreateInterior(ceed, 1, DivDiffusiveFluxVolumeRHS_AdvDif_2D, DivDiffusiveFluxVolumeRHS_AdvDif_2D_loc, &qf_rhs_volume)); 131 break; 132 case 3: 133 PetscCallCeed( 134 ceed, CeedQFunctionCreateInterior(ceed, 1, DivDiffusiveFluxVolumeRHS_AdvDif_3D, DivDiffusiveFluxVolumeRHS_AdvDif_3D_loc, &qf_rhs_volume)); 135 break; 136 } 137 PetscCheck(qf_rhs_volume, honee->comm, PETSC_ERR_SUP, "%s not valid for DM of dimension %" PetscInt_FMT, __func__, dim); 138 139 PetscCallCeed(ceed, CeedQFunctionSetContext(qf_rhs_volume, advection_qfctx)); 140 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_rhs_volume, "Grad_q", num_comp_q * dim, CEED_EVAL_GRAD)); 141 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_rhs_volume, "qdata", q_data_size, CEED_EVAL_NONE)); 142 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_rhs_volume, "diffusive flux RHS", projection->num_comp * dim, CEED_EVAL_GRAD)); 143 144 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_rhs_volume, NULL, NULL, &op_rhs_volume)); 145 PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_volume, "Grad_q", honee->elem_restr_q, honee->basis_q, CEED_VECTOR_ACTIVE)); 146 PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_volume, "qdata", elem_restr_qd, CEED_BASIS_NONE, q_data)); 147 PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_volume, "diffusive flux RHS", elem_restr_diff_flux_volume, basis_diff_flux, CEED_VECTOR_ACTIVE)); 148 149 PetscCallCeed(ceed, CeedCompositeOperatorAddSub(*op_rhs, op_rhs_volume)); 150 151 PetscCallCeed(ceed, CeedVectorDestroy(&q_data)); 152 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd)); 153 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_diff_flux_volume)); 154 PetscCallCeed(ceed, CeedBasisDestroy(&basis_diff_flux)); 155 PetscCallCeed(ceed, CeedOperatorDestroy(&op_rhs_volume)); 156 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_rhs_volume)); 157 } 158 159 { // Add the boundary integral CeedOperator 160 CeedQFunction qf_rhs_boundary; 161 DMLabel face_sets_label; 162 PetscInt num_face_set_values, *face_set_values; 163 CeedInt q_data_size; 164 165 // -- Build RHS operator 166 switch (dim) { 167 case 2: 168 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, DivDiffusiveFluxBoundaryRHS_AdvDif_2D, DivDiffusiveFluxBoundaryRHS_AdvDif_2D_loc, 169 &qf_rhs_boundary)); 170 break; 171 case 3: 172 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, DivDiffusiveFluxBoundaryRHS_AdvDif_3D, DivDiffusiveFluxBoundaryRHS_AdvDif_3D_loc, 173 &qf_rhs_boundary)); 174 break; 175 } 176 177 PetscCall(QDataBoundaryGradientGetNumComponents(honee->dm, &q_data_size)); 178 PetscCallCeed(ceed, CeedQFunctionSetContext(qf_rhs_boundary, advection_qfctx)); 179 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_rhs_boundary, "Grad_q", num_comp_q * dim, CEED_EVAL_GRAD)); 180 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_rhs_boundary, "qdata", q_data_size, CEED_EVAL_NONE)); 181 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_rhs_boundary, "diffusive flux RHS", projection->num_comp, CEED_EVAL_INTERP)); 182 183 PetscCall(DMGetLabel(projection->dm, "Face Sets", &face_sets_label)); 184 PetscCall(DMLabelCreateGlobalValueArray(projection->dm, face_sets_label, &num_face_set_values, &face_set_values)); 185 for (PetscInt f = 0; f < num_face_set_values; f++) { 186 DMLabel face_orientation_label; 187 PetscInt num_orientations_values, *orientation_values; 188 189 { 190 char *face_orientation_label_name; 191 192 PetscCall(DMPlexCreateFaceLabel(projection->dm, face_set_values[f], &face_orientation_label_name)); 193 PetscCall(DMGetLabel(projection->dm, face_orientation_label_name, &face_orientation_label)); 194 PetscCall(PetscFree(face_orientation_label_name)); 195 } 196 PetscCall(DMLabelCreateGlobalValueArray(projection->dm, face_orientation_label, &num_orientations_values, &orientation_values)); 197 for (PetscInt o = 0; o < num_orientations_values; o++) { 198 CeedOperator op_rhs_boundary; 199 CeedBasis basis_q, basis_diff_flux_boundary; 200 CeedElemRestriction elem_restr_qdata, elem_restr_q, elem_restr_diff_flux_boundary; 201 CeedVector q_data; 202 CeedInt q_data_size; 203 PetscInt orientation = orientation_values[o], dm_field_q = 0, height_cell = 0, height_face = 1; 204 205 PetscCall(DMPlexCeedElemRestrictionCreate(ceed, honee->dm, face_orientation_label, orientation, height_cell, dm_field_q, &elem_restr_q)); 206 PetscCall(DMPlexCeedBasisCellToFaceCreate(ceed, honee->dm, face_orientation_label, orientation, orientation, dm_field_q, &basis_q)); 207 PetscCall(DMPlexCeedElemRestrictionCreate(ceed, projection->dm, face_orientation_label, orientation, height_face, 0, 208 &elem_restr_diff_flux_boundary)); 209 PetscCall(CreateBasisFromPlex(ceed, projection->dm, face_orientation_label, orientation, height_face, 0, &basis_diff_flux_boundary)); 210 PetscCall( 211 QDataBoundaryGradientGet(ceed, honee->dm, face_orientation_label, orientation, honee->x_coord, &elem_restr_qdata, &q_data, &q_data_size)); 212 213 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_rhs_boundary, NULL, NULL, &op_rhs_boundary)); 214 PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_boundary, "Grad_q", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE)); 215 PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_boundary, "qdata", elem_restr_qdata, CEED_BASIS_NONE, q_data)); 216 PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_boundary, "diffusive flux RHS", elem_restr_diff_flux_boundary, basis_diff_flux_boundary, 217 CEED_VECTOR_ACTIVE)); 218 219 PetscCallCeed(ceed, CeedCompositeOperatorAddSub(*op_rhs, op_rhs_boundary)); 220 221 PetscCallCeed(ceed, CeedOperatorDestroy(&op_rhs_boundary)); 222 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qdata)); 223 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_q)); 224 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_diff_flux_boundary)); 225 PetscCallCeed(ceed, CeedBasisDestroy(&basis_q)); 226 PetscCallCeed(ceed, CeedBasisDestroy(&basis_diff_flux_boundary)); 227 PetscCallCeed(ceed, CeedVectorDestroy(&q_data)); 228 } 229 PetscCall(PetscFree(orientation_values)); 230 } 231 PetscCall(PetscFree(face_set_values)); 232 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_rhs_boundary)); 233 } 234 235 PetscCallCeed(ceed, CeedQFunctionContextDestroy(&advection_qfctx)); 236 PetscFunctionReturn(PETSC_SUCCESS); 237 } 238 239 /** 240 @brief Create RHS CeedOperator for indirect projection of divergence of diffusive flux 241 242 @param[in] honee `Honee` context 243 @param[in] diff_flux_proj `DivDiffFluxProjectionData` object 244 @param[out] op_rhs Operator to calculate the RHS of the L^2 projection 245 **/ 246 static PetscErrorCode DivDiffFluxProjectionCreateRHS_Indirect_AdvDif(Honee honee, DivDiffFluxProjectionData diff_flux_proj, CeedOperator *op_rhs) { 247 Ceed ceed = honee->ceed; 248 NodalProjectionData projection = diff_flux_proj->projection; 249 CeedBasis basis_diff_flux; 250 CeedElemRestriction elem_restr_diff_flux, elem_restr_qd; 251 CeedVector q_data; 252 CeedInt num_comp_q, q_data_size; 253 PetscInt dim; 254 PetscInt label_value = 0, height = 0, dm_field = 0; 255 DMLabel domain_label = NULL; 256 CeedQFunction qf_rhs = NULL; 257 CeedQFunctionContext advection_qfctx = NULL; 258 259 PetscFunctionBeginUser; 260 PetscCall(DMGetDimension(projection->dm, &dim)); 261 PetscCallCeed(ceed, CeedBasisGetNumComponents(honee->basis_q, &num_comp_q)); 262 263 { // Get advection-diffusion QF context 264 CeedOperator *sub_ops; 265 PetscInt sub_op_index = 0; // will be 0 for the volume op 266 267 if (honee->op_ifunction) PetscCallCeed(ceed, CeedCompositeOperatorGetSubList(honee->op_ifunction, &sub_ops)); 268 else PetscCallCeed(ceed, CeedCompositeOperatorGetSubList(honee->op_rhs_ctx->op, &sub_ops)); 269 PetscCallCeed(ceed, CeedOperatorGetContext(sub_ops[sub_op_index], &advection_qfctx)); 270 } 271 PetscCall(DMPlexCeedElemRestrictionCreate(ceed, projection->dm, domain_label, label_value, height, dm_field, &elem_restr_diff_flux)); 272 PetscCall(CreateBasisFromPlex(ceed, projection->dm, domain_label, label_value, height, dm_field, &basis_diff_flux)); 273 PetscCall(QDataGet(ceed, projection->dm, domain_label, label_value, honee->elem_restr_x, honee->basis_x, honee->x_coord, &elem_restr_qd, &q_data, 274 &q_data_size)); 275 276 switch (dim) { 277 case 2: 278 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, DiffusiveFluxRHS_AdvDif_2D, DiffusiveFluxRHS_AdvDif_2D_loc, &qf_rhs)); 279 break; 280 case 3: 281 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, DiffusiveFluxRHS_AdvDif_3D, DiffusiveFluxRHS_AdvDif_3D_loc, &qf_rhs)); 282 break; 283 } 284 PetscCheck(qf_rhs, honee->comm, PETSC_ERR_SUP, "%s not valid for DM of dimension %" PetscInt_FMT, __func__, dim); 285 286 PetscCallCeed(ceed, CeedQFunctionSetContext(qf_rhs, advection_qfctx)); 287 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_rhs, "Grad_q", num_comp_q * dim, CEED_EVAL_GRAD)); 288 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_rhs, "qdata", q_data_size, CEED_EVAL_NONE)); 289 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_rhs, "F_diff RHS", projection->num_comp, CEED_EVAL_INTERP)); 290 291 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_rhs, NULL, NULL, op_rhs)); 292 PetscCallCeed(ceed, CeedOperatorSetField(*op_rhs, "Grad_q", honee->elem_restr_q, honee->basis_q, CEED_VECTOR_ACTIVE)); 293 PetscCallCeed(ceed, CeedOperatorSetField(*op_rhs, "qdata", elem_restr_qd, CEED_BASIS_NONE, q_data)); 294 PetscCallCeed(ceed, CeedOperatorSetField(*op_rhs, "F_diff RHS", elem_restr_diff_flux, basis_diff_flux, CEED_VECTOR_ACTIVE)); 295 296 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_rhs)); 297 PetscCallCeed(ceed, CeedQFunctionContextDestroy(&advection_qfctx)); 298 PetscCallCeed(ceed, CeedBasisDestroy(&basis_diff_flux)); 299 PetscCallCeed(ceed, CeedVectorDestroy(&q_data)); 300 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd)); 301 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_diff_flux)); 302 PetscFunctionReturn(PETSC_SUCCESS); 303 } 304 305 static PetscErrorCode AdvectionInflowBCSetup_CreateIFunctionQF(BCDefinition bc_def, CeedQFunction *qf) { 306 HoneeBCStruct honee_bc; 307 DM dm; 308 PetscInt dim; 309 310 PetscFunctionBeginUser; 311 PetscCall(BCDefinitionGetContext(bc_def, &honee_bc)); 312 PetscCall(BCDefinitionGetDM(bc_def, &dm)); 313 PetscCall(DMGetDimension(dm, &dim)); 314 switch (dim) { 315 case 2: 316 PetscCall(HoneeBCCreateIFunctionQF(bc_def, Advection2d_InOutFlow, Advection2d_InOutFlow_loc, honee_bc->qfctx, qf)); 317 break; 318 case 3: 319 PetscCall(HoneeBCCreateIFunctionQF(bc_def, Advection_InOutFlow, Advection_InOutFlow_loc, honee_bc->qfctx, qf)); 320 break; 321 } 322 PetscFunctionReturn(PETSC_SUCCESS); 323 } 324 325 PetscErrorCode NS_ADVECTION(ProblemData problem, DM dm, void *ctx) { 326 AdvDifWindType wind_type; 327 AdvDifICType advectionic_type; 328 AdvDifBubbleContinuityType bubble_continuity_type = -1; 329 StabilizationType stab; 330 StabilizationTauType stab_tau; 331 SetupContextAdv setup_context; 332 Honee honee = *(Honee *)ctx; 333 MPI_Comm comm = honee->comm; 334 Ceed ceed = honee->ceed; 335 PetscBool implicit; 336 AdvectionContext advection_ctx; 337 CeedQFunctionContext advection_qfctx; 338 PetscInt dim; 339 340 PetscFunctionBeginUser; 341 PetscCall(PetscCalloc1(1, &setup_context)); 342 PetscCall(PetscCalloc1(1, &advection_ctx)); 343 PetscCall(DMGetDimension(dm, &dim)); 344 345 // ------------------------------------------------------ 346 // SET UP ADVECTION 347 // ------------------------------------------------------ 348 problem->print_info = PRINT_ADVECTION; 349 problem->num_comps_jac_data = 0; 350 switch (dim) { 351 case 2: 352 problem->ics.qf_func_ptr = ICsAdvection2d; 353 problem->ics.qf_loc = ICsAdvection2d_loc; 354 problem->apply_vol_rhs.qf_func_ptr = RHS_Advection2d; 355 problem->apply_vol_rhs.qf_loc = RHS_Advection2d_loc; 356 problem->apply_vol_ifunction.qf_func_ptr = IFunction_Advection2d; 357 problem->apply_vol_ifunction.qf_loc = IFunction_Advection2d_loc; 358 problem->compute_exact_solution_error = PETSC_TRUE; 359 break; 360 case 3: 361 problem->ics.qf_func_ptr = ICsAdvection; 362 problem->ics.qf_loc = ICsAdvection_loc; 363 problem->apply_vol_rhs.qf_func_ptr = RHS_Advection; 364 problem->apply_vol_rhs.qf_loc = RHS_Advection_loc; 365 problem->apply_vol_ifunction.qf_func_ptr = IFunction_Advection; 366 problem->apply_vol_ifunction.qf_loc = IFunction_Advection_loc; 367 problem->compute_exact_solution_error = PETSC_FALSE; 368 break; 369 } 370 371 PetscCall(DivDiffFluxProjectionCreate(honee, 1, &honee->diff_flux_proj)); 372 if (honee->diff_flux_proj) { 373 DivDiffFluxProjectionData diff_flux_proj = honee->diff_flux_proj; 374 NodalProjectionData projection = diff_flux_proj->projection; 375 376 diff_flux_proj->CreateRHSOperator_Direct = DivDiffFluxProjectionCreateRHS_Direct_AdvDif; 377 diff_flux_proj->CreateRHSOperator_Indirect = DivDiffFluxProjectionCreateRHS_Indirect_AdvDif; 378 379 switch (honee->diff_flux_proj->method) { 380 case DIV_DIFF_FLUX_PROJ_DIRECT: { 381 PetscSection section; 382 383 PetscCall(DMGetLocalSection(projection->dm, §ion)); 384 PetscCall(PetscSectionSetFieldName(section, 0, "")); 385 PetscCall(PetscSectionSetComponentName(section, 0, 0, "DivDiffusiveFlux_Scalar")); 386 } break; 387 case DIV_DIFF_FLUX_PROJ_INDIRECT: { 388 PetscSection section; 389 390 PetscCall(DMGetLocalSection(projection->dm, §ion)); 391 PetscCall(PetscSectionSetFieldName(section, 0, "")); 392 PetscCall(PetscSectionSetComponentName(section, 0, 0, "DiffusiveFlux_ScalarX")); 393 PetscCall(PetscSectionSetComponentName(section, 0, 1, "DiffusiveFlux_ScalarY")); 394 if (dim >= 3) PetscCall(PetscSectionSetComponentName(section, 0, 2, "DiffusiveFlux_ScalarZ")); 395 } break; 396 case DIV_DIFF_FLUX_PROJ_NONE: 397 SETERRQ(PetscObjectComm((PetscObject)honee->dm), PETSC_ERR_ARG_WRONG, "Should not reach here with div_diff_flux_projection_method %s", 398 DivDiffFluxProjectionMethods[honee->app_ctx->divFdiffproj_method]); 399 break; 400 } 401 } 402 403 // ------------------------------------------------------ 404 // Create the QFunction context 405 // ------------------------------------------------------ 406 CeedScalar rc = 1000.; // m (Radius of bubble) 407 CeedScalar CtauS = 0.; // dimensionless 408 PetscBool strong_form = PETSC_FALSE; 409 CeedScalar E_wind = 1.e6; // J 410 CeedScalar Ctau_a = PetscPowScalarInt(honee->app_ctx->degree, 2); 411 CeedScalar Ctau_d = PetscPowScalarInt(honee->app_ctx->degree, 4); 412 CeedScalar Ctau_t = 0.; 413 PetscReal wind[3] = {1., 0, 0}; // m/s 414 CeedScalar diffusion_coeff = 0.; 415 CeedScalar wave_frequency = 2 * M_PI; 416 CeedScalar wave_phase = 0; 417 AdvDifWaveType wave_type = -1; 418 PetscScalar bl_height_factor = 1.; 419 PetscReal domain_min[3], domain_max[3], domain_size[3] = {0.}; 420 PetscCall(DMGetBoundingBox(dm, domain_min, domain_max)); 421 for (PetscInt i = 0; i < dim; i++) domain_size[i] = domain_max[i] - domain_min[i]; 422 423 // ------------------------------------------------------ 424 // Create the PETSc context 425 // ------------------------------------------------------ 426 PetscScalar meter = 1e-2; // 1 meter in scaled length units 427 PetscScalar kilogram = 1e-6; // 1 kilogram in scaled mass units 428 PetscScalar second = 1e-2; // 1 second in scaled time units 429 PetscScalar Joule; 430 431 // ------------------------------------------------------ 432 // Command line Options 433 // ------------------------------------------------------ 434 PetscOptionsBegin(comm, NULL, "Options for ADVECTION problem", NULL); 435 // -- Physics 436 PetscBool translation; 437 PetscCall(PetscOptionsEnum("-wind_type", "Wind type in Advection", NULL, AdvDifWindTypes, (PetscEnum)(wind_type = ADVDIF_WIND_ROTATION), 438 (PetscEnum *)&wind_type, &translation)); 439 PetscInt n = dim; 440 PetscBool user_wind; 441 PetscCall(PetscOptionsRealArray("-wind_translation", "Constant wind vector", NULL, wind, &n, &user_wind)); 442 PetscCall(PetscOptionsScalar("-diffusion_coeff", "Diffusion coefficient", NULL, diffusion_coeff, &diffusion_coeff, NULL)); 443 PetscCall(PetscOptionsScalar("-CtauS", "Scale coefficient for tau (nondimensional)", NULL, CtauS, &CtauS, NULL)); 444 PetscCall(PetscOptionsBool("-strong_form", "Strong (true) or weak/integrated by parts (false) advection residual", NULL, strong_form, &strong_form, 445 NULL)); 446 PetscCall(PetscOptionsScalar("-E_wind", "Total energy of inflow wind", NULL, E_wind, &E_wind, NULL)); 447 PetscCall(PetscOptionsEnum("-stab", "Stabilization method", NULL, StabilizationTypes, (PetscEnum)(stab = STAB_NONE), (PetscEnum *)&stab, NULL)); 448 PetscCall(PetscOptionsEnum("-stab_tau", "Stabilization constant, tau", NULL, StabilizationTauTypes, (PetscEnum)(stab_tau = STAB_TAU_CTAU), 449 (PetscEnum *)&stab_tau, NULL)); 450 PetscCall(PetscOptionsScalar("-Ctau_t", "Stabilization time constant", NULL, Ctau_t, &Ctau_t, NULL)); 451 PetscCall(PetscOptionsScalar("-Ctau_a", "Coefficient for the stabilization, advection component", NULL, Ctau_a, &Ctau_a, NULL)); 452 PetscCall(PetscOptionsScalar("-Ctau_d", "Coefficient for the stabilization, diffusion component", NULL, Ctau_d, &Ctau_d, NULL)); 453 PetscCall(PetscOptionsBool("-implicit", "Use implicit (IFunction) formulation", NULL, implicit = PETSC_FALSE, &implicit, NULL)); 454 PetscCall(PetscOptionsEnum("-advection_ic_type", "Initial condition for Advection problem", NULL, AdvDifICTypes, 455 (PetscEnum)(advectionic_type = ADVDIF_IC_BUBBLE_SPHERE), (PetscEnum *)&advectionic_type, NULL)); 456 // IC-specific options 457 switch (advectionic_type) { 458 case ADVDIF_IC_WAVE: 459 PetscCall(PetscOptionsDeprecated("-wave_type", "-advection_ic_wave_type", "HONEE 0.0", NULL)); 460 PetscCall(PetscOptionsDeprecated("-wave_frequency", "-advection_ic_wave_frequency", "HONEE 0.0", NULL)); 461 PetscCall(PetscOptionsDeprecated("-wave_phase", "-advection_ic_wave_phase", "HONEE 0.0", NULL)); 462 PetscCall(PetscOptionsEnum("-advection_ic_wave_type", "Type of wave", NULL, AdvDifWaveTypes, (PetscEnum)(wave_type = ADVDIF_WAVE_SINE), 463 (PetscEnum *)&wave_type, NULL)); 464 PetscCall(PetscOptionsScalar("-advection_ic_wave_frequency", "Frequency of sine wave", NULL, wave_frequency, &wave_frequency, NULL)); 465 PetscCall(PetscOptionsScalar("-advection_ic_wave_phase", "Length correction", NULL, wave_phase, &wave_phase, NULL)); 466 break; 467 case ADVDIF_IC_BOUNDARY_LAYER: 468 PetscCall( 469 PetscOptionsScalar("-advection_ic_bl_height_factor", "Height of boundary layer in IC", NULL, bl_height_factor, &bl_height_factor, NULL)); 470 break; 471 case ADVDIF_IC_BUBBLE_CYLINDER: 472 case ADVDIF_IC_BUBBLE_SPHERE: 473 PetscCall(PetscOptionsDeprecated("-rc", "-advection_ic_bubble_rc", "HONEE 0.0", NULL)); 474 PetscCall(PetscOptionsDeprecated("-bubble_continuity", "-advection_ic_bubble_continuity", "HONEE 0.0", NULL)); 475 PetscCall(PetscOptionsScalar("-advection_ic_bubble_rc", "Characteristic radius of thermal bubble", NULL, rc, &rc, NULL)); 476 bubble_continuity_type = dim == 3 ? ADVDIF_BUBBLE_CONTINUITY_SMOOTH : ADVDIF_BUBBLE_CONTINUITY_COSINE; 477 PetscCall(PetscOptionsEnum("-advection_ic_bubble_continuity", "Smooth, back_sharp, or thick", NULL, AdvDifBubbleContinuityTypes, 478 (PetscEnum)bubble_continuity_type, (PetscEnum *)&bubble_continuity_type, NULL)); 479 break; 480 case ADVDIF_IC_SKEW: 481 case ADVDIF_IC_COSINE_HILL: 482 break; 483 } 484 485 // -- Units 486 PetscCall(PetscOptionsScalar("-units_meter", "1 meter in scaled length units", NULL, meter, &meter, NULL)); 487 meter = fabs(meter); 488 PetscCall(PetscOptionsScalar("-units_kilogram", "1 kilogram in scaled mass units", NULL, kilogram, &kilogram, NULL)); 489 kilogram = fabs(kilogram); 490 PetscCall(PetscOptionsScalar("-units_second", "1 second in scaled time units", NULL, second, &second, NULL)); 491 second = fabs(second); 492 493 // -- Warnings 494 if (wind_type == ADVDIF_WIND_ROTATION && user_wind) { 495 PetscCall(PetscPrintf(comm, "Warning! Use -wind_translation only with -wind_type translation\n")); 496 } 497 if (wind_type == ADVDIF_WIND_TRANSLATION && advectionic_type == ADVDIF_IC_BUBBLE_CYLINDER && wind[2] != 0.) { 498 wind[2] = 0; 499 PetscCall( 500 PetscPrintf(comm, "Warning! Background wind in the z direction should be zero (-wind_translation x,x,0) with -advection_ic_type cylinder\n")); 501 } 502 if (stab == STAB_NONE && CtauS != 0) { 503 PetscCall(PetscPrintf(comm, "Warning! Use -CtauS only with -stab su or -stab supg\n")); 504 } 505 PetscOptionsEnd(); 506 507 if (stab == STAB_SUPG) problem->create_mass_operator = CreateKSPMassOperator_AdvectionStabilized; 508 509 // ------------------------------------------------------ 510 // Set up the PETSc context 511 // ------------------------------------------------------ 512 // -- Define derived units 513 Joule = kilogram * PetscSqr(meter) / PetscSqr(second); 514 515 honee->units->meter = meter; 516 honee->units->kilogram = kilogram; 517 honee->units->second = second; 518 honee->units->Joule = Joule; 519 520 // ------------------------------------------------------ 521 // Set up the QFunction contexts 522 // ------------------------------------------------------ 523 // -- Scale variables to desired units 524 E_wind *= Joule; 525 rc = fabs(rc) * meter; 526 for (PetscInt i = 0; i < dim; i++) { 527 wind[i] *= (meter / second); 528 domain_size[i] *= meter; 529 } 530 531 // -- Setup Context 532 setup_context->rc = rc; 533 setup_context->lx = domain_size[0]; 534 setup_context->ly = domain_size[1]; 535 setup_context->lz = dim == 3 ? domain_size[2] : 0.; 536 setup_context->wind[0] = wind[0]; 537 setup_context->wind[1] = wind[1]; 538 setup_context->wind[2] = dim == 3 ? wind[2] : 0.; 539 setup_context->wind_type = wind_type; 540 setup_context->initial_condition_type = advectionic_type; 541 setup_context->bubble_continuity_type = bubble_continuity_type; 542 setup_context->time = 0; 543 setup_context->wave_frequency = wave_frequency; 544 setup_context->wave_phase = wave_phase; 545 setup_context->wave_type = wave_type; 546 setup_context->bl_height_factor = bl_height_factor; 547 548 // -- QFunction Context 549 honee->phys->implicit = implicit; 550 advection_ctx->CtauS = CtauS; 551 advection_ctx->E_wind = E_wind; 552 advection_ctx->implicit = implicit; 553 advection_ctx->strong_form = strong_form; 554 advection_ctx->stabilization = stab; 555 advection_ctx->stabilization_tau = stab_tau; 556 advection_ctx->Ctau_a = Ctau_a; 557 advection_ctx->Ctau_d = Ctau_d; 558 advection_ctx->Ctau_t = Ctau_t; 559 advection_ctx->diffusion_coeff = diffusion_coeff; 560 advection_ctx->divFdiff_method = honee->app_ctx->divFdiffproj_method; 561 562 PetscCallCeed(ceed, CeedQFunctionContextCreate(honee->ceed, &problem->ics.qfctx)); 563 PetscCallCeed(ceed, CeedQFunctionContextSetData(problem->ics.qfctx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*setup_context), setup_context)); 564 PetscCallCeed(ceed, CeedQFunctionContextSetDataDestroy(problem->ics.qfctx, CEED_MEM_HOST, FreeContextPetsc)); 565 566 PetscCallCeed(ceed, CeedQFunctionContextCreate(honee->ceed, &advection_qfctx)); 567 PetscCallCeed(ceed, CeedQFunctionContextSetData(advection_qfctx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*advection_ctx), advection_ctx)); 568 PetscCallCeed(ceed, CeedQFunctionContextSetDataDestroy(advection_qfctx, CEED_MEM_HOST, FreeContextPetsc)); 569 PetscCallCeed(ceed, CeedQFunctionContextRegisterDouble(advection_qfctx, "timestep size", offsetof(struct AdvectionContext_, dt), 1, 570 "Size of timestep, delta t")); 571 problem->apply_vol_rhs.qfctx = advection_qfctx; 572 PetscCallCeed(ceed, CeedQFunctionContextReferenceCopy(advection_qfctx, &problem->apply_vol_ifunction.qfctx)); 573 574 for (PetscCount b = 0; b < problem->num_bc_defs; b++) { 575 BCDefinition bc_def = problem->bc_defs[b]; 576 const char *name; 577 578 PetscCall(BCDefinitionGetInfo(bc_def, &name, NULL, NULL)); 579 if (!strcmp(name, "inflow")) { 580 HoneeBCStruct honee_bc; 581 582 PetscCall(PetscNew(&honee_bc)); 583 PetscCallCeed(ceed, CeedQFunctionContextReferenceCopy(advection_qfctx, &honee_bc->qfctx)); 584 honee_bc->honee = honee; 585 honee_bc->num_comps_jac_data = 0; 586 PetscCall(BCDefinitionSetContext(bc_def, HoneeBCDestroy, honee_bc)); 587 588 PetscCall(BCDefinitionSetIFunction(bc_def, AdvectionInflowBCSetup_CreateIFunctionQF, HoneeBCAddIFunctionOp)); 589 PetscCall(BCDefinitionSetIJacobian(bc_def, NULL, NULL)); 590 } 591 } 592 PetscFunctionReturn(PETSC_SUCCESS); 593 } 594 595 PetscErrorCode PRINT_ADVECTION(Honee honee, ProblemData problem, AppCtx app_ctx) { 596 MPI_Comm comm = honee->comm; 597 Ceed ceed = honee->ceed; 598 SetupContextAdv setup_ctx; 599 AdvectionContext advection_ctx; 600 PetscInt dim; 601 602 PetscFunctionBeginUser; 603 PetscCall(DMGetDimension(honee->dm, &dim)); 604 PetscCallCeed(ceed, CeedQFunctionContextGetData(problem->ics.qfctx, CEED_MEM_HOST, &setup_ctx)); 605 PetscCallCeed(ceed, CeedQFunctionContextGetData(problem->apply_vol_rhs.qfctx, CEED_MEM_HOST, &advection_ctx)); 606 PetscCall(PetscPrintf(comm, 607 " Problem:\n" 608 " Problem Name : %s\n" 609 " Stabilization : %s\n" 610 " Stabilization Tau : %s\n" 611 " Wind Type : %s\n", 612 app_ctx->problem_name, StabilizationTypes[advection_ctx->stabilization], 613 StabilizationTauTypes[advection_ctx->stabilization_tau], AdvDifWindTypes[setup_ctx->wind_type])); 614 615 if (setup_ctx->wind_type == ADVDIF_WIND_TRANSLATION) { 616 CeedScalar *wind = setup_ctx->wind; 617 switch (dim) { 618 case 2: 619 PetscCall(PetscPrintf(comm, " Background Wind : %f,%f\n", wind[0], wind[1])); 620 break; 621 case 3: 622 PetscCall(PetscPrintf(comm, " Background Wind : %f,%f,%f\n", wind[0], wind[1], wind[2])); 623 break; 624 } 625 } 626 627 PetscCall(PetscPrintf(comm, " Initial Condition Type : %s\n", AdvDifICTypes[setup_ctx->initial_condition_type])); 628 switch (setup_ctx->initial_condition_type) { 629 case ADVDIF_IC_SKEW: 630 case ADVDIF_IC_COSINE_HILL: 631 case ADVDIF_IC_BOUNDARY_LAYER: 632 break; 633 case ADVDIF_IC_BUBBLE_SPHERE: 634 case ADVDIF_IC_BUBBLE_CYLINDER: 635 PetscCall(PetscPrintf(comm, " Bubble Continuity : %s\n", AdvDifBubbleContinuityTypes[setup_ctx->bubble_continuity_type])); 636 break; 637 case ADVDIF_IC_WAVE: 638 PetscCall(PetscPrintf(comm, " Wave Type : %s\n", AdvDifWaveTypes[setup_ctx->wave_type])); 639 break; 640 } 641 642 PetscCallCeed(ceed, CeedQFunctionContextRestoreData(problem->ics.qfctx, &setup_ctx)); 643 PetscCallCeed(ceed, CeedQFunctionContextRestoreData(problem->apply_vol_rhs.qfctx, &advection_ctx)); 644 PetscFunctionReturn(PETSC_SUCCESS); 645 } 646