1 // Copyright (c) 2017-2023, Lawrence Livermore National Security, LLC and other CEED contributors. 2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 3 // 4 // SPDX-License-Identifier: BSD-2-Clause 5 // 6 // This file is part of CEED: http://github.com/ceed 7 /// @file 8 /// Functions for setting up and projecting the velocity gradient 9 10 #include "../qfunctions/velocity_gradient_projection.h" 11 12 #include <petscdmplex.h> 13 14 #include "../navierstokes.h" 15 16 PetscErrorCode VelocityGradientProjectionCreateDM(NodalProjectionData grad_velo_proj, User user, PetscInt degree) { 17 PetscSection section; 18 19 PetscFunctionBeginUser; 20 grad_velo_proj->num_comp = 9; // 9 velocity gradient 21 22 PetscCall(DMClone(user->dm, &grad_velo_proj->dm)); 23 PetscCall(PetscObjectSetName((PetscObject)grad_velo_proj->dm, "Velocity Gradient Projection")); 24 25 PetscCall( 26 DMSetupByOrder_FEM(PETSC_TRUE, PETSC_TRUE, user->app_ctx->degree, 1, user->app_ctx->q_extra, 1, &grad_velo_proj->num_comp, grad_velo_proj->dm)); 27 28 PetscCall(DMGetLocalSection(grad_velo_proj->dm, §ion)); 29 PetscCall(PetscSectionSetFieldName(section, 0, "")); 30 PetscCall(PetscSectionSetComponentName(section, 0, 0, "VelocityGradientXX")); 31 PetscCall(PetscSectionSetComponentName(section, 0, 1, "VelocityGradientXY")); 32 PetscCall(PetscSectionSetComponentName(section, 0, 2, "VelocityGradientXZ")); 33 PetscCall(PetscSectionSetComponentName(section, 0, 3, "VelocityGradientYX")); 34 PetscCall(PetscSectionSetComponentName(section, 0, 4, "VelocityGradientYY")); 35 PetscCall(PetscSectionSetComponentName(section, 0, 5, "VelocityGradientYZ")); 36 PetscCall(PetscSectionSetComponentName(section, 0, 6, "VelocityGradientZX")); 37 PetscCall(PetscSectionSetComponentName(section, 0, 7, "VelocityGradientZY")); 38 PetscCall(PetscSectionSetComponentName(section, 0, 8, "VelocityGradientZZ")); 39 PetscFunctionReturn(PETSC_SUCCESS); 40 }; 41 42 PetscErrorCode VelocityGradientProjectionSetup(Ceed ceed, User user, CeedData ceed_data, ProblemData *problem, StateVariable state_var_input, 43 CeedElemRestriction elem_restr_input, CeedBasis basis_input, NodalProjectionData *pgrad_velo_proj) { 44 NodalProjectionData grad_velo_proj; 45 CeedOperator op_rhs_assemble, op_mass; 46 CeedQFunction qf_rhs_assemble, qf_mass; 47 CeedBasis basis_grad_velo; 48 CeedElemRestriction elem_restr_grad_velo; 49 PetscInt dim, label_value = 0, height = 0, dm_field = 0; 50 CeedInt num_comp_x, num_comp_input, q_data_size; 51 DMLabel domain_label = NULL; 52 53 PetscFunctionBeginUser; 54 PetscCall(PetscNew(&grad_velo_proj)); 55 56 PetscCall(VelocityGradientProjectionCreateDM(grad_velo_proj, user, user->app_ctx->degree)); 57 58 // -- Get Pre-requisite things 59 PetscCall(DMGetDimension(grad_velo_proj->dm, &dim)); 60 PetscCallCeed(ceed, CeedBasisGetNumComponents(ceed_data->basis_x, &num_comp_x)); 61 PetscCallCeed(ceed, CeedBasisGetNumComponents(ceed_data->basis_q, &num_comp_input)); 62 PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_qd_i, &q_data_size)); 63 PetscCall(DMPlexCeedElemRestrictionCreate(ceed, grad_velo_proj->dm, domain_label, label_value, height, dm_field, &elem_restr_grad_velo)); 64 65 PetscCall(CreateBasisFromPlex(ceed, grad_velo_proj->dm, domain_label, label_value, height, dm_field, &basis_grad_velo)); 66 67 // -- Build RHS operator 68 switch (state_var_input) { 69 case STATEVAR_PRIMITIVE: 70 PetscCallCeed( 71 ceed, CeedQFunctionCreateInterior(ceed, 1, VelocityGradientProjectionRHS_Prim, VelocityGradientProjectionRHS_Prim_loc, &qf_rhs_assemble)); 72 break; 73 case STATEVAR_CONSERVATIVE: 74 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, VelocityGradientProjectionRHS_Conserv, VelocityGradientProjectionRHS_Conserv_loc, 75 &qf_rhs_assemble)); 76 break; 77 default: 78 SETERRQ(PetscObjectComm((PetscObject)user->dm), PETSC_ERR_SUP, "No velocity gradient projection QFunction for chosen state variable"); 79 } 80 81 PetscCallCeed(ceed, CeedQFunctionSetContext(qf_rhs_assemble, problem->apply_vol_ifunction.qfunction_context)); 82 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_rhs_assemble, "q", num_comp_input, CEED_EVAL_INTERP)); 83 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_rhs_assemble, "Grad_q", num_comp_input * dim, CEED_EVAL_GRAD)); 84 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_rhs_assemble, "qdata", q_data_size, CEED_EVAL_NONE)); 85 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_rhs_assemble, "x", num_comp_x, CEED_EVAL_INTERP)); 86 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_rhs_assemble, "velocity gradient", grad_velo_proj->num_comp, CEED_EVAL_INTERP)); 87 88 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_rhs_assemble, NULL, NULL, &op_rhs_assemble)); 89 PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_assemble, "q", elem_restr_input, basis_input, CEED_VECTOR_ACTIVE)); 90 PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_assemble, "Grad_q", elem_restr_input, basis_input, CEED_VECTOR_ACTIVE)); 91 PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_assemble, "qdata", ceed_data->elem_restr_qd_i, CEED_BASIS_NONE, ceed_data->q_data)); 92 PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_assemble, "x", ceed_data->elem_restr_x, ceed_data->basis_x, ceed_data->x_coord)); 93 PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_assemble, "velocity gradient", elem_restr_grad_velo, basis_grad_velo, CEED_VECTOR_ACTIVE)); 94 95 PetscCall(OperatorApplyContextCreate(user->dm, grad_velo_proj->dm, ceed, op_rhs_assemble, NULL, NULL, NULL, NULL, &grad_velo_proj->l2_rhs_ctx)); 96 97 // -- Build Mass operator 98 PetscCall(CreateMassQFunction(ceed, grad_velo_proj->num_comp, q_data_size, &qf_mass)); 99 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_mass, NULL, NULL, &op_mass)); 100 PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "u", elem_restr_grad_velo, basis_grad_velo, CEED_VECTOR_ACTIVE)); 101 PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "qdata", ceed_data->elem_restr_qd_i, CEED_BASIS_NONE, ceed_data->q_data)); 102 PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "v", elem_restr_grad_velo, basis_grad_velo, CEED_VECTOR_ACTIVE)); 103 104 { // -- Setup KSP for L^2 projection with lumped mass operator 105 Mat mat_mass = NULL; 106 OperatorApplyContext mass_matop_ctx; 107 MPI_Comm comm = PetscObjectComm((PetscObject)grad_velo_proj->dm); 108 109 PetscCall(OperatorApplyContextCreate(grad_velo_proj->dm, grad_velo_proj->dm, ceed, op_mass, NULL, NULL, NULL, NULL, &mass_matop_ctx)); 110 PetscCall(CreateMatShell_Ceed(mass_matop_ctx, &mat_mass)); 111 112 PetscCall(KSPCreate(comm, &grad_velo_proj->ksp)); 113 PetscCall(KSPSetOptionsPrefix(grad_velo_proj->ksp, "velocity_gradient_projection_")); 114 { 115 PC pc; 116 PetscCall(KSPGetPC(grad_velo_proj->ksp, &pc)); 117 PetscCall(PCSetType(pc, PCJACOBI)); 118 PetscCall(PCJacobiSetType(pc, PC_JACOBI_ROWSUM)); 119 PetscCall(KSPSetType(grad_velo_proj->ksp, KSPPREONLY)); 120 } 121 PetscCall(KSPSetOperators(grad_velo_proj->ksp, mat_mass, mat_mass)); 122 PetscCall(KSPSetFromOptions(grad_velo_proj->ksp)); 123 PetscCall(MatDestroy(&mat_mass)); 124 } 125 126 *pgrad_velo_proj = grad_velo_proj; 127 128 PetscCallCeed(ceed, CeedBasisDestroy(&basis_grad_velo)); 129 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_grad_velo)); 130 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_rhs_assemble)); 131 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_mass)); 132 PetscCallCeed(ceed, CeedOperatorDestroy(&op_rhs_assemble)); 133 PetscCallCeed(ceed, CeedOperatorDestroy(&op_mass)); 134 PetscFunctionReturn(PETSC_SUCCESS); 135 } 136 137 PetscErrorCode VelocityGradientProjectionApply(NodalProjectionData grad_velo_proj, Vec Q_loc, Vec VelocityGradient) { 138 OperatorApplyContext l2_rhs_ctx = grad_velo_proj->l2_rhs_ctx; 139 140 PetscFunctionBeginUser; 141 PetscCall(PetscLogEventBegin(FLUIDS_VelocityGradientProjection, Q_loc, VelocityGradient, 0, 0)); 142 PetscCall(ApplyCeedOperatorLocalToGlobal(Q_loc, VelocityGradient, l2_rhs_ctx)); 143 144 PetscCall(KSPSolve(grad_velo_proj->ksp, VelocityGradient, VelocityGradient)); 145 PetscCall(PetscLogEventEnd(FLUIDS_VelocityGradientProjection, Q_loc, VelocityGradient, 0, 0)); 146 PetscFunctionReturn(PETSC_SUCCESS); 147 } 148