PICurv 0.1.0
A Parallel Particle-In-Cell Solver for Curvilinear LES
Loading...
Searching...
No Matches
BC_Handlers.c
Go to the documentation of this file.
1#include "BC_Handlers.h" // The header that declares this file's "constructor" functions
2
3
4//================================================================================
5// VALIDATORS
6//================================================================================
7
8
9#undef __FUNCT__
10#define __FUNCT__ "Validate_DrivenFlowConfiguration"
11/**
12 * @brief Internal helper implementation: `Validate_DrivenFlowConfiguration()`.
13 * @details Local to this translation unit.
14 */
15PetscErrorCode Validate_DrivenFlowConfiguration(UserCtx *user)
16{
17 PetscFunctionBeginUser;
18
19 // --- CHECK 1: Detect if a driven flow is active. ---
20 PetscBool is_driven_flow_active = PETSC_FALSE;
21 char driven_direction = ' ';
22 const char* first_driven_face_name = "";
23
24 for (int i = 0; i < 6; i++) {
25 BCHandlerType handler_type = user->boundary_faces[i].handler_type;
26 if (handler_type == BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX ||
27 handler_type == BC_HANDLER_PERIODIC_DRIVEN_INITIAL_FLUX)
28 {
29 is_driven_flow_active = PETSC_TRUE;
30 first_driven_face_name = BCFaceToString((BCFace)i);
31
32 if (i <= 1) driven_direction = 'X';
33 else if (i <= 3) driven_direction = 'Y';
34 else driven_direction = 'Z';
35
36 break; // Exit loop once we've confirmed it's active and found the direction.
37 }
38 }
39
40 // If no driven flow handler is found, validation for this rule set is complete.
41 if (!is_driven_flow_active) {
42 PetscFunctionReturn(0);
43 }
44
45 LOG_ALLOW(GLOBAL, LOG_DEBUG, " - Driven Flow Handler detected on face %s. Applying driven flow validation rules...\n", first_driven_face_name);
46
47 // --- CHECK 2: Ensure no conflicting BCs (Inlet/Outlet/Far-field) are present. ---
48 LOG_ALLOW(GLOBAL, LOG_DEBUG, " - Checking for incompatible Inlet/Outlet/Far-field BCs...\n");
49 for (int i = 0; i < 6; i++) {
50 BCType math_type = user->boundary_faces[i].mathematical_type;
51 if (math_type == INLET || math_type == OUTLET || math_type == FARFIELD) {
52 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT,
53 "Configuration Error: A DRIVEN flow handler is active, which is incompatible with the %s boundary condition found on face %s.",
54 BCTypeToString(math_type), BCFaceToString((BCFace)i));
55 }
56 }
57 LOG_ALLOW(GLOBAL, LOG_DEBUG, " ... No conflicting BC types found. OK.\n");
58
59 // --- CHECK 3: Ensure both ends of the driven direction have identical, valid setups. ---
60 LOG_ALLOW(GLOBAL, LOG_DEBUG, " - Validating symmetry and mathematical types for the '%c' direction...\n", driven_direction);
61
62 PetscInt neg_face_idx = 0, pos_face_idx = 0;
63 if (driven_direction == 'X') {
64 neg_face_idx = BC_FACE_NEG_X; pos_face_idx = BC_FACE_POS_X;
65 } else if (driven_direction == 'Y') {
66 neg_face_idx = BC_FACE_NEG_Y; pos_face_idx = BC_FACE_POS_Y;
67 } else { // 'Z'
68 neg_face_idx = BC_FACE_NEG_Z; pos_face_idx = BC_FACE_POS_Z;
69 }
70
71 BoundaryFaceConfig *neg_face_cfg = &user->boundary_faces[neg_face_idx];
72 BoundaryFaceConfig *pos_face_cfg = &user->boundary_faces[pos_face_idx];
73
74 // Rule 3a: Both faces must be PERIODIC.
75 if (neg_face_cfg->mathematical_type != PERIODIC || pos_face_cfg->mathematical_type != PERIODIC) {
76 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT,
77 "Configuration Error: For a driven flow in the '%c' direction, both the %s and %s faces must be of mathematical_type PERIODIC.",
78 driven_direction, BCFaceToString((BCFace)neg_face_idx), BCFaceToString((BCFace)pos_face_idx));
79 }
80
81 // Rule 3b: Both faces must use the exact same handler type.
82 if (neg_face_cfg->handler_type != pos_face_cfg->handler_type) {
83 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT,
84 "Configuration Error: The DRIVEN handlers on the %s and %s faces of the '%c' direction do not match. Both must be the same type (e.g., both CONSTANT_FLUX).",
85 BCFaceToString((BCFace)neg_face_idx), BCFaceToString((BCFace)pos_face_idx), driven_direction);
86 }
87
88 LOG_ALLOW(GLOBAL, LOG_DEBUG, " ... Symmetry and mathematical types are valid. OK.\n");
89
90 PetscFunctionReturn(0);
91}
92
93//================================================================================
94//
95// HANDLER IMPLEMENTATION: NO-SLIP WALL
96// (Corresponds to BC_HANDLER_WALL_NOSLIP)
97//
98// This handler implements a stationary, impenetrable wall where the fluid
99// velocity is zero (no-slip condition).
100//
101//================================================================================
102
103// --- FORWARD DECLARATIONS ---
104static PetscErrorCode Apply_WallNoSlip(BoundaryCondition *self, BCContext *ctx);
105
106#undef __FUNCT__
107#define __FUNCT__ "Create_WallNoSlip"
108/**
109 * @brief Implementation of \ref Create_WallNoSlip().
110 * @details Full API contract (arguments, ownership, side effects) is documented with
111 * the header declaration in `include/BC_Handlers.h`.
112 * @see Create_WallNoSlip()
113 */
114PetscErrorCode Create_WallNoSlip(BoundaryCondition *bc)
115{
116 PetscFunctionBeginUser;
117
118 if (!bc) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
119 "Input BoundaryCondition object is NULL in Create_WallNoSlip");
120
121 // ✅ Set priority
122 bc->priority = BC_PRIORITY_WALL;
123
124 // Assign function pointers
125 bc->Initialize = NULL;
126 bc->PreStep = NULL;
127 bc->Apply = Apply_WallNoSlip;
128 bc->PostStep = NULL;
129 bc->UpdateUbcs = NULL;
130 bc->Destroy = NULL;
131
132 // No private data needed for this simple handler
133 bc->data = NULL;
134
135 PetscFunctionReturn(0);
136}
137
138#undef __FUNCT__
139#define __FUNCT__ "Apply_WallNoSlip"
140/**
141 * @brief Internal helper implementation: `Apply_WallNoSlip()`.
142 * @details Local to this translation unit.
143 */
144static PetscErrorCode Apply_WallNoSlip(BoundaryCondition *self, BCContext *ctx)
145{
146 PetscErrorCode ierr;
147 UserCtx* user = ctx->user;
148 BCFace face_id = ctx->face_id;
149 PetscBool can_service;
150
151 (void)self; // Unused for simple handlers
152
153 PetscFunctionBeginUser;
154 DMDALocalInfo *info = &user->info;
155 Cmpnts ***ubcs, ***ucont;
156 PetscInt IM_nodes_global, JM_nodes_global,KM_nodes_global;
157
158 IM_nodes_global = user->IM;
159 JM_nodes_global = user->JM;
160 KM_nodes_global = user->KM;
161
162 ierr = CanRankServiceFace(info,IM_nodes_global,JM_nodes_global,KM_nodes_global,face_id,&can_service); CHKERRQ(ierr);
163 // Check if this rank owns part of this boundary face
164 if (!can_service) PetscFunctionReturn(0);
165
166 LOG_ALLOW(LOCAL, LOG_DEBUG, "Apply_WallNoSlip: Applying to Face %d (%s).\n",
167 face_id, BCFaceToString(face_id));
168
169 // Get arrays
170
171 ierr = DMDAVecGetArray(user->fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
172 ierr = DMDAVecGetArray(user->fda, user->Ucont, &ucont); CHKERRQ(ierr);
173
174 PetscInt xs = info->xs, xe = info->xs + info->xm;
175 PetscInt ys = info->ys, ye = info->ys + info->ym;
176 PetscInt zs = info->zs, ze = info->zs + info->zm;
177 PetscInt mx = info->mx, my = info->my, mz = info->mz;
178
179 // ✅ Use shrunken loop bounds (avoids edges/corners like inlet handler)
180 PetscInt lxs = xs, lxe = xe, lys = ys, lye = ye, lzs = zs, lze = ze;
181 if (xs == 0) lxs = xs + 1;
182 if (xe == mx) lxe = xe - 1;
183 if (ys == 0) lys = ys + 1;
184 if (ye == my) lye = ye - 1;
185 if (zs == 0) lzs = zs + 1;
186 if (ze == mz) lze = ze - 1;
187
188 switch (face_id) {
189 case BC_FACE_NEG_X: {
190 if (xs == 0){
191 PetscInt i = xs;
192 for (PetscInt k = lzs; k < lze; k++) {
193 for (PetscInt j = lys; j < lye; j++) {
194 // ✅ Set contravariant flux to zero (no penetration)
195 ucont[k][j][i].x = 0.0;
196
197 // ✅ Set boundary velocity to zero (no slip)
198 ubcs[k][j][i].x = 0.0;
199 ubcs[k][j][i].y = 0.0;
200 ubcs[k][j][i].z = 0.0;
201 }
202 }
203 }
204 break;
205 }
206
207 case BC_FACE_POS_X: {
208 if (xe == mx){
209 PetscInt i = xe - 1;
210 for (PetscInt k = lzs; k < lze; k++) {
211 for (PetscInt j = lys; j < lye; j++) {
212 ucont[k][j][i-1].x = 0.0;
213
214 ubcs[k][j][i].x = 0.0;
215 ubcs[k][j][i].y = 0.0;
216 ubcs[k][j][i].z = 0.0;
217 }
218 }
219 }
220 break;
221 }
222 case BC_FACE_NEG_Y: {
223 if (ys == 0){
224 PetscInt j = ys;
225 for (PetscInt k = lzs; k < lze; k++) {
226 for (PetscInt i = lxs; i < lxe; i++) {
227 ucont[k][j][i].y = 0.0;
228
229 ubcs[k][j][i].x = 0.0;
230 ubcs[k][j][i].y = 0.0;
231 ubcs[k][j][i].z = 0.0;
232 }
233 }
234 }
235 } break;
236
237 case BC_FACE_POS_Y: {
238 if (ye == my){
239 PetscInt j = ye - 1;
240 for (PetscInt k = lzs; k < lze; k++) {
241 for (PetscInt i = lxs; i < lxe; i++) {
242 ucont[k][j-1][i].y = 0.0;
243
244 ubcs[k][j][i].x = 0.0;
245 ubcs[k][j][i].y = 0.0;
246 ubcs[k][j][i].z = 0.0;
247 }
248 }
249 }
250 } break;
251
252 case BC_FACE_NEG_Z: {
253 if (zs == 0){
254 PetscInt k = zs;
255 for (PetscInt j = lys; j < lye; j++) {
256 for (PetscInt i = lxs; i < lxe; i++) {
257 ucont[k][j][i].z = 0.0;
258
259 ubcs[k][j][i].x = 0.0;
260 ubcs[k][j][i].y = 0.0;
261 ubcs[k][j][i].z = 0.0;
262 }
263 }
264 }
265 } break;
266
267 case BC_FACE_POS_Z: {
268 if (ze == mz) {
269 PetscInt k = ze - 1;
270 for (PetscInt j = lys; j < lye; j++) {
271 for (PetscInt i = lxs; i < lxe; i++) {
272 ucont[k-1][j][i].z = 0.0;
273
274 ubcs[k][j][i].x = 0.0;
275 ubcs[k][j][i].y = 0.0;
276 ubcs[k][j][i].z = 0.0;
277 }
278 }
279 }
280 } break;
281 }
282
283 // Restore arrays
284 ierr = DMDAVecRestoreArray(user->fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
285 ierr = DMDAVecRestoreArray(user->fda, user->Ucont, &ucont); CHKERRQ(ierr);
286
287 PetscFunctionReturn(0);
288}
289
290////////////////////////////////////////////////////////////////////////
291
292//================================================================================
293//
294// HANDLER IMPLEMENTATION: CONSTANT VELOCITY INLET
295// (Corresponds to BC_HANDLER_INLET_CONSTANT_VELOCITY)
296//
297//================================================================================
298
299// --- FORWARD DECLARATIONS ---
300static PetscErrorCode Initialize_InletConstantVelocity(BoundaryCondition *self, BCContext *ctx);
301static PetscErrorCode PreStep_InletConstantVelocity(BoundaryCondition *self, BCContext *ctx,
302 PetscReal *in, PetscReal *out);
303static PetscErrorCode Apply_InletConstantVelocity(BoundaryCondition *self, BCContext *ctx);
304static PetscErrorCode PostStep_InletConstantVelocity(BoundaryCondition *self, BCContext *ctx,
305 PetscReal *in, PetscReal *out);
306static PetscErrorCode Destroy_InletConstantVelocity(BoundaryCondition *self);
307
308/**
309 * @brief Private data structure for the Constant Velocity Inlet handler.
310 */
311typedef struct{
312 PetscReal normal_velocity; // The desired Cartesian velocity (vx, vy, vz)
313}InletConstantData;
314
315#undef __FUNCT__
316#define __FUNCT__ "Create_InletConstantVelocity"
317/**
318 * @brief Implementation of \ref Create_InletConstantVelocity().
319 * @details Full API contract (arguments, ownership, side effects) is documented with
320 * the header declaration in `include/BC_Handlers.h`.
321 * @see Create_InletConstantVelocity()
322 */
323PetscErrorCode Create_InletConstantVelocity(BoundaryCondition *bc)
324{
325 PetscErrorCode ierr;
326 PetscFunctionBeginUser;
327
328 if (!bc) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "BoundaryCondition is NULL");
329
330 InletConstantData *data = NULL;
331 ierr = PetscMalloc1(1, &data); CHKERRQ(ierr);
332 bc->data = (void*)data;
333
334 bc->priority = BC_PRIORITY_INLET;
335 bc->Initialize = Initialize_InletConstantVelocity;
336 bc->PreStep = PreStep_InletConstantVelocity;
337 bc->Apply = Apply_InletConstantVelocity;
338 bc->PostStep = PostStep_InletConstantVelocity;
339 bc->UpdateUbcs = NULL;
340 bc->Destroy = Destroy_InletConstantVelocity;
341
342 PetscFunctionReturn(0);
343}
344
345
346#undef __FUNCT__
347#define __FUNCT__ "Initialize_InletConstantVelocity"
348/**
349 * @brief Internal helper implementation: `Initialize_InletConstantVelocity()`.
350 * @details Local to this translation unit.
351 */
352static PetscErrorCode Initialize_InletConstantVelocity(BoundaryCondition *self, BCContext *ctx)
353{
354 PetscErrorCode ierr;
355 UserCtx* user = ctx->user;
356 BCFace face_id = ctx->face_id;
357 InletConstantData *data = (InletConstantData*)self->data;
358 PetscBool found;
359
360 PetscFunctionBeginUser;
361 LOG_ALLOW(LOCAL, LOG_DEBUG, "Initialize_InletConstantVelocity: Initializing handler for Face %d. \n", face_id);
362 data->normal_velocity = 0.0;
363
364 switch (face_id) {
365 case BC_FACE_NEG_X:
366 case BC_FACE_POS_X:
367 // For X-faces, read "vx" as normal velocity
368 ierr = GetBCParamReal(user->boundary_faces[face_id].params, "vx",
369 &data->normal_velocity, &found); CHKERRQ(ierr);
370 break;
371
372 case BC_FACE_NEG_Y:
373 case BC_FACE_POS_Y:
374 // For Y-faces, read "vy" as normal velocity
375 ierr = GetBCParamReal(user->boundary_faces[face_id].params, "vy",
376 &data->normal_velocity, &found); CHKERRQ(ierr);
377 break;
378
379 case BC_FACE_NEG_Z:
380 case BC_FACE_POS_Z:
381 // For Z-faces, read "vz" as normal velocity
382 ierr = GetBCParamReal(user->boundary_faces[face_id].params, "vz",
383 &data->normal_velocity, &found); CHKERRQ(ierr);
384 break;
385 }
386
387 LOG_ALLOW(LOCAL, LOG_INFO, " Inlet Face %d: normal velocity = %.4f\n",
388 face_id, data->normal_velocity);
389
390 // Set initial boundary state
391 ierr = Apply_InletConstantVelocity(self, ctx); CHKERRQ(ierr);
392
393 PetscFunctionReturn(0);
394}
395
396#undef __FUNCT__
397#define __FUNCT__ "PreStep_InletConstantVelocity"
398/**
399 * @brief Internal helper implementation: `PreStep_InletConstantVelocity()`.
400 * @details Local to this translation unit.
401 */
402static PetscErrorCode PreStep_InletConstantVelocity(BoundaryCondition *self, BCContext *ctx,
403 PetscReal *local_inflow_contribution,
404 PetscReal *local_outflow_contribution)
405{
406 // No preparation needed for constant velocity inlet.
407 // The velocity is already stored in self->data from Initialize.
408 // Apply will set ucont, and PostStep will measure the actual flux.
409
410 (void)self;
411 (void)ctx;
412 (void)local_inflow_contribution;
413 (void)local_outflow_contribution;
414
415 PetscFunctionBeginUser;
416 PetscFunctionReturn(0);
417}
418
419#undef __FUNCT__
420#define __FUNCT__ "Apply_InletConstantVelocity"
421/**
422 * @brief Internal helper implementation: `Apply_InletConstantVelocity()`.
423 * @details Local to this translation unit.
424 */
425static PetscErrorCode Apply_InletConstantVelocity(BoundaryCondition *self, BCContext *ctx)
426{
427 PetscErrorCode ierr;
428 UserCtx* user = ctx->user;
429 BCFace face_id = ctx->face_id;
430 InletConstantData *data = (InletConstantData*)self->data;
431 PetscBool can_service;
432
433 PetscFunctionBeginUser;
434
435 DMDALocalInfo *info = &user->info;
436 Cmpnts ***ubcs, ***ucont, ***csi, ***eta, ***zet;
437 PetscReal ***nvert;
438 PetscInt IM_nodes_global, JM_nodes_global,KM_nodes_global;
439
440 IM_nodes_global = user->IM;
441 JM_nodes_global = user->JM;
442 KM_nodes_global = user->KM;
443
444 ierr = CanRankServiceFace(info,IM_nodes_global,JM_nodes_global,KM_nodes_global,face_id,&can_service); CHKERRQ(ierr);
445
446 if (!can_service) PetscFunctionReturn(0);
447
448 // Get arrays
449
450 ierr = DMDAVecGetArray(user->fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
451 ierr = DMDAVecGetArray(user->fda, user->Ucont, &ucont); CHKERRQ(ierr);
452 ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
453 ierr = DMDAVecGetArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
454 ierr = DMDAVecGetArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
455 ierr = DMDAVecGetArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
456
457 // Get SCALAR velocity (not vector!)
458 PetscReal uin_this_point = data->normal_velocity;
459
460 PetscInt xs = info->xs, xe = info->xs + info->xm;
461 PetscInt ys = info->ys, ye = info->ys + info->ym;
462 PetscInt zs = info->zs, ze = info->zs + info->zm;
463 PetscInt mx = info->mx, my = info->my, mz = info->mz;
464
465 PetscInt lxs = xs, lxe = xe, lys = ys, lye = ye, lzs = zs, lze = ze;
466 if (xs == 0) lxs = xs + 1;
467 if (xe == mx) lxe = xe - 1;
468 if (ys == 0) lys = ys + 1;
469 if (ye == my) lye = ye - 1;
470 if (zs == 0) lzs = zs + 1;
471 if (ze == mz) lze = ze - 1;
472
473 switch (face_id) {
474 case BC_FACE_NEG_X:
475 case BC_FACE_POS_X: {
476 PetscReal sign = (face_id == BC_FACE_NEG_X) ? 1.0 : -1.0;
477 PetscInt i = (face_id == BC_FACE_NEG_X) ? xs : mx - 2;
478
479 for (PetscInt k = lzs; k < lze; k++) {
480 for (PetscInt j = lys; j < lye; j++) {
481 if ((sign > 0 && nvert[k][j][i+1] > 0.1) ||
482 (sign < 0 && nvert[k][j][i] > 0.1)) continue;
483
484 PetscReal CellArea = sqrt(csi[k][j][i].x * csi[k][j][i].x +
485 csi[k][j][i].y * csi[k][j][i].y +
486 csi[k][j][i].z * csi[k][j][i].z);
487
488 ucont[k][j][i].x = sign * uin_this_point * CellArea;
489
490 ubcs[k][j][i + (sign < 0)].x = sign * uin_this_point * csi[k][j][i].x / CellArea;
491 ubcs[k][j][i + (sign < 0)].y = sign * uin_this_point * csi[k][j][i].y / CellArea;
492 ubcs[k][j][i + (sign < 0)].z = sign * uin_this_point * csi[k][j][i].z / CellArea;
493 }
494 }
495 } break;
496
497 case BC_FACE_NEG_Y:
498 case BC_FACE_POS_Y: {
499 PetscReal sign = (face_id == BC_FACE_NEG_Y) ? 1.0 : -1.0;
500 PetscInt j = (face_id == BC_FACE_NEG_Y) ? ys : my - 2;
501
502 for (PetscInt k = lzs; k < lze; k++) {
503 for (PetscInt i = lxs; i < lxe; i++) {
504 if ((sign > 0 && nvert[k][j+1][i] > 0.1) ||
505 (sign < 0 && nvert[k][j][i] > 0.1)) continue;
506
507 PetscReal CellArea = sqrt(eta[k][j][i].x * eta[k][j][i].x +
508 eta[k][j][i].y * eta[k][j][i].y +
509 eta[k][j][i].z * eta[k][j][i].z);
510
511 ucont[k][j][i].y = sign * uin_this_point * CellArea;
512
513 ubcs[k][j + (sign < 0)][i].x = sign * uin_this_point * eta[k][j][i].x / CellArea;
514 ubcs[k][j + (sign < 0)][i].y = sign * uin_this_point * eta[k][j][i].y / CellArea;
515 ubcs[k][j + (sign < 0)][i].z = sign * uin_this_point * eta[k][j][i].z / CellArea;
516 }
517 }
518 } break;
519
520 case BC_FACE_NEG_Z:
521 case BC_FACE_POS_Z: {
522 PetscReal sign = (face_id == BC_FACE_NEG_Z) ? 1.0 : -1.0;
523 PetscInt k = (face_id == BC_FACE_NEG_Z) ? zs : mz - 2;
524
525 for (PetscInt j = lys; j < lye; j++) {
526 for (PetscInt i = lxs; i < lxe; i++) {
527 if ((sign > 0 && nvert[k+1][j][i] > 0.1) ||
528 (sign < 0 && nvert[k][j][i] > 0.1)) continue;
529
530 PetscReal CellArea = sqrt(zet[k][j][i].x * zet[k][j][i].x +
531 zet[k][j][i].y * zet[k][j][i].y +
532 zet[k][j][i].z * zet[k][j][i].z);
533
534 ucont[k][j][i].z = sign * uin_this_point * CellArea;
535
536 ubcs[k + (sign < 0)][j][i].x = sign * uin_this_point * zet[k][j][i].x / CellArea;
537 ubcs[k + (sign < 0)][j][i].y = sign * uin_this_point * zet[k][j][i].y / CellArea;
538 ubcs[k + (sign < 0)][j][i].z = sign * uin_this_point * zet[k][j][i].z / CellArea;
539 }
540 }
541 } break;
542 }
543
544 // Restore arrays
545 ierr = DMDAVecRestoreArray(user->fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
546 ierr = DMDAVecRestoreArray(user->fda, user->Ucont, &ucont); CHKERRQ(ierr);
547 ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
548 ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
549 ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
550 ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
551
552 PetscFunctionReturn(0);
553}
554
555#undef __FUNCT__
556#define __FUNCT__ "PostStep_InletConstantVelocity"
557/**
558 * @brief Internal helper implementation: `PostStep_InletConstantVelocity()`.
559 * @details Local to this translation unit.
560 */
561static PetscErrorCode PostStep_InletConstantVelocity(BoundaryCondition *self, BCContext *ctx,
562 PetscReal *local_inflow_contribution,
563 PetscReal *local_outflow_contribution)
564{
565 PetscErrorCode ierr;
566 UserCtx* user = ctx->user;
567 BCFace face_id = ctx->face_id;
568 PetscBool can_service;
569
570 (void)self;
571 (void)local_outflow_contribution;
572
573 PetscFunctionBeginUser;
574
575 DMDALocalInfo *info = &user->info;
576 Cmpnts ***ucont;
577
578 PetscInt IM_nodes_global, JM_nodes_global,KM_nodes_global;
579
580 IM_nodes_global = user->IM;
581 JM_nodes_global = user->JM;
582 KM_nodes_global = user->KM;
583
584 ierr = CanRankServiceFace(info,IM_nodes_global,JM_nodes_global,KM_nodes_global,face_id,&can_service); CHKERRQ(ierr);
585
586
587 if (!can_service) PetscFunctionReturn(0);
588
589 ierr = DMDAVecGetArrayRead(user->fda, user->Ucont, (const Cmpnts***)&ucont); CHKERRQ(ierr);
590
591 PetscReal local_flux = 0.0;
592
593 PetscInt xs = info->xs, xe = info->xs + info->xm;
594 PetscInt ys = info->ys, ye = info->ys + info->ym;
595 PetscInt zs = info->zs, ze = info->zs + info->zm;
596 PetscInt mx = info->mx, my = info->my, mz = info->mz;
597
598 PetscInt lxs = xs, lxe = xe, lys = ys, lye = ye, lzs = zs, lze = ze;
599 if (xs == 0) lxs = xs + 1;
600 if (xe == mx) lxe = xe - 1;
601 if (ys == 0) lys = ys + 1;
602 if (ye == my) lye = ye - 1;
603 if (zs == 0) lzs = zs + 1;
604 if (ze == mz) lze = ze - 1;
605
606 // Sum ucont components
607 switch (face_id) {
608 case BC_FACE_NEG_X:
609 case BC_FACE_POS_X: {
610 PetscInt i = (face_id == BC_FACE_NEG_X) ? xs : mx - 2;
611 for (PetscInt k = lzs; k < lze; k++) {
612 for (PetscInt j = lys; j < lye; j++) {
613 local_flux += ucont[k][j][i].x;
614 }
615 }
616 } break;
617
618 case BC_FACE_NEG_Y:
619 case BC_FACE_POS_Y: {
620 PetscInt j = (face_id == BC_FACE_NEG_Y) ? ys : my - 2;
621 for (PetscInt k = lzs; k < lze; k++) {
622 for (PetscInt i = lxs; i < lxe; i++) {
623 local_flux += ucont[k][j][i].y;
624 }
625 }
626 } break;
627
628 case BC_FACE_NEG_Z:
629 case BC_FACE_POS_Z: {
630 PetscInt k = (face_id == BC_FACE_NEG_Z) ? zs : mz - 2;
631 for (PetscInt j = lys; j < lye; j++) {
632 for (PetscInt i = lxs; i < lxe; i++) {
633 local_flux += ucont[k][j][i].z;
634 }
635 }
636 } break;
637 }
638
639 ierr = DMDAVecRestoreArrayRead(user->fda, user->Ucont, (const Cmpnts***)&ucont); CHKERRQ(ierr);
640
641 *local_inflow_contribution += local_flux;
642
643 LOG_ALLOW(LOCAL, LOG_DEBUG, "PostStep_InletConstantVelocity: Face %d, flux = %.6e\n",
644 face_id, local_flux);
645
646 PetscFunctionReturn(0);
647}
648
649
650#undef __FUNCT__
651#define __FUNCT__ "Destroy_InletConstantVelocity"
652/**
653 * @brief Internal helper implementation: `Destroy_InletConstantVelocity()`.
654 * @details Local to this translation unit.
655 */
656static PetscErrorCode Destroy_InletConstantVelocity(BoundaryCondition *self)
657{
658 PetscFunctionBeginUser;
659 if (self && self->data) {
660 PetscFree(self->data);
661 self->data = NULL;
662 }
663 PetscFunctionReturn(0);
664}
665
666//================================================================================
667//
668// HANDLER IMPLEMENTATION: PARABOLIC VELOCITY INLET (POISEUILLE PROFILE)
669// (Corresponds to BC_HANDLER_INLET_PARABOLIC)
670//
671// This handler enforces a fully-developed parabolic (Poiseuille) velocity profile
672// on a rectangular/square inlet face. The profile shape is:
673//
674// V(cs1, cs2) = v_max * (1 - cs1_norm^2) * (1 - cs2_norm^2)
675//
676// where cs1 and cs2 are the two cross-stream index directions for the given face,
677// normalized to [-1, +1] across the interior nodes. The profile is zero at the
678// walls and peaks at v_max at the center.
679//
680// Workflow:
681// Constructor -> Allocate private data, wire function pointers.
682// Initialize -> Parse v_max from params, compute cross-stream geometry, call Apply.
683// PreStep -> No-op (profile is static in time).
684// Apply -> Set ucont and ubcs on the inlet face using the parabolic profile.
685// PostStep -> Measure actual volumetric flux through the face.
686// Destroy -> Free private data.
687//
688//================================================================================
689
690// --- FORWARD DECLARATIONS ---
691static PetscErrorCode Initialize_InletParabolicProfile(BoundaryCondition *self, BCContext *ctx);
692static PetscErrorCode PreStep_InletParabolicProfile(BoundaryCondition *self, BCContext *ctx,
693 PetscReal *in, PetscReal *out);
694static PetscErrorCode Apply_InletParabolicProfile(BoundaryCondition *self, BCContext *ctx);
695static PetscErrorCode PostStep_InletParabolicProfile(BoundaryCondition *self, BCContext *ctx,
696 PetscReal *in, PetscReal *out);
697static PetscErrorCode Destroy_InletParabolicProfile(BoundaryCondition *self);
698
699/**
700 * @brief Private data structure for the Parabolic Velocity Inlet handler.
701 *
702 * Stores the peak velocity and pre-computed cross-stream geometry needed
703 * to evaluate the parabolic profile at each boundary node.
704 */
705typedef struct {
706 PetscReal v_max; /**< Peak centerline velocity (from user params). */
707 PetscReal cs1_center; /**< Center index in cross-stream direction 1. */
708 PetscReal cs2_center; /**< Center index in cross-stream direction 2. */
709 PetscReal cs1_half; /**< Half-width (in index space) in cross-stream direction 1. */
710 PetscReal cs2_half; /**< Half-width (in index space) in cross-stream direction 2. */
711} InletParabolicData;
712
713#undef __FUNCT__
714#define __FUNCT__ "Create_InletParabolicProfile"
715/**
716 * @brief Implementation of \ref Create_InletParabolicProfile().
717 * @details Full API contract (arguments, ownership, side effects) is documented with
718 * the header declaration in `include/BC_Handlers.h`.
719 * @see Create_InletParabolicProfile()
720 */
721PetscErrorCode Create_InletParabolicProfile(BoundaryCondition *bc)
722{
723 PetscErrorCode ierr;
724 PetscFunctionBeginUser;
725
726 if (!bc) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "BoundaryCondition is NULL");
727
728 InletParabolicData *data = NULL;
729 ierr = PetscMalloc1(1, &data); CHKERRQ(ierr);
730 bc->data = (void*)data;
731
732 bc->priority = BC_PRIORITY_INLET;
733 bc->Initialize = Initialize_InletParabolicProfile;
734 bc->PreStep = PreStep_InletParabolicProfile;
735 bc->Apply = Apply_InletParabolicProfile;
736 bc->PostStep = PostStep_InletParabolicProfile;
737 bc->UpdateUbcs = NULL;
738 bc->Destroy = Destroy_InletParabolicProfile;
739
740 PetscFunctionReturn(0);
741}
742
743
744#undef __FUNCT__
745#define __FUNCT__ "Initialize_InletParabolicProfile"
746/**
747 * @brief Internal helper implementation: `Initialize_InletParabolicProfile()`.
748 * @details Local to this translation unit.
749 */
750static PetscErrorCode Initialize_InletParabolicProfile(BoundaryCondition *self, BCContext *ctx)
751{
752 PetscErrorCode ierr;
753 UserCtx* user = ctx->user;
754 BCFace face_id = ctx->face_id;
755 InletParabolicData *data = (InletParabolicData*)self->data;
756 PetscBool found;
757
758 PetscFunctionBeginUser;
759 LOG_ALLOW(LOCAL, LOG_DEBUG, "Initialize_InletParabolicProfile: Initializing handler for Face %d.\n", face_id);
760
761 // --- Parse v_max from boundary condition parameters ---
762 data->v_max = 0.0;
763 ierr = GetBCParamReal(user->boundary_faces[face_id].params, "v_max",
764 &data->v_max, &found); CHKERRQ(ierr);
765 if (!found) {
766 LOG_ALLOW(GLOBAL, LOG_WARNING, "Initialize_InletParabolicProfile: 'v_max' not found in params for face %d. Defaulting to 0.0.\n", face_id);
767 }
768
769 // --- Determine cross-stream dimensions based on face orientation ---
770 PetscReal cs1_dim, cs2_dim;
771 switch (face_id) {
772 case BC_FACE_NEG_X:
773 case BC_FACE_POS_X:
774 cs1_dim = (PetscReal)user->JM; // j-direction
775 cs2_dim = (PetscReal)user->KM; // k-direction
776 break;
777 case BC_FACE_NEG_Y:
778 case BC_FACE_POS_Y:
779 cs1_dim = (PetscReal)user->IM; // i-direction
780 cs2_dim = (PetscReal)user->KM; // k-direction
781 break;
782 case BC_FACE_NEG_Z:
783 case BC_FACE_POS_Z:
784 default:
785 cs1_dim = (PetscReal)user->IM; // i-direction
786 cs2_dim = (PetscReal)user->JM; // j-direction
787 break;
788 }
789
790 // Interior width = dim - 2 (nodes 1 through dim-2 are interior)
791 PetscReal cs1_width = cs1_dim - 2.0;
792 PetscReal cs2_width = cs2_dim - 2.0;
793
794 data->cs1_center = 1.0 + cs1_width / 2.0;
795 data->cs2_center = 1.0 + cs2_width / 2.0;
796 data->cs1_half = cs1_width / 2.0;
797 data->cs2_half = cs2_width / 2.0;
798
799 LOG_ALLOW(LOCAL, LOG_INFO, " Inlet Face %d (Parabolic): v_max = %.4f\n", face_id, data->v_max);
800 LOG_ALLOW(LOCAL, LOG_DEBUG, " Cross-stream 1: center=%.1f, half=%.1f\n", data->cs1_center, data->cs1_half);
801 LOG_ALLOW(LOCAL, LOG_DEBUG, " Cross-stream 2: center=%.1f, half=%.1f\n", data->cs2_center, data->cs2_half);
802
803 // Set initial boundary state
804 ierr = Apply_InletParabolicProfile(self, ctx); CHKERRQ(ierr);
805
806 PetscFunctionReturn(0);
807}
808
809
810#undef __FUNCT__
811#define __FUNCT__ "PreStep_InletParabolicProfile"
812/**
813 * @brief Internal helper implementation: `PreStep_InletParabolicProfile()`.
814 * @details Local to this translation unit.
815 */
816static PetscErrorCode PreStep_InletParabolicProfile(BoundaryCondition *self, BCContext *ctx,
817 PetscReal *local_inflow_contribution,
818 PetscReal *local_outflow_contribution)
819{
820 (void)self;
821 (void)ctx;
822 (void)local_inflow_contribution;
823 (void)local_outflow_contribution;
824
825 PetscFunctionBeginUser;
826 PetscFunctionReturn(0);
827}
828
829
830#undef __FUNCT__
831#define __FUNCT__ "Apply_InletParabolicProfile"
832/**
833 * @brief Internal helper implementation: `Apply_InletParabolicProfile()`.
834 * @details Local to this translation unit.
835 */
836static PetscErrorCode Apply_InletParabolicProfile(BoundaryCondition *self, BCContext *ctx)
837{
838 PetscErrorCode ierr;
839 UserCtx* user = ctx->user;
840 BCFace face_id = ctx->face_id;
841 InletParabolicData *data = (InletParabolicData*)self->data;
842 PetscBool can_service;
843
844 PetscFunctionBeginUser;
845
846 DMDALocalInfo *info = &user->info;
847 Cmpnts ***ubcs, ***ucont, ***csi, ***eta, ***zet;
848 PetscReal ***nvert;
849 PetscInt IM_nodes_global, JM_nodes_global, KM_nodes_global;
850
851 IM_nodes_global = user->IM;
852 JM_nodes_global = user->JM;
853 KM_nodes_global = user->KM;
854
855 ierr = CanRankServiceFace(info, IM_nodes_global, JM_nodes_global, KM_nodes_global,
856 face_id, &can_service); CHKERRQ(ierr);
857
858 if (!can_service) PetscFunctionReturn(0);
859
860 // Get arrays
861 ierr = DMDAVecGetArray(user->fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
862 ierr = DMDAVecGetArray(user->fda, user->Ucont, &ucont); CHKERRQ(ierr);
863 ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
864 ierr = DMDAVecGetArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
865 ierr = DMDAVecGetArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
866 ierr = DMDAVecGetArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
867
868 PetscInt xs = info->xs, xe = info->xs + info->xm;
869 PetscInt ys = info->ys, ye = info->ys + info->ym;
870 PetscInt zs = info->zs, ze = info->zs + info->zm;
871 PetscInt mx = info->mx, my = info->my, mz = info->mz;
872
873 PetscInt lxs = xs, lxe = xe, lys = ys, lye = ye, lzs = zs, lze = ze;
874 if (xs == 0) lxs = xs + 1;
875 if (xe == mx) lxe = xe - 1;
876 if (ys == 0) lys = ys + 1;
877 if (ye == my) lye = ye - 1;
878 if (zs == 0) lzs = zs + 1;
879 if (ze == mz) lze = ze - 1;
880
881 switch (face_id) {
882 case BC_FACE_NEG_X:
883 case BC_FACE_POS_X: {
884 // X-faces: normal = i, cross-stream = (j, k)
885 PetscReal sign = (face_id == BC_FACE_NEG_X) ? 1.0 : -1.0;
886 PetscInt i = (face_id == BC_FACE_NEG_X) ? xs : mx - 2;
887
888 for (PetscInt k = lzs; k < lze; k++) {
889 for (PetscInt j = lys; j < lye; j++) {
890 if ((sign > 0 && nvert[k][j][i+1] > 0.1) ||
891 (sign < 0 && nvert[k][j][i] > 0.1)) continue;
892
893 // Evaluate parabolic profile: cs1 = j, cs2 = k
894 PetscReal cs1_norm = ((PetscReal)j - data->cs1_center) / data->cs1_half;
895 PetscReal cs2_norm = ((PetscReal)k - data->cs2_center) / data->cs2_half;
896 PetscReal profile = PetscMax(0.0, 1.0 - cs1_norm * cs1_norm)
897 * PetscMax(0.0, 1.0 - cs2_norm * cs2_norm);
898 PetscReal uin_local = data->v_max * profile;
899
900 PetscReal CellArea = sqrt(csi[k][j][i].x * csi[k][j][i].x +
901 csi[k][j][i].y * csi[k][j][i].y +
902 csi[k][j][i].z * csi[k][j][i].z);
903
904 ucont[k][j][i].x = sign * uin_local * CellArea;
905
906 ubcs[k][j][i + (sign < 0)].x = sign * uin_local * csi[k][j][i].x / CellArea;
907 ubcs[k][j][i + (sign < 0)].y = sign * uin_local * csi[k][j][i].y / CellArea;
908 ubcs[k][j][i + (sign < 0)].z = sign * uin_local * csi[k][j][i].z / CellArea;
909 }
910 }
911 } break;
912
913 case BC_FACE_NEG_Y:
914 case BC_FACE_POS_Y: {
915 // Y-faces: normal = j, cross-stream = (i, k)
916 PetscReal sign = (face_id == BC_FACE_NEG_Y) ? 1.0 : -1.0;
917 PetscInt j = (face_id == BC_FACE_NEG_Y) ? ys : my - 2;
918
919 for (PetscInt k = lzs; k < lze; k++) {
920 for (PetscInt i = lxs; i < lxe; i++) {
921 if ((sign > 0 && nvert[k][j+1][i] > 0.1) ||
922 (sign < 0 && nvert[k][j][i] > 0.1)) continue;
923
924 // Evaluate parabolic profile: cs1 = i, cs2 = k
925 PetscReal cs1_norm = ((PetscReal)i - data->cs1_center) / data->cs1_half;
926 PetscReal cs2_norm = ((PetscReal)k - data->cs2_center) / data->cs2_half;
927 PetscReal profile = PetscMax(0.0, 1.0 - cs1_norm * cs1_norm)
928 * PetscMax(0.0, 1.0 - cs2_norm * cs2_norm);
929 PetscReal uin_local = data->v_max * profile;
930
931 PetscReal CellArea = sqrt(eta[k][j][i].x * eta[k][j][i].x +
932 eta[k][j][i].y * eta[k][j][i].y +
933 eta[k][j][i].z * eta[k][j][i].z);
934
935 ucont[k][j][i].y = sign * uin_local * CellArea;
936
937 ubcs[k][j + (sign < 0)][i].x = sign * uin_local * eta[k][j][i].x / CellArea;
938 ubcs[k][j + (sign < 0)][i].y = sign * uin_local * eta[k][j][i].y / CellArea;
939 ubcs[k][j + (sign < 0)][i].z = sign * uin_local * eta[k][j][i].z / CellArea;
940 }
941 }
942 } break;
943
944 case BC_FACE_NEG_Z:
945 case BC_FACE_POS_Z: {
946 // Z-faces: normal = k, cross-stream = (i, j)
947 PetscReal sign = (face_id == BC_FACE_NEG_Z) ? 1.0 : -1.0;
948 PetscInt k = (face_id == BC_FACE_NEG_Z) ? zs : mz - 2;
949
950 for (PetscInt j = lys; j < lye; j++) {
951 for (PetscInt i = lxs; i < lxe; i++) {
952 if ((sign > 0 && nvert[k+1][j][i] > 0.1) ||
953 (sign < 0 && nvert[k][j][i] > 0.1)) continue;
954
955 // Evaluate parabolic profile: cs1 = i, cs2 = j
956 PetscReal cs1_norm = ((PetscReal)i - data->cs1_center) / data->cs1_half;
957 PetscReal cs2_norm = ((PetscReal)j - data->cs2_center) / data->cs2_half;
958 PetscReal profile = PetscMax(0.0, 1.0 - cs1_norm * cs1_norm)
959 * PetscMax(0.0, 1.0 - cs2_norm * cs2_norm);
960 PetscReal uin_local = data->v_max * profile;
961
962 PetscReal CellArea = sqrt(zet[k][j][i].x * zet[k][j][i].x +
963 zet[k][j][i].y * zet[k][j][i].y +
964 zet[k][j][i].z * zet[k][j][i].z);
965
966 ucont[k][j][i].z = sign * uin_local * CellArea;
967
968 ubcs[k + (sign < 0)][j][i].x = sign * uin_local * zet[k][j][i].x / CellArea;
969 ubcs[k + (sign < 0)][j][i].y = sign * uin_local * zet[k][j][i].y / CellArea;
970 ubcs[k + (sign < 0)][j][i].z = sign * uin_local * zet[k][j][i].z / CellArea;
971 }
972 }
973 } break;
974 }
975
976 // Restore arrays
977 ierr = DMDAVecRestoreArray(user->fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
978 ierr = DMDAVecRestoreArray(user->fda, user->Ucont, &ucont); CHKERRQ(ierr);
979 ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
980 ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
981 ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
982 ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
983
984 PetscFunctionReturn(0);
985}
986
987
988#undef __FUNCT__
989#define __FUNCT__ "PostStep_InletParabolicProfile"
990/**
991 * @brief Internal helper implementation: `PostStep_InletParabolicProfile()`.
992 * @details Local to this translation unit.
993 */
994static PetscErrorCode PostStep_InletParabolicProfile(BoundaryCondition *self, BCContext *ctx,
995 PetscReal *local_inflow_contribution,
996 PetscReal *local_outflow_contribution)
997{
998 PetscErrorCode ierr;
999 UserCtx* user = ctx->user;
1000 BCFace face_id = ctx->face_id;
1001 PetscBool can_service;
1002
1003 (void)self;
1004 (void)local_outflow_contribution;
1005
1006 PetscFunctionBeginUser;
1007
1008 DMDALocalInfo *info = &user->info;
1009 Cmpnts ***ucont;
1010
1011 PetscInt IM_nodes_global, JM_nodes_global, KM_nodes_global;
1012
1013 IM_nodes_global = user->IM;
1014 JM_nodes_global = user->JM;
1015 KM_nodes_global = user->KM;
1016
1017 ierr = CanRankServiceFace(info, IM_nodes_global, JM_nodes_global, KM_nodes_global,
1018 face_id, &can_service); CHKERRQ(ierr);
1019
1020 if (!can_service) PetscFunctionReturn(0);
1021
1022 ierr = DMDAVecGetArrayRead(user->fda, user->Ucont, (const Cmpnts***)&ucont); CHKERRQ(ierr);
1023
1024 PetscReal local_flux = 0.0;
1025
1026 PetscInt xs = info->xs, xe = info->xs + info->xm;
1027 PetscInt ys = info->ys, ye = info->ys + info->ym;
1028 PetscInt zs = info->zs, ze = info->zs + info->zm;
1029 PetscInt mx = info->mx, my = info->my, mz = info->mz;
1030
1031 PetscInt lxs = xs, lxe = xe, lys = ys, lye = ye, lzs = zs, lze = ze;
1032 if (xs == 0) lxs = xs + 1;
1033 if (xe == mx) lxe = xe - 1;
1034 if (ys == 0) lys = ys + 1;
1035 if (ye == my) lye = ye - 1;
1036 if (zs == 0) lzs = zs + 1;
1037 if (ze == mz) lze = ze - 1;
1038
1039 switch (face_id) {
1040 case BC_FACE_NEG_X:
1041 case BC_FACE_POS_X: {
1042 PetscInt i = (face_id == BC_FACE_NEG_X) ? xs : mx - 2;
1043 for (PetscInt k = lzs; k < lze; k++) {
1044 for (PetscInt j = lys; j < lye; j++) {
1045 local_flux += ucont[k][j][i].x;
1046 }
1047 }
1048 } break;
1049
1050 case BC_FACE_NEG_Y:
1051 case BC_FACE_POS_Y: {
1052 PetscInt j = (face_id == BC_FACE_NEG_Y) ? ys : my - 2;
1053 for (PetscInt k = lzs; k < lze; k++) {
1054 for (PetscInt i = lxs; i < lxe; i++) {
1055 local_flux += ucont[k][j][i].y;
1056 }
1057 }
1058 } break;
1059
1060 case BC_FACE_NEG_Z:
1061 case BC_FACE_POS_Z: {
1062 PetscInt k = (face_id == BC_FACE_NEG_Z) ? zs : mz - 2;
1063 for (PetscInt j = lys; j < lye; j++) {
1064 for (PetscInt i = lxs; i < lxe; i++) {
1065 local_flux += ucont[k][j][i].z;
1066 }
1067 }
1068 } break;
1069 }
1070
1071 ierr = DMDAVecRestoreArrayRead(user->fda, user->Ucont, (const Cmpnts***)&ucont); CHKERRQ(ierr);
1072
1073 *local_inflow_contribution += local_flux;
1074
1075 LOG_ALLOW(LOCAL, LOG_DEBUG, "PostStep_InletParabolicProfile: Face %d, flux = %.6e\n",
1076 face_id, local_flux);
1077
1078 PetscFunctionReturn(0);
1079}
1080
1081
1082#undef __FUNCT__
1083#define __FUNCT__ "Destroy_InletParabolicProfile"
1084/**
1085 * @brief Internal helper implementation: `Destroy_InletParabolicProfile()`.
1086 * @details Local to this translation unit.
1087 */
1088static PetscErrorCode Destroy_InletParabolicProfile(BoundaryCondition *self)
1089{
1090 PetscFunctionBeginUser;
1091 if (self && self->data) {
1092 PetscFree(self->data);
1093 self->data = NULL;
1094 }
1095 PetscFunctionReturn(0);
1096}
1097
1098//================================================================================
1099//
1100// HANDLER IMPLEMENTATION: OUTLET WITH MASS CONSERVATION
1101// (Corresponds to BC_HANDLER_OUTLET_CONSERVATION)
1102//
1103// This handler ensures that the total flux leaving through outlet boundaries
1104// balances the total flux entering through inlet and far-field boundaries.
1105//
1106// Workflow:
1107// 1. PreStep: Measures the *uncorrected* flux based on interior velocities.
1108// 2. Apply: Calculates a global correction factor based on the flux imbalance
1109// and applies it to the contravariant velocity (ucont) on the outlet face.
1110// 3. PostStep: Measures the *corrected* flux for verification and logging.
1111//
1112//================================================================================
1113
1114// --- 1. FORWARD DECLARATIONS ---
1115static PetscErrorCode PreStep_OutletConservation(BoundaryCondition *self, BCContext *ctx,
1116 PetscReal *local_inflow_contribution, PetscReal *local_outflow_contribution);
1117static PetscErrorCode Apply_OutletConservation(BoundaryCondition *self, BCContext *ctx);
1118static PetscErrorCode PostStep_OutletConservation(BoundaryCondition *self, BCContext *ctx,
1119 PetscReal *in, PetscReal *out);
1120
1121#undef __FUNCT__
1122#define __FUNCT__ "Create_OutletConservation"
1123/**
1124 * @brief Implementation of \ref Create_OutletConservation().
1125 * @details Full API contract (arguments, ownership, side effects) is documented with
1126 * the header declaration in `include/BC_Handlers.h`.
1127 * @see Create_OutletConservation()
1128 */
1129PetscErrorCode Create_OutletConservation(BoundaryCondition *bc)
1130{
1131 PetscFunctionBeginUser;
1132
1133 if (!bc) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Input BoundaryCondition is NULL");
1134
1135 // This handler has the highest priority to ensure it runs after
1136 // all inflow fluxes have been calculated.
1137 bc->priority = BC_PRIORITY_OUTLET;
1138
1139 // Assign function pointers
1140 bc->Initialize = NULL; // No initialization needed
1141 bc->PreStep = PreStep_OutletConservation;
1142 bc->Apply = Apply_OutletConservation;
1143 bc->PostStep = PostStep_OutletConservation;
1144 bc->UpdateUbcs = NULL;
1145 bc->Destroy = NULL; // No private data to destroy
1146
1147 bc->data = NULL;
1148
1149 PetscFunctionReturn(0);
1150}
1151
1152#undef __FUNCT__
1153#define __FUNCT__ "PreStep_OutletConservation"
1154/**
1155 * @brief Internal helper implementation: `PreStep_OutletConservation()`.
1156 * @details Local to this translation unit.
1157 */
1158static PetscErrorCode PreStep_OutletConservation(BoundaryCondition *self, BCContext *ctx,
1159 PetscReal *local_inflow_contribution, PetscReal *local_outflow_contribution)
1160{
1161 PetscErrorCode ierr;
1162 UserCtx* user = ctx->user;
1163 BCFace face_id = ctx->face_id;
1164 DMDALocalInfo* info = &user->info;
1165 PetscBool can_service;
1166
1167 // Suppress unused parameter warnings for clarity.
1168 (void)self;
1169 (void)local_inflow_contribution;
1170
1171 PetscFunctionBeginUser;
1172
1173 // Step 1: Use the robust utility function to determine if this MPI rank owns a computable
1174 // portion of the specified boundary face. If not, there is no work to do, so we exit immediately.
1175 const PetscInt IM_nodes_global = user->IM;
1176 const PetscInt JM_nodes_global = user->JM;
1177 const PetscInt KM_nodes_global = user->KM;
1178 ierr = CanRankServiceFace(info, IM_nodes_global, JM_nodes_global, KM_nodes_global, face_id, &can_service); CHKERRQ(ierr);
1179
1180 if (!can_service) {
1181 PetscFunctionReturn(0);
1182 }
1183
1184 // Step 2: Get read-only access to the necessary PETSc arrays.
1185 // We use the local versions (`lUcat`, `lNvert`) which include ghost cell data,
1186 // ensuring we have the correct interior values adjacent to the boundary.
1187 Cmpnts ***ucat, ***csi, ***eta, ***zet;
1188 PetscReal ***nvert;
1189 ierr = DMDAVecGetArrayRead(user->fda, user->lUcat, (const Cmpnts***)&ucat); CHKERRQ(ierr);
1190 ierr = DMDAVecGetArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
1191 ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
1192 ierr = DMDAVecGetArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
1193 ierr = DMDAVecGetArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
1194
1195 PetscReal local_flux_out = 0.0;
1196 const PetscInt xs=info->xs, xe=info->xs+info->xm;
1197 const PetscInt ys=info->ys, ye=info->ys+info->ym;
1198 const PetscInt zs=info->zs, ze=info->zs+info->zm;
1199 const PetscInt mx=info->mx, my=info->my, mz=info->mz;
1200
1201 // Step 3: Replicate the legacy shrunk loop bounds to exclude corners and edges.
1202 PetscInt lxs = xs; if (xs == 0) lxs = xs + 1;
1203 PetscInt lxe = xe; if (xe == mx) lxe = xe - 1;
1204 PetscInt lys = ys; if (ys == 0) lys = ys + 1;
1205 PetscInt lye = ye; if (ye == my) lye = ye - 1;
1206 PetscInt lzs = zs; if (zs == 0) lzs = zs + 1;
1207 PetscInt lze = ze; if (ze == mz) lze = ze - 1;
1208
1209 // Step 4: Loop over the specified face using the corrected bounds and indexing to calculate flux.
1210 switch (face_id) {
1211 case BC_FACE_NEG_X: {
1212 const PetscInt i_cell = xs + 1; // Index for first interior cell-centered data
1213 const PetscInt i_face = xs; // Index for the -X face of that cell
1214 for (int k=lzs; k<lze; k++) for (int j=lys; j<lye; j++) {
1215 if (nvert[k][j][i_cell] < 0.1) {
1216 local_flux_out += (ucat[k][j][i_cell].x * csi[k][j][i_face].x + ucat[k][j][i_cell].y * csi[k][j][i_face].y + ucat[k][j][i_cell].z * csi[k][j][i_face].z);
1217 }
1218 }
1219 break;
1220 }
1221 case BC_FACE_POS_X: {
1222 const PetscInt i_cell = xe - 2; // Index for last interior cell-centered data
1223 const PetscInt i_face = xe - 2; // Index for the +X face of that cell
1224 for (int k=lzs; k<lze; k++) for (int j=lys; j<lye; j++) {
1225 if (nvert[k][j][i_cell] < 0.1) {
1226 local_flux_out += (ucat[k][j][i_cell].x * csi[k][j][i_face].x + ucat[k][j][i_cell].y * csi[k][j][i_face].y + ucat[k][j][i_cell].z * csi[k][j][i_face].z);
1227 }
1228 }
1229 break;
1230 }
1231 case BC_FACE_NEG_Y: {
1232 const PetscInt j_cell = ys + 1;
1233 const PetscInt j_face = ys;
1234 for (int k=lzs; k<lze; k++) for (int i=lxs; i<lxe; i++) {
1235 if (nvert[k][j_cell][i] < 0.1) {
1236 local_flux_out += (ucat[k][j_cell][i].x * eta[k][j_face][i].x + ucat[k][j_cell][i].y * eta[k][j_face][i].y + ucat[k][j_cell][i].z * eta[k][j_face][i].z);
1237 }
1238 }
1239 break;
1240 }
1241 case BC_FACE_POS_Y: {
1242 const PetscInt j_cell = ye - 2;
1243 const PetscInt j_face = ye - 2;
1244 for (int k=lzs; k<lze; k++) for (int i=lxs; i<lxe; i++) {
1245 if (nvert[k][j_cell][i] < 0.1) {
1246 local_flux_out += (ucat[k][j_cell][i].x * eta[k][j_face][i].x + ucat[k][j_cell][i].y * eta[k][j_face][i].y + ucat[k][j_cell][i].z * eta[k][j_face][i].z);
1247 }
1248 }
1249 break;
1250 }
1251 case BC_FACE_NEG_Z: {
1252 const PetscInt k_cell = zs + 1;
1253 const PetscInt k_face = zs;
1254 for (int j=lys; j<lye; j++) for (int i=lxs; i<lxe; i++) {
1255 if (nvert[k_cell][j][i] < 0.1) {
1256 local_flux_out += (ucat[k_cell][j][i].x * zet[k_face][j][i].x + ucat[k_cell][j][i].y * zet[k_face][j][i].y + ucat[k_cell][j][i].z * zet[k_face][j][i].z);
1257 }
1258 }
1259 break;
1260 }
1261 case BC_FACE_POS_Z: {
1262 const PetscInt k_cell = ze - 2;
1263 const PetscInt k_face = ze - 2;
1264 for (int j=lys; j<lye; j++) for (int i=lxs; i<lxe; i++) {
1265 if (nvert[k_cell][j][i] < 0.1) {
1266 local_flux_out += (ucat[k_cell][j][i].x * zet[k_face][j][i].x + ucat[k_cell][j][i].y * zet[k_face][j][i].y + ucat[k_cell][j][i].z * zet[k_face][j][i].z);
1267 }
1268 }
1269 break;
1270 }
1271 }
1272
1273 // Step 5: Restore the PETSc arrays.
1274 ierr = DMDAVecRestoreArrayRead(user->fda, user->lUcat, (const Cmpnts***)&ucat); CHKERRQ(ierr);
1275 ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
1276 ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
1277 ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
1278 ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
1279
1280 // Step 6: Add this face's calculated flux to the accumulator for this rank.
1281 *local_outflow_contribution += local_flux_out;
1282
1283 PetscFunctionReturn(0);
1284}
1285
1286#undef __FUNCT__
1287#define __FUNCT__ "Apply_OutletConservation"
1288/**
1289 * @brief (Handler Action) Applies mass conservation correction to the outlet face.
1290 *
1291 * This function calculates a global correction factor based on the total inflow and outflow fluxes
1292 * and applies it to the contravariant velocity (`ucont`) on the outlet face to ensure mass conservation.
1293 */
1294static PetscErrorCode Apply_OutletConservation(BoundaryCondition *self, BCContext *ctx)
1295{
1296 PetscErrorCode ierr;
1297 (void)self;
1298 UserCtx* user = ctx->user;
1299 BCFace face_id = ctx->face_id;
1300 DMDALocalInfo* info = &user->info;
1301 PetscBool can_service;
1302
1303 PetscFunctionBeginUser;
1304 PROFILE_FUNCTION_BEGIN;
1305
1306 const PetscInt IM_nodes_global = user->IM;
1307 const PetscInt JM_nodes_global = user->JM;
1308 const PetscInt KM_nodes_global = user->KM;
1309 ierr = CanRankServiceFace(info, IM_nodes_global, JM_nodes_global, KM_nodes_global, face_id, &can_service); CHKERRQ(ierr);
1310
1311 if (!can_service) {
1312 PROFILE_FUNCTION_END;
1313 PetscFunctionReturn(0);
1314 }
1315
1316 // --- STEP 1: Calculate the correction factor using pre-calculated area ---
1317 PetscReal total_inflow = *ctx->global_inflow_sum + *ctx->global_farfield_inflow_sum;
1318 PetscReal flux_imbalance = total_inflow - *ctx->global_outflow_sum;
1319
1320 // Directly use the pre-calculated area from the simulation context.
1321 PetscReal velocity_correction = (PetscAbsReal(user->simCtx->AreaOutSum) > 1e-12)
1322 ? flux_imbalance / user->simCtx->AreaOutSum
1323 : 0.0;
1324
1325 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Outlet Correction on Face %d: Imbalance=%.4e, Pre-calc Area=%.4e, V_corr=%.4e\n",
1326 face_id, flux_imbalance, user->simCtx->AreaOutSum, velocity_correction);
1327
1328 // --- STEP 2: Apply the correction to ucont on the outlet face ---
1329
1330 // Get read/write access to necessary arrays
1331
1332 Cmpnts ***ubcs, ***ucont, ***csi, ***eta, ***zet, ***ucat;
1333 PetscReal ***nvert;
1334 ierr = DMDAVecGetArray(user->fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
1335 ierr = DMDAVecGetArray(user->fda, user->Ucont, &ucont); CHKERRQ(ierr);
1336 ierr = DMDAVecGetArrayRead(user->fda,user->lUcat, (const Cmpnts***)&ucat); CHKERRQ(ierr);
1337 ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
1338 ierr = DMDAVecGetArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
1339 ierr = DMDAVecGetArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
1340 ierr = DMDAVecGetArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
1341
1342 // Get local grid bounds to exclude corners/edges
1343 PetscInt xs = info->xs, xe = info->xs + info->xm;
1344 PetscInt ys = info->ys, ye = info->ys + info->ym;
1345 PetscInt zs = info->zs, ze = info->zs + info->zm;
1346 PetscInt mx = info->mx, my = info->my, mz = info->mz;
1347 PetscInt lxs = xs, lxe = xe, lys = ys, lye = ye, lzs = zs, lze = ze;
1348
1349 if (xs == 0) lxs = xs + 1;
1350 if (xe == mx) lxe = xe - 1;
1351 if (ys == 0) lys = ys + 1;
1352 if (ye == my) lye = ye - 1;
1353 if (zs == 0) lzs = zs + 1;
1354 if (ze == mz) lze = ze - 1;
1355
1356 // Loop over faces and apply correction
1357 switch(face_id){
1358 case BC_FACE_NEG_X:{
1359 const PetscInt i_cell = xs + 1;
1360 const PetscInt i_face = xs;
1361 const PetscInt i_dummy = xs;
1362 for (PetscInt k = lzs; k < lze; k++) {
1363 for (PetscInt j = lys; j < lye; j++) {
1364 if (nvert[k][j][i_cell] < 0.1) {
1365 // Set ubcs
1366 ubcs[k][j][i_dummy] = ucat[k][j][i_cell];
1367
1368 // Calculate Local uncorrected original flux
1369 PetscReal Uncorrected_local_flux = (ubcs[k][j][i_dummy].x * csi[k][j][i_face].x) + (ubcs[k][j][i_dummy].y * csi[k][j][i_face].y) + (ubcs[k][j][i_dummy].z * csi[k][j][i_face].z);
1370
1371 PetscReal Cell_Area = sqrt((csi[k][j][i_face].x*csi[k][j][i_face].x) + (csi[k][j][i_face].y*csi[k][j][i_face].y) + (csi[k][j][i_face].z*csi[k][j][i_face].z));
1372
1373 PetscReal Correction_flux = velocity_correction*Cell_Area;
1374
1375 ucont[k][j][i_face].x = Uncorrected_local_flux + Correction_flux;
1376 }
1377 }
1378 }
1379 break;
1380 }
1381 case BC_FACE_POS_X:{
1382 const PetscInt i_cell = xe - 2;
1383 const PetscInt i_face = xe - 2;
1384 const PetscInt i_dummy = xe - 1;
1385 for(PetscInt k = lzs; k < lze; k++) for (PetscInt j = lys; j < lye; j++){
1386 if(nvert[k][j][i_cell]<0.1){
1387 // Set ubcs
1388 ubcs[k][j][i_dummy] = ucat[k][j][i_cell];
1389
1390 // Calculate Local uncorrected original flux
1391 PetscReal Uncorrected_local_flux = (ubcs[k][j][i_dummy].x * csi[k][j][i_face].x) + (ubcs[k][j][i_dummy].y * csi[k][j][i_face].y) + (ubcs[k][j][i_dummy].z * csi[k][j][i_face].z);
1392
1393 PetscReal Cell_Area = sqrt((csi[k][j][i_face].x*csi[k][j][i_face].x) + (csi[k][j][i_face].y*csi[k][j][i_face].y) + (csi[k][j][i_face].z*csi[k][j][i_face].z));
1394
1395 PetscReal Correction_flux = velocity_correction*Cell_Area;
1396
1397 ucont[k][j][i_face].x = Uncorrected_local_flux + Correction_flux;
1398 }
1399 }
1400 break;
1401 }
1402 case BC_FACE_NEG_Y:{
1403 const PetscInt j_cell = ys + 1;
1404 const PetscInt j_face = ys;
1405 const PetscInt j_dummy = ys;
1406 for(PetscInt k = lzs; k < lze; k++) for (PetscInt i = lxs; i < lxe; i++){
1407 if(nvert[k][j_cell][i]<0.1){
1408 // Set ubcs
1409 ubcs[k][j_dummy][i] = ucat[k][j_cell][i];
1410
1411 // Calculate Local uncorrected original flux
1412 PetscReal Uncorrected_local_flux = (ubcs[k][j_dummy][i].x*eta[k][j_face][i].x) + (ubcs[k][j_dummy][i].y*eta[k][j_face][i].y) + (ubcs[k][j_dummy][i].z*eta[k][j_face][i].z);
1413
1414 PetscReal Cell_Area = sqrt((eta[k][j_face][i].x*eta[k][j_face][i].x)+(eta[k][j_face][i].y*eta[k][j_face][i].y)+(eta[k][j_face][i].z*eta[k][j_face][i].z));
1415
1416 PetscReal Correction_flux = velocity_correction*Cell_Area;
1417
1418 ucont[k][j_face][i].y = Uncorrected_local_flux + Correction_flux;
1419 }
1420 }
1421 break;
1422 }
1423 case BC_FACE_POS_Y:{
1424 const PetscInt j_cell = ye - 2;
1425 const PetscInt j_face = ye - 2;
1426 const PetscInt j_dummy = ye - 1;
1427 for(PetscInt k = lzs; k < lze; k++) for (PetscInt i = lxs; i < lxe; i++){
1428 if(nvert[k][j_cell][i]<0.1){
1429 // Set ubcs
1430 ubcs[k][j_dummy][i] = ucat[k][j_cell][i];
1431
1432 // Calculate Local uncorrected original flux
1433 PetscReal Uncorrected_local_flux = (ubcs[k][j_dummy][i].x*eta[k][j_face][i].x) + (ubcs[k][j_dummy][i].y*eta[k][j_face][i].y) + (ubcs[k][j_dummy][i].z*eta[k][j_face][i].z);
1434
1435 PetscReal Cell_Area = sqrt((eta[k][j_face][i].x*eta[k][j_face][i].x)+(eta[k][j_face][i].y*eta[k][j_face][i].y)+(eta[k][j_face][i].z*eta[k][j_face][i].z));
1436
1437 PetscReal Correction_flux = velocity_correction*Cell_Area;
1438
1439 ucont[k][j_face][i].y = Uncorrected_local_flux + Correction_flux;
1440 }
1441 }
1442 break;
1443 }
1444 case BC_FACE_NEG_Z:{
1445 const PetscInt k_cell = zs + 1;
1446 const PetscInt k_face = zs;
1447 const PetscInt k_dummy = zs;
1448 for(PetscInt j = lys; j < lye; j++) for (PetscInt i = lxs; i < lxe; i++){
1449 if(nvert[k_cell][j][i]<0.1){
1450 // Set ubcs
1451 ubcs[k_dummy][j][i] = ucat[k_cell][j][i];
1452
1453 // Calculate Local uncorrected original flux
1454 PetscReal Uncorrected_local_flux = ((ubcs[k_dummy][j][i].x*zet[k_face][j][i].x) + (ubcs[k_dummy][j][i].y*zet[k_face][j][i].y) + (ubcs[k_dummy][j][i].z*zet[k_face][j][i].z));
1455
1456 PetscReal Cell_Area = sqrt((zet[k_face][j][i].x*zet[k_face][j][i].x)+(zet[k_face][j][i].y*zet[k_face][j][i].y)+(zet[k_face][j][i].z*zet[k_face][j][i].z));
1457
1458 PetscReal Correction_flux = velocity_correction*Cell_Area;
1459
1460 ucont[k_face][j][i].z = Uncorrected_local_flux + Correction_flux;
1461 }
1462 }
1463 break;
1464 }
1465 case BC_FACE_POS_Z:{
1466 const PetscInt k_cell = ze - 2;
1467 const PetscInt k_face = ze - 2;
1468 const PetscInt k_dummy = ze - 1;
1469 for(PetscInt j = lys; j < lye; j++) for (PetscInt i = lxs; i < lxe; i++){
1470 if(nvert[k_cell][j][i]<0.1){
1471 // Set ubcs
1472 ubcs[k_dummy][j][i] = ucat[k_cell][j][i];
1473
1474 // Calculate Local uncorrected original flux
1475 PetscReal Uncorrected_local_flux = ((ubcs[k_dummy][j][i].x*zet[k_face][j][i].x) + (ubcs[k_dummy][j][i].y*zet[k_face][j][i].y) + (ubcs[k_dummy][j][i].z*zet[k_face][j][i].z));
1476
1477 PetscReal Cell_Area = sqrt((zet[k_face][j][i].x*zet[k_face][j][i].x)+(zet[k_face][j][i].y*zet[k_face][j][i].y)+(zet[k_face][j][i].z*zet[k_face][j][i].z));
1478
1479 PetscReal Correction_flux = velocity_correction*Cell_Area;
1480
1481 ucont[k_face][j][i].z = Uncorrected_local_flux + Correction_flux;
1482 }
1483 }
1484 break;
1485 }
1486 }
1487
1488 // Restore all arrays
1489 ierr = DMDAVecRestoreArray(user->fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
1490 ierr = DMDAVecRestoreArray(user->fda, user->Ucont, &ucont); CHKERRQ(ierr);
1491 ierr = DMDAVecRestoreArrayRead(user->fda,user->lUcat, (const Cmpnts***)&ucat); CHKERRQ(ierr);
1492 ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
1493 ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
1494 ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
1495 ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
1496
1497 PROFILE_FUNCTION_END;
1498 PetscFunctionReturn(0);
1499}
1500
1501#undef __FUNCT__
1502#define __FUNCT__ "PostStep_OutletConservation"
1503/**
1504 * @brief Internal helper implementation: `PostStep_OutletConservation()`.
1505 * @details Local to this translation unit.
1506 */
1507static PetscErrorCode PostStep_OutletConservation(BoundaryCondition *self, BCContext *ctx,
1508 PetscReal *local_inflow_contribution,
1509 PetscReal *local_outflow_contribution)
1510{
1511 PetscErrorCode ierr;
1512 UserCtx* user = ctx->user;
1513 BCFace face_id = ctx->face_id;
1514 DMDALocalInfo* info = &user->info;
1515 PetscBool can_service;
1516
1517 (void)self;
1518 (void)local_inflow_contribution;
1519
1520 PetscFunctionBeginUser;
1521 const PetscInt IM_nodes_global = user->IM;
1522 const PetscInt JM_nodes_global = user->JM;
1523 const PetscInt KM_nodes_global = user->KM;
1524 ierr = CanRankServiceFace(info, IM_nodes_global, JM_nodes_global, KM_nodes_global, face_id, &can_service); CHKERRQ(ierr);
1525
1526 if (!can_service) PetscFunctionReturn(0);
1527
1528 // Get arrays (need both ucont and nvert)
1529 Cmpnts ***ucont;
1530 PetscReal ***nvert; // ✅ ADD nvert
1531
1532 ierr = DMDAVecGetArrayRead(user->fda, user->Ucont, (const Cmpnts***)&ucont); CHKERRQ(ierr);
1533 ierr = DMDAVecGetArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr); // ✅ ADD
1534
1535 PetscReal local_flux = 0.0;
1536
1537 PetscInt xs = info->xs, xe = info->xs + info->xm;
1538 PetscInt ys = info->ys, ye = info->ys + info->ym;
1539 PetscInt zs = info->zs, ze = info->zs + info->zm;
1540 PetscInt mx = info->mx, my = info->my, mz = info->mz;
1541
1542 PetscInt lxs = xs, lxe = xe, lys = ys, lye = ye, lzs = zs, lze = ze;
1543 if (xs == 0) lxs = xs + 1;
1544 if (xe == mx) lxe = xe - 1;
1545 if (ys == 0) lys = ys + 1;
1546 if (ye == my) lye = ye - 1;
1547 if (zs == 0) lzs = zs + 1;
1548 if (ze == mz) lze = ze - 1;
1549
1550 // Sum ucont components, skipping solid cells (same indices as PreStep)
1551 switch (face_id) {
1552 case BC_FACE_NEG_X: {
1553 const PetscInt i_cell = xs + 1; // ✅ Match PreStep
1554 const PetscInt i_face = xs;
1555 for (PetscInt k = lzs; k < lze; k++) {
1556 for (PetscInt j = lys; j < lye; j++) {
1557 if (nvert[k][j][i_cell] < 0.1) { // ✅ Skip solid cells
1558 local_flux += ucont[k][j][i_face].x;
1559 }
1560 }
1561 }
1562 } break;
1563
1564 case BC_FACE_POS_X: {
1565 const PetscInt i_cell = xe - 2; // ✅ Match PreStep
1566 const PetscInt i_face = xe - 2;
1567 for (PetscInt k = lzs; k < lze; k++) {
1568 for (PetscInt j = lys; j < lye; j++) {
1569 if (nvert[k][j][i_cell] < 0.1) { // ✅ Skip solid cells
1570 local_flux += ucont[k][j][i_face].x;
1571 }
1572 }
1573 }
1574 } break;
1575
1576 case BC_FACE_NEG_Y: {
1577 const PetscInt j_cell = ys + 1; // ✅ Match PreStep
1578 const PetscInt j_face = ys;
1579 for (PetscInt k = lzs; k < lze; k++) {
1580 for (PetscInt i = lxs; i < lxe; i++) {
1581 if (nvert[k][j_cell][i] < 0.1) { // ✅ Skip solid cells
1582 local_flux += ucont[k][j_face][i].y;
1583 }
1584 }
1585 }
1586 } break;
1587
1588 case BC_FACE_POS_Y: {
1589 const PetscInt j_cell = ye - 2; // ✅ Match PreStep
1590 const PetscInt j_face = ye - 2;
1591 for (PetscInt k = lzs; k < lze; k++) {
1592 for (PetscInt i = lxs; i < lxe; i++) {
1593 if (nvert[k][j_cell][i] < 0.1) { // ✅ Skip solid cells
1594 local_flux += ucont[k][j_face][i].y;
1595 }
1596 }
1597 }
1598 } break;
1599
1600 case BC_FACE_NEG_Z: {
1601 const PetscInt k_cell = zs + 1; // ✅ Match PreStep
1602 const PetscInt k_face = zs;
1603 for (PetscInt j = lys; j < lye; j++) {
1604 for (PetscInt i = lxs; i < lxe; i++) {
1605 if (nvert[k_cell][j][i] < 0.1) { // ✅ Skip solid cells
1606 local_flux += ucont[k_face][j][i].z;
1607 }
1608 }
1609 }
1610 } break;
1611
1612 case BC_FACE_POS_Z: {
1613 const PetscInt k_cell = ze - 2; // ✅ Match PreStep
1614 const PetscInt k_face = ze - 2;
1615 for (PetscInt j = lys; j < lye; j++) {
1616 for (PetscInt i = lxs; i < lxe; i++) {
1617 if (nvert[k_cell][j][i] < 0.1) { // ✅ Skip solid cells
1618 local_flux += ucont[k_face][j][i].z;
1619 }
1620 }
1621 }
1622 } break;
1623 }
1624
1625 // Restore arrays
1626 ierr = DMDAVecRestoreArrayRead(user->fda, user->Ucont, (const Cmpnts***)&ucont); CHKERRQ(ierr);
1627 ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr); // ✅ ADD
1628
1629 // Add to accumulator
1630 *local_outflow_contribution += local_flux;
1631
1632 LOG_ALLOW(LOCAL, LOG_DEBUG, "PostStep_OutletConservation: Face %d, corrected flux = %.6e\n",
1633 face_id, local_flux);
1634
1635 PetscFunctionReturn(0);
1636}
1637
1638
1639/**
1640 * @brief Implementation of \ref Create_PeriodicGeometric().
1641 * @details Full API contract (arguments, ownership, side effects) is documented with
1642 * the header declaration in `include/BC_Handlers.h`.
1643 * @see Create_PeriodicGeometric()
1644 */
1645
1646PetscErrorCode Create_PeriodicGeometric(BoundaryCondition *bc){
1647 PetscFunctionBeginUser;
1648
1649 if (!bc) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Input BoundaryCondition is NULL");
1650 bc->priority = BC_PRIORITY_WALL;
1651
1652 // Assign function pointers
1653 bc->Initialize = NULL; // No initialization needed
1654 bc->PreStep = NULL;
1655 bc->Apply = NULL;
1656 bc->PostStep = NULL;
1657 bc->UpdateUbcs = NULL;
1658 bc->Destroy = NULL; // No private data to destroy
1659
1660 bc->data = NULL;
1661
1662 PetscFunctionReturn(0);
1663}
1664
1665
1666// ===============================================================================
1667//
1668// HANDLER IMPLEMENTATION: PERIODIC DRIVEN CONSTANT FLUX
1669// (Corresponds to BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX)
1670//
1671// ===============================================================================
1672
1673// --- 1. FORWARD DECLARATIONS & PRIVATE DATA ---
1674
1675// Forward declarations for the static functions that implement this handler's behavior.
1676static PetscErrorCode Initialize_PeriodicDrivenConstant(BoundaryCondition *self, BCContext *ctx);
1677static PetscErrorCode PreStep_PeriodicDrivenConstant(BoundaryCondition *self, BCContext *ctx, PetscReal *in, PetscReal *out);
1678static PetscErrorCode Apply_PeriodicDrivenConstant(BoundaryCondition *self, BCContext *ctx);
1679static PetscErrorCode Destroy_PeriodicDrivenConstant(BoundaryCondition *self);
1680
1681/** @brief Private data structure for the handler. */
1682typedef struct {
1683 char direction; // 'X', 'Y', or 'Z', determined at initialization.
1684 PetscReal targetVolumetricFlux; // The constant target flux, parsed from parameters.
1685 PetscReal boundaryVelocityCorrection; // The "Boundary Trim" value for the Apply() step.
1686 PetscBool isMasterController; // Flag: PETSC_TRUE only for the handler on the negative face.
1687 PetscBool applyBoundaryTrim; // Flag: PETSC_TRUE if applying Boundary trim on ucont.
1688} DrivenConstantData;
1689
1690
1691// --- 2. HANDLER CONSTRUCTOR ---
1692
1693#undef __FUNCT__
1694#define __FUNCT__ "Create_PeriodicDrivenConstant"
1695/**
1696 * @brief Internal helper implementation: `Create_PeriodicDrivenConstant()`.
1697 * @details Local to this translation unit.
1698 */
1699PetscErrorCode Create_PeriodicDrivenConstant(BoundaryCondition *bc)
1700{
1701 PetscErrorCode ierr;
1702 PetscFunctionBeginUser;
1703
1704 if (!bc) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Input BoundaryCondition object is NULL in Create_PeriodicDrivenConstantFlux");
1705
1706 // --- Allocate the private data structure ---
1707 DrivenConstantData *data = NULL;
1708 ierr = PetscNew(&data); CHKERRQ(ierr);
1709 // Initialize fields to safe default values
1710 data->direction = ' ';
1711 data->targetVolumetricFlux = 0.0;
1712 data->boundaryVelocityCorrection = 0.0;
1713 data->isMasterController = PETSC_FALSE;
1714 data->applyBoundaryTrim = PETSC_FALSE;
1715
1716 // Attach the private data to the generic handler object
1717 bc->data = (void*)data;
1718
1719 // --- Configure the handler's properties and methods ---
1720
1721 // Set priority: Using BC_PRIORITY_INLET ensures this handler's PreStep runs
1722 // before other handlers (like outlets) that might depend on its calculations.
1723 // It is the caller's responsibility that there are no Inlets called along with driven periodic to avoid clash.
1724 bc->priority = BC_PRIORITY_INLET;
1725
1726 // Assign the function pointers to the implementations in this file.
1727 bc->Initialize = Initialize_PeriodicDrivenConstant;
1728 bc->PreStep = PreStep_PeriodicDrivenConstant;
1729 bc->Apply = Apply_PeriodicDrivenConstant;
1730 bc->PostStep = NULL; // This handler has no action after the main solver step.
1731 bc->UpdateUbcs = NULL; // The boundary value is not flow-dependent (it's periodic).
1732 bc->Destroy = Destroy_PeriodicDrivenConstant;
1733
1734 PetscFunctionReturn(0);
1735}
1736
1737#undef __FUNCT__
1738#define __FUNCT__ "Initialize_PeriodicDrivenConstant"
1739/**
1740 * @brief Internal helper implementation: `Initialize_PeriodicDrivenConstant()`.
1741 * @details Local to this translation unit.
1742 */
1743static PetscErrorCode Initialize_PeriodicDrivenConstant(BoundaryCondition *self, BCContext *ctx)
1744{
1745 PetscErrorCode ierr;
1746 DrivenConstantData *data = (DrivenConstantData*)self->data;
1747 BCFace face_id = ctx->face_id;
1748 UserCtx* user = ctx->user;
1749
1750 PetscFunctionBeginUser;
1751
1752 LOG_ALLOW(LOCAL, LOG_DEBUG, "Initializing PERIODIC_DRIVEN_CONSTANT_FLUX handler on Face %s...\n", BCFaceToString(face_id));
1753
1754 // --- 1. Validation: Ensure the mathematical type is PERIODIC ---
1755 if (user->boundary_faces[face_id].mathematical_type != PERIODIC) {
1756 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT,
1757 "Configuration Error: Handler PERIODIC_DRIVEN_CONSTANT_FLUX on Face %s must be applied to a face with mathematical_type PERIODIC.",
1758 BCFaceToString(face_id));
1759 }
1760
1761 // --- 2. Role Assignment: Determine direction and master status ---
1762 data->isMasterController = PETSC_FALSE;
1763 switch (face_id) {
1764 case BC_FACE_NEG_X: data->direction = 'X'; data->isMasterController = PETSC_TRUE; break;
1765 case BC_FACE_POS_X: data->direction = 'X'; break;
1766 case BC_FACE_NEG_Y: data->direction = 'Y'; data->isMasterController = PETSC_TRUE; break;
1767 case BC_FACE_POS_Y: data->direction = 'Y'; break;
1768 case BC_FACE_NEG_Z: data->direction = 'Z'; data->isMasterController = PETSC_TRUE; break;
1769 case BC_FACE_POS_Z: data->direction = 'Z'; break;
1770 }
1771
1772 // --- 3. Parameter Parsing (Master Controller only) ---
1773 if (data->isMasterController) {
1774 PetscBool found;
1775
1776 // Attempt to read the 'target_flux' parameter from the bcs.run file.
1777 ierr = GetBCParamReal(user->boundary_faces[face_id].params, "target_flux",
1778 &data->targetVolumetricFlux, &found); CHKERRQ(ierr);
1779
1780 // If the required parameter is not found, halt with an informative error.
1781 if (!found) {
1782 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT,
1783 "Configuration Error: Handler PERIODIC_DRIVEN_CONSTANT_FLUX on Face %s requires a 'target_flux' parameter in the bcs file (e.g., target_flux=10.0).",
1784 BCFaceToString(face_id));
1785 }
1786
1787 LOG_ALLOW(GLOBAL, LOG_INFO, "Driven Flow (Dir %c): Constant target volumetric flux set to %le.\n",
1788 data->direction, data->targetVolumetricFlux);
1789
1790 // Store the target flux in the UserCtx. This makes it globally accessible
1791 // to other parts of the solver, such as the `CorrectChannelFluxProfile` enforcer function.
1792 user->simCtx->targetVolumetricFlux = data->targetVolumetricFlux;
1793 }
1794
1795 PetscBool trimfound;
1796 // Attempt to read Trim flag.
1797 ierr = GetBCParamBool(user->boundary_faces[face_id].params, "apply_trim",
1798 &data->applyBoundaryTrim, &trimfound); CHKERRQ(ierr);
1799
1800 if(!trimfound) LOG_ALLOW(GLOBAL,LOG_DEBUG,"Trim Application not found,defaults to %s.\n",data->applyBoundaryTrim? "True":"False");
1801
1802 PetscFunctionReturn(0);
1803}
1804
1805#undef __FUNCT__
1806#define __FUNCT__ "PreStep_PeriodicDrivenConstant"
1807/**
1808 * @brief Internal helper implementation: `PreStep_PeriodicDrivenConstant()`.
1809 * @details Local to this translation unit.
1810 */
1811static PetscErrorCode PreStep_PeriodicDrivenConstant(BoundaryCondition *self, BCContext *ctx,
1812 PetscReal *local_inflow_contribution,
1813 PetscReal *local_outflow_contribution)
1814{
1815 PetscErrorCode ierr;
1816 DrivenConstantData *data = (DrivenConstantData*)self->data;
1817 UserCtx* user = ctx->user;
1818 SimCtx* simCtx = user->simCtx;
1819
1820 PetscFunctionBeginUser;
1821
1822 // --- Master Check: Only the handler on the negative face performs calculations ---
1823 if (!data->isMasterController) {
1824 PetscFunctionReturn(0);
1825 }
1826
1827 // This function is the generalized implementation of the old InitializeChannelFlowController.
1828 char direction = data->direction;
1829 DMDALocalInfo info = user->info;
1830 PetscInt i, j, k;
1831
1832 // --- Get read-only access to necessary field data ---
1833 Cmpnts ***ucont;
1834 PetscReal ***nvert;
1835 ierr = DMDAVecGetArrayRead(user->fda, user->lUcont, (const Cmpnts***)&ucont); CHKERRQ(ierr);
1836 ierr = DMDAVecGetArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
1837
1838 // --- Define local loop bounds ---
1839 PetscInt lxs = (info.xs == 0) ? 1 : info.xs;
1840 PetscInt lys = (info.ys == 0) ? 1 : info.ys;
1841 PetscInt lzs = (info.zs == 0) ? 1 : info.zs;
1842 PetscInt lxe = (info.xs + info.xm == info.mx) ? info.mx - 1 : info.xs + info.xm;
1843 PetscInt lye = (info.ys + info.ym == info.my) ? info.my - 1 : info.ys + info.ym;
1844 PetscInt lze = (info.zs + info.zm == info.mz) ? info.mz - 1 : info.zs + info.zm;
1845
1846 // ===================================================================================
1847 // CONTROLLER SENSORS: DUAL-FLUX MEASUREMENT STRATEGY
1848 //
1849 // To create a control system that is both stable and responsive, we measure the
1850 // volumetric flux in two distinct ways. One provides a fast, local signal for
1851 // immediate corrections, while the other provides a slow, stable signal for
1852 // strategic, domain-wide adjustments.
1853 //
1854 // -----------------------------------------------------------------------------------
1855 // ** 1. Fast/Local Sensor: `globalCurrentBoundaryFlux` **
1856 // -----------------------------------------------------------------------------------
1857 //
1858 // - WHAT IT IS: The total volumetric flux (m³/s) passing through the single,
1859 // representative periodic boundary plane at k=0.
1860 //
1861 // - PHYSICAL MEANING: An instantaneous "snapshot" of the flow rate at the
1862 // critical "seam" where the domain wraps around.
1863 //
1864 // - CHARACTERISTICS:
1865 // - Fast & Responsive: It immediately reflects any local flow changes or
1866 // errors occurring at the periodic interface.
1867 // - Noisy: It is susceptible to local fluctuations from turbulence or
1868 // numerical artifacts, causing its value to jitter from step to step.
1869 //
1870 // - CONTROLLER'S USE: This measurement is used to compute the
1871 // `boundaryVelocityCorrection`. This is a TACTICAL, fast-acting "trim" that is
1872 // applied immediately and directly to the boundary velocities to ensure perfect
1873 // continuity at the seam.
1874 //
1875 // -----------------------------------------------------------------------------------
1876 // ** 2. Stable/Global Sensor: `globalAveragePlanarVolumetricFlux` **
1877 // -----------------------------------------------------------------------------------
1878 //
1879 // - WHAT IT IS: The average of the volumetric fluxes across ALL cross-sectional
1880 // k-planes in the domain.
1881 // (i.e., [Flux(k=0) + Flux(k=1) + ... + Flux(k=N)] / N)
1882 //
1883 // - PHYSICAL MEANING: It represents the overall, bulk momentum of the fluid,
1884 // effectively filtering out local, transient fluctuations.
1885 //
1886 // - CHARACTERISTICS:
1887 // - Stable & Robust: Local noise at one plane is averaged out by all other
1888 // planes, providing a very smooth signal.
1889 // - Inertial: It changes more slowly, reflecting the inertia of the entire
1890 // fluid volume.
1891 //
1892 // - CONTROLLER'S USE: This measurement is used to compute the
1893 // `bulkVelocityCorrection`. This is a STRATEGIC, long-term adjustment used to
1894 // scale the main momentum source (body force). Using this stable signal prevents
1895 // the main driving force from oscillating wildly, ensuring simulation stability.
1896 //
1897 // ===================================================================================
1898
1899 // --- Initialize local accumulators ---
1900 PetscReal localCurrentBoundaryFlux = 0.0;
1901 PetscReal localAveragePlanarVolumetricFluxTerm = 0.0;
1902
1903 // --- Measure local contributions to the two flux types, generalized by direction ---
1904 switch (direction) {
1905 case 'X':
1906 if (info.xs == 0) { // Only the rank on the negative face contributes to boundary flux
1907 i = 0;
1908 for (k = lzs; k < lze; k++) for (j = lys; j < lye; j++) {
1909 if (nvert[k][j][i + 1] < 0.1) localCurrentBoundaryFlux += ucont[k][j][i].x;
1910 }
1911 }
1912 for (i = info.xs; i < lxe; i++) {
1913 for (k = lzs; k < lze; k++) for (j = lys; j < lye; j++) {
1914 if (nvert[k][j][i + 1] < 0.1) localAveragePlanarVolumetricFluxTerm += ucont[k][j][i].x / (PetscReal)(info.mx - 1);
1915 }
1916 }
1917 break;
1918 case 'Y':
1919 if (info.ys == 0) {
1920 j = 0;
1921 for (k = lzs; k < lze; k++) for (i = lxs; i < lxe; i++) {
1922 if (nvert[k][j + 1][i] < 0.1) localCurrentBoundaryFlux += ucont[k][j][i].y;
1923 }
1924 }
1925 for (j = info.ys; j < lye; j++) {
1926 for (k = lzs; k < lze; k++) for (i = lxs; i < lxe; i++) {
1927 if (nvert[k][j + 1][i] < 0.1) localAveragePlanarVolumetricFluxTerm += ucont[k][j][i].y / (PetscReal)(info.my - 1);
1928 }
1929 }
1930 break;
1931 case 'Z':
1932 if (info.zs == 0) {
1933 k = 0;
1934 for (j = lys; j < lye; j++) for (i = lxs; i < lxe; i++) {
1935 if (nvert[k + 1][j][i] < 0.1) localCurrentBoundaryFlux += ucont[k][j][i].z;
1936 }
1937 }
1938 for (k = info.zs; k < lze; k++) {
1939 for (j = lys; j < lye; j++) for (i = lxs; i < lxe; i++) {
1940 if (nvert[k + 1][j][i] < 0.1) localAveragePlanarVolumetricFluxTerm += ucont[k][j][i].z / (PetscReal)(info.mz - 1);
1941 }
1942 }
1943 break;
1944 }
1945
1946 // --- Release array access as soon as possible ---
1947 ierr = DMDAVecRestoreArrayRead(user->fda, user->lUcont, (const Cmpnts***)&ucont); CHKERRQ(ierr);
1948 ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
1949
1950 // --- Get cross-sectional area using the dedicated geometry function ---
1951 Cmpnts ignored_center;
1952 PetscReal globalBoundaryArea;
1953 BCFace neg_face_id = (direction == 'X') ? BC_FACE_NEG_X : (direction == 'Y') ? BC_FACE_NEG_Y : BC_FACE_NEG_Z;
1954 ierr = CalculateFaceCenterAndArea(user, neg_face_id, &ignored_center, &globalBoundaryArea); CHKERRQ(ierr);
1955
1956 // --- Perform global reductions to get the final flux values ---
1957 PetscReal globalCurrentBoundaryFlux, globalAveragePlanarVolumetricFlux;
1958 ierr = MPI_Allreduce(&localCurrentBoundaryFlux, &globalCurrentBoundaryFlux, 1, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
1959 ierr = MPI_Allreduce(&localAveragePlanarVolumetricFluxTerm, &globalAveragePlanarVolumetricFlux, 1, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
1960
1961 // --- Calculate the two correction terms ---
1962 if (globalBoundaryArea > 1.0e-12) {
1963 // The main correction for the body force is based on the STABLE domain-averaged flux.
1964 simCtx->bulkVelocityCorrection = (data->targetVolumetricFlux - globalAveragePlanarVolumetricFlux) / globalBoundaryArea;
1965
1966 // The immediate correction for the boundary trim is based on the FAST boundary-specific flux.
1967 data->boundaryVelocityCorrection = (data->targetVolumetricFlux - globalCurrentBoundaryFlux) / globalBoundaryArea;
1968 } else {
1969 simCtx->bulkVelocityCorrection = 0.0;
1970 data->boundaryVelocityCorrection = 0.0;
1971 }
1972
1973 LOG_ALLOW(GLOBAL, LOG_INFO, "Driven Flow Controller Update (Dir %c):\n", data->direction);
1974 LOG_ALLOW(GLOBAL, LOG_INFO, " - Target Volumetric Flux: %.6e\n", data->targetVolumetricFlux);
1975 LOG_ALLOW(GLOBAL, LOG_INFO, " - Avg Planar Volumetric Flux (Stable): %.6e\n", globalAveragePlanarVolumetricFlux);
1976 LOG_ALLOW(GLOBAL, LOG_INFO, " - Boundary Flux (Fast): %.6e\n", globalCurrentBoundaryFlux);
1977 LOG_ALLOW(GLOBAL, LOG_INFO, " - Bulk Velocity Correction: %.6e (For Momentum Source)\n", simCtx->bulkVelocityCorrection);
1978 LOG_ALLOW(GLOBAL, LOG_INFO, " - Boundary Velocity Correction: %.6e (For Boundary Trim)\n", data->boundaryVelocityCorrection);
1979
1980 // Suppress unused parameter warnings for this handler
1981 (void)local_inflow_contribution;
1982 (void)local_outflow_contribution;
1983
1984 PetscFunctionReturn(0);
1985}
1986
1987#undef __FUNCT__
1988#define __FUNCT__ "Apply_PeriodicDrivenConstant"
1989/**
1990 * @brief Internal helper implementation: `Apply_PeriodicDrivenConstant()`.
1991 * @details Local to this translation unit.
1992 */
1993static PetscErrorCode Apply_PeriodicDrivenConstant(BoundaryCondition *self, BCContext *ctx)
1994{
1995 PetscErrorCode ierr;
1996 DrivenConstantData *data = (DrivenConstantData*)self->data;
1997 UserCtx* user = ctx->user;
1998 BCFace face_id = ctx->face_id;
1999 PetscBool can_service;
2000
2001 PetscFunctionBeginUser;
2002
2003 // Check if this rank owns part of this boundary face
2004 ierr = CanRankServiceFace(&user->info, user->IM, user->JM, user->KM, face_id, &can_service); CHKERRQ(ierr);
2005 if (!can_service) {
2006 PetscFunctionReturn(0);
2007 }
2008
2009 // If the correction is negligible, no work is needed.
2010 if (PetscAbsReal(data->boundaryVelocityCorrection) < 1e-12) {
2011 PetscFunctionReturn(0);
2012 }
2013
2014 LOG_ALLOW(LOCAL, LOG_TRACE, "Apply_PeriodicDrivenConstant: Applying boundary trim on Face %s...\n", BCFaceToString(face_id));
2015
2016 // --- Get read/write access to necessary arrays ---
2017 DMDALocalInfo info = user->info;
2018 Cmpnts ***ucont, ***uch, ***csi, ***eta, ***zet;
2019 PetscReal ***nvert;
2020
2021 ierr = DMDAVecGetArray(user->fda, user->lUcont, &ucont); CHKERRQ(ierr);
2022 ierr = DMDAVecGetArray(user->fda, user->Bcs.Uch, &uch); CHKERRQ(ierr);
2023 ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
2024 ierr = DMDAVecGetArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
2025 ierr = DMDAVecGetArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
2026 ierr = DMDAVecGetArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
2027
2028 PetscInt lxs = (info.xs == 0) ? 1 : info.xs;
2029 PetscInt lys = (info.ys == 0) ? 1 : info.ys;
2030 PetscInt lzs = (info.zs == 0) ? 1 : info.zs;
2031 PetscInt lxe = (info.xs + info.xm == info.mx) ? info.mx - 1 : info.xs + info.xm;
2032 PetscInt lye = (info.ys + info.ym == info.my) ? info.my - 1 : info.ys + info.ym;
2033 PetscInt lze = (info.zs + info.zm == info.mz) ? info.mz - 1 : info.zs + info.zm;
2034
2035 // --- Apply correction to the appropriate face and velocity component ---
2036 switch (face_id) {
2037 case BC_FACE_NEG_X: case BC_FACE_POS_X: {
2038 PetscInt i_face = (face_id == BC_FACE_NEG_X) ? info.xs : info.mx - 2;
2039 PetscInt i_nvert = (face_id == BC_FACE_NEG_X) ? info.xs + 1 : info.mx - 2;
2040
2041 for (PetscInt k = lzs; k < lze; k++) for (PetscInt j = lys; j < lye; j++) {
2042 if (nvert[k][j][i_nvert] < 0.1) {
2043 PetscReal faceArea = sqrt(csi[k][j][i_face].x*csi[k][j][i_nvert].x + csi[k][j][i_nvert].y*csi[k][j][i_nvert].y + csi[k][j][i_face].z*csi[k][j][i_face].z);
2044 PetscReal fluxTrim = data->boundaryVelocityCorrection * faceArea;
2045 if(data->applyBoundaryTrim) ucont[k][j][i_face].x += fluxTrim;
2046 uch[k][j][i_face].x = fluxTrim; // Store correction for diagnostics
2047 }
2048 }
2049 } break;
2050
2051 case BC_FACE_NEG_Y: case BC_FACE_POS_Y: {
2052 PetscInt j_face = (face_id == BC_FACE_NEG_Y) ? info.ys : info.my - 2;
2053 PetscInt j_nvert = (face_id == BC_FACE_NEG_Y) ? info.ys + 1 : info.my - 2;
2054
2055 for (PetscInt k = lzs; k < lze; k++) for (PetscInt i = lxs; i < lxe; i++) {
2056 if (nvert[k][j_nvert][i] < 0.1) {
2057 PetscReal faceArea = sqrt(eta[k][j_face][i].x*eta[k][j_face][i].x + eta[k][j_face][i].y*eta[k][j_face][i].y + eta[k][j_face][i].z*eta[k][j_face][i].z);
2058 PetscReal fluxTrim = data->boundaryVelocityCorrection * faceArea;
2059 if(data->applyBoundaryTrim) ucont[k][j_face][i].y += fluxTrim;
2060 uch[k][j_face][i].y = fluxTrim;
2061 }
2062 }
2063 } break;
2064
2065 case BC_FACE_NEG_Z: case BC_FACE_POS_Z: {
2066 PetscInt k_face = (face_id == BC_FACE_NEG_Z) ? info.zs : info.mz - 2;
2067 PetscInt k_nvert = (face_id == BC_FACE_NEG_Z) ? info.zs + 1 : info.mz - 2;
2068
2069 for (PetscInt j = lys; j < lye; j++) for (PetscInt i = lxs; i < lxe; i++) {
2070 if (nvert[k_nvert][j][i] < 0.1) {
2071 PetscReal faceArea = sqrt(zet[k_nvert][j][i].x*zet[k_nvert][j][i].x + zet[k_nvert][j][i].y*zet[k_nvert][j][i].y + zet[k_nvert][j][i].z*zet[k_nvert][j][i].z);
2072 PetscReal fluxTrim = data->boundaryVelocityCorrection * faceArea;
2073 if(data->applyBoundaryTrim) ucont[k_face][j][i].z += fluxTrim;
2074 uch[k_face][j][i].z = fluxTrim;
2075 }
2076 }
2077 } break;
2078 }
2079
2080 // --- Restore arrays ---
2081 ierr = DMDAVecRestoreArray(user->fda, user->lUcont, &ucont); CHKERRQ(ierr);
2082 ierr = DMDAVecRestoreArray(user->fda, user->Bcs.Uch, &uch); CHKERRQ(ierr);
2083 ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
2084 ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
2085 ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
2086 ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
2087
2088 PetscFunctionReturn(0);
2089}
2090
2091#undef __FUNCT__
2092#define __FUNCT__ "Destroy_PeriodicDrivenConstant"
2093/**
2094 * @brief Internal helper implementation: `Destroy_PeriodicDrivenConstant()`.
2095 * @details Local to this translation unit.
2096 */
2097static PetscErrorCode Destroy_PeriodicDrivenConstant(BoundaryCondition *self)
2098{
2099 PetscFunctionBeginUser;
2100
2101 // Check that the handler object and its private data pointer are valid before trying to free.
2102 if (self && self->data) {
2103 // The private data was allocated with PetscNew(), so it must be freed with PetscFree().
2104 PetscFree(self->data);
2105
2106 // It is good practice to nullify the pointer after freeing to prevent
2107 // any accidental use of the dangling pointer (use-after-free).
2108 self->data = NULL;
2109
2110 LOG_ALLOW(LOCAL, LOG_TRACE, "Destroy_PeriodicDrivenConstant: Private data freed successfully.\n");
2111 }
2112
2113 PetscFunctionReturn(0);
2114}
InletParabolicData::cs2_half
PetscReal cs2_half
Half-width (in index space) in cross-stream direction 2.
Definition BC_Handlers.c:710
Create_InletConstantVelocity
PetscErrorCode Create_InletConstantVelocity(BoundaryCondition *bc)
Implementation of Create_InletConstantVelocity().
Definition BC_Handlers.c:323
Apply_WallNoSlip
static PetscErrorCode Apply_WallNoSlip(BoundaryCondition *self, BCContext *ctx)
Internal helper implementation: Apply_WallNoSlip().
Definition BC_Handlers.c:144
DrivenConstantData::applyBoundaryTrim
PetscBool applyBoundaryTrim
Definition BC_Handlers.c:1687
Initialize_InletParabolicProfile
static PetscErrorCode Initialize_InletParabolicProfile(BoundaryCondition *self, BCContext *ctx)
Internal helper implementation: Initialize_InletParabolicProfile().
Definition BC_Handlers.c:750
PreStep_InletConstantVelocity
static PetscErrorCode PreStep_InletConstantVelocity(BoundaryCondition *self, BCContext *ctx, PetscReal *in, PetscReal *out)
Internal helper implementation: PreStep_InletConstantVelocity().
Definition BC_Handlers.c:402
Apply_InletConstantVelocity
static PetscErrorCode Apply_InletConstantVelocity(BoundaryCondition *self, BCContext *ctx)
Internal helper implementation: Apply_InletConstantVelocity().
Definition BC_Handlers.c:425
InletParabolicData::cs2_center
PetscReal cs2_center
Center index in cross-stream direction 2.
Definition BC_Handlers.c:708
Apply_OutletConservation
static PetscErrorCode Apply_OutletConservation(BoundaryCondition *self, BCContext *ctx)
(Handler Action) Applies mass conservation correction to the outlet face.
Definition BC_Handlers.c:1294
Apply_InletParabolicProfile
static PetscErrorCode Apply_InletParabolicProfile(BoundaryCondition *self, BCContext *ctx)
Internal helper implementation: Apply_InletParabolicProfile().
Definition BC_Handlers.c:836
PostStep_InletParabolicProfile
static PetscErrorCode PostStep_InletParabolicProfile(BoundaryCondition *self, BCContext *ctx, PetscReal *in, PetscReal *out)
Internal helper implementation: PostStep_InletParabolicProfile().
Definition BC_Handlers.c:994
DrivenConstantData::targetVolumetricFlux
PetscReal targetVolumetricFlux
Definition BC_Handlers.c:1684
DrivenConstantData::direction
char direction
Definition BC_Handlers.c:1683
DrivenConstantData::isMasterController
PetscBool isMasterController
Definition BC_Handlers.c:1686
Apply_PeriodicDrivenConstant
static PetscErrorCode Apply_PeriodicDrivenConstant(BoundaryCondition *self, BCContext *ctx)
Internal helper implementation: Apply_PeriodicDrivenConstant().
Definition BC_Handlers.c:1993
Initialize_PeriodicDrivenConstant
static PetscErrorCode Initialize_PeriodicDrivenConstant(BoundaryCondition *self, BCContext *ctx)
Internal helper implementation: Initialize_PeriodicDrivenConstant().
Definition BC_Handlers.c:1743
Create_PeriodicGeometric
PetscErrorCode Create_PeriodicGeometric(BoundaryCondition *bc)
Implementation of Create_PeriodicGeometric().
Definition BC_Handlers.c:1646
InletParabolicData::v_max
PetscReal v_max
Peak centerline velocity (from user params).
Definition BC_Handlers.c:706
Validate_DrivenFlowConfiguration
PetscErrorCode Validate_DrivenFlowConfiguration(UserCtx *user)
Internal helper implementation: Validate_DrivenFlowConfiguration().
Definition BC_Handlers.c:15
Create_InletParabolicProfile
PetscErrorCode Create_InletParabolicProfile(BoundaryCondition *bc)
Implementation of Create_InletParabolicProfile().
Definition BC_Handlers.c:721
Destroy_InletConstantVelocity
static PetscErrorCode Destroy_InletConstantVelocity(BoundaryCondition *self)
Internal helper implementation: Destroy_InletConstantVelocity().
Definition BC_Handlers.c:656
Destroy_PeriodicDrivenConstant
static PetscErrorCode Destroy_PeriodicDrivenConstant(BoundaryCondition *self)
Internal helper implementation: Destroy_PeriodicDrivenConstant().
Definition BC_Handlers.c:2097
Create_PeriodicDrivenConstant
PetscErrorCode Create_PeriodicDrivenConstant(BoundaryCondition *bc)
Internal helper implementation: Create_PeriodicDrivenConstant().
Definition BC_Handlers.c:1699
PreStep_InletParabolicProfile
static PetscErrorCode PreStep_InletParabolicProfile(BoundaryCondition *self, BCContext *ctx, PetscReal *in, PetscReal *out)
Internal helper implementation: PreStep_InletParabolicProfile().
Definition BC_Handlers.c:816
PostStep_OutletConservation
static PetscErrorCode PostStep_OutletConservation(BoundaryCondition *self, BCContext *ctx, PetscReal *in, PetscReal *out)
Internal helper implementation: PostStep_OutletConservation().
Definition BC_Handlers.c:1507
InletConstantData::normal_velocity
PetscReal normal_velocity
Definition BC_Handlers.c:312
Destroy_InletParabolicProfile
static PetscErrorCode Destroy_InletParabolicProfile(BoundaryCondition *self)
Internal helper implementation: Destroy_InletParabolicProfile().
Definition BC_Handlers.c:1088
Initialize_InletConstantVelocity
static PetscErrorCode Initialize_InletConstantVelocity(BoundaryCondition *self, BCContext *ctx)
Internal helper implementation: Initialize_InletConstantVelocity().
Definition BC_Handlers.c:352
Create_WallNoSlip
PetscErrorCode Create_WallNoSlip(BoundaryCondition *bc)
Implementation of Create_WallNoSlip().
Definition BC_Handlers.c:114
PreStep_PeriodicDrivenConstant
static PetscErrorCode PreStep_PeriodicDrivenConstant(BoundaryCondition *self, BCContext *ctx, PetscReal *in, PetscReal *out)
Internal helper implementation: PreStep_PeriodicDrivenConstant().
Definition BC_Handlers.c:1811
PostStep_InletConstantVelocity
static PetscErrorCode PostStep_InletConstantVelocity(BoundaryCondition *self, BCContext *ctx, PetscReal *in, PetscReal *out)
Internal helper implementation: PostStep_InletConstantVelocity().
Definition BC_Handlers.c:561
PreStep_OutletConservation
static PetscErrorCode PreStep_OutletConservation(BoundaryCondition *self, BCContext *ctx, PetscReal *local_inflow_contribution, PetscReal *local_outflow_contribution)
Internal helper implementation: PreStep_OutletConservation().
Definition BC_Handlers.c:1158
DrivenConstantData::boundaryVelocityCorrection
PetscReal boundaryVelocityCorrection
Definition BC_Handlers.c:1685
Create_OutletConservation
PetscErrorCode Create_OutletConservation(BoundaryCondition *bc)
Implementation of Create_OutletConservation().
Definition BC_Handlers.c:1129
InletParabolicData::cs1_half
PetscReal cs1_half
Half-width (in index space) in cross-stream direction 1.
Definition BC_Handlers.c:709
InletParabolicData::cs1_center
PetscReal cs1_center
Center index in cross-stream direction 1.
Definition BC_Handlers.c:707
DrivenConstantData
Private data structure for the handler.
Definition BC_Handlers.c:1682
InletConstantData
Private data structure for the Constant Velocity Inlet handler.
Definition BC_Handlers.c:311
InletParabolicData
Private data structure for the Parabolic Velocity Inlet handler.
Definition BC_Handlers.c:705
BC_Handlers.h
CanRankServiceFace
PetscErrorCode CanRankServiceFace(const DMDALocalInfo *info, PetscInt IM_nodes_global, PetscInt JM_nodes_global, PetscInt KM_nodes_global, BCFace face_id, PetscBool *can_service_out)
Determines if the current MPI rank owns any part of a specified global face.
Definition Boundaries.c:126
CalculateFaceCenterAndArea
PetscErrorCode CalculateFaceCenterAndArea(UserCtx *user, BCFace face_id, Cmpnts *face_center, PetscReal *face_area)
Calculates the geometric center and total area of a specified boundary face.
Definition grid.c:1029
GetBCParamReal
PetscErrorCode GetBCParamReal(BC_Param *params, const char *key, PetscReal *value_out, PetscBool *found)
Searches a BC_Param linked list for a key and returns its value as a double.
Definition io.c:387
GetBCParamBool
PetscErrorCode GetBCParamBool(BC_Param *params, const char *key, PetscBool *value_out, PetscBool *found)
Searches a BC_Param linked list for a key and returns its value as a bool.
Definition io.c:408
LOCAL
#define LOCAL
Logging scope definitions for controlling message output.
Definition logging.h:44
GLOBAL
#define GLOBAL
Scope for global logging across all processes.
Definition logging.h:45
BCFaceToString
const char * BCFaceToString(BCFace face)
Helper function to convert BCFace enum to a string representation.
Definition logging.c:669
LOG_ALLOW
#define LOG_ALLOW(scope, level, fmt,...)
Logging macro that checks both the log level and whether the calling function is in the allowed-funct...
Definition logging.h:199
PROFILE_FUNCTION_END
#define PROFILE_FUNCTION_END
Marks the end of a profiled code block.
Definition logging.h:779
BCTypeToString
const char * BCTypeToString(BCType type)
Helper function to convert BCType enum to a string representation.
Definition logging.c:751
LOG_TRACE
@ LOG_TRACE
Very fine-grained tracing information for in-depth debugging.
Definition logging.h:32
LOG_INFO
@ LOG_INFO
Informational messages about program execution.
Definition logging.h:30
LOG_WARNING
@ LOG_WARNING
Non-critical issues that warrant attention.
Definition logging.h:29
LOG_DEBUG
@ LOG_DEBUG
Detailed debugging information.
Definition logging.h:31
PROFILE_FUNCTION_BEGIN
#define PROFILE_FUNCTION_BEGIN
Marks the beginning of a profiled code block (typically a function).
Definition logging.h:770
BoundaryCondition
The "virtual table" struct for a boundary condition handler object.
Definition variables.h:321
BoundaryCondition::PostStep
PetscErrorCode(* PostStep)(BoundaryCondition *self, BCContext *ctx, PetscReal *local_inflow, PetscReal *local_outflow)
Definition variables.h:328
BoundaryCondition::PreStep
PetscErrorCode(* PreStep)(BoundaryCondition *self, BCContext *ctx, PetscReal *local_inflow, PetscReal *local_outflow)
Definition variables.h:326
BoundaryCondition::Destroy
PetscErrorCode(* Destroy)(BoundaryCondition *self)
Definition variables.h:330
BoundaryCondition::Initialize
PetscErrorCode(* Initialize)(BoundaryCondition *self, BCContext *ctx)
Definition variables.h:325
BoundaryCondition::UpdateUbcs
PetscErrorCode(* UpdateUbcs)(BoundaryCondition *self, BCContext *ctx)
Definition variables.h:329
BoundaryCondition::Apply
PetscErrorCode(* Apply)(BoundaryCondition *self, BCContext *ctx)
Definition variables.h:327
BoundaryCondition::priority
BCPriorityType priority
Definition variables.h:323
BCContext::global_outflow_sum
const PetscReal * global_outflow_sum
Definition variables.h:317
BCType
BCType
Defines the general mathematical/physical Category of a boundary.
Definition variables.h:251
INLET
@ INLET
Definition variables.h:258
FARFIELD
@ FARFIELD
Definition variables.h:259
OUTLET
@ OUTLET
Definition variables.h:257
PERIODIC
@ PERIODIC
Definition variables.h:260
UserCtx::boundary_faces
BoundaryFaceConfig boundary_faces[6]
Definition variables.h:829
SimCtx::targetVolumetricFlux
PetscReal targetVolumetricFlux
Definition variables.h:724
UserCtx::lNvert
Vec lNvert
Definition variables.h:837
UserCtx::simCtx
SimCtx * simCtx
Back-pointer to the master simulation context.
Definition variables.h:814
UserCtx::KM
PetscInt KM
Definition variables.h:820
UserCtx::lZet
Vec lZet
Definition variables.h:858
BCHandlerType
BCHandlerType
Defines the specific computational "strategy" for a boundary handler.
Definition variables.h:271
BC_HANDLER_PERIODIC_DRIVEN_INITIAL_FLUX
@ BC_HANDLER_PERIODIC_DRIVEN_INITIAL_FLUX
Definition variables.h:287
BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX
@ BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX
Definition variables.h:286
BoundaryFaceConfig::handler_type
BCHandlerType handler_type
Definition variables.h:337
SimCtx::bulkVelocityCorrection
PetscReal bulkVelocityCorrection
Definition variables.h:725
BCContext::face_id
BCFace face_id
Definition variables.h:313
UserCtx::Ucont
Vec Ucont
Definition variables.h:837
BCS::Ubcs
Vec Ubcs
Physical Cartesian velocity at boundary faces. Full 3D array but only boundary-face entries are meani...
Definition variables.h:121
Cmpnts::x
PetscScalar x
Definition variables.h:101
UserCtx::Bcs
BCS Bcs
Definition variables.h:832
BCContext::user
UserCtx * user
Definition variables.h:312
BCContext::global_inflow_sum
const PetscReal * global_inflow_sum
Definition variables.h:314
UserCtx::lCsi
Vec lCsi
Definition variables.h:858
BoundaryFaceConfig::params
BC_Param * params
Definition variables.h:338
Cmpnts::z
PetscScalar z
Definition variables.h:101
UserCtx::JM
PetscInt JM
Definition variables.h:820
BCContext::global_farfield_inflow_sum
const PetscReal * global_farfield_inflow_sum
Definition variables.h:315
UserCtx::lUcont
Vec lUcont
Definition variables.h:837
SimCtx::AreaOutSum
PetscReal AreaOutSum
Definition variables.h:726
UserCtx::info
DMDALocalInfo info
Definition variables.h:818
BC_PRIORITY_OUTLET
@ BC_PRIORITY_OUTLET
Definition variables.h:296
BC_PRIORITY_WALL
@ BC_PRIORITY_WALL
Definition variables.h:295
BC_PRIORITY_INLET
@ BC_PRIORITY_INLET
Definition variables.h:293
UserCtx::lUcat
Vec lUcat
Definition variables.h:837
Cmpnts::y
PetscScalar y
Definition variables.h:101
UserCtx::IM
PetscInt IM
Definition variables.h:820
UserCtx::lEta
Vec lEta
Definition variables.h:858
BoundaryFaceConfig::mathematical_type
BCType mathematical_type
Definition variables.h:336
BCS::Uch
Vec Uch
Characteristic velocity for boundary conditions.
Definition variables.h:122
BCFace
BCFace
Identifies the six logical faces of a structured computational block.
Definition variables.h:244
BC_FACE_NEG_X
@ BC_FACE_NEG_X
Definition variables.h:245
BC_FACE_POS_Z
@ BC_FACE_POS_Z
Definition variables.h:247
BC_FACE_POS_Y
@ BC_FACE_POS_Y
Definition variables.h:246
BC_FACE_NEG_Z
@ BC_FACE_NEG_Z
Definition variables.h:247
BC_FACE_POS_X
@ BC_FACE_POS_X
Definition variables.h:245
BC_FACE_NEG_Y
@ BC_FACE_NEG_Y
Definition variables.h:246
BCContext
Provides execution context for a boundary condition handler.
Definition variables.h:311
BoundaryFaceConfig
Holds the complete configuration for one of the six boundary faces.
Definition variables.h:334
Cmpnts
A 3D point or vector with PetscScalar components.
Definition variables.h:100
SimCtx
The master context for the entire simulation.
Definition variables.h:643
UserCtx
User-defined context containing data specific to a single computational grid level.
Definition variables.h:811
Generated by doxygen 1.9.8