xref: /honee/src/boundary_condition.c (revision d510ce3ce7d7a59b0a32036570b9db0675ee3035) 
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 **/
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 **/
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 **/
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 **/
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, &section));
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 **/
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, &section));
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 **/
219 PetscErrorCode HoneeBCAddIFunctionOp(BCDefinition bc_def, DMLabel domain_label, PetscInt label_value, CeedQFunction qf_ifunc, CeedOperator op_ifunc,
220                                      CeedOperator *sub_op_ifunc) {
221   Honee               honee;
222   Ceed                ceed;
223   DM                  dm;
224   HoneeBCStruct       honee_bc;
225   PetscInt            dim, height = 1, num_comp_x, num_comp_q;
226   CeedInt             q_data_size, num_comps_jac_data;
227   CeedVector          q_data, jac_data;
228   CeedOperator        sub_op_ifunc_;
229   CeedElemRestriction elem_restr_qd, elem_restr_q, elem_restr_x, elem_restr_jac;
230   CeedBasis           basis_x, basis_q;
231 
232   PetscFunctionBeginUser;
233   PetscCall(BCDefinitionGetDM(bc_def, &dm));
234   PetscCall(BCDefinitionGetContext(bc_def, &honee_bc));
235   honee              = honee_bc->honee;
236   ceed               = honee_bc->honee->ceed;
237   num_comps_jac_data = honee_bc->num_comps_jac_data;
238 
239   PetscCall(DMGetDimension(dm, &dim));
240   PetscCall(DMGetCoordinateNumComps(dm, &num_comp_x));
241   PetscCall(DMGetFieldNumComps(dm, bc_def->dm_field, &num_comp_q));
242 
243   PetscCall(DMPlexCeedBasisCreate(ceed, dm, domain_label, label_value, height, bc_def->dm_field, &basis_q));
244   PetscCall(DMPlexCeedBasisCoordinateCreate(ceed, dm, domain_label, label_value, height, &basis_x));
245   PetscCall(DMPlexCeedElemRestrictionCreate(ceed, dm, domain_label, label_value, height, bc_def->dm_field, &elem_restr_q));
246   PetscCall(DMPlexCeedElemRestrictionCoordinateCreate(ceed, dm, domain_label, label_value, height, &elem_restr_x));
247   PetscCall(QDataBoundaryGet(ceed, dm, domain_label, label_value, &elem_restr_qd, &q_data, &q_data_size));
248   if (num_comps_jac_data > 0) {
249     PetscCall(DMPlexCeedElemRestrictionQDataCreate(ceed, dm, domain_label, label_value, height, num_comps_jac_data, &elem_restr_jac));
250     PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_jac, &jac_data, NULL));
251   }
252 
253   PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_ifunc, NULL, NULL, &sub_op_ifunc_));
254   PetscCallCeed(ceed, CeedOperatorSetField(sub_op_ifunc_, "q", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE));
255   PetscCallCeed(ceed, CeedOperatorSetField(sub_op_ifunc_, "Grad_q", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE));
256   PetscCallCeed(ceed, CeedOperatorSetField(sub_op_ifunc_, "surface qdata", elem_restr_qd, CEED_BASIS_NONE, q_data));
257   PetscCallCeed(ceed, CeedOperatorSetField(sub_op_ifunc_, "x", elem_restr_x, basis_x, honee->x_coord));
258   PetscCallCeed(ceed, CeedOperatorSetField(sub_op_ifunc_, "v", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE));
259   if (num_comps_jac_data > 0)
260     PetscCallCeed(ceed, CeedOperatorSetField(sub_op_ifunc_, "surface jacobian data", elem_restr_jac, CEED_BASIS_NONE, jac_data));
261 
262   PetscCallCeed(ceed, CeedOperatorCompositeAddSub(op_ifunc, sub_op_ifunc_));
263   *sub_op_ifunc = sub_op_ifunc_;
264 
265   PetscCallCeed(ceed, CeedVectorDestroy(&q_data));
266   PetscCallCeed(ceed, CeedBasisDestroy(&basis_q));
267   PetscCallCeed(ceed, CeedBasisDestroy(&basis_x));
268   PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd));
269   PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_q));
270   PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_x));
271   if (num_comps_jac_data > 0) {
272     PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_jac));
273     PetscCallCeed(ceed, CeedVectorDestroy(&jac_data));
274   }
275   PetscFunctionReturn(PETSC_SUCCESS);
276 }
277 
278 /**
279   @brief Setups and adds IJacobian operator to given composite operator
280 
281   The field data (element restriction, vector, etc) is read from `sub_op_ifunc` and used for the IJacobian.
282 
283   @param[in]     bc_def       `BCDefinition` for the boundary condition IJacobian
284   @param[out]    sub_op_ifunc IFunction operator corresponding to this IJacobian operator.
285   @param[in]     domain_label `DMLabel` for the face orientation label
286   @param[in]     label_value  Orientation value
287   @param[in]     qf_ijac      QFunction for the IJacobian
288   @param[in,out] op_ijac      Composite operator to be added to
289 **/
290 PetscErrorCode HoneeBCAddIJacobianOp(BCDefinition bc_def, CeedOperator sub_op_ifunc, DMLabel domain_label, PetscInt label_value,
291                                      CeedQFunction qf_ijac, CeedOperator op_ijac) {
292   Honee               honee;
293   Ceed                ceed;
294   DM                  dm;
295   HoneeBCStruct       honee_bc;
296   PetscInt            dim, height = 1, num_comp_x, num_comp_q;
297   CeedInt             q_data_size;
298   CeedVector          q_data, jac_data;
299   CeedElemRestriction elem_restr_qd, elem_restr_q, elem_restr_x, elem_restr_jac;
300   CeedOperator        sub_op_ijac;
301   CeedBasis           basis_x, basis_q;
302 
303   PetscFunctionBeginUser;
304   PetscCall(BCDefinitionGetContext(bc_def, &honee_bc));
305   PetscBool use_jac_data = honee_bc->num_comps_jac_data > 0;
306   PetscCall(BCDefinitionGetDM(bc_def, &dm));
307   honee = honee_bc->honee;
308   ceed  = honee_bc->honee->ceed;
309 
310   PetscCall(DMGetDimension(dm, &dim));
311   PetscCall(DMGetCoordinateNumComps(dm, &num_comp_x));
312   PetscCall(DMGetFieldNumComps(dm, bc_def->dm_field, &num_comp_q));
313 
314   PetscCall(DMPlexCeedBasisCreate(ceed, dm, domain_label, label_value, height, bc_def->dm_field, &basis_q));
315   PetscCall(DMPlexCeedBasisCoordinateCreate(ceed, dm, domain_label, label_value, height, &basis_x));
316   PetscCall(QDataBoundaryGet(ceed, dm, domain_label, label_value, &elem_restr_qd, &q_data, &q_data_size));
317 
318   {  // Get restriction and basis from the IFunction function
319     CeedOperatorField op_field;
320 
321     PetscCallCeed(ceed, CeedOperatorGetFieldByName(sub_op_ifunc, "q", &op_field));
322     PetscCallCeed(ceed, CeedOperatorFieldGetElemRestriction(op_field, &elem_restr_q));
323 
324     PetscCallCeed(ceed, CeedOperatorGetFieldByName(sub_op_ifunc, "x", &op_field));
325     PetscCallCeed(ceed, CeedOperatorFieldGetElemRestriction(op_field, &elem_restr_x));
326 
327     if (use_jac_data) {
328       PetscCallCeed(ceed, CeedOperatorGetFieldByName(sub_op_ifunc, "surface jacobian data", &op_field));
329       PetscCallCeed(ceed, CeedOperatorFieldGetData(op_field, NULL, &elem_restr_jac, NULL, &jac_data));
330     }
331   }
332 
333   PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_ijac, NULL, NULL, &sub_op_ijac));
334   PetscCallCeed(ceed, CeedOperatorSetField(sub_op_ijac, "dq", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE));
335   PetscCallCeed(ceed, CeedOperatorSetField(sub_op_ijac, "Grad_dq", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE));
336   PetscCallCeed(ceed, CeedOperatorSetField(sub_op_ijac, "surface qdata", elem_restr_qd, CEED_BASIS_NONE, q_data));
337   PetscCallCeed(ceed, CeedOperatorSetField(sub_op_ijac, "x", elem_restr_x, basis_x, honee->x_coord));
338   if (use_jac_data) PetscCallCeed(ceed, CeedOperatorSetField(sub_op_ijac, "surface jacobian data", elem_restr_jac, CEED_BASIS_NONE, jac_data));
339   PetscCallCeed(ceed, CeedOperatorSetField(sub_op_ijac, "v", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE));
340 
341   PetscCallCeed(ceed, CeedOperatorCompositeAddSub(op_ijac, sub_op_ijac));
342 
343   PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_q));
344   PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_x));
345   if (use_jac_data) {
346     PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_jac));
347     PetscCallCeed(ceed, CeedVectorDestroy(&jac_data));
348   }
349   PetscCallCeed(ceed, CeedVectorDestroy(&q_data));
350   PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd));
351   PetscCallCeed(ceed, CeedBasisDestroy(&basis_q));
352   PetscCallCeed(ceed, CeedBasisDestroy(&basis_x));
353   PetscCallCeed(ceed, CeedOperatorDestroy(&sub_op_ijac));
354   PetscFunctionReturn(PETSC_SUCCESS);
355 }
356