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
VelocityGradientProjectionCreateDM(NodalProjectionData grad_velo_proj,Honee honee,PetscInt degree)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
VelocityGradientProjectionSetup(Ceed ceed,Honee honee,ProblemData problem,StateVariable state_var_input,CeedElemRestriction elem_restr_input,CeedBasis basis_input,NodalProjectionData * pgrad_velo_proj)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
VelocityGradientProjectionApply(NodalProjectionData grad_velo_proj,Vec Q_loc,Vec VelocityGradient)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