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