1 // SPDX-FileCopyrightText: Copyright (c) 2017-2024, HONEE contributors.
2 // SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause
3
4 #include <navierstokes.h>
5 #include "petscerror.h"
6
7 /**
8 @brief Add `BCDefinition` to a `PetscSegBuffer`
9
10 @param[in] bc_def `BCDefinition` to add
11 @param[in,out] bc_defs_seg `PetscSegBuffer` to add to
12 **/
AddBCDefinitionToSegBuffer(BCDefinition bc_def,PetscSegBuffer bc_defs_seg)13 static PetscErrorCode AddBCDefinitionToSegBuffer(BCDefinition bc_def, PetscSegBuffer bc_defs_seg) {
14 BCDefinition *bc_def_ptr;
15
16 PetscFunctionBeginUser;
17 if (bc_def == NULL) PetscFunctionReturn(PETSC_SUCCESS);
18 PetscCall(PetscSegBufferGet(bc_defs_seg, 1, &bc_def_ptr));
19 *bc_def_ptr = bc_def;
20 PetscFunctionReturn(PETSC_SUCCESS);
21 }
22
23 /**
24 @brief Create and setup `BCDefinition`s from commandline options
25
26 @param[in] honee `Honee`
27 @param[in,out] problem `ProblemData`
28 @param[in] app_ctx `AppCtx`
29 **/
BoundaryConditionSetUp(Honee honee,ProblemData problem,AppCtx app_ctx)30 PetscErrorCode BoundaryConditionSetUp(Honee honee, ProblemData problem, AppCtx app_ctx) {
31 PetscSegBuffer bc_defs_seg;
32 PetscBool flg;
33 BCDefinition bc_def;
34
35 PetscFunctionBeginUser;
36 PetscCall(PetscSegBufferCreate(sizeof(BCDefinition), 4, &bc_defs_seg));
37
38 PetscOptionsBegin(honee->comm, NULL, "Boundary Condition Options", NULL);
39
40 PetscCall(PetscOptionsBCDefinition("-bc_wall", "Face IDs to apply wall BC", NULL, "wall", &bc_def, NULL));
41 PetscCall(AddBCDefinitionToSegBuffer(bc_def, bc_defs_seg));
42 if (bc_def) {
43 PetscInt num_essential_comps = 16, essential_comps[16];
44
45 PetscCall(PetscOptionsIntArray("-wall_comps", "An array of constrained component numbers", NULL, essential_comps, &num_essential_comps, &flg));
46 PetscCall(BCDefinitionSetEssential(bc_def, num_essential_comps, essential_comps));
47
48 app_ctx->wall_forces.num_wall = bc_def->num_label_values;
49 PetscCall(PetscMalloc1(bc_def->num_label_values, &app_ctx->wall_forces.walls));
50 PetscCall(PetscArraycpy(app_ctx->wall_forces.walls, bc_def->label_values, bc_def->num_label_values));
51 }
52
53 { // Symmetry Boundary Conditions
54 const char *deprecated[3] = {"-bc_slip_x", "-bc_slip_y", "-bc_slip_z"};
55 const char *flags[3] = {"-bc_symmetry_x", "-bc_symmetry_y", "-bc_symmetry_z"};
56
57 for (PetscInt j = 0; j < 3; j++) {
58 PetscCall(PetscOptionsDeprecated(deprecated[j], flags[j], "libCEED 0.12.0",
59 "Use -bc_symmetry_[x,y,z] for direct equivalency, or -bc_slip for weak, Riemann-based, direction-invariant "
60 "slip/no-penatration boundary conditions"));
61 PetscCall(PetscOptionsBCDefinition(flags[j], "Face IDs to apply symmetry BC", NULL, "symmetry", &bc_def, NULL));
62 if (!bc_def) {
63 PetscCall(PetscOptionsBCDefinition(deprecated[j], "Face IDs to apply symmetry BC", NULL, "symmetry", &bc_def, NULL));
64 }
65 PetscCall(AddBCDefinitionToSegBuffer(bc_def, bc_defs_seg));
66 if (bc_def) {
67 PetscInt essential_comps[1] = {j + 1};
68
69 PetscCall(BCDefinitionSetEssential(bc_def, 1, essential_comps));
70 }
71 }
72 }
73
74 PetscCall(PetscOptionsBCDefinition("-bc_inflow", "Face IDs to apply inflow BC", NULL, "inflow", &bc_def, NULL));
75 PetscCall(AddBCDefinitionToSegBuffer(bc_def, bc_defs_seg));
76
77 PetscCall(PetscOptionsBCDefinition("-bc_outflow", "Face IDs to apply outflow BC", NULL, "outflow", &bc_def, NULL));
78 PetscCall(AddBCDefinitionToSegBuffer(bc_def, bc_defs_seg));
79
80 PetscCall(PetscOptionsBCDefinition("-bc_freestream", "Face IDs to apply freestream BC", NULL, "freestream", &bc_def, NULL));
81 PetscCall(AddBCDefinitionToSegBuffer(bc_def, bc_defs_seg));
82
83 PetscCall(PetscOptionsBCDefinition("-bc_slip", "Face IDs to apply slip BC", NULL, "slip", &bc_def, NULL));
84 PetscCall(AddBCDefinitionToSegBuffer(bc_def, bc_defs_seg));
85
86 PetscOptionsEnd();
87
88 PetscCall(PetscSegBufferGetSize(bc_defs_seg, &problem->num_bc_defs));
89 PetscCall(PetscSegBufferExtractAlloc(bc_defs_seg, &problem->bc_defs));
90 PetscCall(PetscSegBufferDestroy(&bc_defs_seg));
91
92 //TODO: Verify that the BCDefinition don't have overlapping claims to boundary faces
93 PetscFunctionReturn(PETSC_SUCCESS);
94 }
95
96 /**
97 @brief Destroy `HoneeBCStruct` object
98
99 @param[in] ctx Pointer to `HoneeBCStruct`
100 **/
HoneeBCDestroy(void ** ctx)101 PetscErrorCode HoneeBCDestroy(void **ctx) {
102 HoneeBCStruct honee_bc = *(HoneeBCStruct *)ctx;
103 Ceed ceed = honee_bc->honee->ceed;
104
105 PetscFunctionBeginUser;
106 if (honee_bc->qfctx) PetscCallCeed(ceed, CeedQFunctionContextDestroy(&honee_bc->qfctx));
107 if (honee_bc->DestroyCtx) PetscCall((*honee_bc->DestroyCtx)(&honee_bc->ctx));
108 PetscCall(PetscFree(honee_bc));
109 *ctx = NULL;
110 PetscFunctionReturn(PETSC_SUCCESS);
111 }
112
113 /**
114 @brief Create a QFunction matching the "standard" HONEE inputs for IFunctions
115
116 This assumes the `bc_def` context is a `HoneeBCStruct` and inputs/outputs of the IFunction QFunction are in the correct order.
117
118 @param[in] bc_def `BCDefinition` that hosts the IFunction
119 @param[in] qf_func_ptr `CeedQFunctionUser` for the IFunction
120 @param[in] qf_loc Absolute path to source of `CeedQFunctionUser`
121 @param[in] qfctx `CeedQFunctionContext` for the IFunction (also shared with the IJacobian, if applicable)
122 @param[out] qf_ifunc QFunction for the IFunction
123 **/
HoneeBCCreateIFunctionQF(BCDefinition bc_def,CeedQFunctionUser qf_func_ptr,const char * qf_loc,CeedQFunctionContext qfctx,CeedQFunction * qf_ifunc)124 PetscErrorCode HoneeBCCreateIFunctionQF(BCDefinition bc_def, CeedQFunctionUser qf_func_ptr, const char *qf_loc, CeedQFunctionContext qfctx,
125 CeedQFunction *qf_ifunc) {
126 Ceed ceed;
127 DM dm;
128 PetscInt dim, dim_sur, height = 1, num_comp_x, num_comp_q;
129 CeedInt q_data_size_sur, num_comps_jac_data;
130 HoneeBCStruct honee_bc;
131
132 PetscFunctionBeginUser;
133 PetscCall(BCDefinitionGetDM(bc_def, &dm));
134 PetscCall(BCDefinitionGetContext(bc_def, &honee_bc));
135 ceed = honee_bc->honee->ceed;
136 num_comps_jac_data = honee_bc->num_comps_jac_data;
137
138 PetscCall(DMGetDimension(dm, &dim));
139 dim_sur = dim - height;
140 PetscCall(QDataBoundaryGetNumComponents(dm, &q_data_size_sur));
141 PetscCall(DMGetCoordinateNumComps(dm, &num_comp_x));
142 {
143 PetscSection section;
144
145 PetscCall(DMGetLocalSection(dm, §ion));
146 PetscCall(PetscSectionGetFieldComponents(section, bc_def->dm_field, &num_comp_q));
147 }
148
149 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, qf_func_ptr, qf_loc, qf_ifunc));
150 PetscCallCeed(ceed, CeedQFunctionSetContext(*qf_ifunc, qfctx));
151 PetscCallCeed(ceed, CeedQFunctionSetUserFlopsEstimate(*qf_ifunc, 0));
152 PetscCallCeed(ceed, CeedQFunctionAddInput(*qf_ifunc, "q", num_comp_q, CEED_EVAL_INTERP));
153 PetscCallCeed(ceed, CeedQFunctionAddInput(*qf_ifunc, "Grad_q", num_comp_q * dim_sur, CEED_EVAL_GRAD));
154 PetscCallCeed(ceed, CeedQFunctionAddInput(*qf_ifunc, "surface qdata", q_data_size_sur, CEED_EVAL_NONE));
155 PetscCallCeed(ceed, CeedQFunctionAddInput(*qf_ifunc, "x", num_comp_x, CEED_EVAL_INTERP));
156 PetscCallCeed(ceed, CeedQFunctionAddOutput(*qf_ifunc, "v", num_comp_q, CEED_EVAL_INTERP));
157 if (num_comps_jac_data) PetscCallCeed(ceed, CeedQFunctionAddOutput(*qf_ifunc, "surface jacobian data", num_comps_jac_data, CEED_EVAL_NONE));
158 PetscFunctionReturn(PETSC_SUCCESS);
159 }
160
161 /**
162 @brief Create a QFunction matching the "standard" HONEE inputs for IJacobian
163
164 This assumes the `bc_def` context is a `HoneeBCStruct` and inputs/outputs of the IJacobian QFunction are in the correct order.
165
166 @param[in] bc_def `BCDefinition` that hosts the IJacobian
167 @param[in] qf_func_ptr `CeedQFunctionUser` for the IJacobian
168 @param[in] qf_loc Absolute path to source of `CeedQFunctionUser`
169 @param[in] qfctx `CeedQFunctionContext` for the IJacobian (also shared with the IFunction)
170 @param[out] qf_ijac QFunction for the IJacobian
171 **/
HoneeBCCreateIJacobianQF(BCDefinition bc_def,CeedQFunctionUser qf_func_ptr,const char * qf_loc,CeedQFunctionContext qfctx,CeedQFunction * qf_ijac)172 PetscErrorCode HoneeBCCreateIJacobianQF(BCDefinition bc_def, CeedQFunctionUser qf_func_ptr, const char *qf_loc, CeedQFunctionContext qfctx,
173 CeedQFunction *qf_ijac) {
174 Ceed ceed;
175 DM dm;
176 PetscInt dim, dim_sur, height = 1, num_comp_x, num_comp_q;
177 CeedInt q_data_size_sur, num_comps_jac_data;
178 HoneeBCStruct honee_bc;
179
180 PetscFunctionBeginUser;
181 PetscCall(BCDefinitionGetDM(bc_def, &dm));
182 PetscCall(BCDefinitionGetContext(bc_def, &honee_bc));
183 ceed = honee_bc->honee->ceed;
184 num_comps_jac_data = honee_bc->num_comps_jac_data;
185
186 PetscCall(DMGetDimension(dm, &dim));
187 dim_sur = dim - height;
188 PetscCall(QDataBoundaryGetNumComponents(dm, &q_data_size_sur));
189 PetscCall(DMGetCoordinateNumComps(dm, &num_comp_x));
190 {
191 PetscSection section;
192
193 PetscCall(DMGetLocalSection(dm, §ion));
194 PetscCall(PetscSectionGetFieldComponents(section, bc_def->dm_field, &num_comp_q));
195 }
196
197 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, qf_func_ptr, qf_loc, qf_ijac));
198 PetscCallCeed(ceed, CeedQFunctionSetContext(*qf_ijac, qfctx));
199 PetscCallCeed(ceed, CeedQFunctionSetUserFlopsEstimate(*qf_ijac, 0));
200 PetscCallCeed(ceed, CeedQFunctionAddInput(*qf_ijac, "dq", num_comp_q, CEED_EVAL_INTERP));
201 PetscCallCeed(ceed, CeedQFunctionAddInput(*qf_ijac, "Grad_dq", num_comp_q * dim_sur, CEED_EVAL_GRAD));
202 PetscCallCeed(ceed, CeedQFunctionAddInput(*qf_ijac, "surface qdata", q_data_size_sur, CEED_EVAL_NONE));
203 PetscCallCeed(ceed, CeedQFunctionAddInput(*qf_ijac, "x", num_comp_x, CEED_EVAL_INTERP));
204 if (num_comps_jac_data) PetscCallCeed(ceed, CeedQFunctionAddInput(*qf_ijac, "surface jacobian data", num_comps_jac_data, CEED_EVAL_NONE));
205 PetscCallCeed(ceed, CeedQFunctionAddOutput(*qf_ijac, "v", num_comp_q, CEED_EVAL_INTERP));
206 PetscFunctionReturn(PETSC_SUCCESS);
207 }
208
209 /**
210 @brief Setups and adds IFunction operator to given composite operator
211
212 @param[in] bc_def `BCDefinition` for the boundary condition IFunction
213 @param[in] domain_label `DMLabel` for the face orientation label
214 @param[in] label_value Orientation value
215 @param[in] qf_ifunc QFunction for the IFunction
216 @param[in,out] op_ifunc Composite operator to be added to
217 @param[out] sub_op_ifunc Non-composite operator that is added in this function. To be passed to `HoneeBCAddIJacobianOp()`
218 **/
HoneeBCAddIFunctionOp(BCDefinition bc_def,DMLabel domain_label,PetscInt label_value,CeedQFunction qf_ifunc,CeedOperator op_ifunc,CeedOperator * sub_op_ifunc)219 PetscErrorCode HoneeBCAddIFunctionOp(BCDefinition bc_def, DMLabel domain_label, PetscInt label_value, CeedQFunction qf_ifunc, CeedOperator op_ifunc,
220 CeedOperator *sub_op_ifunc) {
221 Ceed ceed;
222 DM dm;
223 HoneeBCStruct honee_bc;
224 PetscInt dim, height = 1, num_comp_x, num_comp_q;
225 CeedInt q_data_size, num_comps_jac_data;
226 CeedVector q_data, jac_data, x_coord;
227 CeedOperator sub_op_ifunc_;
228 CeedElemRestriction elem_restr_qd, elem_restr_q, elem_restr_x, elem_restr_jac;
229 CeedBasis basis_x, basis_q;
230
231 PetscFunctionBeginUser;
232 PetscCall(BCDefinitionGetDM(bc_def, &dm));
233 PetscCall(BCDefinitionGetContext(bc_def, &honee_bc));
234 ceed = honee_bc->honee->ceed;
235 num_comps_jac_data = honee_bc->num_comps_jac_data;
236
237 PetscCall(DMGetDimension(dm, &dim));
238 PetscCall(DMGetCoordinateNumComps(dm, &num_comp_x));
239 PetscCall(DMGetFieldNumComps(dm, bc_def->dm_field, &num_comp_q));
240
241 PetscCall(DMPlexCeedBasisCreate(ceed, dm, domain_label, label_value, height, bc_def->dm_field, &basis_q));
242 PetscCall(DMPlexCeedElemRestrictionCreate(ceed, dm, domain_label, label_value, height, bc_def->dm_field, &elem_restr_q));
243 PetscCall(DMPlexCeedCoordinateCreateField(ceed, dm, domain_label, label_value, height, &elem_restr_x, &basis_x, &x_coord));
244 PetscCall(QDataBoundaryGet(ceed, dm, domain_label, label_value, &elem_restr_qd, &q_data, &q_data_size));
245 if (num_comps_jac_data > 0) {
246 PetscCall(DMPlexCeedElemRestrictionQDataCreate(ceed, dm, domain_label, label_value, height, num_comps_jac_data, &elem_restr_jac));
247 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_jac, &jac_data, NULL));
248 }
249
250 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_ifunc, NULL, NULL, &sub_op_ifunc_));
251 PetscCallCeed(ceed, CeedOperatorSetField(sub_op_ifunc_, "q", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE));
252 PetscCallCeed(ceed, CeedOperatorSetField(sub_op_ifunc_, "Grad_q", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE));
253 PetscCallCeed(ceed, CeedOperatorSetField(sub_op_ifunc_, "surface qdata", elem_restr_qd, CEED_BASIS_NONE, q_data));
254 PetscCallCeed(ceed, CeedOperatorSetField(sub_op_ifunc_, "x", elem_restr_x, basis_x, x_coord));
255 PetscCallCeed(ceed, CeedOperatorSetField(sub_op_ifunc_, "v", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE));
256 if (num_comps_jac_data > 0)
257 PetscCallCeed(ceed, CeedOperatorSetField(sub_op_ifunc_, "surface jacobian data", elem_restr_jac, CEED_BASIS_NONE, jac_data));
258
259 PetscCallCeed(ceed, CeedOperatorCompositeAddSub(op_ifunc, sub_op_ifunc_));
260 *sub_op_ifunc = sub_op_ifunc_;
261
262 PetscCallCeed(ceed, CeedVectorDestroy(&q_data));
263 PetscCallCeed(ceed, CeedVectorDestroy(&x_coord));
264 PetscCallCeed(ceed, CeedBasisDestroy(&basis_q));
265 PetscCallCeed(ceed, CeedBasisDestroy(&basis_x));
266 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd));
267 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_q));
268 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_x));
269 if (num_comps_jac_data > 0) {
270 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_jac));
271 PetscCallCeed(ceed, CeedVectorDestroy(&jac_data));
272 }
273 PetscFunctionReturn(PETSC_SUCCESS);
274 }
275
276 /**
277 @brief Setups and adds IJacobian operator to given composite operator
278
279 The field data (element restriction, vector, etc) is read from `sub_op_ifunc` and used for the IJacobian.
280
281 @param[in] bc_def `BCDefinition` for the boundary condition IJacobian
282 @param[out] sub_op_ifunc IFunction operator corresponding to this IJacobian operator.
283 @param[in] domain_label `DMLabel` for the face orientation label
284 @param[in] label_value Orientation value
285 @param[in] qf_ijac QFunction for the IJacobian
286 @param[in,out] op_ijac Composite operator to be added to
287 **/
HoneeBCAddIJacobianOp(BCDefinition bc_def,CeedOperator sub_op_ifunc,DMLabel domain_label,PetscInt label_value,CeedQFunction qf_ijac,CeedOperator op_ijac)288 PetscErrorCode HoneeBCAddIJacobianOp(BCDefinition bc_def, CeedOperator sub_op_ifunc, DMLabel domain_label, PetscInt label_value,
289 CeedQFunction qf_ijac, CeedOperator op_ijac) {
290 Ceed ceed;
291 DM dm;
292 HoneeBCStruct honee_bc;
293 PetscInt dim, height = 1, num_comp_x, num_comp_q;
294 CeedInt q_data_size;
295 CeedVector q_data, jac_data, x_coord;
296 CeedElemRestriction elem_restr_qd, elem_restr_q, elem_restr_x, elem_restr_jac;
297 CeedOperator sub_op_ijac;
298 CeedBasis basis_x, basis_q;
299
300 PetscFunctionBeginUser;
301 PetscCall(BCDefinitionGetContext(bc_def, &honee_bc));
302 PetscBool use_jac_data = honee_bc->num_comps_jac_data > 0;
303 PetscCall(BCDefinitionGetDM(bc_def, &dm));
304 ceed = honee_bc->honee->ceed;
305
306 PetscCall(DMGetDimension(dm, &dim));
307 PetscCall(DMGetCoordinateNumComps(dm, &num_comp_x));
308 PetscCall(DMGetFieldNumComps(dm, bc_def->dm_field, &num_comp_q));
309
310 PetscCall(DMPlexCeedCoordinateCreateField(ceed, dm, domain_label, label_value, height, &elem_restr_x, &basis_x, &x_coord));
311 PetscCall(QDataBoundaryGet(ceed, dm, domain_label, label_value, &elem_restr_qd, &q_data, &q_data_size));
312
313 { // Get restriction and basis from the IFunction function
314 CeedOperatorField op_field;
315
316 PetscCallCeed(ceed, CeedOperatorGetFieldByName(sub_op_ifunc, "q", &op_field));
317 PetscCallCeed(ceed, CeedOperatorFieldGetData(op_field, NULL, &elem_restr_q, &basis_q, NULL));
318 if (use_jac_data) {
319 PetscCallCeed(ceed, CeedOperatorGetFieldByName(sub_op_ifunc, "surface jacobian data", &op_field));
320 PetscCallCeed(ceed, CeedOperatorFieldGetData(op_field, NULL, &elem_restr_jac, NULL, &jac_data));
321 }
322 }
323
324 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_ijac, NULL, NULL, &sub_op_ijac));
325 PetscCallCeed(ceed, CeedOperatorSetField(sub_op_ijac, "dq", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE));
326 PetscCallCeed(ceed, CeedOperatorSetField(sub_op_ijac, "Grad_dq", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE));
327 PetscCallCeed(ceed, CeedOperatorSetField(sub_op_ijac, "surface qdata", elem_restr_qd, CEED_BASIS_NONE, q_data));
328 PetscCallCeed(ceed, CeedOperatorSetField(sub_op_ijac, "x", elem_restr_x, basis_x, x_coord));
329 if (use_jac_data) PetscCallCeed(ceed, CeedOperatorSetField(sub_op_ijac, "surface jacobian data", elem_restr_jac, CEED_BASIS_NONE, jac_data));
330 PetscCallCeed(ceed, CeedOperatorSetField(sub_op_ijac, "v", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE));
331
332 PetscCallCeed(ceed, CeedOperatorCompositeAddSub(op_ijac, sub_op_ijac));
333
334 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_q));
335 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_x));
336 if (use_jac_data) {
337 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_jac));
338 PetscCallCeed(ceed, CeedVectorDestroy(&jac_data));
339 }
340 PetscCallCeed(ceed, CeedVectorDestroy(&q_data));
341 PetscCallCeed(ceed, CeedVectorDestroy(&x_coord));
342 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd));
343 PetscCallCeed(ceed, CeedBasisDestroy(&basis_q));
344 PetscCallCeed(ceed, CeedBasisDestroy(&basis_x));
345 PetscCallCeed(ceed, CeedOperatorDestroy(&sub_op_ijac));
346 PetscFunctionReturn(PETSC_SUCCESS);
347 }
348