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 velocity gradient 5 6 #include "../qfunctions/velocity_gradient_projection.h" 7 8 #include <petscdmplex.h> 9 10 #include <navierstokes.h> 11 12 static PetscErrorCode VelocityGradientProjectionCreateDM(NodalProjectionData grad_velo_proj, Honee honee, PetscInt degree) { 13 PetscSection section; 14 15 PetscFunctionBeginUser; 16 grad_velo_proj->num_comp = 9; // 9 velocity gradient 17 18 PetscCall(DMClone(honee->dm, &grad_velo_proj->dm)); 19 PetscCall(DMSetMatrixPreallocateSkip(grad_velo_proj->dm, PETSC_TRUE)); 20 PetscCall(PetscObjectSetName((PetscObject)grad_velo_proj->dm, "Velocity Gradient Projection")); 21 22 PetscCall(DMSetupByOrder_FEM(PETSC_TRUE, PETSC_TRUE, honee->app_ctx->degree, 1, honee->app_ctx->q_extra, 1, &grad_velo_proj->num_comp, 23 grad_velo_proj->dm)); 24 25 PetscCall(DMGetLocalSection(grad_velo_proj->dm, §ion)); 26 PetscCall(PetscSectionSetFieldName(section, 0, "")); 27 PetscCall(PetscSectionSetComponentName(section, 0, 0, "VelocityGradientXX")); 28 PetscCall(PetscSectionSetComponentName(section, 0, 1, "VelocityGradientXY")); 29 PetscCall(PetscSectionSetComponentName(section, 0, 2, "VelocityGradientXZ")); 30 PetscCall(PetscSectionSetComponentName(section, 0, 3, "VelocityGradientYX")); 31 PetscCall(PetscSectionSetComponentName(section, 0, 4, "VelocityGradientYY")); 32 PetscCall(PetscSectionSetComponentName(section, 0, 5, "VelocityGradientYZ")); 33 PetscCall(PetscSectionSetComponentName(section, 0, 6, "VelocityGradientZX")); 34 PetscCall(PetscSectionSetComponentName(section, 0, 7, "VelocityGradientZY")); 35 PetscCall(PetscSectionSetComponentName(section, 0, 8, "VelocityGradientZZ")); 36 PetscFunctionReturn(PETSC_SUCCESS); 37 }; 38 39 PetscErrorCode VelocityGradientProjectionSetup(Ceed ceed, Honee honee, ProblemData problem, StateVariable state_var_input, 40 CeedElemRestriction elem_restr_input, CeedBasis basis_input, NodalProjectionData *pgrad_velo_proj) { 41 NodalProjectionData grad_velo_proj; 42 CeedBasis basis_grad_velo; 43 CeedElemRestriction elem_restr_grad_velo, elem_restr_qd; 44 CeedVector q_data; 45 PetscInt dim, height = 0, dm_field = 0, num_comp_input; 46 CeedInt q_data_size; 47 48 PetscFunctionBeginUser; 49 PetscCall(PetscNew(&grad_velo_proj)); 50 PetscCall(VelocityGradientProjectionCreateDM(grad_velo_proj, honee, honee->app_ctx->degree)); 51 52 // -- Get Pre-requisite things 53 PetscCall(DMGetDimension(grad_velo_proj->dm, &dim)); 54 PetscCall(DMGetFieldNumComps(honee->dm, 0, &num_comp_input)); 55 PetscCall(DMPlexCeedElemRestrictionCreate(ceed, grad_velo_proj->dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, height, dm_field, 56 &elem_restr_grad_velo)); 57 PetscCall(DMPlexCeedBasisCreate(ceed, grad_velo_proj->dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, height, dm_field, &basis_grad_velo)); 58 PetscCall(QDataGet(ceed, grad_velo_proj->dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, &elem_restr_qd, &q_data, &q_data_size)); 59 60 { // -- Build RHS operator 61 CeedOperator op_rhs_assemble; 62 CeedQFunction qf_rhs_assemble; 63 CeedElemRestriction elem_restr_x; 64 CeedBasis basis_x; 65 CeedVector x_coord; 66 PetscInt num_comp_x; 67 68 PetscCall(DMGetCoordinateNumComps(honee->dm, &num_comp_x)); 69 PetscCall(DMPlexCeedCoordinateCreateField(ceed, honee->dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, height, &elem_restr_x, &basis_x, &x_coord)); 70 switch (state_var_input) { 71 case STATEVAR_PRIMITIVE: 72 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, VelocityGradientProjectionRHS_Prim, VelocityGradientProjectionRHS_Prim_loc, 73 &qf_rhs_assemble)); 74 break; 75 case STATEVAR_CONSERVATIVE: 76 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, VelocityGradientProjectionRHS_Conserv, VelocityGradientProjectionRHS_Conserv_loc, 77 &qf_rhs_assemble)); 78 break; 79 case STATEVAR_ENTROPY: 80 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, VelocityGradientProjectionRHS_Entropy, VelocityGradientProjectionRHS_Entropy_loc, 81 &qf_rhs_assemble)); 82 break; 83 } 84 85 PetscCallCeed(ceed, CeedQFunctionSetContext(qf_rhs_assemble, problem->apply_vol_ifunction.qfctx)); 86 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_rhs_assemble, "q", num_comp_input, CEED_EVAL_INTERP)); 87 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_rhs_assemble, "Grad_q", num_comp_input * dim, CEED_EVAL_GRAD)); 88 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_rhs_assemble, "qdata", q_data_size, CEED_EVAL_NONE)); 89 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_rhs_assemble, "x", num_comp_x, CEED_EVAL_INTERP)); 90 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_rhs_assemble, "velocity gradient", grad_velo_proj->num_comp, CEED_EVAL_INTERP)); 91 92 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_rhs_assemble, NULL, NULL, &op_rhs_assemble)); 93 PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_assemble, "q", elem_restr_input, basis_input, CEED_VECTOR_ACTIVE)); 94 PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_assemble, "Grad_q", elem_restr_input, basis_input, CEED_VECTOR_ACTIVE)); 95 PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_assemble, "qdata", elem_restr_qd, CEED_BASIS_NONE, q_data)); 96 PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_assemble, "x", elem_restr_x, basis_x, x_coord)); 97 PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_assemble, "velocity gradient", elem_restr_grad_velo, basis_grad_velo, CEED_VECTOR_ACTIVE)); 98 99 PetscCall(OperatorApplyContextCreate(honee->dm, grad_velo_proj->dm, ceed, op_rhs_assemble, NULL, NULL, NULL, NULL, &grad_velo_proj->l2_rhs_ctx)); 100 101 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_x)); 102 PetscCallCeed(ceed, CeedBasisDestroy(&basis_x)); 103 PetscCallCeed(ceed, CeedVectorDestroy(&x_coord)); 104 PetscCallCeed(ceed, CeedOperatorDestroy(&op_rhs_assemble)); 105 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_rhs_assemble)); 106 } 107 108 { // -- Build Mass operator 109 CeedOperator op_mass; 110 CeedQFunction qf_mass; 111 Mat mat_mass; 112 113 PetscCall(HoneeMassQFunctionCreate(ceed, grad_velo_proj->num_comp, q_data_size, &qf_mass)); 114 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_mass, NULL, NULL, &op_mass)); 115 PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "u", elem_restr_grad_velo, basis_grad_velo, CEED_VECTOR_ACTIVE)); 116 PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "qdata", elem_restr_qd, CEED_BASIS_NONE, q_data)); 117 PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "v", elem_restr_grad_velo, basis_grad_velo, CEED_VECTOR_ACTIVE)); 118 119 PetscCall(MatCreateCeed(grad_velo_proj->dm, grad_velo_proj->dm, op_mass, NULL, &mat_mass)); 120 121 // Setup KSP for L^2 projection with lumped mass operator 122 PetscCall(KSPCreate(PetscObjectComm((PetscObject)grad_velo_proj->dm), &grad_velo_proj->ksp)); 123 PetscCall(KSPSetOptionsPrefix(grad_velo_proj->ksp, "velocity_gradient_projection_")); 124 { 125 PC pc; 126 PetscCall(KSPGetPC(grad_velo_proj->ksp, &pc)); 127 PetscCall(PCSetType(pc, PCJACOBI)); 128 PetscCall(PCJacobiSetType(pc, PC_JACOBI_ROWSUM)); 129 PetscCall(KSPSetType(grad_velo_proj->ksp, KSPPREONLY)); 130 } 131 PetscCall(KSPSetFromOptions_WithMatCeed(grad_velo_proj->ksp, mat_mass)); 132 133 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_mass)); 134 PetscCallCeed(ceed, CeedOperatorDestroy(&op_mass)); 135 PetscCall(MatDestroy(&mat_mass)); 136 } 137 138 *pgrad_velo_proj = grad_velo_proj; 139 140 PetscCallCeed(ceed, CeedBasisDestroy(&basis_grad_velo)); 141 PetscCallCeed(ceed, CeedVectorDestroy(&q_data)); 142 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd)); 143 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_grad_velo)); 144 PetscFunctionReturn(PETSC_SUCCESS); 145 } 146 147 PetscErrorCode VelocityGradientProjectionApply(NodalProjectionData grad_velo_proj, Vec Q_loc, Vec VelocityGradient) { 148 OperatorApplyContext l2_rhs_ctx = grad_velo_proj->l2_rhs_ctx; 149 150 PetscFunctionBeginUser; 151 PetscCall(PetscLogEventBegin(HONEE_VelocityGradientProjection, Q_loc, VelocityGradient, 0, 0)); 152 PetscCall(ApplyCeedOperatorLocalToGlobal(Q_loc, VelocityGradient, l2_rhs_ctx)); 153 154 PetscCall(KSPSolve(grad_velo_proj->ksp, VelocityGradient, VelocityGradient)); 155 PetscCall(PetscLogEventEnd(HONEE_VelocityGradientProjection, Q_loc, VelocityGradient, 0, 0)); 156 PetscFunctionReturn(PETSC_SUCCESS); 157 } 158