1 // Copyright (c) 2017-2026, 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
VelocityGradientProjectionCreateDM(NodalProjectionData grad_velo_proj,User user,PetscInt degree)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(DMSetupByOrder_FEM(PETSC_TRUE, PETSC_TRUE, user->app_ctx->degree, 1, user->app_ctx->q_extra, 1, &grad_velo_proj->num_comp,
26 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
VelocityGradientProjectionSetup(Ceed ceed,User user,CeedData ceed_data,ProblemData problem,StateVariable state_var_input,CeedElemRestriction elem_restr_input,CeedBasis basis_input,NodalProjectionData * pgrad_velo_proj)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(ceed, CeedQFunctionCreateInterior(ceed, 1, VelocityGradientProjectionRHS_Prim, VelocityGradientProjectionRHS_Prim_loc,
71 &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 case STATEVAR_ENTROPY:
78 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, VelocityGradientProjectionRHS_Entropy, VelocityGradientProjectionRHS_Entropy_loc,
79 &qf_rhs_assemble));
80 break;
81 }
82
83 PetscCallCeed(ceed, CeedQFunctionSetContext(qf_rhs_assemble, problem->apply_vol_ifunction.qfunction_context));
84 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_rhs_assemble, "q", num_comp_input, CEED_EVAL_INTERP));
85 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_rhs_assemble, "Grad_q", num_comp_input * dim, CEED_EVAL_GRAD));
86 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_rhs_assemble, "qdata", q_data_size, CEED_EVAL_NONE));
87 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_rhs_assemble, "x", num_comp_x, CEED_EVAL_INTERP));
88 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_rhs_assemble, "velocity gradient", grad_velo_proj->num_comp, CEED_EVAL_INTERP));
89
90 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_rhs_assemble, NULL, NULL, &op_rhs_assemble));
91 PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_assemble, "q", elem_restr_input, basis_input, CEED_VECTOR_ACTIVE));
92 PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_assemble, "Grad_q", elem_restr_input, basis_input, CEED_VECTOR_ACTIVE));
93 PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_assemble, "qdata", ceed_data->elem_restr_qd_i, CEED_BASIS_NONE, ceed_data->q_data));
94 PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_assemble, "x", ceed_data->elem_restr_x, ceed_data->basis_x, ceed_data->x_coord));
95 PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_assemble, "velocity gradient", elem_restr_grad_velo, basis_grad_velo, CEED_VECTOR_ACTIVE));
96
97 PetscCall(OperatorApplyContextCreate(user->dm, grad_velo_proj->dm, ceed, op_rhs_assemble, NULL, NULL, NULL, NULL, &grad_velo_proj->l2_rhs_ctx));
98
99 // -- Build Mass operator
100 PetscCall(CreateMassQFunction(ceed, grad_velo_proj->num_comp, q_data_size, &qf_mass));
101 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_mass, NULL, NULL, &op_mass));
102 PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "u", elem_restr_grad_velo, basis_grad_velo, CEED_VECTOR_ACTIVE));
103 PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "qdata", ceed_data->elem_restr_qd_i, CEED_BASIS_NONE, ceed_data->q_data));
104 PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "v", elem_restr_grad_velo, basis_grad_velo, CEED_VECTOR_ACTIVE));
105
106 { // -- Setup KSP for L^2 projection with lumped mass operator
107 Mat mat_mass;
108 MPI_Comm comm = PetscObjectComm((PetscObject)grad_velo_proj->dm);
109
110 PetscCall(MatCreateCeed(grad_velo_proj->dm, grad_velo_proj->dm, op_mass, NULL, &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(KSPSetFromOptions_WithMatCeed(grad_velo_proj->ksp, mat_mass));
122 PetscCall(MatDestroy(&mat_mass));
123 }
124
125 *pgrad_velo_proj = grad_velo_proj;
126
127 PetscCallCeed(ceed, CeedBasisDestroy(&basis_grad_velo));
128 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_grad_velo));
129 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_rhs_assemble));
130 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_mass));
131 PetscCallCeed(ceed, CeedOperatorDestroy(&op_rhs_assemble));
132 PetscCallCeed(ceed, CeedOperatorDestroy(&op_mass));
133 PetscFunctionReturn(PETSC_SUCCESS);
134 }
135
VelocityGradientProjectionApply(NodalProjectionData grad_velo_proj,Vec Q_loc,Vec VelocityGradient)136 PetscErrorCode VelocityGradientProjectionApply(NodalProjectionData grad_velo_proj, Vec Q_loc, Vec VelocityGradient) {
137 OperatorApplyContext l2_rhs_ctx = grad_velo_proj->l2_rhs_ctx;
138
139 PetscFunctionBeginUser;
140 PetscCall(PetscLogEventBegin(FLUIDS_VelocityGradientProjection, Q_loc, VelocityGradient, 0, 0));
141 PetscCall(ApplyCeedOperatorLocalToGlobal(Q_loc, VelocityGradient, l2_rhs_ctx));
142
143 PetscCall(KSPSolve(grad_velo_proj->ksp, VelocityGradient, VelocityGradient));
144 PetscCall(PetscLogEventEnd(FLUIDS_VelocityGradientProjection, Q_loc, VelocityGradient, 0, 0));
145 PetscFunctionReturn(PETSC_SUCCESS);
146 }
147