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