1 // SPDX-FileCopyrightText: Copyright (c) 2017-2024, HONEE contributors.
2 // SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause
3
4 #include "../qfunctions/strong_boundary_conditions.h"
5
6 #include <ceed.h>
7 #include <petscdmplex.h>
8
9 #include <navierstokes.h>
10 #include "../problems/stg_shur14.h"
11
SetupStrongSTG_Ceed(Ceed ceed,Honee honee,DM dm,ProblemData problem,Physics phys,CeedOperator op_strong_bc)12 PetscErrorCode SetupStrongSTG_Ceed(Ceed ceed, Honee honee, DM dm, ProblemData problem, Physics phys, CeedOperator op_strong_bc) {
13 CeedInt num_comp_x, num_comp_q = 5, stg_data_size = 1, dXdx_size;
14 CeedVector multiplicity, x_stored, scale_stored, stg_data, dXdx, x_coord;
15 CeedBasis basis_x_to_q_face;
16 CeedElemRestriction elem_restr_x_face, elem_restr_q_face, elem_restr_x_stored, elem_restr_scale, elem_restr_stgdata, elem_restr_dXdx;
17 CeedQFunction qf_setup, qf_strongbc, qf_stgdata;
18 CeedOperator op_setup, op_strong_bc_sub, op_stgdata;
19 DMLabel domain_label;
20 const PetscInt dm_field = 0, height_face = 1, height_cell = 0;
21 PetscInt dim;
22 MPI_Comm comm = honee->comm;
23
24 PetscFunctionBeginUser;
25 PetscCall(DMGetLabel(dm, "Face Sets", &domain_label));
26 PetscCall(DMGetDimension(dm, &dim));
27 num_comp_x = dim;
28 dXdx_size = num_comp_x * (dim - height_cell);
29
30 { // Basis
31 CeedBasis basis_x_face, basis_q_face;
32 DM dm_coord;
33
34 PetscCall(DMGetCoordinateDM(dm, &dm_coord));
35 PetscCall(DMPlexCeedBasisCreate(ceed, dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, height_face, dm_field, &basis_q_face));
36 PetscCall(DMPlexCeedBasisCreate(ceed, dm_coord, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, height_face, dm_field, &basis_x_face));
37
38 PetscCallCeed(ceed, CeedBasisCreateProjection(basis_x_face, basis_q_face, &basis_x_to_q_face));
39
40 PetscCallCeed(ceed, CeedBasisDestroy(&basis_q_face));
41 PetscCallCeed(ceed, CeedBasisDestroy(&basis_x_face));
42 }
43
44 // Setup QFunction
45 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, SetupStrongBC, SetupStrongBC_loc, &qf_setup));
46 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup, "x", num_comp_x, CEED_EVAL_INTERP));
47 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup, "dxdX", dXdx_size, CEED_EVAL_GRAD));
48 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup, "multiplicity", num_comp_q, CEED_EVAL_NONE));
49 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_setup, "x stored", num_comp_x, CEED_EVAL_NONE));
50 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_setup, "scale", 1, CEED_EVAL_NONE));
51 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_setup, "dXdx", dXdx_size, CEED_EVAL_NONE));
52
53 // Setup STG QFunctions
54 PetscCall(SetupStrongStg_PreProcessing(ceed, problem, num_comp_x, stg_data_size, dXdx_size, &qf_stgdata));
55 PetscCall(SetupStrongStg_QF(ceed, problem, num_comp_x, num_comp_q, stg_data_size, dXdx_size, &qf_strongbc));
56
57 // Get face label values for the inflow BC
58 PetscInt num_inflow_faces;
59 const PetscInt *inflow_faces;
60 {
61 BCDefinition bc_def = NULL;
62 for (PetscCount b = 0; b < problem->num_bc_defs; b++) {
63 BCDefinition bc_def_ = problem->bc_defs[b];
64 const char *name;
65
66 PetscCall(BCDefinitionGetInfo(bc_def_, &name, NULL, NULL));
67 if (!strcmp(name, "inflow")) {
68 bc_def = bc_def_;
69 break;
70 }
71 }
72 PetscCheck(bc_def, comm, PETSC_ERR_SUP, "Could not find BCDefintion named 'inflow' for STG setup.");
73 PetscCall(BCDefinitionGetInfo(bc_def, NULL, &num_inflow_faces, &inflow_faces));
74 }
75
76 // Compute contribution on each boundary face
77 for (CeedInt i = 0; i < num_inflow_faces; i++) {
78 DMLabel face_orientation_label;
79 PetscInt num_orientations_values, *orientation_values;
80
81 {
82 char *face_orientation_label_name;
83
84 PetscCall(DMPlexCreateFaceLabel(dm, inflow_faces[i], &face_orientation_label_name));
85 PetscCall(DMGetLabel(dm, face_orientation_label_name, &face_orientation_label));
86 PetscCall(PetscFree(face_orientation_label_name));
87 }
88 PetscCall(DMLabelCreateGlobalValueArray(dm, face_orientation_label, &num_orientations_values, &orientation_values));
89 for (PetscInt o = 0; o < num_orientations_values; o++) {
90 CeedBasis basis_x_to_q_cell;
91 CeedElemRestriction elem_restr_x_cell;
92 PetscInt orientation = orientation_values[o];
93
94 {
95 CeedBasis basis_x_cell_to_face, basis_q_face;
96
97 PetscCall(DMPlexCeedBasisCellToFaceCoordinateCreate(ceed, dm, face_orientation_label, orientation, orientation, &basis_x_cell_to_face));
98 PetscCall(DMPlexCeedBasisCreate(ceed, dm, face_orientation_label, orientation, height_face, dm_field, &basis_q_face));
99 PetscCallCeed(ceed, CeedBasisCreateProjection(basis_x_cell_to_face, basis_q_face, &basis_x_to_q_cell));
100 PetscCallCeed(ceed, CeedBasisDestroy(&basis_x_cell_to_face));
101 PetscCallCeed(ceed, CeedBasisDestroy(&basis_q_face));
102 }
103
104 PetscCall(DMPlexCeedElemRestrictionCreate(ceed, dm, face_orientation_label, orientation, height_face, dm_field, &elem_restr_q_face));
105 PetscCall(DMPlexCeedElemRestrictionCoordinateCreate(ceed, dm, face_orientation_label, orientation, height_face, &elem_restr_x_face));
106 PetscCall(DMPlexCeedCoordinateCreateField(ceed, dm, face_orientation_label, orientation, height_cell, &elem_restr_x_cell, NULL, &x_coord));
107 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_q_face, &multiplicity, NULL));
108 PetscCallCeed(ceed, CeedElemRestrictionGetMultiplicity(elem_restr_q_face, multiplicity));
109
110 PetscCall(DMPlexCeedElemRestrictionCollocatedCreate(ceed, dm, face_orientation_label, orientation, height_face, num_comp_x,
111 &elem_restr_x_stored));
112 PetscCall(DMPlexCeedElemRestrictionCollocatedCreate(ceed, dm, face_orientation_label, orientation, height_face, 1, &elem_restr_scale));
113 PetscCall(DMPlexCeedElemRestrictionCollocatedCreate(ceed, dm, face_orientation_label, orientation, height_face, stg_data_size,
114 &elem_restr_stgdata));
115 PetscCall(DMPlexCeedElemRestrictionCollocatedCreate(ceed, dm, face_orientation_label, orientation, height_face, dXdx_size, &elem_restr_dXdx));
116 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_x_stored, &x_stored, NULL));
117 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_scale, &scale_stored, NULL));
118 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_stgdata, &stg_data, NULL));
119 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_dXdx, &dXdx, NULL));
120
121 // -- Setup Operator
122 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_setup, NULL, NULL, &op_setup));
123 PetscCallCeed(ceed, CeedOperatorSetName(op_setup, "Precomputed data for strong boundary conditions"));
124 PetscCallCeed(ceed, CeedOperatorSetField(op_setup, "x", elem_restr_x_face, basis_x_to_q_face, CEED_VECTOR_ACTIVE));
125 PetscCallCeed(ceed, CeedOperatorSetField(op_setup, "dxdX", elem_restr_x_cell, basis_x_to_q_cell, CEED_VECTOR_ACTIVE));
126 PetscCallCeed(ceed, CeedOperatorSetField(op_setup, "multiplicity", elem_restr_q_face, CEED_BASIS_NONE, multiplicity));
127 PetscCallCeed(ceed, CeedOperatorSetField(op_setup, "x stored", elem_restr_x_stored, CEED_BASIS_NONE, x_stored));
128 PetscCallCeed(ceed, CeedOperatorSetField(op_setup, "scale", elem_restr_scale, CEED_BASIS_NONE, scale_stored));
129 PetscCallCeed(ceed, CeedOperatorSetField(op_setup, "dXdx", elem_restr_dXdx, CEED_BASIS_NONE, dXdx));
130
131 // -- Compute geometric factors
132 PetscCallCeed(ceed, CeedOperatorApply(op_setup, x_coord, CEED_VECTOR_NONE, CEED_REQUEST_IMMEDIATE));
133
134 // -- Compute STGData
135 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_stgdata, NULL, NULL, &op_stgdata));
136 PetscCallCeed(ceed, CeedOperatorSetField(op_stgdata, "dXdx", elem_restr_dXdx, CEED_BASIS_NONE, dXdx));
137 PetscCallCeed(ceed, CeedOperatorSetField(op_stgdata, "x", elem_restr_x_stored, CEED_BASIS_NONE, x_stored));
138 PetscCallCeed(ceed, CeedOperatorSetField(op_stgdata, "stg data", elem_restr_stgdata, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE));
139
140 PetscCallCeed(ceed, CeedOperatorApply(op_stgdata, CEED_VECTOR_NONE, stg_data, CEED_REQUEST_IMMEDIATE));
141
142 // -- Setup BC QFunctions
143 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_strongbc, NULL, NULL, &op_strong_bc_sub));
144 PetscCallCeed(ceed, CeedOperatorSetName(op_strong_bc_sub, "Strong STG"));
145
146 PetscCallCeed(ceed, CeedOperatorSetField(op_strong_bc_sub, "dXdx", elem_restr_dXdx, CEED_BASIS_NONE, dXdx));
147 PetscCallCeed(ceed, CeedOperatorSetField(op_strong_bc_sub, "x", elem_restr_x_stored, CEED_BASIS_NONE, x_stored));
148 PetscCallCeed(ceed, CeedOperatorSetField(op_strong_bc_sub, "scale", elem_restr_scale, CEED_BASIS_NONE, scale_stored));
149 PetscCallCeed(ceed, CeedOperatorSetField(op_strong_bc_sub, "stg data", elem_restr_stgdata, CEED_BASIS_NONE, stg_data));
150 PetscCallCeed(ceed, CeedOperatorSetField(op_strong_bc_sub, "q", elem_restr_q_face, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE));
151
152 // -- Add to composite operator
153 PetscCallCeed(ceed, CeedOperatorCompositeAddSub(op_strong_bc, op_strong_bc_sub));
154
155 PetscCallCeed(ceed, CeedVectorDestroy(&multiplicity));
156 PetscCallCeed(ceed, CeedVectorDestroy(&x_stored));
157 PetscCallCeed(ceed, CeedVectorDestroy(&scale_stored));
158 PetscCallCeed(ceed, CeedVectorDestroy(&stg_data));
159 PetscCallCeed(ceed, CeedVectorDestroy(&dXdx));
160 PetscCallCeed(ceed, CeedVectorDestroy(&x_coord));
161 PetscCallCeed(ceed, CeedBasisDestroy(&basis_x_to_q_cell));
162 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_x_face));
163 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_x_cell));
164 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_q_face));
165 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_x_stored));
166 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_scale));
167 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_stgdata));
168 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_dXdx));
169 PetscCallCeed(ceed, CeedOperatorDestroy(&op_strong_bc_sub));
170 PetscCallCeed(ceed, CeedOperatorDestroy(&op_setup));
171 PetscCallCeed(ceed, CeedOperatorDestroy(&op_stgdata));
172 }
173 PetscCall(PetscFree(orientation_values));
174 }
175
176 PetscCallCeed(ceed, CeedOperatorGetContextFieldLabel(op_strong_bc, "solution time", &phys->stg_solution_time_label));
177
178 PetscCallCeed(ceed, CeedBasisDestroy(&basis_x_to_q_face));
179 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_strongbc));
180 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_stgdata));
181 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_setup));
182 PetscFunctionReturn(PETSC_SUCCESS);
183 }
184
DMPlexInsertBoundaryValues_StrongBCCeed(DM dm,PetscBool insert_essential,Vec Q_loc,PetscReal time,Vec face_geom_FVM,Vec cell_geom_FVM,Vec grad_FVM)185 PetscErrorCode DMPlexInsertBoundaryValues_StrongBCCeed(DM dm, PetscBool insert_essential, Vec Q_loc, PetscReal time, Vec face_geom_FVM,
186 Vec cell_geom_FVM, Vec grad_FVM) {
187 Vec boundary_mask;
188 Honee honee;
189
190 PetscFunctionBeginUser;
191 PetscCall(PetscLogEventBegin(HONEE_StrongBCInsert, dm, Q_loc, 0, 0));
192 PetscCall(DMGetApplicationContext(dm, &honee));
193
194 if (honee->phys->stg_solution_time_label) {
195 PetscCallCeed(honee->ceed, CeedOperatorSetContextDouble(honee->op_strong_bc_ctx->op, honee->phys->stg_solution_time_label, &time));
196 }
197
198 // Mask Strong BC entries
199 PetscCall(DMGetNamedLocalVector(dm, "boundary mask", &boundary_mask));
200 PetscCall(VecPointwiseMult(Q_loc, Q_loc, boundary_mask));
201 PetscCall(DMRestoreNamedLocalVector(dm, "boundary mask", &boundary_mask));
202
203 PetscCall(PetscLogEventBegin(HONEE_StrongBCCeed, dm, Q_loc, 0, 0));
204 PetscCall(ApplyAddCeedOperatorLocalToLocal(NULL, Q_loc, honee->op_strong_bc_ctx));
205 PetscCall(PetscLogEventEnd(HONEE_StrongBCCeed, dm, Q_loc, 0, 0));
206 PetscCall(PetscLogEventEnd(HONEE_StrongBCInsert, dm, Q_loc, 0, 0));
207 PetscFunctionReturn(PETSC_SUCCESS);
208 }
209
SetupStrongBC_Ceed(Ceed ceed,DM dm,Honee honee,ProblemData problem)210 PetscErrorCode SetupStrongBC_Ceed(Ceed ceed, DM dm, Honee honee, ProblemData problem) {
211 CeedOperator op_strong_bc;
212
213 PetscFunctionBeginUser;
214 {
215 Vec boundary_mask, global_vec;
216
217 PetscCall(DMGetNamedLocalVector(dm, "boundary mask", &boundary_mask));
218 PetscCall(DMGetGlobalVector(dm, &global_vec));
219 PetscCall(VecZeroEntries(boundary_mask));
220 PetscCall(VecSet(global_vec, 1.0));
221 PetscCall(DMGlobalToLocal(dm, global_vec, INSERT_VALUES, boundary_mask));
222 PetscCall(DMRestoreNamedLocalVector(dm, "boundary mask", &boundary_mask));
223 PetscCall(DMRestoreGlobalVector(dm, &global_vec));
224 }
225
226 PetscCallCeed(ceed, CeedOperatorCreateComposite(ceed, &op_strong_bc));
227 {
228 PetscBool use_strongstg = PETSC_FALSE;
229 PetscCall(PetscOptionsGetBool(NULL, NULL, "-stg_strong", &use_strongstg, NULL));
230
231 if (use_strongstg) {
232 PetscCall(SetupStrongSTG_Ceed(ceed, honee, dm, problem, honee->phys, op_strong_bc));
233 }
234 }
235
236 PetscCall(OperatorApplyContextCreate(NULL, NULL, ceed, op_strong_bc, CEED_VECTOR_NONE, NULL, NULL, NULL, &honee->op_strong_bc_ctx));
237
238 PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMPlexInsertBoundaryValues_C", DMPlexInsertBoundaryValues_StrongBCCeed));
239 PetscCallCeed(ceed, CeedOperatorDestroy(&op_strong_bc));
240 PetscFunctionReturn(PETSC_SUCCESS);
241 }
242