1 // SPDX-FileCopyrightText: Copyright (c) 2017-2024, HONEE contributors. 2 // SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause 3 /// @file 4 /// Functions for setting up and projecting the divergence of the diffusive flux 5 6 #include "../qfunctions/diff_flux_projection.h" 7 8 #include <petscdmplex.h> 9 10 #include <navierstokes.h> 11 12 /** 13 @brief Create `DivDiffFluxProjectionData` for solution DM in `user` 14 15 @param[in] user `User` context 16 @param[in] num_diff_flux_comps Number of components that makes up the diffusive flux (e.g. 1 for scalar advection-diffusion) 17 @param[out] pdiff_flux_proj The `DivDiffFluxProjectionData` object created, or `NULL` if not created 18 **/ 19 PetscErrorCode DivDiffFluxProjectionCreate(User user, PetscInt num_diff_flux_comps, DivDiffFluxProjectionData *pdiff_flux_proj) { 20 PetscInt label_value = 0, height = 0, dm_field = 0, dim, degree = user->app_ctx->degree, q_extra = user->app_ctx->q_extra; 21 DMLabel domain_label = NULL; 22 DivDiffFluxProjectionData diff_flux_proj; 23 NodalProjectionData projection; 24 25 PetscFunctionBeginUser; 26 if (user->app_ctx->divFdiffproj_method == DIV_DIFF_FLUX_PROJ_NONE) { 27 *pdiff_flux_proj = NULL; 28 PetscFunctionReturn(PETSC_SUCCESS); 29 } 30 PetscCall(PetscNew(pdiff_flux_proj)); 31 diff_flux_proj = *pdiff_flux_proj; 32 PetscCall(PetscNew(&user->diff_flux_proj->projection)); 33 projection = user->diff_flux_proj->projection; 34 diff_flux_proj->method = user->app_ctx->divFdiffproj_method; 35 diff_flux_proj->num_diff_flux_comps = num_diff_flux_comps; 36 37 PetscCall(DMClone(user->dm, &projection->dm)); 38 PetscCall(DMSetMatrixPreallocateSkip(projection->dm, PETSC_TRUE)); 39 PetscCall(DMGetDimension(projection->dm, &dim)); 40 switch (diff_flux_proj->method) { 41 case DIV_DIFF_FLUX_PROJ_DIRECT: { 42 projection->num_comp = diff_flux_proj->num_diff_flux_comps; 43 PetscCall(PetscObjectSetName((PetscObject)projection->dm, "DivDiffFluxProj")); 44 PetscCall(DMSetupByOrder_FEM(PETSC_TRUE, PETSC_TRUE, degree, 1, q_extra, 1, &projection->num_comp, projection->dm)); 45 46 PetscCall(DMPlexCeedElemRestrictionCreate(user->ceed, projection->dm, domain_label, label_value, height, dm_field, 47 &diff_flux_proj->elem_restr_div_diff_flux)); 48 PetscCallCeed(user->ceed, CeedElemRestrictionCreateVector(diff_flux_proj->elem_restr_div_diff_flux, &diff_flux_proj->div_diff_flux_ceed, NULL)); 49 PetscCall(CreateBasisFromPlex(user->ceed, projection->dm, domain_label, label_value, height, dm_field, &diff_flux_proj->basis_div_diff_flux)); 50 diff_flux_proj->eval_mode_div_diff_flux = CEED_EVAL_INTERP; 51 52 { // Create face labels on projection->dm for boundary integrals 53 DMLabel face_sets_label; 54 PetscInt num_face_set_values, *face_set_values; 55 56 PetscCall(DMGetLabel(user->dm, "Face Sets", &face_sets_label)); 57 PetscCall(DMLabelCreateGlobalValueArray(user->dm, face_sets_label, &num_face_set_values, &face_set_values)); 58 for (PetscInt f = 0; f < num_face_set_values; f++) { 59 DMLabel face_orientation_label; 60 char *face_orientation_label_name; 61 62 PetscCall(DMPlexCreateFaceLabel(user->dm, face_set_values[f], &face_orientation_label_name)); 63 PetscCall(DMGetLabel(user->dm, face_orientation_label_name, &face_orientation_label)); 64 PetscCall(DMAddLabel(projection->dm, face_orientation_label)); 65 PetscCall(PetscFree(face_orientation_label_name)); 66 } 67 PetscCall(PetscFree(face_set_values)); 68 } 69 } break; 70 case DIV_DIFF_FLUX_PROJ_INDIRECT: { 71 projection->num_comp = diff_flux_proj->num_diff_flux_comps * dim; 72 PetscCall(PetscObjectSetName((PetscObject)projection->dm, "DiffFluxProj")); 73 PetscCall(DMSetupByOrder_FEM(PETSC_TRUE, PETSC_TRUE, degree, 1, q_extra, 1, &projection->num_comp, projection->dm)); 74 75 PetscCall(DMPlexCeedElemRestrictionQDataCreate(user->ceed, projection->dm, domain_label, label_value, height, 76 diff_flux_proj->num_diff_flux_comps, &diff_flux_proj->elem_restr_div_diff_flux)); 77 PetscCallCeed(user->ceed, CeedElemRestrictionCreateVector(diff_flux_proj->elem_restr_div_diff_flux, &diff_flux_proj->div_diff_flux_ceed, NULL)); 78 diff_flux_proj->basis_div_diff_flux = CEED_BASIS_NONE; 79 diff_flux_proj->eval_mode_div_diff_flux = CEED_EVAL_NONE; 80 } break; 81 case DIV_DIFF_FLUX_PROJ_NONE: 82 SETERRQ(PetscObjectComm((PetscObject)user->dm), PETSC_ERR_ARG_WRONG, "Should not reach here with div_diff_flux_projection_method %s", 83 DivDiffFluxProjectionMethods[user->app_ctx->divFdiffproj_method]); 84 break; 85 } 86 PetscFunctionReturn(PETSC_SUCCESS); 87 }; 88 89 /** 90 @brief Return the objects required for the Divergence of Diffusive flux to be read by a `CeedOperator` 91 92 @param[in] diff_flux_proj Projection object 93 @param[out] elem_restr Element restriction for the divergence of diffusive flux, or `NULL` 94 @param[out] basis Basis for the divergence of diffusive flux, or `NULL` 95 @param[out] vector Vector for the divergence of diffusive flux, or `NULL` 96 @param[out] eval_mode Eval mode for the divergence of diffusive flux, or `NULL` 97 **/ 98 PetscErrorCode DivDiffFluxProjectionGetOperatorFieldData(DivDiffFluxProjectionData diff_flux_proj, CeedElemRestriction *elem_restr, CeedBasis *basis, 99 CeedVector *vector, CeedEvalMode *eval_mode) { 100 Ceed ceed = CeedVectorReturnCeed(diff_flux_proj->div_diff_flux_ceed); 101 102 PetscFunctionBeginUser; 103 if (elem_restr) PetscCallCeed(ceed, CeedElemRestrictionReferenceCopy(diff_flux_proj->elem_restr_div_diff_flux, elem_restr)); 104 if (basis) PetscCallCeed(ceed, CeedBasisReferenceCopy(diff_flux_proj->basis_div_diff_flux, basis)); 105 if (vector) PetscCallCeed(ceed, CeedVectorReferenceCopy(diff_flux_proj->div_diff_flux_ceed, vector)); 106 if (eval_mode) *eval_mode = diff_flux_proj->eval_mode_div_diff_flux; 107 PetscFunctionReturn(PETSC_SUCCESS); 108 } 109 110 /** 111 @brief Setup direct projection of divergence of diffusive flux 112 113 @param[in] user `User` context 114 @param[in] ceed_data `CeedData` context 115 @param[in,out] diff_flux_proj Flux projection object to setup 116 **/ 117 static PetscErrorCode DivDiffFluxProjectionSetup_Direct(User user, CeedData ceed_data, DivDiffFluxProjectionData diff_flux_proj) { 118 Ceed ceed = user->ceed; 119 NodalProjectionData projection = diff_flux_proj->projection; 120 MPI_Comm comm = PetscObjectComm((PetscObject)projection->dm); 121 122 PetscFunctionBeginUser; 123 { // Create Projection RHS OperatorApplyContext 124 CeedOperator op_rhs; 125 126 PetscCheck(diff_flux_proj->CreateRHSOperator_Direct, comm, PETSC_ERR_ARG_WRONGSTATE, 127 "Must define CreateRHSOperator_Direct to use indirect div_diff_flux projection"); 128 PetscCall(diff_flux_proj->CreateRHSOperator_Direct(user, ceed_data, diff_flux_proj, &op_rhs)); 129 PetscCall(DMCreateLocalVector(projection->dm, &diff_flux_proj->DivDiffFlux_loc)); 130 diff_flux_proj->ceed_vec_has_array = PETSC_FALSE; 131 PetscCall(OperatorApplyContextCreate(user->dm, projection->dm, ceed, op_rhs, NULL, NULL, NULL, diff_flux_proj->DivDiffFlux_loc, 132 &projection->l2_rhs_ctx)); 133 PetscCallCeed(ceed, CeedOperatorDestroy(&op_rhs)); 134 } 135 136 { // -- Build Mass operator 137 CeedQFunction qf_mass; 138 CeedOperator op_mass; 139 CeedBasis basis_div_diff_flux = NULL; 140 CeedElemRestriction elem_restr_div_diff_flux_volume = NULL, elem_restr_qd; 141 CeedVector q_data; 142 CeedInt q_data_size; 143 PetscInt label_value = 0; 144 DMLabel domain_label = NULL; 145 146 PetscCall(DivDiffFluxProjectionGetOperatorFieldData(diff_flux_proj, &elem_restr_div_diff_flux_volume, &basis_div_diff_flux, NULL, NULL)); 147 PetscCall(QDataGet(ceed, projection->dm, domain_label, label_value, ceed_data->elem_restr_x, ceed_data->basis_x, ceed_data->x_coord, 148 &elem_restr_qd, &q_data, &q_data_size)); 149 150 PetscCall(CreateMassQFunction(ceed, projection->num_comp, q_data_size, &qf_mass)); 151 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_mass, NULL, NULL, &op_mass)); 152 PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "u", elem_restr_div_diff_flux_volume, basis_div_diff_flux, CEED_VECTOR_ACTIVE)); 153 PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "qdata", elem_restr_qd, CEED_BASIS_NONE, q_data)); 154 PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "v", elem_restr_div_diff_flux_volume, basis_div_diff_flux, CEED_VECTOR_ACTIVE)); 155 156 { // -- Setup KSP for L^2 projection 157 Mat mat_mass; 158 159 PetscCall(MatCreateCeed(projection->dm, projection->dm, op_mass, NULL, &mat_mass)); 160 161 PetscCall(KSPCreate(comm, &projection->ksp)); 162 PetscCall(KSPSetOptionsPrefix(projection->ksp, "div_diff_flux_projection_")); 163 { // lumped by default 164 PC pc; 165 PetscCall(KSPGetPC(projection->ksp, &pc)); 166 PetscCall(PCSetType(pc, PCJACOBI)); 167 PetscCall(PCJacobiSetType(pc, PC_JACOBI_ROWSUM)); 168 PetscCall(KSPSetType(projection->ksp, KSPPREONLY)); 169 } 170 PetscCall(KSPSetFromOptions_WithMatCeed(projection->ksp, mat_mass)); 171 PetscCall(MatDestroy(&mat_mass)); 172 } 173 174 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_div_diff_flux_volume)); 175 PetscCallCeed(ceed, CeedBasisDestroy(&basis_div_diff_flux)); 176 PetscCallCeed(ceed, CeedVectorDestroy(&q_data)); 177 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd)); 178 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_mass)); 179 PetscCallCeed(ceed, CeedOperatorDestroy(&op_mass)); 180 } 181 PetscFunctionReturn(PETSC_SUCCESS); 182 } 183 184 /** 185 @brief Setup indirect projection of divergence of diffusive flux 186 187 @param[in] user `User` context 188 @param[in] ceed_data `CeedData` context 189 @param[in,out] diff_flux_proj Flux projection object to setup 190 **/ 191 static PetscErrorCode DivDiffFluxProjectionSetup_Indirect(User user, CeedData ceed_data, DivDiffFluxProjectionData diff_flux_proj) { 192 Ceed ceed = user->ceed; 193 NodalProjectionData projection = diff_flux_proj->projection; 194 CeedBasis basis_diff_flux; 195 CeedElemRestriction elem_restr_diff_flux, elem_restr_qd; 196 CeedVector q_data; 197 CeedInt q_data_size; 198 MPI_Comm comm = PetscObjectComm((PetscObject)projection->dm); 199 200 PetscFunctionBeginUser; 201 { 202 PetscInt label_value = 0, height = 0, dm_field = 0; 203 DMLabel domain_label = NULL; 204 205 PetscCall(DMPlexCeedElemRestrictionCreate(ceed, projection->dm, domain_label, label_value, height, dm_field, &elem_restr_diff_flux)); 206 PetscCall(CreateBasisFromPlex(ceed, projection->dm, domain_label, label_value, height, dm_field, &basis_diff_flux)); 207 PetscCall(QDataGet(ceed, projection->dm, domain_label, label_value, ceed_data->elem_restr_x, ceed_data->basis_x, ceed_data->x_coord, 208 &elem_restr_qd, &q_data, &q_data_size)); 209 } 210 211 { 212 CeedOperator op_rhs; 213 214 PetscCheck(diff_flux_proj->CreateRHSOperator_Indirect, comm, PETSC_ERR_ARG_WRONGSTATE, 215 "Must define CreateRHSOperator_Indirect to use indirect div_diff_flux projection"); 216 PetscCall(diff_flux_proj->CreateRHSOperator_Indirect(user, ceed_data, diff_flux_proj, &op_rhs)); 217 PetscCall(OperatorApplyContextCreate(user->dm, projection->dm, ceed, op_rhs, NULL, NULL, NULL, NULL, &projection->l2_rhs_ctx)); 218 PetscCallCeed(ceed, CeedOperatorDestroy(&op_rhs)); 219 } 220 221 { // -- Build Mass operator 222 CeedQFunction qf_mass; 223 CeedOperator op_mass; 224 225 PetscCall(CreateMassQFunction(ceed, projection->num_comp, q_data_size, &qf_mass)); 226 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_mass, NULL, NULL, &op_mass)); 227 PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "u", elem_restr_diff_flux, basis_diff_flux, CEED_VECTOR_ACTIVE)); 228 PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "qdata", elem_restr_qd, CEED_BASIS_NONE, q_data)); 229 PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "v", elem_restr_diff_flux, basis_diff_flux, CEED_VECTOR_ACTIVE)); 230 231 { // -- Setup KSP for L^2 projection 232 Mat mat_mass; 233 234 PetscCall(MatCreateCeed(projection->dm, projection->dm, op_mass, NULL, &mat_mass)); 235 236 PetscCall(KSPCreate(comm, &projection->ksp)); 237 PetscCall(KSPSetOptionsPrefix(projection->ksp, "div_diff_flux_projection_")); 238 { // lumped by default 239 PC pc; 240 PetscCall(KSPGetPC(projection->ksp, &pc)); 241 PetscCall(PCSetType(pc, PCJACOBI)); 242 PetscCall(PCJacobiSetType(pc, PC_JACOBI_ROWSUM)); 243 PetscCall(KSPSetType(projection->ksp, KSPPREONLY)); 244 } 245 PetscCall(KSPSetFromOptions_WithMatCeed(projection->ksp, mat_mass)); 246 PetscCall(MatDestroy(&mat_mass)); 247 } 248 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_mass)); 249 PetscCallCeed(ceed, CeedOperatorDestroy(&op_mass)); 250 } 251 252 { // Create OperatorApplyContext to calculate divergence at quadrature points 253 CeedQFunction qf_calc_divergence = NULL; 254 CeedOperator op_calc_divergence; 255 CeedElemRestriction elem_restr_div_diff_flux = NULL; 256 PetscInt dim; 257 258 PetscCall(DMGetDimension(projection->dm, &dim)); 259 PetscCall(DivDiffFluxProjectionGetOperatorFieldData(diff_flux_proj, &elem_restr_div_diff_flux, NULL, NULL, NULL)); 260 261 switch (dim) { 262 case 2: 263 switch (diff_flux_proj->num_diff_flux_comps) { 264 case 1: 265 PetscCallCeed(ceed, 266 CeedQFunctionCreateInterior(ceed, 1, ComputeDivDiffusiveFlux2D_1, ComputeDivDiffusiveFlux2D_1_loc, &qf_calc_divergence)); 267 break; 268 } 269 break; 270 case 3: 271 switch (diff_flux_proj->num_diff_flux_comps) { 272 case 1: 273 PetscCallCeed(ceed, 274 CeedQFunctionCreateInterior(ceed, 1, ComputeDivDiffusiveFlux3D_1, ComputeDivDiffusiveFlux3D_1_loc, &qf_calc_divergence)); 275 break; 276 case 4: 277 PetscCallCeed(ceed, 278 CeedQFunctionCreateInterior(ceed, 1, ComputeDivDiffusiveFlux3D_4, ComputeDivDiffusiveFlux3D_4_loc, &qf_calc_divergence)); 279 break; 280 } 281 break; 282 } 283 PetscCheck(qf_calc_divergence, comm, PETSC_ERR_SUP, 284 "QFunction for calculating divergence of diffusive flux does not exist for" 285 " %" PetscInt_FMT " dimensional grid and %" PetscInt_FMT 286 " number of components.\nA new qfunction can be easily added; see source code for pattern.", 287 dim, diff_flux_proj->num_diff_flux_comps); 288 289 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_calc_divergence, "Grad F_diff", projection->num_comp * dim, CEED_EVAL_GRAD)); 290 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_calc_divergence, "qdata", q_data_size, CEED_EVAL_NONE)); 291 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_calc_divergence, "Div F_diff", diff_flux_proj->num_diff_flux_comps, CEED_EVAL_NONE)); 292 293 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_calc_divergence, NULL, NULL, &op_calc_divergence)); 294 PetscCallCeed(ceed, CeedOperatorSetField(op_calc_divergence, "Grad F_diff", elem_restr_diff_flux, basis_diff_flux, CEED_VECTOR_ACTIVE)); 295 PetscCallCeed(ceed, CeedOperatorSetField(op_calc_divergence, "qdata", elem_restr_qd, CEED_BASIS_NONE, q_data)); 296 PetscCallCeed( 297 ceed, CeedOperatorSetField(op_calc_divergence, "Div F_diff", elem_restr_div_diff_flux, CEED_BASIS_NONE, diff_flux_proj->div_diff_flux_ceed)); 298 299 PetscCall( 300 OperatorApplyContextCreate(projection->dm, NULL, ceed, op_calc_divergence, NULL, NULL, NULL, NULL, &diff_flux_proj->calc_div_diff_flux)); 301 302 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_div_diff_flux)); 303 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_calc_divergence)); 304 PetscCallCeed(ceed, CeedOperatorDestroy(&op_calc_divergence)); 305 } 306 PetscCallCeed(ceed, CeedBasisDestroy(&basis_diff_flux)); 307 PetscCallCeed(ceed, CeedVectorDestroy(&q_data)); 308 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd)); 309 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_diff_flux)); 310 PetscFunctionReturn(PETSC_SUCCESS); 311 } 312 313 /** 314 @brief Setup projection of divergence of diffusive flux 315 316 @param[in] user `User` context 317 @param[in] ceed_data `CeedData` context 318 @param[in,out] diff_flux_proj Flux projection object to setup 319 **/ 320 PetscErrorCode DivDiffFluxProjectionSetup(User user, CeedData ceed_data, DivDiffFluxProjectionData diff_flux_proj) { 321 PetscFunctionBeginUser; 322 switch (user->app_ctx->divFdiffproj_method) { 323 case DIV_DIFF_FLUX_PROJ_DIRECT: 324 PetscCall(DivDiffFluxProjectionSetup_Direct(user, ceed_data, diff_flux_proj)); 325 break; 326 case DIV_DIFF_FLUX_PROJ_INDIRECT: 327 PetscCall(DivDiffFluxProjectionSetup_Indirect(user, ceed_data, diff_flux_proj)); 328 break; 329 case DIV_DIFF_FLUX_PROJ_NONE: 330 SETERRQ(PetscObjectComm((PetscObject)user->dm), PETSC_ERR_ARG_WRONG, "Should not reach here with div_diff_flux_projection_method %s", 331 DivDiffFluxProjectionMethods[user->app_ctx->divFdiffproj_method]); 332 break; 333 } 334 PetscFunctionReturn(PETSC_SUCCESS); 335 } 336 337 /** 338 @brief Project the divergence of diffusive flux 339 340 This implicitly sets the `CeedVector` input (`div_diff_flux_ceed`) to the divergence of diffusive flux. 341 342 @param[in] diff_flux_proj `NodalProjectionData` for the projection 343 @param[in] Q_loc Localized solution vector 344 **/ 345 PetscErrorCode DivDiffFluxProjectionApply(DivDiffFluxProjectionData diff_flux_proj, Vec Q_loc) { 346 NodalProjectionData projection = diff_flux_proj->projection; 347 348 PetscFunctionBeginUser; 349 PetscCall(PetscLogEventBegin(FLUIDS_DivDiffFluxProjection, Q_loc, 0, 0, 0)); 350 switch (diff_flux_proj->method) { 351 case DIV_DIFF_FLUX_PROJ_DIRECT: { 352 Vec DivDiffFlux; 353 354 PetscCall(DMGetGlobalVector(projection->dm, &DivDiffFlux)); 355 if (diff_flux_proj->ceed_vec_has_array) { 356 PetscCall(VecReadCeedToPetsc(diff_flux_proj->div_diff_flux_ceed, diff_flux_proj->DivDiffFlux_memtype, diff_flux_proj->DivDiffFlux_loc)); 357 diff_flux_proj->ceed_vec_has_array = PETSC_FALSE; 358 } 359 PetscCall(ApplyCeedOperatorLocalToGlobal(Q_loc, DivDiffFlux, projection->l2_rhs_ctx)); 360 PetscCall(VecViewFromOptions(DivDiffFlux, NULL, "-div_diff_flux_projection_rhs_view")); 361 362 PetscCall(KSPSolve(projection->ksp, DivDiffFlux, DivDiffFlux)); 363 PetscCall(VecViewFromOptions(DivDiffFlux, NULL, "-div_diff_flux_projection_view")); 364 365 PetscCall(DMGlobalToLocal(projection->dm, DivDiffFlux, INSERT_VALUES, diff_flux_proj->DivDiffFlux_loc)); 366 PetscCall(VecReadPetscToCeed(diff_flux_proj->DivDiffFlux_loc, &diff_flux_proj->DivDiffFlux_memtype, diff_flux_proj->div_diff_flux_ceed)); 367 diff_flux_proj->ceed_vec_has_array = PETSC_TRUE; 368 369 PetscCall(DMRestoreGlobalVector(projection->dm, &DivDiffFlux)); 370 break; 371 } 372 case DIV_DIFF_FLUX_PROJ_INDIRECT: { 373 Vec DiffFlux; 374 375 PetscCall(DMGetGlobalVector(projection->dm, &DiffFlux)); 376 PetscCall(ApplyCeedOperatorLocalToGlobal(Q_loc, DiffFlux, projection->l2_rhs_ctx)); 377 PetscCall(VecViewFromOptions(DiffFlux, NULL, "-div_diff_flux_projection_rhs_view")); 378 379 PetscCall(KSPSolve(projection->ksp, DiffFlux, DiffFlux)); 380 PetscCall(VecViewFromOptions(DiffFlux, NULL, "-div_diff_flux_projection_view")); 381 382 PetscCall(ApplyCeedOperatorGlobalToLocal(DiffFlux, NULL, diff_flux_proj->calc_div_diff_flux)); 383 PetscCall(DMRestoreGlobalVector(projection->dm, &DiffFlux)); 384 } break; 385 case DIV_DIFF_FLUX_PROJ_NONE: 386 SETERRQ(PetscObjectComm((PetscObject)projection->dm), PETSC_ERR_ARG_WRONG, "Should not reach here with div_diff_flux_projection_method %s", 387 DivDiffFluxProjectionMethods[diff_flux_proj->method]); 388 break; 389 } 390 PetscCall(PetscLogEventEnd(FLUIDS_DivDiffFluxProjection, Q_loc, 0, 0, 0)); 391 PetscFunctionReturn(PETSC_SUCCESS); 392 } 393 394 /** 395 @brief Destroy `DivDiffFluxProjectionData` object 396 397 @param[in,out] diff_flux_proj Object to destroy 398 **/ 399 PetscErrorCode DivDiffFluxProjectionDataDestroy(DivDiffFluxProjectionData diff_flux_proj) { 400 PetscFunctionBeginUser; 401 if (diff_flux_proj == NULL) PetscFunctionReturn(PETSC_SUCCESS); 402 Ceed ceed = CeedVectorReturnCeed(diff_flux_proj->div_diff_flux_ceed); 403 404 PetscCall(NodalProjectionDataDestroy(diff_flux_proj->projection)); 405 PetscCall(OperatorApplyContextDestroy(diff_flux_proj->calc_div_diff_flux)); 406 if (diff_flux_proj->ceed_vec_has_array) { 407 PetscCall(VecReadCeedToPetsc(diff_flux_proj->div_diff_flux_ceed, diff_flux_proj->DivDiffFlux_memtype, diff_flux_proj->DivDiffFlux_loc)); 408 diff_flux_proj->ceed_vec_has_array = PETSC_FALSE; 409 } 410 PetscCallCeed(ceed, CeedVectorDestroy(&diff_flux_proj->div_diff_flux_ceed)); 411 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&diff_flux_proj->elem_restr_div_diff_flux)); 412 PetscCallCeed(ceed, CeedBasisDestroy(&diff_flux_proj->basis_div_diff_flux)); 413 PetscCall(VecDestroy(&diff_flux_proj->DivDiffFlux_loc)); 414 PetscCall(PetscFree(diff_flux_proj)); 415 PetscFunctionReturn(PETSC_SUCCESS); 416 } 417