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 PetscFE fe; 18 PetscSection section; 19 PetscInt dim; 20 21 PetscFunctionBeginUser; 22 grad_velo_proj->num_comp = 9; // 9 velocity gradient 23 24 PetscCall(DMClone(user->dm, &grad_velo_proj->dm)); 25 PetscCall(DMGetDimension(grad_velo_proj->dm, &dim)); 26 PetscCall(PetscObjectSetName((PetscObject)grad_velo_proj->dm, "Velocity Gradient Projection")); 27 28 PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, dim, grad_velo_proj->num_comp, PETSC_FALSE, degree, PETSC_DECIDE, &fe)); 29 PetscCall(PetscObjectSetName((PetscObject)fe, "Velocity Gradient Projection")); 30 PetscCall(DMAddField(grad_velo_proj->dm, NULL, (PetscObject)fe)); 31 PetscCall(DMCreateDS(grad_velo_proj->dm)); 32 PetscCall(DMPlexSetClosurePermutationTensor(grad_velo_proj->dm, PETSC_DETERMINE, NULL)); 33 34 PetscCall(DMGetLocalSection(grad_velo_proj->dm, §ion)); 35 PetscCall(PetscSectionSetFieldName(section, 0, "")); 36 PetscCall(PetscSectionSetComponentName(section, 0, 0, "VelocityGradientXX")); 37 PetscCall(PetscSectionSetComponentName(section, 0, 1, "VelocityGradientXY")); 38 PetscCall(PetscSectionSetComponentName(section, 0, 2, "VelocityGradientXZ")); 39 PetscCall(PetscSectionSetComponentName(section, 0, 3, "VelocityGradientYX")); 40 PetscCall(PetscSectionSetComponentName(section, 0, 4, "VelocityGradientYY")); 41 PetscCall(PetscSectionSetComponentName(section, 0, 5, "VelocityGradientYZ")); 42 PetscCall(PetscSectionSetComponentName(section, 0, 6, "VelocityGradientZX")); 43 PetscCall(PetscSectionSetComponentName(section, 0, 7, "VelocityGradientZY")); 44 PetscCall(PetscSectionSetComponentName(section, 0, 8, "VelocityGradientZZ")); 45 46 PetscCall(PetscFEDestroy(&fe)); 47 PetscFunctionReturn(0); 48 }; 49 50 PetscErrorCode VelocityGradientProjectionSetup(Ceed ceed, User user, CeedData ceed_data, ProblemData *problem) { 51 OperatorApplyContext mass_matop_ctx; 52 CeedOperator op_rhs_assemble, op_mass; 53 CeedQFunction qf_rhs_assemble, qf_mass; 54 CeedBasis basis_grad_velo; 55 CeedVector q_ceed, rhs_ceed, mass_output; 56 CeedElemRestriction elem_restr_grad_velo; 57 PetscInt dim, num_comp_x, num_comp_q, q_data_size, num_qpts_1d, num_nodes_1d; 58 59 PetscFunctionBeginUser; 60 PetscCall(PetscNew(&user->grad_velo_proj)); 61 NodalProjectionData grad_velo_proj = user->grad_velo_proj; 62 63 PetscCall(VelocityGradientProjectionCreateDM(grad_velo_proj, user, user->app_ctx->degree)); 64 65 // -- Get Pre-requisite things 66 PetscCall(DMGetDimension(grad_velo_proj->dm, &dim)); 67 CeedBasisGetNumComponents(ceed_data->basis_x, &num_comp_x); 68 CeedBasisGetNumComponents(ceed_data->basis_q, &num_comp_q); 69 CeedBasisGetNumQuadraturePoints1D(ceed_data->basis_q, &num_qpts_1d); 70 CeedBasisGetNumNodes1D(ceed_data->basis_q, &num_nodes_1d); 71 CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_qd_i, &q_data_size); 72 73 PetscCall(GetRestrictionForDomain(ceed, grad_velo_proj->dm, 0, 0, 0, num_qpts_1d, q_data_size, &elem_restr_grad_velo, NULL, NULL)); 74 75 CeedBasisCreateTensorH1Lagrange(ceed, dim, grad_velo_proj->num_comp, num_nodes_1d, num_qpts_1d, CEED_GAUSS, &basis_grad_velo); 76 77 // -- Build RHS operator 78 switch (user->phys->state_var) { 79 case STATEVAR_PRIMITIVE: 80 CeedQFunctionCreateInterior(ceed, 1, VelocityGradientProjectionRHS_Prim, VelocityGradientProjectionRHS_Prim_loc, &qf_rhs_assemble); 81 break; 82 case STATEVAR_CONSERVATIVE: 83 CeedQFunctionCreateInterior(ceed, 1, VelocityGradientProjectionRHS_Conserv, VelocityGradientProjectionRHS_Conserv_loc, &qf_rhs_assemble); 84 break; 85 default: 86 SETERRQ(PetscObjectComm((PetscObject)user->dm), PETSC_ERR_SUP, "No velocity gradient projection QFunction for chosen state variable"); 87 } 88 89 CeedQFunctionSetContext(qf_rhs_assemble, problem->apply_vol_ifunction.qfunction_context); 90 CeedQFunctionAddInput(qf_rhs_assemble, "q", num_comp_q, CEED_EVAL_INTERP); 91 CeedQFunctionAddInput(qf_rhs_assemble, "Grad_q", num_comp_q * dim, CEED_EVAL_GRAD); 92 CeedQFunctionAddInput(qf_rhs_assemble, "qdata", q_data_size, CEED_EVAL_NONE); 93 CeedQFunctionAddInput(qf_rhs_assemble, "x", num_comp_x, CEED_EVAL_INTERP); 94 CeedQFunctionAddOutput(qf_rhs_assemble, "velocity gradient", grad_velo_proj->num_comp, CEED_EVAL_INTERP); 95 96 CeedOperatorCreate(ceed, qf_rhs_assemble, NULL, NULL, &op_rhs_assemble); 97 CeedOperatorSetField(op_rhs_assemble, "q", ceed_data->elem_restr_q, ceed_data->basis_q, CEED_VECTOR_ACTIVE); 98 CeedOperatorSetField(op_rhs_assemble, "Grad_q", ceed_data->elem_restr_q, ceed_data->basis_q, CEED_VECTOR_ACTIVE); 99 CeedOperatorSetField(op_rhs_assemble, "qdata", ceed_data->elem_restr_qd_i, CEED_BASIS_COLLOCATED, ceed_data->q_data); 100 CeedOperatorSetField(op_rhs_assemble, "x", ceed_data->elem_restr_x, ceed_data->basis_x, ceed_data->x_coord); 101 CeedOperatorSetField(op_rhs_assemble, "velocity gradient", elem_restr_grad_velo, basis_grad_velo, CEED_VECTOR_ACTIVE); 102 103 CeedElemRestrictionCreateVector(ceed_data->elem_restr_q, &q_ceed, NULL); 104 CeedElemRestrictionCreateVector(elem_restr_grad_velo, &rhs_ceed, NULL); 105 106 PetscCall( 107 OperatorApplyContextCreate(user->dm, grad_velo_proj->dm, ceed, op_rhs_assemble, q_ceed, rhs_ceed, NULL, NULL, &grad_velo_proj->l2_rhs_ctx)); 108 109 // -- Build Mass operator 110 PetscCall(CreateMassQFunction(ceed, grad_velo_proj->num_comp, q_data_size, &qf_mass)); 111 CeedOperatorCreate(ceed, qf_mass, NULL, NULL, &op_mass); 112 CeedOperatorSetField(op_mass, "u", elem_restr_grad_velo, basis_grad_velo, CEED_VECTOR_ACTIVE); 113 CeedOperatorSetField(op_mass, "qdata", ceed_data->elem_restr_qd_i, CEED_BASIS_COLLOCATED, ceed_data->q_data); 114 CeedOperatorSetField(op_mass, "v", elem_restr_grad_velo, basis_grad_velo, CEED_VECTOR_ACTIVE); 115 116 { // -- Setup KSP for L^2 projection with lumped mass operator 117 PetscInt l_size, g_size; 118 Mat mat_mass; 119 VecType vec_type; 120 Vec M_inv; 121 MPI_Comm comm = PetscObjectComm((PetscObject)grad_velo_proj->dm); 122 123 PetscCall(DMGetGlobalVector(grad_velo_proj->dm, &M_inv)); 124 PetscCall(VecGetLocalSize(M_inv, &l_size)); 125 PetscCall(VecGetSize(M_inv, &g_size)); 126 PetscCall(VecGetType(M_inv, &vec_type)); 127 PetscCall(DMRestoreGlobalVector(grad_velo_proj->dm, &M_inv)); 128 129 CeedElemRestrictionCreateVector(elem_restr_grad_velo, &mass_output, NULL); 130 PetscCall(OperatorApplyContextCreate(grad_velo_proj->dm, grad_velo_proj->dm, ceed, op_mass, rhs_ceed, mass_output, NULL, NULL, &mass_matop_ctx)); 131 CeedVectorDestroy(&mass_output); 132 133 PetscCall(MatCreateShell(comm, l_size, l_size, g_size, g_size, mass_matop_ctx, &mat_mass)); 134 PetscCall(MatShellSetContextDestroy(mat_mass, (PetscErrorCode(*)(void *))OperatorApplyContextDestroy)); 135 PetscCall(MatShellSetOperation(mat_mass, MATOP_MULT, (void (*)(void))MatMult_Ceed)); 136 PetscCall(MatShellSetOperation(mat_mass, MATOP_GET_DIAGONAL, (void (*)(void))MatGetDiag_Ceed)); 137 PetscCall(MatShellSetVecType(mat_mass, vec_type)); 138 139 PetscCall(KSPCreate(comm, &grad_velo_proj->ksp)); 140 PetscCall(KSPSetOptionsPrefix(grad_velo_proj->ksp, "velocity_gradient_projection_")); 141 { 142 PC pc; 143 PetscCall(KSPGetPC(grad_velo_proj->ksp, &pc)); 144 PetscCall(PCSetType(pc, PCJACOBI)); 145 PetscCall(PCJacobiSetType(pc, PC_JACOBI_ROWSUM)); 146 PetscCall(KSPSetType(grad_velo_proj->ksp, KSPPREONLY)); 147 // TODO Not sure if the option below are necessary 148 PetscCall(KSPSetConvergenceTest(grad_velo_proj->ksp, KSPConvergedSkip, NULL, NULL)); 149 } 150 PetscCall(KSPSetOperators(grad_velo_proj->ksp, mat_mass, mat_mass)); 151 PetscCall(KSPSetFromOptions(grad_velo_proj->ksp)); 152 } 153 154 CeedVectorDestroy(&q_ceed); 155 CeedVectorDestroy(&rhs_ceed); 156 CeedBasisDestroy(&basis_grad_velo); 157 CeedElemRestrictionDestroy(&elem_restr_grad_velo); 158 CeedQFunctionDestroy(&qf_rhs_assemble); 159 CeedQFunctionDestroy(&qf_mass); 160 CeedOperatorDestroy(&op_rhs_assemble); 161 CeedOperatorDestroy(&op_mass); 162 PetscFunctionReturn(0); 163 } 164 165 PetscErrorCode VelocityGradientProjectionApply(User user, Vec Q_loc, Vec VelocityGradient) { 166 NodalProjectionData grad_velo_proj = user->grad_velo_proj; 167 OperatorApplyContext l2_rhs_ctx = grad_velo_proj->l2_rhs_ctx; 168 169 PetscFunctionBeginUser; 170 PetscCall(ApplyCeedOperatorLocalToGlobal(Q_loc, VelocityGradient, l2_rhs_ctx)); 171 172 PetscCall(KSPSolve(grad_velo_proj->ksp, VelocityGradient, VelocityGradient)); 173 174 PetscFunctionReturn(0); 175 } 176