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 Initialize projection of divergence of diffusive flux 14 15 Creates underlying `DM` for the projection operation and creates the restriction and basis to use with the CeedOperator 16 17 @param[in] user `User` context 18 @param[out] elem_restr_div_diff_flux `CeedElemRestriction` of the divergence of diffusive flux vector 19 @param[out] basis_div_diff_flux `CeedBasis` of the divergence of diffusive flux vector 20 @param[out] eval_mode_diff_flux `CeedEvalMode` for the divergence of the diffusive flux 21 **/ 22 PetscErrorCode DiffFluxProjectionInitialize(User user, CeedElemRestriction *elem_restr_div_diff_flux, CeedBasis *basis_div_diff_flux, 23 CeedEvalMode *eval_mode_diff_flux) { 24 PetscSection section; 25 PetscInt label_value = 0, height = 0, dm_field = 0, dim; 26 DMLabel domain_label = NULL; 27 DivDiffFluxProjectionData diff_flux_proj; 28 NodalProjectionData projection; 29 30 PetscFunctionBeginUser; 31 PetscCall(PetscNew(&user->diff_flux_proj)); 32 diff_flux_proj = user->diff_flux_proj; 33 PetscCall(PetscNew(&user->diff_flux_proj->projection)); 34 projection = user->diff_flux_proj->projection; 35 diff_flux_proj->method = user->app_ctx->divFdiffproj_method; 36 diff_flux_proj->num_diff_flux_comps = 4; 37 38 PetscCall(DMClone(user->dm, &projection->dm)); 39 PetscCall(DMSetMatrixPreallocateSkip(projection->dm, PETSC_TRUE)); 40 PetscCall(DMGetDimension(projection->dm, &dim)); 41 switch (diff_flux_proj->method) { 42 case DIV_DIFF_FLUX_PROJ_DIRECT: { 43 projection->num_comp = diff_flux_proj->num_diff_flux_comps; 44 PetscCall(PetscObjectSetName((PetscObject)projection->dm, "DivDiffFluxProj")); 45 PetscCall( 46 DMSetupByOrder_FEM(PETSC_TRUE, PETSC_TRUE, user->app_ctx->degree, 1, user->app_ctx->q_extra, 1, &projection->num_comp, projection->dm)); 47 48 PetscCall(DMGetLocalSection(projection->dm, §ion)); 49 PetscCall(PetscSectionSetFieldName(section, 0, "")); 50 PetscCall(PetscSectionSetComponentName(section, 0, 0, "DivDiffusiveFlux_MomentumX")); 51 PetscCall(PetscSectionSetComponentName(section, 0, 1, "DivDiffusiveFlux_MomentumY")); 52 PetscCall(PetscSectionSetComponentName(section, 0, 2, "DivDiffusiveFlux_MomentumZ")); 53 PetscCall(PetscSectionSetComponentName(section, 0, 3, "DivDiffusiveFlux_Energy")); 54 55 PetscCall(DMPlexCeedElemRestrictionCreate(user->ceed, projection->dm, domain_label, label_value, height, dm_field, elem_restr_div_diff_flux)); 56 PetscCallCeed(user->ceed, CeedElemRestrictionCreateVector(*elem_restr_div_diff_flux, &diff_flux_proj->div_diff_flux_ceed, NULL)); 57 PetscCall(CreateBasisFromPlex(user->ceed, projection->dm, domain_label, label_value, height, dm_field, basis_div_diff_flux)); 58 *eval_mode_diff_flux = CEED_EVAL_INTERP; 59 60 { // Create face labels on projection->dm for boundary integrals 61 DMLabel face_sets_label; 62 PetscInt num_face_set_values, *face_set_values; 63 64 PetscCall(DMGetLabel(user->dm, "Face Sets", &face_sets_label)); 65 PetscCall(DMLabelCreateGlobalValueArray(user->dm, face_sets_label, &num_face_set_values, &face_set_values)); 66 for (PetscInt f = 0; f < num_face_set_values; f++) { 67 DMLabel face_orientation_label; 68 char *face_orientation_label_name; 69 70 PetscCall(DMPlexCreateFaceLabel(user->dm, face_set_values[f], &face_orientation_label_name)); 71 PetscCall(DMGetLabel(user->dm, face_orientation_label_name, &face_orientation_label)); 72 PetscCall(DMAddLabel(projection->dm, face_orientation_label)); 73 PetscCall(PetscFree(face_orientation_label_name)); 74 } 75 } 76 } break; 77 case DIV_DIFF_FLUX_PROJ_INDIRECT: { 78 projection->num_comp = diff_flux_proj->num_diff_flux_comps * dim; 79 PetscCall(PetscObjectSetName((PetscObject)projection->dm, "DiffFluxProj")); 80 PetscCall( 81 DMSetupByOrder_FEM(PETSC_TRUE, PETSC_TRUE, user->app_ctx->degree, 1, user->app_ctx->q_extra, 1, &projection->num_comp, projection->dm)); 82 83 PetscCall(DMGetLocalSection(projection->dm, §ion)); 84 PetscCall(PetscSectionSetFieldName(section, 0, "")); 85 PetscCall(PetscSectionSetComponentName(section, 0, 0, "DiffusiveFlux_MomentumXX")); 86 PetscCall(PetscSectionSetComponentName(section, 0, 1, "DiffusiveFlux_MomentumXY")); 87 PetscCall(PetscSectionSetComponentName(section, 0, 2, "DiffusiveFlux_MomentumXZ")); 88 PetscCall(PetscSectionSetComponentName(section, 0, 3, "DiffusiveFlux_MomentumYX")); 89 PetscCall(PetscSectionSetComponentName(section, 0, 4, "DiffusiveFlux_MomentumYY")); 90 PetscCall(PetscSectionSetComponentName(section, 0, 5, "DiffusiveFlux_MomentumYZ")); 91 PetscCall(PetscSectionSetComponentName(section, 0, 6, "DiffusiveFlux_MomentumZX")); 92 PetscCall(PetscSectionSetComponentName(section, 0, 7, "DiffusiveFlux_MomentumZY")); 93 PetscCall(PetscSectionSetComponentName(section, 0, 8, "DiffusiveFlux_MomentumZZ")); 94 PetscCall(PetscSectionSetComponentName(section, 0, 9, "DiffusiveFlux_EnergyX")); 95 PetscCall(PetscSectionSetComponentName(section, 0, 10, "DiffusiveFlux_EnergyY")); 96 PetscCall(PetscSectionSetComponentName(section, 0, 11, "DiffusiveFlux_EnergyZ")); 97 98 PetscCall(DMPlexCeedElemRestrictionQDataCreate(user->ceed, projection->dm, domain_label, label_value, height, 99 diff_flux_proj->num_diff_flux_comps, elem_restr_div_diff_flux)); 100 PetscCallCeed(user->ceed, CeedElemRestrictionCreateVector(*elem_restr_div_diff_flux, &diff_flux_proj->div_diff_flux_ceed, NULL)); 101 *basis_div_diff_flux = CEED_BASIS_NONE; 102 *eval_mode_diff_flux = CEED_EVAL_NONE; 103 } break; 104 case DIV_DIFF_FLUX_PROJ_NONE: 105 SETERRQ(PetscObjectComm((PetscObject)user->dm), PETSC_ERR_ARG_WRONG, "Should not reach here with div_diff_flux_projection_method %s", 106 DivDiffFluxProjectionMethods[user->app_ctx->divFdiffproj_method]); 107 break; 108 } 109 PetscFunctionReturn(PETSC_SUCCESS); 110 }; 111 112 /** 113 @brief Setup direct projection of divergence of diffusive flux 114 115 @param[in] ceed `Ceed` context 116 @param[in] user `User` context 117 @param[in] ceed_data `CeedData` context 118 @param[in] problem `ProblemData` context 119 **/ 120 static PetscErrorCode DivDiffFluxProjectionSetup_Direct(Ceed ceed, User user, CeedData ceed_data, ProblemData problem) { 121 DivDiffFluxProjectionData diff_flux_proj = user->diff_flux_proj; 122 NodalProjectionData projection = diff_flux_proj->projection; 123 CeedOperator op_rhs; 124 CeedBasis basis_diff_flux; 125 CeedElemRestriction elem_restr_diff_flux_volume, elem_restr_qd; 126 CeedVector q_data; 127 CeedInt num_comp_q, q_data_size; 128 PetscInt dim, label_value = 0; 129 DMLabel domain_label = NULL; 130 131 PetscFunctionBeginUser; 132 // -- Get Pre-requisite things 133 PetscCall(DMGetDimension(projection->dm, &dim)); 134 PetscCallCeed(ceed, CeedBasisGetNumComponents(ceed_data->basis_q, &num_comp_q)); 135 136 { // Get elem_restr_diff_flux and basis_diff_flux 137 CeedOperator *sub_ops; 138 CeedOperatorField op_field; 139 PetscInt sub_op_index = 0; // will be 0 for the volume op 140 141 PetscCallCeed(ceed, CeedCompositeOperatorGetSubList(user->op_ifunction, &sub_ops)); 142 PetscCallCeed(ceed, CeedOperatorGetFieldByName(sub_ops[sub_op_index], "div F_diff", &op_field)); 143 PetscCallCeed(ceed, CeedOperatorFieldGetElemRestriction(op_field, &elem_restr_diff_flux_volume)); 144 PetscCallCeed(ceed, CeedOperatorFieldGetBasis(op_field, &basis_diff_flux)); 145 } 146 PetscCall(QDataGet(ceed, projection->dm, domain_label, label_value, ceed_data->elem_restr_x, ceed_data->basis_x, ceed_data->x_coord, &elem_restr_qd, 147 &q_data, &q_data_size)); 148 149 PetscCallCeed(ceed, CeedCompositeOperatorCreate(ceed, &op_rhs)); 150 { // Add the volume integral CeedOperator 151 CeedQFunction qf_rhs_volume; 152 CeedOperator op_rhs_volume; 153 154 switch (user->phys->state_var) { 155 case STATEVAR_PRIMITIVE: 156 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, DivDiffusiveFluxVolumeRHS_Prim, DivDiffusiveFluxVolumeRHS_Prim_loc, &qf_rhs_volume)); 157 break; 158 case STATEVAR_CONSERVATIVE: 159 PetscCallCeed(ceed, 160 CeedQFunctionCreateInterior(ceed, 1, DivDiffusiveFluxVolumeRHS_Conserv, DivDiffusiveFluxVolumeRHS_Conserv_loc, &qf_rhs_volume)); 161 break; 162 case STATEVAR_ENTROPY: 163 PetscCallCeed(ceed, 164 CeedQFunctionCreateInterior(ceed, 1, DivDiffusiveFluxVolumeRHS_Entropy, DivDiffusiveFluxVolumeRHS_Entropy_loc, &qf_rhs_volume)); 165 break; 166 } 167 168 PetscCallCeed(ceed, CeedQFunctionSetContext(qf_rhs_volume, problem->apply_vol_ifunction.qfctx)); 169 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_rhs_volume, "q", num_comp_q, CEED_EVAL_INTERP)); 170 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_rhs_volume, "Grad_q", num_comp_q * dim, CEED_EVAL_GRAD)); 171 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_rhs_volume, "qdata", q_data_size, CEED_EVAL_NONE)); 172 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_rhs_volume, "diffusive flux RHS", projection->num_comp * dim, CEED_EVAL_GRAD)); 173 174 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_rhs_volume, NULL, NULL, &op_rhs_volume)); 175 PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_volume, "q", ceed_data->elem_restr_q, ceed_data->basis_q, CEED_VECTOR_ACTIVE)); 176 PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_volume, "Grad_q", ceed_data->elem_restr_q, ceed_data->basis_q, CEED_VECTOR_ACTIVE)); 177 PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_volume, "qdata", elem_restr_qd, CEED_BASIS_NONE, q_data)); 178 PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_volume, "diffusive flux RHS", elem_restr_diff_flux_volume, basis_diff_flux, CEED_VECTOR_ACTIVE)); 179 180 PetscCallCeed(ceed, CeedCompositeOperatorAddSub(op_rhs, op_rhs_volume)); 181 182 PetscCallCeed(ceed, CeedOperatorDestroy(&op_rhs_volume)); 183 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_rhs_volume)); 184 } 185 186 { // Add the boundary integral CeedOperator 187 CeedQFunction qf_rhs_boundary; 188 DMLabel face_sets_label; 189 PetscInt num_face_set_values, *face_set_values; 190 CeedInt q_data_size; 191 192 // -- Build RHS operator 193 switch (user->phys->state_var) { 194 case STATEVAR_PRIMITIVE: 195 PetscCallCeed(ceed, 196 CeedQFunctionCreateInterior(ceed, 1, DivDiffusiveFluxBoundaryRHS_Prim, DivDiffusiveFluxBoundaryRHS_Prim_loc, &qf_rhs_boundary)); 197 break; 198 case STATEVAR_CONSERVATIVE: 199 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, DivDiffusiveFluxBoundaryRHS_Conserv, DivDiffusiveFluxBoundaryRHS_Conserv_loc, 200 &qf_rhs_boundary)); 201 break; 202 case STATEVAR_ENTROPY: 203 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, DivDiffusiveFluxBoundaryRHS_Entropy, DivDiffusiveFluxBoundaryRHS_Entropy_loc, 204 &qf_rhs_boundary)); 205 break; 206 } 207 208 PetscCall(QDataBoundaryGradientGetNumComponents(user->dm, &q_data_size)); 209 PetscCallCeed(ceed, CeedQFunctionSetContext(qf_rhs_boundary, problem->apply_vol_ifunction.qfctx)); 210 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_rhs_boundary, "q", num_comp_q, CEED_EVAL_INTERP)); 211 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_rhs_boundary, "Grad_q", num_comp_q * dim, CEED_EVAL_GRAD)); 212 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_rhs_boundary, "qdata", q_data_size, CEED_EVAL_NONE)); 213 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_rhs_boundary, "diffusive flux RHS", projection->num_comp, CEED_EVAL_INTERP)); 214 215 PetscCall(DMGetLabel(projection->dm, "Face Sets", &face_sets_label)); 216 PetscCall(DMLabelCreateGlobalValueArray(projection->dm, face_sets_label, &num_face_set_values, &face_set_values)); 217 for (PetscInt f = 0; f < num_face_set_values; f++) { 218 DMLabel face_orientation_label; 219 PetscInt num_orientations_values, *orientation_values; 220 221 { 222 char *face_orientation_label_name; 223 224 PetscCall(DMPlexCreateFaceLabel(projection->dm, face_set_values[f], &face_orientation_label_name)); 225 PetscCall(DMGetLabel(projection->dm, face_orientation_label_name, &face_orientation_label)); 226 PetscCall(PetscFree(face_orientation_label_name)); 227 } 228 PetscCall(DMLabelCreateGlobalValueArray(projection->dm, face_orientation_label, &num_orientations_values, &orientation_values)); 229 for (PetscInt o = 0; o < num_orientations_values; o++) { 230 CeedOperator op_rhs_boundary; 231 CeedBasis basis_q, basis_diff_flux_boundary; 232 CeedElemRestriction elem_restr_qdata, elem_restr_q, elem_restr_diff_flux_boundary; 233 CeedVector q_data; 234 CeedInt q_data_size; 235 PetscInt orientation = orientation_values[o], dm_field_q = 0, height_cell = 0, height_face = 1; 236 237 PetscCall(DMPlexCeedElemRestrictionCreate(ceed, user->dm, face_orientation_label, orientation, height_cell, dm_field_q, &elem_restr_q)); 238 PetscCall(DMPlexCeedElemRestrictionCreate(ceed, projection->dm, face_orientation_label, orientation, height_face, 0, 239 &elem_restr_diff_flux_boundary)); 240 PetscCall(QDataBoundaryGradientGet(ceed, user->dm, face_orientation_label, orientation, ceed_data->x_coord, &elem_restr_qdata, &q_data, 241 &q_data_size)); 242 PetscCall(DMPlexCeedBasisCellToFaceCreate(ceed, user->dm, face_orientation_label, orientation, orientation, dm_field_q, &basis_q)); 243 PetscCall(CreateBasisFromPlex(ceed, projection->dm, face_orientation_label, orientation, height_face, 0, &basis_diff_flux_boundary)); 244 245 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_rhs_boundary, NULL, NULL, &op_rhs_boundary)); 246 PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_boundary, "q", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE)); 247 PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_boundary, "Grad_q", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE)); 248 PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_boundary, "qdata", elem_restr_qdata, CEED_BASIS_NONE, q_data)); 249 PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_boundary, "diffusive flux RHS", elem_restr_diff_flux_boundary, basis_diff_flux_boundary, 250 CEED_VECTOR_ACTIVE)); 251 252 PetscCallCeed(ceed, CeedCompositeOperatorAddSub(op_rhs, op_rhs_boundary)); 253 254 PetscCallCeed(ceed, CeedOperatorDestroy(&op_rhs_boundary)); 255 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qdata)); 256 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_q)); 257 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_diff_flux_boundary)); 258 PetscCallCeed(ceed, CeedBasisDestroy(&basis_q)); 259 PetscCallCeed(ceed, CeedBasisDestroy(&basis_diff_flux_boundary)); 260 PetscCallCeed(ceed, CeedVectorDestroy(&q_data)); 261 } 262 PetscCall(PetscFree(orientation_values)); 263 } 264 PetscCall(PetscFree(face_set_values)); 265 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_rhs_boundary)); 266 } 267 268 PetscCall(DMCreateLocalVector(projection->dm, &diff_flux_proj->DivDiffFlux_loc)); 269 diff_flux_proj->ceed_vec_has_array = PETSC_FALSE; 270 PetscCall( 271 OperatorApplyContextCreate(user->dm, projection->dm, ceed, op_rhs, NULL, NULL, NULL, diff_flux_proj->DivDiffFlux_loc, &projection->l2_rhs_ctx)); 272 273 { // -- Build Mass operator 274 CeedQFunction qf_mass; 275 CeedOperator op_mass; 276 277 PetscCall(CreateMassQFunction(ceed, projection->num_comp, q_data_size, &qf_mass)); 278 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_mass, NULL, NULL, &op_mass)); 279 PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "u", elem_restr_diff_flux_volume, basis_diff_flux, CEED_VECTOR_ACTIVE)); 280 PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "qdata", elem_restr_qd, CEED_BASIS_NONE, q_data)); 281 PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "v", elem_restr_diff_flux_volume, basis_diff_flux, CEED_VECTOR_ACTIVE)); 282 283 { // -- Setup KSP for L^2 projection 284 Mat mat_mass; 285 MPI_Comm comm = PetscObjectComm((PetscObject)projection->dm); 286 287 PetscCall(MatCreateCeed(projection->dm, projection->dm, op_mass, NULL, &mat_mass)); 288 289 PetscCall(KSPCreate(comm, &projection->ksp)); 290 PetscCall(KSPSetOptionsPrefix(projection->ksp, "div_diff_flux_projection_")); 291 { // lumped by default 292 PC pc; 293 PetscCall(KSPGetPC(projection->ksp, &pc)); 294 PetscCall(PCSetType(pc, PCJACOBI)); 295 PetscCall(PCJacobiSetType(pc, PC_JACOBI_ROWSUM)); 296 PetscCall(KSPSetType(projection->ksp, KSPPREONLY)); 297 } 298 PetscCall(KSPSetFromOptions_WithMatCeed(projection->ksp, mat_mass)); 299 PetscCall(MatDestroy(&mat_mass)); 300 } 301 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_mass)); 302 PetscCallCeed(ceed, CeedOperatorDestroy(&op_mass)); 303 } 304 PetscCallCeed(ceed, CeedVectorDestroy(&q_data)); 305 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd)); 306 PetscCallCeed(ceed, CeedOperatorDestroy(&op_rhs)); 307 PetscFunctionReturn(PETSC_SUCCESS); 308 } 309 310 /** 311 @brief Setup indirect projection of divergence of diffusive flux 312 313 @param[in] ceed `Ceed` context 314 @param[in,out] user `User` context 315 @param[in] ceed_data `CeedData` context 316 @param[in] problem `ProblemData` context 317 **/ 318 static PetscErrorCode DivDiffFluxProjectionSetup_Indirect(Ceed ceed, User user, CeedData ceed_data, ProblemData problem) { 319 DivDiffFluxProjectionData diff_flux_proj = user->diff_flux_proj; 320 NodalProjectionData projection = diff_flux_proj->projection; 321 CeedBasis basis_diff_flux; 322 CeedElemRestriction elem_restr_diff_flux, elem_restr_qd; 323 CeedVector q_data; 324 CeedInt num_comp_q, q_data_size; 325 PetscInt dim; 326 PetscInt label_value = 0, height = 0, dm_field = 0; 327 DMLabel domain_label = NULL; 328 329 PetscFunctionBeginUser; 330 PetscCall(DMGetDimension(projection->dm, &dim)); 331 PetscCallCeed(ceed, CeedBasisGetNumComponents(ceed_data->basis_q, &num_comp_q)); 332 333 PetscCall(DMPlexCeedElemRestrictionCreate(ceed, projection->dm, domain_label, label_value, height, dm_field, &elem_restr_diff_flux)); 334 PetscCall(CreateBasisFromPlex(ceed, projection->dm, domain_label, label_value, height, dm_field, &basis_diff_flux)); 335 PetscCall(QDataGet(ceed, projection->dm, domain_label, label_value, ceed_data->elem_restr_x, ceed_data->basis_x, ceed_data->x_coord, &elem_restr_qd, 336 &q_data, &q_data_size)); 337 338 { // Create RHS CeedOperator for L^2 projection 339 CeedQFunction qf_rhs; 340 CeedOperator op_rhs; 341 342 switch (user->phys->state_var) { 343 case STATEVAR_PRIMITIVE: 344 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, DiffusiveFluxRHS_Prim, DiffusiveFluxRHS_Prim_loc, &qf_rhs)); 345 break; 346 case STATEVAR_CONSERVATIVE: 347 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, DiffusiveFluxRHS_Conserv, DiffusiveFluxRHS_Conserv_loc, &qf_rhs)); 348 break; 349 case STATEVAR_ENTROPY: 350 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, DiffusiveFluxRHS_Entropy, DiffusiveFluxRHS_Entropy_loc, &qf_rhs)); 351 break; 352 } 353 354 PetscCallCeed(ceed, CeedQFunctionSetContext(qf_rhs, problem->apply_vol_ifunction.qfctx)); 355 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_rhs, "q", num_comp_q, CEED_EVAL_INTERP)); 356 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_rhs, "Grad_q", num_comp_q * dim, CEED_EVAL_GRAD)); 357 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_rhs, "qdata", q_data_size, CEED_EVAL_NONE)); 358 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_rhs, "F_diff RHS", projection->num_comp, CEED_EVAL_INTERP)); 359 360 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_rhs, NULL, NULL, &op_rhs)); 361 PetscCallCeed(ceed, CeedOperatorSetField(op_rhs, "q", ceed_data->elem_restr_q, ceed_data->basis_q, CEED_VECTOR_ACTIVE)); 362 PetscCallCeed(ceed, CeedOperatorSetField(op_rhs, "Grad_q", ceed_data->elem_restr_q, ceed_data->basis_q, CEED_VECTOR_ACTIVE)); 363 PetscCallCeed(ceed, CeedOperatorSetField(op_rhs, "qdata", elem_restr_qd, CEED_BASIS_NONE, q_data)); 364 PetscCallCeed(ceed, CeedOperatorSetField(op_rhs, "F_diff RHS", elem_restr_diff_flux, basis_diff_flux, CEED_VECTOR_ACTIVE)); 365 366 PetscCall(OperatorApplyContextCreate(user->dm, projection->dm, ceed, op_rhs, NULL, NULL, NULL, NULL, &projection->l2_rhs_ctx)); 367 368 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_rhs)); 369 PetscCallCeed(ceed, CeedOperatorDestroy(&op_rhs)); 370 } 371 372 { // -- Build Mass operator 373 CeedQFunction qf_mass; 374 CeedOperator op_mass; 375 376 PetscCall(CreateMassQFunction(ceed, projection->num_comp, q_data_size, &qf_mass)); 377 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_mass, NULL, NULL, &op_mass)); 378 PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "u", elem_restr_diff_flux, basis_diff_flux, CEED_VECTOR_ACTIVE)); 379 PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "qdata", elem_restr_qd, CEED_BASIS_NONE, q_data)); 380 PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "v", elem_restr_diff_flux, basis_diff_flux, CEED_VECTOR_ACTIVE)); 381 382 { // -- Setup KSP for L^2 projection 383 Mat mat_mass; 384 MPI_Comm comm = PetscObjectComm((PetscObject)projection->dm); 385 386 PetscCall(MatCreateCeed(projection->dm, projection->dm, op_mass, NULL, &mat_mass)); 387 388 PetscCall(KSPCreate(comm, &projection->ksp)); 389 PetscCall(KSPSetOptionsPrefix(projection->ksp, "div_diff_flux_projection_")); 390 { // lumped by default 391 PC pc; 392 PetscCall(KSPGetPC(projection->ksp, &pc)); 393 PetscCall(PCSetType(pc, PCJACOBI)); 394 PetscCall(PCJacobiSetType(pc, PC_JACOBI_ROWSUM)); 395 PetscCall(KSPSetType(projection->ksp, KSPPREONLY)); 396 } 397 PetscCall(KSPSetFromOptions_WithMatCeed(projection->ksp, mat_mass)); 398 PetscCall(MatDestroy(&mat_mass)); 399 } 400 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_mass)); 401 PetscCallCeed(ceed, CeedOperatorDestroy(&op_mass)); 402 } 403 404 { // Create OperatorApplyContext to calculate divergence at quadrature points 405 CeedQFunction qf_calc_divergence; 406 CeedOperator op_calc_divergence; 407 CeedElemRestriction elem_restr_div_diff_flux; 408 409 { // Get elem_restr_div_diff_flux 410 CeedOperator *sub_ops; 411 CeedOperatorField op_field; 412 PetscInt sub_op_index = 0; // will be 0 for the volume op 413 414 PetscCallCeed(ceed, CeedCompositeOperatorGetSubList(user->op_ifunction, &sub_ops)); 415 PetscCallCeed(ceed, CeedOperatorGetFieldByName(sub_ops[sub_op_index], "div F_diff", &op_field)); 416 PetscCallCeed(ceed, CeedOperatorFieldGetElemRestriction(op_field, &elem_restr_div_diff_flux)); 417 } 418 419 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, ComputeDivDiffusiveFlux3D_4, ComputeDivDiffusiveFlux3D_4_loc, &qf_calc_divergence)); 420 421 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_calc_divergence, "Grad F_diff", projection->num_comp * dim, CEED_EVAL_GRAD)); 422 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_calc_divergence, "qdata", q_data_size, CEED_EVAL_NONE)); 423 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_calc_divergence, "Div F_diff", 4, CEED_EVAL_NONE)); 424 425 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_calc_divergence, NULL, NULL, &op_calc_divergence)); 426 PetscCallCeed(ceed, CeedOperatorSetField(op_calc_divergence, "Grad F_diff", elem_restr_diff_flux, basis_diff_flux, CEED_VECTOR_ACTIVE)); 427 PetscCallCeed(ceed, CeedOperatorSetField(op_calc_divergence, "qdata", elem_restr_qd, CEED_BASIS_NONE, q_data)); 428 PetscCallCeed( 429 ceed, CeedOperatorSetField(op_calc_divergence, "Div F_diff", elem_restr_div_diff_flux, CEED_BASIS_NONE, diff_flux_proj->div_diff_flux_ceed)); 430 431 PetscCall(OperatorApplyContextCreate(projection->dm, NULL, ceed, op_calc_divergence, NULL, NULL, NULL, NULL, 432 &user->diff_flux_proj->calc_div_diff_flux)); 433 434 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_calc_divergence)); 435 PetscCallCeed(ceed, CeedOperatorDestroy(&op_calc_divergence)); 436 } 437 PetscCallCeed(ceed, CeedBasisDestroy(&basis_diff_flux)); 438 PetscCallCeed(ceed, CeedVectorDestroy(&q_data)); 439 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd)); 440 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_diff_flux)); 441 PetscFunctionReturn(PETSC_SUCCESS); 442 } 443 444 /** 445 @brief Setup projection of divergence of diffusive flux 446 447 @param[in] ceed `Ceed` context 448 @param[in,out] user `User` context 449 @param[in] ceed_data `CeedData` context 450 @param[in] problem `ProblemData` context 451 **/ 452 PetscErrorCode DivDiffFluxProjectionSetup(Ceed ceed, User user, CeedData ceed_data, ProblemData problem) { 453 PetscFunctionBeginUser; 454 switch (user->app_ctx->divFdiffproj_method) { 455 case DIV_DIFF_FLUX_PROJ_DIRECT: 456 PetscCall(DivDiffFluxProjectionSetup_Direct(ceed, user, ceed_data, problem)); 457 break; 458 case DIV_DIFF_FLUX_PROJ_INDIRECT: 459 PetscCall(DivDiffFluxProjectionSetup_Indirect(ceed, user, ceed_data, problem)); 460 break; 461 case DIV_DIFF_FLUX_PROJ_NONE: 462 SETERRQ(PetscObjectComm((PetscObject)user->dm), PETSC_ERR_ARG_WRONG, "Should not reach here with div_diff_flux_projection_method %s", 463 DivDiffFluxProjectionMethods[user->app_ctx->divFdiffproj_method]); 464 break; 465 } 466 PetscFunctionReturn(PETSC_SUCCESS); 467 } 468 469 /** 470 @brief Project the divergence of diffusive flux 471 472 This implicitly sets the `CeedVector` input (`div_diff_flux_ceed`) to the divergence of diffusive flux. 473 474 @param[in] diff_flux_proj `NodalProjectionData` for the projection 475 @param[in] Q_loc Localized solution vector 476 **/ 477 PetscErrorCode DiffFluxProjectionApply(DivDiffFluxProjectionData diff_flux_proj, Vec Q_loc) { 478 NodalProjectionData projection = diff_flux_proj->projection; 479 480 PetscFunctionBeginUser; 481 PetscCall(PetscLogEventBegin(FLUIDS_DivDiffFluxProjection, Q_loc, 0, 0, 0)); 482 switch (diff_flux_proj->method) { 483 case DIV_DIFF_FLUX_PROJ_DIRECT: { 484 Vec DivDiffFlux; 485 486 PetscCall(DMGetGlobalVector(projection->dm, &DivDiffFlux)); 487 if (diff_flux_proj->ceed_vec_has_array) { 488 PetscCall(VecReadCeedToPetsc(diff_flux_proj->div_diff_flux_ceed, diff_flux_proj->DivDiffFlux_memtype, diff_flux_proj->DivDiffFlux_loc)); 489 diff_flux_proj->ceed_vec_has_array = PETSC_FALSE; 490 } 491 PetscCall(ApplyCeedOperatorLocalToGlobal(Q_loc, DivDiffFlux, projection->l2_rhs_ctx)); 492 PetscCall(VecViewFromOptions(DivDiffFlux, NULL, "-div_diff_flux_projection_rhs_view")); 493 494 PetscCall(KSPSolve(projection->ksp, DivDiffFlux, DivDiffFlux)); 495 PetscCall(VecViewFromOptions(DivDiffFlux, NULL, "-div_diff_flux_projection_view")); 496 497 PetscCall(DMGlobalToLocal(projection->dm, DivDiffFlux, INSERT_VALUES, diff_flux_proj->DivDiffFlux_loc)); 498 PetscCall(VecReadPetscToCeed(diff_flux_proj->DivDiffFlux_loc, &diff_flux_proj->DivDiffFlux_memtype, diff_flux_proj->div_diff_flux_ceed)); 499 diff_flux_proj->ceed_vec_has_array = PETSC_TRUE; 500 501 PetscCall(DMRestoreGlobalVector(projection->dm, &DivDiffFlux)); 502 break; 503 } 504 case DIV_DIFF_FLUX_PROJ_INDIRECT: { 505 Vec DiffFlux; 506 507 PetscCall(DMGetGlobalVector(projection->dm, &DiffFlux)); 508 PetscCall(ApplyCeedOperatorLocalToGlobal(Q_loc, DiffFlux, projection->l2_rhs_ctx)); 509 PetscCall(VecViewFromOptions(DiffFlux, NULL, "-div_diff_flux_projection_rhs_view")); 510 511 PetscCall(KSPSolve(projection->ksp, DiffFlux, DiffFlux)); 512 PetscCall(VecViewFromOptions(DiffFlux, NULL, "-div_diff_flux_projection_view")); 513 514 PetscCall(ApplyCeedOperatorGlobalToLocal(DiffFlux, NULL, diff_flux_proj->calc_div_diff_flux)); 515 PetscCall(DMRestoreGlobalVector(projection->dm, &DiffFlux)); 516 } break; 517 case DIV_DIFF_FLUX_PROJ_NONE: 518 SETERRQ(PetscObjectComm((PetscObject)projection->dm), PETSC_ERR_ARG_WRONG, "Should not reach here with div_diff_flux_projection_method %s", 519 DivDiffFluxProjectionMethods[diff_flux_proj->method]); 520 break; 521 } 522 PetscCall(PetscLogEventEnd(FLUIDS_DivDiffFluxProjection, Q_loc, 0, 0, 0)); 523 PetscFunctionReturn(PETSC_SUCCESS); 524 } 525 526 /** 527 @brief Destroy `DivDiffFluxProjectionData` object 528 529 @param[in,out] diff_flux_proj Object to destroy 530 **/ 531 PetscErrorCode DivDiffFluxProjectionDataDestroy(DivDiffFluxProjectionData diff_flux_proj) { 532 PetscFunctionBeginUser; 533 if (diff_flux_proj == NULL) PetscFunctionReturn(PETSC_SUCCESS); 534 PetscCall(NodalProjectionDataDestroy(diff_flux_proj->projection)); 535 PetscCall(OperatorApplyContextDestroy(diff_flux_proj->calc_div_diff_flux)); 536 if (diff_flux_proj->ceed_vec_has_array) { 537 PetscCall(VecReadCeedToPetsc(diff_flux_proj->div_diff_flux_ceed, diff_flux_proj->DivDiffFlux_memtype, diff_flux_proj->DivDiffFlux_loc)); 538 diff_flux_proj->ceed_vec_has_array = PETSC_FALSE; 539 } 540 PetscCall(CeedVectorDestroy(&diff_flux_proj->div_diff_flux_ceed)); 541 PetscCall(VecDestroy(&diff_flux_proj->DivDiffFlux_loc)); 542 PetscCall(PetscFree(diff_flux_proj)); 543 PetscFunctionReturn(PETSC_SUCCESS); 544 } 545