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: PRESCRIBED INLET PROFILE FROM FILE
1101// (Corresponds to BC_HANDLER_INLET_PROFILE_FROM_FILE)
1102//
1103// This handler reads positive scalar normal speeds from a canonical PICSLICE file.
1104// It then applies those speeds through the same face sign and metric conversion
1105// used by the constant and parabolic inlet handlers.
1106//
1107//================================================================================
1108
1109static PetscErrorCode Initialize_InletProfileFromFile(BoundaryCondition *self, BCContext *ctx);
1110static PetscErrorCode PreStep_InletProfileFromFile(BoundaryCondition *self, BCContext *ctx,
1111 PetscReal *in, PetscReal *out);
1112static PetscErrorCode Apply_InletProfileFromFile(BoundaryCondition *self, BCContext *ctx);
1113static PetscErrorCode PostStep_InletProfileFromFile(BoundaryCondition *self, BCContext *ctx,
1114 PetscReal *in, PetscReal *out);
1115static PetscErrorCode Destroy_InletProfileFromFile(BoundaryCondition *self);
1116
1117typedef struct {
1118 PetscInt n1;
1119 PetscInt n2;
1120 PetscReal *profile;
1121 PetscReal min_speed;
1122 PetscReal max_speed;
1123 char *source_file;
1124} InletProfileFileData;
1125
1126/**
1127 * @brief Looks up a string-valued boundary-condition parameter in a BC_Param list.
1128 *
1129 * @param params Head of the boundary-condition parameter linked list.
1130 * @param key Case-insensitive key to search for.
1131 * @param[out] value_out Borrowed pointer to the matching value string, or NULL if absent.
1132 * @param[out] found PETSC_TRUE when the key is present, PETSC_FALSE otherwise.
1133 * @return PetscErrorCode 0 on success.
1134 */
1135static PetscErrorCode GetBCParamStringLocal(BC_Param *params, const char *key,
1136 const char **value_out, PetscBool *found)
1137{
1138 PetscFunctionBeginUser;
1139 *found = PETSC_FALSE;
1140 *value_out = NULL;
1141 for (BC_Param *param = params; param; param = param->next) {
1142 if (strcasecmp(param->key, key) == 0) {
1143 *value_out = param->value;
1144 *found = PETSC_TRUE;
1145 PetscFunctionReturn(0);
1146 }
1147 }
1148 PetscFunctionReturn(0);
1149}
1150
1151/**
1152 * @brief Computes the expected PICSLICE dimensions for an inlet face.
1153 *
1154 * @param user User context containing global grid node counts.
1155 * @param face_id Boundary face whose tangential profile dimensions are requested.
1156 * @param[out] n1 First PICSLICE dimension in handler storage order.
1157 * @param[out] n2 Second PICSLICE dimension in handler storage order.
1158 * @return PetscErrorCode 0 on success, or a PETSc error for unsupported faces or invalid dimensions.
1159 */
1160static PetscErrorCode GetProfileFileExpectedDims(UserCtx *user, BCFace face_id,
1161 PetscInt *n1, PetscInt *n2)
1162{
1163 PetscFunctionBeginUser;
1164 switch (face_id) {
1165 case BC_FACE_NEG_X:
1166 case BC_FACE_POS_X:
1167 *n1 = user->KM - 1;
1168 *n2 = user->JM - 1;
1169 break;
1170 case BC_FACE_NEG_Y:
1171 case BC_FACE_POS_Y:
1172 *n1 = user->KM - 1;
1173 *n2 = user->IM - 1;
1174 break;
1175 case BC_FACE_NEG_Z:
1176 case BC_FACE_POS_Z:
1177 *n1 = user->JM - 1;
1178 *n2 = user->IM - 1;
1179 break;
1180 default:
1181 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE,
1182 "Unsupported face id %d for inlet profile dimensions.", face_id);
1183 }
1184 if (*n1 <= 0 || *n2 <= 0) {
1185 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
1186 "Invalid inlet profile dimensions (%d, %d) for grid (%d, %d, %d).",
1187 *n1, *n2, user->IM, user->JM, user->KM);
1188 }
1189 PetscFunctionReturn(0);
1190}
1191
1192/**
1193 * @brief Reads and validates a static scalar inlet profile from a canonical PICSLICE file.
1194 *
1195 * @details The file must contain magic token `PICSLICE`, frame count 1, the expected
1196 * two-dimensional face shape, and exactly one finite nonnegative scalar speed
1197 * value per interior face slot. Values are stored row-major in `data->profile`.
1198 *
1199 * @param source_file Path to the PICSLICE profile file.
1200 * @param expected_n1 Expected first slice dimension.
1201 * @param expected_n2 Expected second slice dimension.
1202 * @param[in,out] data Handler-private storage that receives dimensions, profile values, and min/max speeds.
1203 * @return PetscErrorCode 0 on success, or a PETSc file/validation error on malformed input.
1204 */
1205static PetscErrorCode ReadPicSliceProfile(const char *source_file, PetscInt expected_n1,
1206 PetscInt expected_n2, InletProfileFileData *data)
1207{
1208 PetscErrorCode ierr;
1209 FILE *fd = NULL;
1210 char magic[32] = {0};
1211 PetscInt frame_count = 0, n1 = 0, n2 = 0;
1212
1213 PetscFunctionBeginUser;
1214 fd = fopen(source_file, "r");
1215 if (!fd) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN,
1216 "Cannot open PICSLICE inlet profile file: %s", source_file);
1217
1218 if (fscanf(fd, "%31s", magic) != 1 || strcmp(magic, "PICSLICE") != 0) {
1219 fclose(fd);
1220 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_READ,
1221 "PICSLICE inlet profile file %s must begin with PICSLICE header.", source_file);
1222 }
1223 if (fscanf(fd, "%d", &frame_count) != 1) {
1224 fclose(fd);
1225 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_READ,
1226 "PICSLICE inlet profile file %s missing frame count.", source_file);
1227 }
1228 if (frame_count != 1) {
1229 fclose(fd);
1230 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED,
1231 "PICSLICE inlet profile file %s has %d frames; static handler requires 1.",
1232 source_file, frame_count);
1233 }
1234 if (fscanf(fd, "%d %d", &n1, &n2) != 2) {
1235 fclose(fd);
1236 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_READ,
1237 "PICSLICE inlet profile file %s missing slice dimensions.", source_file);
1238 }
1239 if (n1 != expected_n1 || n2 != expected_n2) {
1240 fclose(fd);
1241 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED,
1242 "PICSLICE inlet profile dimensions mismatch for %s: expected (%d, %d), found (%d, %d).",
1243 source_file, expected_n1, expected_n2, n1, n2);
1244 }
1245
1246 data->n1 = n1;
1247 data->n2 = n2;
1248 ierr = PetscMalloc1(n1 * n2, &data->profile); CHKERRQ(ierr);
1249 data->min_speed = PETSC_MAX_REAL;
1250 data->max_speed = -PETSC_MAX_REAL;
1251
1252 for (PetscInt idx = 0; idx < n1 * n2; idx++) {
1253 PetscReal value = 0.0;
1254 if (fscanf(fd, "%le", &value) != 1) {
1255 fclose(fd);
1256 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_READ,
1257 "PICSLICE inlet profile file %s ended early: expected %d values.",
1258 source_file, n1 * n2);
1259 }
1260 if (PetscIsInfOrNanReal(value) || value < 0.0) {
1261 fclose(fd);
1262 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED,
1263 "PICSLICE inlet profile file %s contains invalid speed %.6e at flat index %d.",
1264 source_file, (double)value, idx);
1265 }
1266 data->profile[idx] = value;
1267 data->min_speed = PetscMin(data->min_speed, value);
1268 data->max_speed = PetscMax(data->max_speed, value);
1269 }
1270
1271 char extra[64];
1272 if (fscanf(fd, "%63s", extra) == 1) {
1273 fclose(fd);
1274 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED,
1275 "PICSLICE inlet profile file %s has extra token after %d values: %s",
1276 source_file, n1 * n2, extra);
1277 }
1278 fclose(fd);
1279 PetscFunctionReturn(0);
1280}
1281
1282/**
1283 * @brief Returns one scalar speed from the flattened PICSLICE profile.
1284 *
1285 * @param data Handler-private profile storage.
1286 * @param a First profile index in face-specific storage order.
1287 * @param b Second profile index in face-specific storage order.
1288 * @return Scalar normal speed at `(a, b)`.
1289 */
1290static inline PetscReal ProfileSpeedAt(const InletProfileFileData *data, PetscInt a, PetscInt b)
1291{
1292 return data->profile[a * data->n2 + b];
1293}
1294
1295#undef __FUNCT__
1296#define __FUNCT__ "Create_InletProfileFromFile"
1297/**
1298 * @brief Implementation of \ref Create_InletProfileFromFile().
1299 * @details Full API contract (arguments, ownership, side effects) is documented with
1300 * the header declaration in `include/BC_Handlers.h`.
1301 * @see Create_InletProfileFromFile()
1302 */
1303PetscErrorCode Create_InletProfileFromFile(BoundaryCondition *bc)
1304{
1305 PetscErrorCode ierr;
1306 PetscFunctionBeginUser;
1307
1308 if (!bc) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "BoundaryCondition is NULL");
1309
1310 InletProfileFileData *data = NULL;
1311 ierr = PetscMalloc1(1, &data); CHKERRQ(ierr);
1312 data->n1 = 0;
1313 data->n2 = 0;
1314 data->profile = NULL;
1315 data->min_speed = 0.0;
1316 data->max_speed = 0.0;
1317 data->source_file = NULL;
1318 bc->data = (void*)data;
1319
1320 bc->priority = BC_PRIORITY_INLET;
1321 bc->Initialize = Initialize_InletProfileFromFile;
1322 bc->PreStep = PreStep_InletProfileFromFile;
1323 bc->Apply = Apply_InletProfileFromFile;
1324 bc->PostStep = PostStep_InletProfileFromFile;
1325 bc->UpdateUbcs = NULL;
1326 bc->Destroy = Destroy_InletProfileFromFile;
1327
1328 PetscFunctionReturn(0);
1329}
1330
1331#undef __FUNCT__
1332#define __FUNCT__ "Initialize_InletProfileFromFile"
1333/**
1334 * @brief Initializes a file-prescribed inlet profile handler for one boundary face.
1335 *
1336 * @details Reads the `source_file` BC parameter, validates the target face dimensions,
1337 * loads the PICSLICE scalar speed profile, records summary statistics for logging,
1338 * and applies the profile once to initialize boundary state.
1339 *
1340 * @param self BoundaryCondition object configured by Create_InletProfileFromFile().
1341 * @param ctx Runtime boundary context containing the UserCtx and face id.
1342 * @return PetscErrorCode 0 on success, or a PETSc error for missing parameters or malformed files.
1343 */
1344static PetscErrorCode Initialize_InletProfileFromFile(BoundaryCondition *self, BCContext *ctx)
1345{
1346 PetscErrorCode ierr;
1347 UserCtx *user = ctx->user;
1348 BCFace face_id = ctx->face_id;
1349 InletProfileFileData *data = (InletProfileFileData*)self->data;
1350 PetscBool found = PETSC_FALSE;
1351 const char *source_file = NULL;
1352 PetscInt expected_n1 = 0, expected_n2 = 0;
1353
1354 PetscFunctionBeginUser;
1355 ierr = GetBCParamStringLocal(user->boundary_faces[face_id].params, "source_file",
1356 &source_file, &found); CHKERRQ(ierr);
1357 if (!found || !source_file || source_file[0] == '\0') {
1358 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
1359 "InletProfileFromFile requires source_file parameter for face %d.", face_id);
1360 }
1361
1362 ierr = PetscStrallocpy(source_file, &data->source_file); CHKERRQ(ierr);
1363 ierr = GetProfileFileExpectedDims(user, face_id, &expected_n1, &expected_n2); CHKERRQ(ierr);
1364 ierr = ReadPicSliceProfile(source_file, expected_n1, expected_n2, data); CHKERRQ(ierr);
1365
1366 LOG_ALLOW(LOCAL, LOG_INFO,
1367 " Inlet Face %d (Prescribed Flow): source=%s dims=(%d,%d) speed[min,max]=[%.6e, %.6e]\n",
1368 face_id, data->source_file, data->n1, data->n2,
1369 (double)data->min_speed, (double)data->max_speed);
1370
1371 ierr = Apply_InletProfileFromFile(self, ctx); CHKERRQ(ierr);
1372 PetscFunctionReturn(0);
1373}
1374
1375#undef __FUNCT__
1376#define __FUNCT__ "PreStep_InletProfileFromFile"
1377/**
1378 * @brief Pre-step hook for the static file-prescribed inlet profile handler.
1379 *
1380 * @details Static profiles require no per-step preparation. The hook is implemented
1381 * so future time-varying profile support can reuse the same handler lifecycle.
1382 *
1383 * @param self BoundaryCondition object for this inlet handler.
1384 * @param ctx Runtime boundary context.
1385 * @param local_inflow_contribution Inflow accumulator, intentionally unchanged.
1386 * @param local_outflow_contribution Outflow accumulator, intentionally unchanged.
1387 * @return PetscErrorCode 0 on success.
1388 */
1389static PetscErrorCode PreStep_InletProfileFromFile(BoundaryCondition *self, BCContext *ctx,
1390 PetscReal *local_inflow_contribution,
1391 PetscReal *local_outflow_contribution)
1392{
1393 (void)self;
1394 (void)ctx;
1395 (void)local_inflow_contribution;
1396 (void)local_outflow_contribution;
1397 PetscFunctionBeginUser;
1398 PetscFunctionReturn(0);
1399}
1400
1401#undef __FUNCT__
1402#define __FUNCT__ "Apply_InletProfileFromFile"
1403/**
1404 * @brief Applies the loaded PICSLICE scalar profile to Ucont and Ubcs on an inlet face.
1405 *
1406 * @details Each stored scalar is treated as a positive normal speed magnitude. The routine
1407 * uses the same negative/positive face sign convention, immersed-boundary skip
1408 * checks, metric vectors, and CellArea conversion used by the constant and
1409 * parabolic inlet handlers.
1410 *
1411 * @param self BoundaryCondition object with InletProfileFileData storage.
1412 * @param ctx Runtime boundary context containing arrays and face id.
1413 * @return PetscErrorCode 0 on success, or a PETSc error from DMDA array access.
1414 */
1415static PetscErrorCode Apply_InletProfileFromFile(BoundaryCondition *self, BCContext *ctx)
1416{
1417 PetscErrorCode ierr;
1418 UserCtx *user = ctx->user;
1419 BCFace face_id = ctx->face_id;
1420 InletProfileFileData *data = (InletProfileFileData*)self->data;
1421 PetscBool can_service;
1422
1423 PetscFunctionBeginUser;
1424 DMDALocalInfo *info = &user->info;
1425 Cmpnts ***ubcs, ***ucont, ***csi, ***eta, ***zet;
1426 PetscReal ***nvert;
1427
1428 ierr = CanRankServiceFace(info, user->IM, user->JM, user->KM, face_id, &can_service); CHKERRQ(ierr);
1429 if (!can_service) PetscFunctionReturn(0);
1430
1431 ierr = DMDAVecGetArray(user->fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
1432 ierr = DMDAVecGetArray(user->fda, user->Ucont, &ucont); CHKERRQ(ierr);
1433 ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
1434 ierr = DMDAVecGetArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
1435 ierr = DMDAVecGetArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
1436 ierr = DMDAVecGetArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
1437
1438 PetscInt xs = info->xs, xe = info->xs + info->xm;
1439 PetscInt ys = info->ys, ye = info->ys + info->ym;
1440 PetscInt zs = info->zs, ze = info->zs + info->zm;
1441 PetscInt mx = info->mx, my = info->my, mz = info->mz;
1442
1443 PetscInt lxs = xs, lxe = xe, lys = ys, lye = ye, lzs = zs, lze = ze;
1444 if (xs == 0) lxs = xs + 1;
1445 if (xe == mx) lxe = xe - 1;
1446 if (ys == 0) lys = ys + 1;
1447 if (ye == my) lye = ye - 1;
1448 if (zs == 0) lzs = zs + 1;
1449 if (ze == mz) lze = ze - 1;
1450
1451 switch (face_id) {
1452 case BC_FACE_NEG_X:
1453 case BC_FACE_POS_X: {
1454 PetscReal sign = (face_id == BC_FACE_NEG_X) ? 1.0 : -1.0;
1455 PetscInt i = (face_id == BC_FACE_NEG_X) ? xs : mx - 2;
1456 for (PetscInt k = lzs; k < lze; k++) {
1457 for (PetscInt j = lys; j < lye; j++) {
1458 if ((sign > 0 && nvert[k][j][i+1] > 0.1) ||
1459 (sign < 0 && nvert[k][j][i] > 0.1)) continue;
1460 PetscReal uin_local = ProfileSpeedAt(data, k - 1, j - 1);
1461 PetscReal CellArea = sqrt(csi[k][j][i].x * csi[k][j][i].x +
1462 csi[k][j][i].y * csi[k][j][i].y +
1463 csi[k][j][i].z * csi[k][j][i].z);
1464 ucont[k][j][i].x = sign * uin_local * CellArea;
1465 ubcs[k][j][i + (sign < 0)].x = sign * uin_local * csi[k][j][i].x / CellArea;
1466 ubcs[k][j][i + (sign < 0)].y = sign * uin_local * csi[k][j][i].y / CellArea;
1467 ubcs[k][j][i + (sign < 0)].z = sign * uin_local * csi[k][j][i].z / CellArea;
1468 }
1469 }
1470 } break;
1471 case BC_FACE_NEG_Y:
1472 case BC_FACE_POS_Y: {
1473 PetscReal sign = (face_id == BC_FACE_NEG_Y) ? 1.0 : -1.0;
1474 PetscInt j = (face_id == BC_FACE_NEG_Y) ? ys : my - 2;
1475 for (PetscInt k = lzs; k < lze; k++) {
1476 for (PetscInt i = lxs; i < lxe; i++) {
1477 if ((sign > 0 && nvert[k][j+1][i] > 0.1) ||
1478 (sign < 0 && nvert[k][j][i] > 0.1)) continue;
1479 PetscReal uin_local = ProfileSpeedAt(data, k - 1, i - 1);
1480 PetscReal CellArea = sqrt(eta[k][j][i].x * eta[k][j][i].x +
1481 eta[k][j][i].y * eta[k][j][i].y +
1482 eta[k][j][i].z * eta[k][j][i].z);
1483 ucont[k][j][i].y = sign * uin_local * CellArea;
1484 ubcs[k][j + (sign < 0)][i].x = sign * uin_local * eta[k][j][i].x / CellArea;
1485 ubcs[k][j + (sign < 0)][i].y = sign * uin_local * eta[k][j][i].y / CellArea;
1486 ubcs[k][j + (sign < 0)][i].z = sign * uin_local * eta[k][j][i].z / CellArea;
1487 }
1488 }
1489 } break;
1490 case BC_FACE_NEG_Z:
1491 case BC_FACE_POS_Z: {
1492 PetscReal sign = (face_id == BC_FACE_NEG_Z) ? 1.0 : -1.0;
1493 PetscInt k = (face_id == BC_FACE_NEG_Z) ? zs : mz - 2;
1494 for (PetscInt j = lys; j < lye; j++) {
1495 for (PetscInt i = lxs; i < lxe; i++) {
1496 if ((sign > 0 && nvert[k+1][j][i] > 0.1) ||
1497 (sign < 0 && nvert[k][j][i] > 0.1)) continue;
1498 PetscReal uin_local = ProfileSpeedAt(data, j - 1, i - 1);
1499 PetscReal CellArea = sqrt(zet[k][j][i].x * zet[k][j][i].x +
1500 zet[k][j][i].y * zet[k][j][i].y +
1501 zet[k][j][i].z * zet[k][j][i].z);
1502 ucont[k][j][i].z = sign * uin_local * CellArea;
1503 ubcs[k + (sign < 0)][j][i].x = sign * uin_local * zet[k][j][i].x / CellArea;
1504 ubcs[k + (sign < 0)][j][i].y = sign * uin_local * zet[k][j][i].y / CellArea;
1505 ubcs[k + (sign < 0)][j][i].z = sign * uin_local * zet[k][j][i].z / CellArea;
1506 }
1507 }
1508 } break;
1509 }
1510
1511 ierr = DMDAVecRestoreArray(user->fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
1512 ierr = DMDAVecRestoreArray(user->fda, user->Ucont, &ucont); CHKERRQ(ierr);
1513 ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
1514 ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
1515 ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
1516 ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
1517
1518 PetscFunctionReturn(0);
1519}
1520
1521#undef __FUNCT__
1522#define __FUNCT__ "PostStep_InletProfileFromFile"
1523/**
1524 * @brief Accumulates the applied inlet flux for a file-prescribed profile.
1525 *
1526 * @details Sums the face-normal Ucont component over the same interior face slots
1527 * populated by Apply_InletProfileFromFile().
1528 *
1529 * @param self BoundaryCondition object for this inlet handler.
1530 * @param ctx Runtime boundary context containing the UserCtx and face id.
1531 * @param local_inflow_contribution Accumulator incremented by the measured inlet flux.
1532 * @param local_outflow_contribution Outflow accumulator, intentionally unchanged.
1533 * @return PetscErrorCode 0 on success, or a PETSc error from DMDA array access.
1534 */
1535static PetscErrorCode PostStep_InletProfileFromFile(BoundaryCondition *self, BCContext *ctx,
1536 PetscReal *local_inflow_contribution,
1537 PetscReal *local_outflow_contribution)
1538{
1539 PetscErrorCode ierr;
1540 UserCtx *user = ctx->user;
1541 BCFace face_id = ctx->face_id;
1542 PetscBool can_service;
1543
1544 (void)self;
1545 (void)local_outflow_contribution;
1546
1547 PetscFunctionBeginUser;
1548 DMDALocalInfo *info = &user->info;
1549 Cmpnts ***ucont;
1550
1551 ierr = CanRankServiceFace(info, user->IM, user->JM, user->KM, face_id, &can_service); CHKERRQ(ierr);
1552 if (!can_service) PetscFunctionReturn(0);
1553
1554 ierr = DMDAVecGetArrayRead(user->fda, user->Ucont, (const Cmpnts***)&ucont); CHKERRQ(ierr);
1555 PetscReal local_flux = 0.0;
1556
1557 PetscInt xs = info->xs, xe = info->xs + info->xm;
1558 PetscInt ys = info->ys, ye = info->ys + info->ym;
1559 PetscInt zs = info->zs, ze = info->zs + info->zm;
1560 PetscInt mx = info->mx, my = info->my, mz = info->mz;
1561
1562 PetscInt lxs = xs, lxe = xe, lys = ys, lye = ye, lzs = zs, lze = ze;
1563 if (xs == 0) lxs = xs + 1;
1564 if (xe == mx) lxe = xe - 1;
1565 if (ys == 0) lys = ys + 1;
1566 if (ye == my) lye = ye - 1;
1567 if (zs == 0) lzs = zs + 1;
1568 if (ze == mz) lze = ze - 1;
1569
1570 switch (face_id) {
1571 case BC_FACE_NEG_X:
1572 case BC_FACE_POS_X: {
1573 PetscInt i = (face_id == BC_FACE_NEG_X) ? xs : mx - 2;
1574 for (PetscInt k = lzs; k < lze; k++)
1575 for (PetscInt j = lys; j < lye; j++)
1576 local_flux += ucont[k][j][i].x;
1577 } break;
1578 case BC_FACE_NEG_Y:
1579 case BC_FACE_POS_Y: {
1580 PetscInt j = (face_id == BC_FACE_NEG_Y) ? ys : my - 2;
1581 for (PetscInt k = lzs; k < lze; k++)
1582 for (PetscInt i = lxs; i < lxe; i++)
1583 local_flux += ucont[k][j][i].y;
1584 } break;
1585 case BC_FACE_NEG_Z:
1586 case BC_FACE_POS_Z: {
1587 PetscInt k = (face_id == BC_FACE_NEG_Z) ? zs : mz - 2;
1588 for (PetscInt j = lys; j < lye; j++)
1589 for (PetscInt i = lxs; i < lxe; i++)
1590 local_flux += ucont[k][j][i].z;
1591 } break;
1592 }
1593
1594 ierr = DMDAVecRestoreArrayRead(user->fda, user->Ucont, (const Cmpnts***)&ucont); CHKERRQ(ierr);
1595 *local_inflow_contribution += local_flux;
1596
1597 LOG_ALLOW(LOCAL, LOG_DEBUG, "PostStep_InletProfileFromFile: Face %d, flux = %.6e\n",
1598 face_id, local_flux);
1599
1600 PetscFunctionReturn(0);
1601}
1602
1603#undef __FUNCT__
1604#define __FUNCT__ "Destroy_InletProfileFromFile"
1605/**
1606 * @brief Releases private storage owned by a file-prescribed inlet profile handler.
1607 *
1608 * @param self BoundaryCondition object whose `data` field stores InletProfileFileData.
1609 * @return PetscErrorCode 0 on success.
1610 */
1611static PetscErrorCode Destroy_InletProfileFromFile(BoundaryCondition *self)
1612{
1613 PetscFunctionBeginUser;
1614 if (self && self->data) {
1615 InletProfileFileData *data = (InletProfileFileData*)self->data;
1616 PetscFree(data->profile);
1617 PetscFree(data->source_file);
1618 PetscFree(self->data);
1619 self->data = NULL;
1620 }
1621 PetscFunctionReturn(0);
1622}
1623
1624//================================================================================
1625//
1626// HANDLER IMPLEMENTATION: OUTLET WITH MASS CONSERVATION
1627// (Corresponds to BC_HANDLER_OUTLET_CONSERVATION)
1628//
1629// This handler ensures that the total flux leaving through outlet boundaries
1630// balances the total flux entering through inlet and far-field boundaries.
1631//
1632// Workflow:
1633// 1. PreStep: Measures the *uncorrected* flux based on interior velocities.
1634// 2. Apply: Calculates a global correction factor based on the flux imbalance
1635// and applies it to the contravariant velocity (ucont) on the outlet face.
1636// 3. PostStep: Measures the *corrected* flux for verification and logging.
1637//
1638//================================================================================
1639
1640// --- 1. FORWARD DECLARATIONS ---
1641static PetscErrorCode PreStep_OutletConservation(BoundaryCondition *self, BCContext *ctx,
1642 PetscReal *local_inflow_contribution, PetscReal *local_outflow_contribution);
1643static PetscErrorCode Apply_OutletConservation(BoundaryCondition *self, BCContext *ctx);
1644static PetscErrorCode PostStep_OutletConservation(BoundaryCondition *self, BCContext *ctx,
1645 PetscReal *in, PetscReal *out);
1646
1647#undef __FUNCT__
1648#define __FUNCT__ "Create_OutletConservation"
1649/**
1650 * @brief Implementation of \ref Create_OutletConservation().
1651 * @details Full API contract (arguments, ownership, side effects) is documented with
1652 * the header declaration in `include/BC_Handlers.h`.
1653 * @see Create_OutletConservation()
1654 */
1655PetscErrorCode Create_OutletConservation(BoundaryCondition *bc)
1656{
1657 PetscFunctionBeginUser;
1658
1659 if (!bc) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Input BoundaryCondition is NULL");
1660
1661 // This handler has the highest priority to ensure it runs after
1662 // all inflow fluxes have been calculated.
1663 bc->priority = BC_PRIORITY_OUTLET;
1664
1665 // Assign function pointers
1666 bc->Initialize = NULL; // No initialization needed
1667 bc->PreStep = PreStep_OutletConservation;
1668 bc->Apply = Apply_OutletConservation;
1669 bc->PostStep = PostStep_OutletConservation;
1670 bc->UpdateUbcs = NULL;
1671 bc->Destroy = NULL; // No private data to destroy
1672
1673 bc->data = NULL;
1674
1675 PetscFunctionReturn(0);
1676}
1677
1678#undef __FUNCT__
1679#define __FUNCT__ "PreStep_OutletConservation"
1680/**
1681 * @brief Internal helper implementation: `PreStep_OutletConservation()`.
1682 * @details Local to this translation unit.
1683 */
1684static PetscErrorCode PreStep_OutletConservation(BoundaryCondition *self, BCContext *ctx,
1685 PetscReal *local_inflow_contribution, PetscReal *local_outflow_contribution)
1686{
1687 PetscErrorCode ierr;
1688 UserCtx* user = ctx->user;
1689 BCFace face_id = ctx->face_id;
1690 DMDALocalInfo* info = &user->info;
1691 PetscBool can_service;
1692
1693 // Suppress unused parameter warnings for clarity.
1694 (void)self;
1695 (void)local_inflow_contribution;
1696
1697 PetscFunctionBeginUser;
1698
1699 // Step 1: Use the robust utility function to determine if this MPI rank owns a computable
1700 // portion of the specified boundary face. If not, there is no work to do, so we exit immediately.
1701 const PetscInt IM_nodes_global = user->IM;
1702 const PetscInt JM_nodes_global = user->JM;
1703 const PetscInt KM_nodes_global = user->KM;
1704 ierr = CanRankServiceFace(info, IM_nodes_global, JM_nodes_global, KM_nodes_global, face_id, &can_service); CHKERRQ(ierr);
1705
1706 if (!can_service) {
1707 PetscFunctionReturn(0);
1708 }
1709
1710 // Step 2: Get read-only access to the necessary PETSc arrays.
1711 // We use the local versions (`lUcat`, `lNvert`) which include ghost cell data,
1712 // ensuring we have the correct interior values adjacent to the boundary.
1713 Cmpnts ***ucat, ***csi, ***eta, ***zet;
1714 PetscReal ***nvert;
1715 ierr = DMDAVecGetArrayRead(user->fda, user->lUcat, (const Cmpnts***)&ucat); CHKERRQ(ierr);
1716 ierr = DMDAVecGetArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
1717 ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
1718 ierr = DMDAVecGetArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
1719 ierr = DMDAVecGetArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
1720
1721 PetscReal local_flux_out = 0.0;
1722 const PetscInt xs=info->xs, xe=info->xs+info->xm;
1723 const PetscInt ys=info->ys, ye=info->ys+info->ym;
1724 const PetscInt zs=info->zs, ze=info->zs+info->zm;
1725 const PetscInt mx=info->mx, my=info->my, mz=info->mz;
1726
1727 // Step 3: Replicate the legacy shrunk loop bounds to exclude corners and edges.
1728 PetscInt lxs = xs; if (xs == 0) lxs = xs + 1;
1729 PetscInt lxe = xe; if (xe == mx) lxe = xe - 1;
1730 PetscInt lys = ys; if (ys == 0) lys = ys + 1;
1731 PetscInt lye = ye; if (ye == my) lye = ye - 1;
1732 PetscInt lzs = zs; if (zs == 0) lzs = zs + 1;
1733 PetscInt lze = ze; if (ze == mz) lze = ze - 1;
1734
1735 // Step 4: Loop over the specified face using the corrected bounds and indexing to calculate flux.
1736 switch (face_id) {
1737 case BC_FACE_NEG_X: {
1738 const PetscInt i_cell = xs + 1; // Index for first interior cell-centered data
1739 const PetscInt i_face = xs; // Index for the -X face of that cell
1740 for (int k=lzs; k<lze; k++) for (int j=lys; j<lye; j++) {
1741 if (nvert[k][j][i_cell] < 0.1) {
1742 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);
1743 }
1744 }
1745 break;
1746 }
1747 case BC_FACE_POS_X: {
1748 const PetscInt i_cell = xe - 2; // Index for last interior cell-centered data
1749 const PetscInt i_face = xe - 2; // Index for the +X face of that cell
1750 for (int k=lzs; k<lze; k++) for (int j=lys; j<lye; j++) {
1751 if (nvert[k][j][i_cell] < 0.1) {
1752 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);
1753 }
1754 }
1755 break;
1756 }
1757 case BC_FACE_NEG_Y: {
1758 const PetscInt j_cell = ys + 1;
1759 const PetscInt j_face = ys;
1760 for (int k=lzs; k<lze; k++) for (int i=lxs; i<lxe; i++) {
1761 if (nvert[k][j_cell][i] < 0.1) {
1762 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);
1763 }
1764 }
1765 break;
1766 }
1767 case BC_FACE_POS_Y: {
1768 const PetscInt j_cell = ye - 2;
1769 const PetscInt j_face = ye - 2;
1770 for (int k=lzs; k<lze; k++) for (int i=lxs; i<lxe; i++) {
1771 if (nvert[k][j_cell][i] < 0.1) {
1772 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);
1773 }
1774 }
1775 break;
1776 }
1777 case BC_FACE_NEG_Z: {
1778 const PetscInt k_cell = zs + 1;
1779 const PetscInt k_face = zs;
1780 for (int j=lys; j<lye; j++) for (int i=lxs; i<lxe; i++) {
1781 if (nvert[k_cell][j][i] < 0.1) {
1782 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);
1783 }
1784 }
1785 break;
1786 }
1787 case BC_FACE_POS_Z: {
1788 const PetscInt k_cell = ze - 2;
1789 const PetscInt k_face = ze - 2;
1790 for (int j=lys; j<lye; j++) for (int i=lxs; i<lxe; i++) {
1791 if (nvert[k_cell][j][i] < 0.1) {
1792 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);
1793 }
1794 }
1795 break;
1796 }
1797 }
1798
1799 // Step 5: Restore the PETSc arrays.
1800 ierr = DMDAVecRestoreArrayRead(user->fda, user->lUcat, (const Cmpnts***)&ucat); CHKERRQ(ierr);
1801 ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
1802 ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
1803 ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
1804 ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
1805
1806 // Step 6: Add this face's calculated flux to the accumulator for this rank.
1807 *local_outflow_contribution += local_flux_out;
1808
1809 PetscFunctionReturn(0);
1810}
1811
1812#undef __FUNCT__
1813#define __FUNCT__ "Apply_OutletConservation"
1814/**
1815 * @brief (Handler Action) Applies mass conservation correction to the outlet face.
1816 *
1817 * This function calculates a global correction factor based on the total inflow and outflow fluxes
1818 * and applies it to the contravariant velocity (`ucont`) on the outlet face to ensure mass conservation.
1819 */
1820static PetscErrorCode Apply_OutletConservation(BoundaryCondition *self, BCContext *ctx)
1821{
1822 PetscErrorCode ierr;
1823 (void)self;
1824 UserCtx* user = ctx->user;
1825 BCFace face_id = ctx->face_id;
1826 DMDALocalInfo* info = &user->info;
1827 PetscBool can_service;
1828
1829 PetscFunctionBeginUser;
1830 PROFILE_FUNCTION_BEGIN;
1831
1832 const PetscInt IM_nodes_global = user->IM;
1833 const PetscInt JM_nodes_global = user->JM;
1834 const PetscInt KM_nodes_global = user->KM;
1835 ierr = CanRankServiceFace(info, IM_nodes_global, JM_nodes_global, KM_nodes_global, face_id, &can_service); CHKERRQ(ierr);
1836
1837 if (!can_service) {
1838 PROFILE_FUNCTION_END;
1839 PetscFunctionReturn(0);
1840 }
1841
1842 // --- STEP 1: Calculate the correction factor using pre-calculated area ---
1843 PetscReal total_inflow = *ctx->global_inflow_sum + *ctx->global_farfield_inflow_sum;
1844 PetscReal flux_imbalance = total_inflow - *ctx->global_outflow_sum;
1845
1846 // Directly use the pre-calculated area from the simulation context.
1847 PetscReal velocity_correction = (PetscAbsReal(user->simCtx->AreaOutSum) > 1e-12)
1848 ? flux_imbalance / user->simCtx->AreaOutSum
1849 : 0.0;
1850
1851 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Outlet Correction on Face %d: Imbalance=%.4e, Pre-calc Area=%.4e, V_corr=%.4e\n",
1852 face_id, flux_imbalance, user->simCtx->AreaOutSum, velocity_correction);
1853
1854 // --- STEP 2: Apply the correction to ucont on the outlet face ---
1855
1856 // Get read/write access to necessary arrays
1857
1858 Cmpnts ***ubcs, ***ucont, ***csi, ***eta, ***zet, ***ucat;
1859 PetscReal ***nvert;
1860 ierr = DMDAVecGetArray(user->fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
1861 ierr = DMDAVecGetArray(user->fda, user->Ucont, &ucont); CHKERRQ(ierr);
1862 ierr = DMDAVecGetArrayRead(user->fda,user->lUcat, (const Cmpnts***)&ucat); CHKERRQ(ierr);
1863 ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
1864 ierr = DMDAVecGetArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
1865 ierr = DMDAVecGetArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
1866 ierr = DMDAVecGetArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
1867
1868 // Get local grid bounds to exclude corners/edges
1869 PetscInt xs = info->xs, xe = info->xs + info->xm;
1870 PetscInt ys = info->ys, ye = info->ys + info->ym;
1871 PetscInt zs = info->zs, ze = info->zs + info->zm;
1872 PetscInt mx = info->mx, my = info->my, mz = info->mz;
1873 PetscInt lxs = xs, lxe = xe, lys = ys, lye = ye, lzs = zs, lze = ze;
1874
1875 if (xs == 0) lxs = xs + 1;
1876 if (xe == mx) lxe = xe - 1;
1877 if (ys == 0) lys = ys + 1;
1878 if (ye == my) lye = ye - 1;
1879 if (zs == 0) lzs = zs + 1;
1880 if (ze == mz) lze = ze - 1;
1881
1882 // Loop over faces and apply correction
1883 switch(face_id){
1884 case BC_FACE_NEG_X:{
1885 const PetscInt i_cell = xs + 1;
1886 const PetscInt i_face = xs;
1887 const PetscInt i_dummy = xs;
1888 for (PetscInt k = lzs; k < lze; k++) {
1889 for (PetscInt j = lys; j < lye; j++) {
1890 if (nvert[k][j][i_cell] < 0.1) {
1891 // Set ubcs
1892 ubcs[k][j][i_dummy] = ucat[k][j][i_cell];
1893
1894 // Calculate Local uncorrected original flux
1895 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);
1896
1897 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));
1898
1899 PetscReal Correction_flux = velocity_correction*Cell_Area;
1900
1901 ucont[k][j][i_face].x = Uncorrected_local_flux + Correction_flux;
1902 }
1903 }
1904 }
1905 break;
1906 }
1907 case BC_FACE_POS_X:{
1908 const PetscInt i_cell = xe - 2;
1909 const PetscInt i_face = xe - 2;
1910 const PetscInt i_dummy = xe - 1;
1911 for(PetscInt k = lzs; k < lze; k++) for (PetscInt j = lys; j < lye; j++){
1912 if(nvert[k][j][i_cell]<0.1){
1913 // Set ubcs
1914 ubcs[k][j][i_dummy] = ucat[k][j][i_cell];
1915
1916 // Calculate Local uncorrected original flux
1917 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);
1918
1919 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));
1920
1921 PetscReal Correction_flux = velocity_correction*Cell_Area;
1922
1923 ucont[k][j][i_face].x = Uncorrected_local_flux + Correction_flux;
1924 }
1925 }
1926 break;
1927 }
1928 case BC_FACE_NEG_Y:{
1929 const PetscInt j_cell = ys + 1;
1930 const PetscInt j_face = ys;
1931 const PetscInt j_dummy = ys;
1932 for(PetscInt k = lzs; k < lze; k++) for (PetscInt i = lxs; i < lxe; i++){
1933 if(nvert[k][j_cell][i]<0.1){
1934 // Set ubcs
1935 ubcs[k][j_dummy][i] = ucat[k][j_cell][i];
1936
1937 // Calculate Local uncorrected original flux
1938 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);
1939
1940 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));
1941
1942 PetscReal Correction_flux = velocity_correction*Cell_Area;
1943
1944 ucont[k][j_face][i].y = Uncorrected_local_flux + Correction_flux;
1945 }
1946 }
1947 break;
1948 }
1949 case BC_FACE_POS_Y:{
1950 const PetscInt j_cell = ye - 2;
1951 const PetscInt j_face = ye - 2;
1952 const PetscInt j_dummy = ye - 1;
1953 for(PetscInt k = lzs; k < lze; k++) for (PetscInt i = lxs; i < lxe; i++){
1954 if(nvert[k][j_cell][i]<0.1){
1955 // Set ubcs
1956 ubcs[k][j_dummy][i] = ucat[k][j_cell][i];
1957
1958 // Calculate Local uncorrected original flux
1959 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);
1960
1961 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));
1962
1963 PetscReal Correction_flux = velocity_correction*Cell_Area;
1964
1965 ucont[k][j_face][i].y = Uncorrected_local_flux + Correction_flux;
1966 }
1967 }
1968 break;
1969 }
1970 case BC_FACE_NEG_Z:{
1971 const PetscInt k_cell = zs + 1;
1972 const PetscInt k_face = zs;
1973 const PetscInt k_dummy = zs;
1974 for(PetscInt j = lys; j < lye; j++) for (PetscInt i = lxs; i < lxe; i++){
1975 if(nvert[k_cell][j][i]<0.1){
1976 // Set ubcs
1977 ubcs[k_dummy][j][i] = ucat[k_cell][j][i];
1978
1979 // Calculate Local uncorrected original flux
1980 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));
1981
1982 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));
1983
1984 PetscReal Correction_flux = velocity_correction*Cell_Area;
1985
1986 ucont[k_face][j][i].z = Uncorrected_local_flux + Correction_flux;
1987 }
1988 }
1989 break;
1990 }
1991 case BC_FACE_POS_Z:{
1992 const PetscInt k_cell = ze - 2;
1993 const PetscInt k_face = ze - 2;
1994 const PetscInt k_dummy = ze - 1;
1995 for(PetscInt j = lys; j < lye; j++) for (PetscInt i = lxs; i < lxe; i++){
1996 if(nvert[k_cell][j][i]<0.1){
1997 // Set ubcs
1998 ubcs[k_dummy][j][i] = ucat[k_cell][j][i];
1999
2000 // Calculate Local uncorrected original flux
2001 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));
2002
2003 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));
2004
2005 PetscReal Correction_flux = velocity_correction*Cell_Area;
2006
2007 ucont[k_face][j][i].z = Uncorrected_local_flux + Correction_flux;
2008 }
2009 }
2010 break;
2011 }
2012 }
2013
2014 // Restore all arrays
2015 ierr = DMDAVecRestoreArray(user->fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
2016 ierr = DMDAVecRestoreArray(user->fda, user->Ucont, &ucont); CHKERRQ(ierr);
2017 ierr = DMDAVecRestoreArrayRead(user->fda,user->lUcat, (const Cmpnts***)&ucat); CHKERRQ(ierr);
2018 ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
2019 ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
2020 ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
2021 ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
2022
2023 PROFILE_FUNCTION_END;
2024 PetscFunctionReturn(0);
2025}
2026
2027#undef __FUNCT__
2028#define __FUNCT__ "PostStep_OutletConservation"
2029/**
2030 * @brief Internal helper implementation: `PostStep_OutletConservation()`.
2031 * @details Local to this translation unit.
2032 */
2033static PetscErrorCode PostStep_OutletConservation(BoundaryCondition *self, BCContext *ctx,
2034 PetscReal *local_inflow_contribution,
2035 PetscReal *local_outflow_contribution)
2036{
2037 PetscErrorCode ierr;
2038 UserCtx* user = ctx->user;
2039 BCFace face_id = ctx->face_id;
2040 DMDALocalInfo* info = &user->info;
2041 PetscBool can_service;
2042
2043 (void)self;
2044 (void)local_inflow_contribution;
2045
2046 PetscFunctionBeginUser;
2047 const PetscInt IM_nodes_global = user->IM;
2048 const PetscInt JM_nodes_global = user->JM;
2049 const PetscInt KM_nodes_global = user->KM;
2050 ierr = CanRankServiceFace(info, IM_nodes_global, JM_nodes_global, KM_nodes_global, face_id, &can_service); CHKERRQ(ierr);
2051
2052 if (!can_service) PetscFunctionReturn(0);
2053
2054 // Get arrays (need both ucont and nvert)
2055 Cmpnts ***ucont;
2056 PetscReal ***nvert; // ✅ ADD nvert
2057
2058 ierr = DMDAVecGetArrayRead(user->fda, user->Ucont, (const Cmpnts***)&ucont); CHKERRQ(ierr);
2059 ierr = DMDAVecGetArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr); // ✅ ADD
2060
2061 PetscReal local_flux = 0.0;
2062
2063 PetscInt xs = info->xs, xe = info->xs + info->xm;
2064 PetscInt ys = info->ys, ye = info->ys + info->ym;
2065 PetscInt zs = info->zs, ze = info->zs + info->zm;
2066 PetscInt mx = info->mx, my = info->my, mz = info->mz;
2067
2068 PetscInt lxs = xs, lxe = xe, lys = ys, lye = ye, lzs = zs, lze = ze;
2069 if (xs == 0) lxs = xs + 1;
2070 if (xe == mx) lxe = xe - 1;
2071 if (ys == 0) lys = ys + 1;
2072 if (ye == my) lye = ye - 1;
2073 if (zs == 0) lzs = zs + 1;
2074 if (ze == mz) lze = ze - 1;
2075
2076 // Sum ucont components, skipping solid cells (same indices as PreStep)
2077 switch (face_id) {
2078 case BC_FACE_NEG_X: {
2079 const PetscInt i_cell = xs + 1; // ✅ Match PreStep
2080 const PetscInt i_face = xs;
2081 for (PetscInt k = lzs; k < lze; k++) {
2082 for (PetscInt j = lys; j < lye; j++) {
2083 if (nvert[k][j][i_cell] < 0.1) { // ✅ Skip solid cells
2084 local_flux += ucont[k][j][i_face].x;
2085 }
2086 }
2087 }
2088 } break;
2089
2090 case BC_FACE_POS_X: {
2091 const PetscInt i_cell = xe - 2; // ✅ Match PreStep
2092 const PetscInt i_face = xe - 2;
2093 for (PetscInt k = lzs; k < lze; k++) {
2094 for (PetscInt j = lys; j < lye; j++) {
2095 if (nvert[k][j][i_cell] < 0.1) { // ✅ Skip solid cells
2096 local_flux += ucont[k][j][i_face].x;
2097 }
2098 }
2099 }
2100 } break;
2101
2102 case BC_FACE_NEG_Y: {
2103 const PetscInt j_cell = ys + 1; // ✅ Match PreStep
2104 const PetscInt j_face = ys;
2105 for (PetscInt k = lzs; k < lze; k++) {
2106 for (PetscInt i = lxs; i < lxe; i++) {
2107 if (nvert[k][j_cell][i] < 0.1) { // ✅ Skip solid cells
2108 local_flux += ucont[k][j_face][i].y;
2109 }
2110 }
2111 }
2112 } break;
2113
2114 case BC_FACE_POS_Y: {
2115 const PetscInt j_cell = ye - 2; // ✅ Match PreStep
2116 const PetscInt j_face = ye - 2;
2117 for (PetscInt k = lzs; k < lze; k++) {
2118 for (PetscInt i = lxs; i < lxe; i++) {
2119 if (nvert[k][j_cell][i] < 0.1) { // ✅ Skip solid cells
2120 local_flux += ucont[k][j_face][i].y;
2121 }
2122 }
2123 }
2124 } break;
2125
2126 case BC_FACE_NEG_Z: {
2127 const PetscInt k_cell = zs + 1; // ✅ Match PreStep
2128 const PetscInt k_face = zs;
2129 for (PetscInt j = lys; j < lye; j++) {
2130 for (PetscInt i = lxs; i < lxe; i++) {
2131 if (nvert[k_cell][j][i] < 0.1) { // ✅ Skip solid cells
2132 local_flux += ucont[k_face][j][i].z;
2133 }
2134 }
2135 }
2136 } break;
2137
2138 case BC_FACE_POS_Z: {
2139 const PetscInt k_cell = ze - 2; // ✅ Match PreStep
2140 const PetscInt k_face = ze - 2;
2141 for (PetscInt j = lys; j < lye; j++) {
2142 for (PetscInt i = lxs; i < lxe; i++) {
2143 if (nvert[k_cell][j][i] < 0.1) { // ✅ Skip solid cells
2144 local_flux += ucont[k_face][j][i].z;
2145 }
2146 }
2147 }
2148 } break;
2149 }
2150
2151 // Restore arrays
2152 ierr = DMDAVecRestoreArrayRead(user->fda, user->Ucont, (const Cmpnts***)&ucont); CHKERRQ(ierr);
2153 ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr); // ✅ ADD
2154
2155 // Add to accumulator
2156 *local_outflow_contribution += local_flux;
2157
2158 LOG_ALLOW(LOCAL, LOG_DEBUG, "PostStep_OutletConservation: Face %d, corrected flux = %.6e\n",
2159 face_id, local_flux);
2160
2161 PetscFunctionReturn(0);
2162}
2163
2164
2165/**
2166 * @brief Implementation of \ref Create_PeriodicGeometric().
2167 * @details Full API contract (arguments, ownership, side effects) is documented with
2168 * the header declaration in `include/BC_Handlers.h`.
2169 * @see Create_PeriodicGeometric()
2170 */
2171
2172PetscErrorCode Create_PeriodicGeometric(BoundaryCondition *bc){
2173 PetscFunctionBeginUser;
2174
2175 if (!bc) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Input BoundaryCondition is NULL");
2176 bc->priority = BC_PRIORITY_WALL;
2177
2178 // Assign function pointers
2179 bc->Initialize = NULL; // No initialization needed
2180 bc->PreStep = NULL;
2181 bc->Apply = NULL;
2182 bc->PostStep = NULL;
2183 bc->UpdateUbcs = NULL;
2184 bc->Destroy = NULL; // No private data to destroy
2185
2186 bc->data = NULL;
2187
2188 PetscFunctionReturn(0);
2189}
2190
2191
2192// ===============================================================================
2193//
2194// HANDLER IMPLEMENTATION: PERIODIC DRIVEN CONSTANT FLUX
2195// (Corresponds to BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX)
2196//
2197// ===============================================================================
2198
2199// --- 1. FORWARD DECLARATIONS & PRIVATE DATA ---
2200
2201// Forward declarations for the static functions that implement this handler's behavior.
2202static PetscErrorCode Initialize_PeriodicDrivenConstant(BoundaryCondition *self, BCContext *ctx);
2203static PetscErrorCode PreStep_PeriodicDrivenConstant(BoundaryCondition *self, BCContext *ctx, PetscReal *in, PetscReal *out);
2204static PetscErrorCode Apply_PeriodicDrivenConstant(BoundaryCondition *self, BCContext *ctx);
2205static PetscErrorCode Destroy_PeriodicDrivenConstant(BoundaryCondition *self);
2206
2207/** @brief Private data structure for the handler. */
2208typedef struct {
2209 char direction; // 'X', 'Y', or 'Z', determined at initialization.
2210 PetscReal targetVolumetricFlux; // The constant target flux, parsed from parameters.
2211 PetscReal boundaryVelocityCorrection; // The "Boundary Trim" value for the Apply() step.
2212 PetscBool isMasterController; // Flag: PETSC_TRUE only for the handler on the negative face.
2213 PetscBool applyBoundaryTrim; // Flag: PETSC_TRUE if applying Boundary trim on ucont.
2214} DrivenConstantData;
2215
2216
2217// --- 2. HANDLER CONSTRUCTOR ---
2218
2219#undef __FUNCT__
2220#define __FUNCT__ "Create_PeriodicDrivenConstant"
2221/**
2222 * @brief Internal helper implementation: `Create_PeriodicDrivenConstant()`.
2223 * @details Local to this translation unit.
2224 */
2225PetscErrorCode Create_PeriodicDrivenConstant(BoundaryCondition *bc)
2226{
2227 PetscErrorCode ierr;
2228 PetscFunctionBeginUser;
2229
2230 if (!bc) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Input BoundaryCondition object is NULL in Create_PeriodicDrivenConstantFlux");
2231
2232 // --- Allocate the private data structure ---
2233 DrivenConstantData *data = NULL;
2234 ierr = PetscNew(&data); CHKERRQ(ierr);
2235 // Initialize fields to safe default values
2236 data->direction = ' ';
2237 data->targetVolumetricFlux = 0.0;
2238 data->boundaryVelocityCorrection = 0.0;
2239 data->isMasterController = PETSC_FALSE;
2240 data->applyBoundaryTrim = PETSC_FALSE;
2241
2242 // Attach the private data to the generic handler object
2243 bc->data = (void*)data;
2244
2245 // --- Configure the handler's properties and methods ---
2246
2247 // Set priority: Using BC_PRIORITY_INLET ensures this handler's PreStep runs
2248 // before other handlers (like outlets) that might depend on its calculations.
2249 // It is the caller's responsibility that there are no Inlets called along with driven periodic to avoid clash.
2250 bc->priority = BC_PRIORITY_INLET;
2251
2252 // Assign the function pointers to the implementations in this file.
2253 bc->Initialize = Initialize_PeriodicDrivenConstant;
2254 bc->PreStep = PreStep_PeriodicDrivenConstant;
2255 bc->Apply = Apply_PeriodicDrivenConstant;
2256 bc->PostStep = NULL; // This handler has no action after the main solver step.
2257 bc->UpdateUbcs = NULL; // The boundary value is not flow-dependent (it's periodic).
2258 bc->Destroy = Destroy_PeriodicDrivenConstant;
2259
2260 PetscFunctionReturn(0);
2261}
2262
2263#undef __FUNCT__
2264#define __FUNCT__ "Initialize_PeriodicDrivenConstant"
2265/**
2266 * @brief Internal helper implementation: `Initialize_PeriodicDrivenConstant()`.
2267 * @details Local to this translation unit.
2268 */
2269static PetscErrorCode Initialize_PeriodicDrivenConstant(BoundaryCondition *self, BCContext *ctx)
2270{
2271 PetscErrorCode ierr;
2272 DrivenConstantData *data = (DrivenConstantData*)self->data;
2273 BCFace face_id = ctx->face_id;
2274 UserCtx* user = ctx->user;
2275
2276 PetscFunctionBeginUser;
2277
2278 LOG_ALLOW(LOCAL, LOG_DEBUG, "Initializing PERIODIC_DRIVEN_CONSTANT_FLUX handler on Face %s...\n", BCFaceToString(face_id));
2279
2280 // --- 1. Validation: Ensure the mathematical type is PERIODIC ---
2281 if (user->boundary_faces[face_id].mathematical_type != PERIODIC) {
2282 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT,
2283 "Configuration Error: Handler PERIODIC_DRIVEN_CONSTANT_FLUX on Face %s must be applied to a face with mathematical_type PERIODIC.",
2284 BCFaceToString(face_id));
2285 }
2286
2287 // --- 2. Role Assignment: Determine direction and master status ---
2288 data->isMasterController = PETSC_FALSE;
2289 switch (face_id) {
2290 case BC_FACE_NEG_X: data->direction = 'X'; data->isMasterController = PETSC_TRUE; break;
2291 case BC_FACE_POS_X: data->direction = 'X'; break;
2292 case BC_FACE_NEG_Y: data->direction = 'Y'; data->isMasterController = PETSC_TRUE; break;
2293 case BC_FACE_POS_Y: data->direction = 'Y'; break;
2294 case BC_FACE_NEG_Z: data->direction = 'Z'; data->isMasterController = PETSC_TRUE; break;
2295 case BC_FACE_POS_Z: data->direction = 'Z'; break;
2296 }
2297
2298 // --- 3. Parameter Parsing (Master Controller only) ---
2299 if (data->isMasterController) {
2300 PetscBool found;
2301
2302 // Attempt to read the 'target_flux' parameter from the bcs.run file.
2303 ierr = GetBCParamReal(user->boundary_faces[face_id].params, "target_flux",
2304 &data->targetVolumetricFlux, &found); CHKERRQ(ierr);
2305
2306 // If the required parameter is not found, halt with an informative error.
2307 if (!found) {
2308 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT,
2309 "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).",
2310 BCFaceToString(face_id));
2311 }
2312
2313 LOG_ALLOW(GLOBAL, LOG_INFO, "Driven Flow (Dir %c): Constant target volumetric flux set to %le.\n",
2314 data->direction, data->targetVolumetricFlux);
2315
2316 // Store the target flux in the UserCtx. This makes it globally accessible
2317 // to other parts of the solver, such as the `CorrectChannelFluxProfile` enforcer function.
2318 user->simCtx->targetVolumetricFlux = data->targetVolumetricFlux;
2319 }
2320
2321 PetscBool trimfound;
2322 // Attempt to read Trim flag.
2323 ierr = GetBCParamBool(user->boundary_faces[face_id].params, "apply_trim",
2324 &data->applyBoundaryTrim, &trimfound); CHKERRQ(ierr);
2325
2326 if(!trimfound) LOG_ALLOW(GLOBAL,LOG_DEBUG,"Trim Application not found,defaults to %s.\n",data->applyBoundaryTrim? "True":"False");
2327
2328 PetscFunctionReturn(0);
2329}
2330
2331#undef __FUNCT__
2332#define __FUNCT__ "PreStep_PeriodicDrivenConstant"
2333/**
2334 * @brief Internal helper implementation: `PreStep_PeriodicDrivenConstant()`.
2335 * @details Local to this translation unit.
2336 */
2337static PetscErrorCode PreStep_PeriodicDrivenConstant(BoundaryCondition *self, BCContext *ctx,
2338 PetscReal *local_inflow_contribution,
2339 PetscReal *local_outflow_contribution)
2340{
2341 PetscErrorCode ierr;
2342 DrivenConstantData *data = (DrivenConstantData*)self->data;
2343 UserCtx* user = ctx->user;
2344 SimCtx* simCtx = user->simCtx;
2345
2346 PetscFunctionBeginUser;
2347
2348 // --- Master Check: Only the handler on the negative face performs calculations ---
2349 if (!data->isMasterController) {
2350 PetscFunctionReturn(0);
2351 }
2352
2353 // This function is the generalized implementation of the old InitializeChannelFlowController.
2354 char direction = data->direction;
2355 DMDALocalInfo info = user->info;
2356 PetscInt i, j, k;
2357
2358 // --- Get read-only access to necessary field data ---
2359 Cmpnts ***ucont;
2360 PetscReal ***nvert;
2361 ierr = DMDAVecGetArrayRead(user->fda, user->lUcont, (const Cmpnts***)&ucont); CHKERRQ(ierr);
2362 ierr = DMDAVecGetArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
2363
2364 // --- Define local loop bounds ---
2365 PetscInt lxs = (info.xs == 0) ? 1 : info.xs;
2366 PetscInt lys = (info.ys == 0) ? 1 : info.ys;
2367 PetscInt lzs = (info.zs == 0) ? 1 : info.zs;
2368 PetscInt lxe = (info.xs + info.xm == info.mx) ? info.mx - 1 : info.xs + info.xm;
2369 PetscInt lye = (info.ys + info.ym == info.my) ? info.my - 1 : info.ys + info.ym;
2370 PetscInt lze = (info.zs + info.zm == info.mz) ? info.mz - 1 : info.zs + info.zm;
2371
2372 // ===================================================================================
2373 // CONTROLLER SENSORS: DUAL-FLUX MEASUREMENT STRATEGY
2374 //
2375 // To create a control system that is both stable and responsive, we measure the
2376 // volumetric flux in two distinct ways. One provides a fast, local signal for
2377 // immediate corrections, while the other provides a slow, stable signal for
2378 // strategic, domain-wide adjustments.
2379 //
2380 // -----------------------------------------------------------------------------------
2381 // ** 1. Fast/Local Sensor: `globalCurrentBoundaryFlux` **
2382 // -----------------------------------------------------------------------------------
2383 //
2384 // - WHAT IT IS: The total volumetric flux (m³/s) passing through the single,
2385 // representative periodic boundary plane at k=0.
2386 //
2387 // - PHYSICAL MEANING: An instantaneous "snapshot" of the flow rate at the
2388 // critical "seam" where the domain wraps around.
2389 //
2390 // - CHARACTERISTICS:
2391 // - Fast & Responsive: It immediately reflects any local flow changes or
2392 // errors occurring at the periodic interface.
2393 // - Noisy: It is susceptible to local fluctuations from turbulence or
2394 // numerical artifacts, causing its value to jitter from step to step.
2395 //
2396 // - CONTROLLER'S USE: This measurement is used to compute the
2397 // `boundaryVelocityCorrection`. This is a TACTICAL, fast-acting "trim" that is
2398 // applied immediately and directly to the boundary velocities to ensure perfect
2399 // continuity at the seam.
2400 //
2401 // -----------------------------------------------------------------------------------
2402 // ** 2. Stable/Global Sensor: `globalAveragePlanarVolumetricFlux` **
2403 // -----------------------------------------------------------------------------------
2404 //
2405 // - WHAT IT IS: The average of the volumetric fluxes across ALL cross-sectional
2406 // k-planes in the domain.
2407 // (i.e., [Flux(k=0) + Flux(k=1) + ... + Flux(k=N)] / N)
2408 //
2409 // - PHYSICAL MEANING: It represents the overall, bulk momentum of the fluid,
2410 // effectively filtering out local, transient fluctuations.
2411 //
2412 // - CHARACTERISTICS:
2413 // - Stable & Robust: Local noise at one plane is averaged out by all other
2414 // planes, providing a very smooth signal.
2415 // - Inertial: It changes more slowly, reflecting the inertia of the entire
2416 // fluid volume.
2417 //
2418 // - CONTROLLER'S USE: This measurement is used to compute the
2419 // `bulkVelocityCorrection`. This is a STRATEGIC, long-term adjustment used to
2420 // scale the main momentum source (body force). Using this stable signal prevents
2421 // the main driving force from oscillating wildly, ensuring simulation stability.
2422 //
2423 // ===================================================================================
2424
2425 // --- Initialize local accumulators ---
2426 PetscReal localCurrentBoundaryFlux = 0.0;
2427 PetscReal localAveragePlanarVolumetricFluxTerm = 0.0;
2428
2429 // --- Measure local contributions to the two flux types, generalized by direction ---
2430 switch (direction) {
2431 case 'X':
2432 if (info.xs == 0) { // Only the rank on the negative face contributes to boundary flux
2433 i = 0;
2434 for (k = lzs; k < lze; k++) for (j = lys; j < lye; j++) {
2435 if (nvert[k][j][i + 1] < 0.1) localCurrentBoundaryFlux += ucont[k][j][i].x;
2436 }
2437 }
2438 for (i = info.xs; i < lxe; i++) {
2439 for (k = lzs; k < lze; k++) for (j = lys; j < lye; j++) {
2440 if (nvert[k][j][i + 1] < 0.1) localAveragePlanarVolumetricFluxTerm += ucont[k][j][i].x / (PetscReal)(info.mx - 1);
2441 }
2442 }
2443 break;
2444 case 'Y':
2445 if (info.ys == 0) {
2446 j = 0;
2447 for (k = lzs; k < lze; k++) for (i = lxs; i < lxe; i++) {
2448 if (nvert[k][j + 1][i] < 0.1) localCurrentBoundaryFlux += ucont[k][j][i].y;
2449 }
2450 }
2451 for (j = info.ys; j < lye; j++) {
2452 for (k = lzs; k < lze; k++) for (i = lxs; i < lxe; i++) {
2453 if (nvert[k][j + 1][i] < 0.1) localAveragePlanarVolumetricFluxTerm += ucont[k][j][i].y / (PetscReal)(info.my - 1);
2454 }
2455 }
2456 break;
2457 case 'Z':
2458 if (info.zs == 0) {
2459 k = 0;
2460 for (j = lys; j < lye; j++) for (i = lxs; i < lxe; i++) {
2461 if (nvert[k + 1][j][i] < 0.1) localCurrentBoundaryFlux += ucont[k][j][i].z;
2462 }
2463 }
2464 for (k = info.zs; k < lze; k++) {
2465 for (j = lys; j < lye; j++) for (i = lxs; i < lxe; i++) {
2466 if (nvert[k + 1][j][i] < 0.1) localAveragePlanarVolumetricFluxTerm += ucont[k][j][i].z / (PetscReal)(info.mz - 1);
2467 }
2468 }
2469 break;
2470 }
2471
2472 // --- Release array access as soon as possible ---
2473 ierr = DMDAVecRestoreArrayRead(user->fda, user->lUcont, (const Cmpnts***)&ucont); CHKERRQ(ierr);
2474 ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
2475
2476 // --- Get cross-sectional area using the dedicated geometry function ---
2477 Cmpnts ignored_center;
2478 PetscReal globalBoundaryArea;
2479 BCFace neg_face_id = (direction == 'X') ? BC_FACE_NEG_X : (direction == 'Y') ? BC_FACE_NEG_Y : BC_FACE_NEG_Z;
2480 ierr = CalculateFaceCenterAndArea(user, neg_face_id, &ignored_center, &globalBoundaryArea); CHKERRQ(ierr);
2481
2482 // --- Perform global reductions to get the final flux values ---
2483 PetscReal globalCurrentBoundaryFlux, globalAveragePlanarVolumetricFlux;
2484 ierr = MPI_Allreduce(&localCurrentBoundaryFlux, &globalCurrentBoundaryFlux, 1, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
2485 ierr = MPI_Allreduce(&localAveragePlanarVolumetricFluxTerm, &globalAveragePlanarVolumetricFlux, 1, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
2486
2487 // --- Calculate the two correction terms ---
2488 if (globalBoundaryArea > 1.0e-12) {
2489 // The main correction for the body force is based on the STABLE domain-averaged flux.
2490 simCtx->bulkVelocityCorrection = (data->targetVolumetricFlux - globalAveragePlanarVolumetricFlux) / globalBoundaryArea;
2491
2492 // The immediate correction for the boundary trim is based on the FAST boundary-specific flux.
2493 data->boundaryVelocityCorrection = (data->targetVolumetricFlux - globalCurrentBoundaryFlux) / globalBoundaryArea;
2494 } else {
2495 simCtx->bulkVelocityCorrection = 0.0;
2496 data->boundaryVelocityCorrection = 0.0;
2497 }
2498
2499 LOG_ALLOW(GLOBAL, LOG_INFO, "Driven Flow Controller Update (Dir %c):\n", data->direction);
2500 LOG_ALLOW(GLOBAL, LOG_INFO, " - Target Volumetric Flux: %.6e\n", data->targetVolumetricFlux);
2501 LOG_ALLOW(GLOBAL, LOG_INFO, " - Avg Planar Volumetric Flux (Stable): %.6e\n", globalAveragePlanarVolumetricFlux);
2502 LOG_ALLOW(GLOBAL, LOG_INFO, " - Boundary Flux (Fast): %.6e\n", globalCurrentBoundaryFlux);
2503 LOG_ALLOW(GLOBAL, LOG_INFO, " - Bulk Velocity Correction: %.6e (For Momentum Source)\n", simCtx->bulkVelocityCorrection);
2504 LOG_ALLOW(GLOBAL, LOG_INFO, " - Boundary Velocity Correction: %.6e (For Boundary Trim)\n", data->boundaryVelocityCorrection);
2505
2506 // Suppress unused parameter warnings for this handler
2507 (void)local_inflow_contribution;
2508 (void)local_outflow_contribution;
2509
2510 PetscFunctionReturn(0);
2511}
2512
2513#undef __FUNCT__
2514#define __FUNCT__ "Apply_PeriodicDrivenConstant"
2515/**
2516 * @brief Internal helper implementation: `Apply_PeriodicDrivenConstant()`.
2517 * @details Local to this translation unit.
2518 */
2519static PetscErrorCode Apply_PeriodicDrivenConstant(BoundaryCondition *self, BCContext *ctx)
2520{
2521 PetscErrorCode ierr;
2522 DrivenConstantData *data = (DrivenConstantData*)self->data;
2523 UserCtx* user = ctx->user;
2524 BCFace face_id = ctx->face_id;
2525 PetscBool can_service;
2526
2527 PetscFunctionBeginUser;
2528
2529 // Check if this rank owns part of this boundary face
2530 ierr = CanRankServiceFace(&user->info, user->IM, user->JM, user->KM, face_id, &can_service); CHKERRQ(ierr);
2531 if (!can_service) {
2532 PetscFunctionReturn(0);
2533 }
2534
2535 // If the correction is negligible, no work is needed.
2536 if (PetscAbsReal(data->boundaryVelocityCorrection) < 1e-12) {
2537 PetscFunctionReturn(0);
2538 }
2539
2540 LOG_ALLOW(LOCAL, LOG_TRACE, "Apply_PeriodicDrivenConstant: Applying boundary trim on Face %s...\n", BCFaceToString(face_id));
2541
2542 // --- Get read/write access to necessary arrays ---
2543 DMDALocalInfo info = user->info;
2544 Cmpnts ***ucont, ***uch, ***csi, ***eta, ***zet;
2545 PetscReal ***nvert;
2546
2547 ierr = DMDAVecGetArray(user->fda, user->lUcont, &ucont); CHKERRQ(ierr);
2548 ierr = DMDAVecGetArray(user->fda, user->Bcs.Uch, &uch); CHKERRQ(ierr);
2549 ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
2550 ierr = DMDAVecGetArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
2551 ierr = DMDAVecGetArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
2552 ierr = DMDAVecGetArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
2553
2554 PetscInt lxs = (info.xs == 0) ? 1 : info.xs;
2555 PetscInt lys = (info.ys == 0) ? 1 : info.ys;
2556 PetscInt lzs = (info.zs == 0) ? 1 : info.zs;
2557 PetscInt lxe = (info.xs + info.xm == info.mx) ? info.mx - 1 : info.xs + info.xm;
2558 PetscInt lye = (info.ys + info.ym == info.my) ? info.my - 1 : info.ys + info.ym;
2559 PetscInt lze = (info.zs + info.zm == info.mz) ? info.mz - 1 : info.zs + info.zm;
2560
2561 // --- Apply correction to the appropriate face and velocity component ---
2562 switch (face_id) {
2563 case BC_FACE_NEG_X: case BC_FACE_POS_X: {
2564 PetscInt i_face = (face_id == BC_FACE_NEG_X) ? info.xs : info.mx - 2;
2565 PetscInt i_nvert = (face_id == BC_FACE_NEG_X) ? info.xs + 1 : info.mx - 2;
2566
2567 for (PetscInt k = lzs; k < lze; k++) for (PetscInt j = lys; j < lye; j++) {
2568 if (nvert[k][j][i_nvert] < 0.1) {
2569 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);
2570 PetscReal fluxTrim = data->boundaryVelocityCorrection * faceArea;
2571 if(data->applyBoundaryTrim) ucont[k][j][i_face].x += fluxTrim;
2572 uch[k][j][i_face].x = fluxTrim; // Store correction for diagnostics
2573 }
2574 }
2575 } break;
2576
2577 case BC_FACE_NEG_Y: case BC_FACE_POS_Y: {
2578 PetscInt j_face = (face_id == BC_FACE_NEG_Y) ? info.ys : info.my - 2;
2579 PetscInt j_nvert = (face_id == BC_FACE_NEG_Y) ? info.ys + 1 : info.my - 2;
2580
2581 for (PetscInt k = lzs; k < lze; k++) for (PetscInt i = lxs; i < lxe; i++) {
2582 if (nvert[k][j_nvert][i] < 0.1) {
2583 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);
2584 PetscReal fluxTrim = data->boundaryVelocityCorrection * faceArea;
2585 if(data->applyBoundaryTrim) ucont[k][j_face][i].y += fluxTrim;
2586 uch[k][j_face][i].y = fluxTrim;
2587 }
2588 }
2589 } break;
2590
2591 case BC_FACE_NEG_Z: case BC_FACE_POS_Z: {
2592 PetscInt k_face = (face_id == BC_FACE_NEG_Z) ? info.zs : info.mz - 2;
2593 PetscInt k_nvert = (face_id == BC_FACE_NEG_Z) ? info.zs + 1 : info.mz - 2;
2594
2595 for (PetscInt j = lys; j < lye; j++) for (PetscInt i = lxs; i < lxe; i++) {
2596 if (nvert[k_nvert][j][i] < 0.1) {
2597 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);
2598 PetscReal fluxTrim = data->boundaryVelocityCorrection * faceArea;
2599 if(data->applyBoundaryTrim) ucont[k_face][j][i].z += fluxTrim;
2600 uch[k_face][j][i].z = fluxTrim;
2601 }
2602 }
2603 } break;
2604 }
2605
2606 // --- Restore arrays ---
2607 ierr = DMDAVecRestoreArray(user->fda, user->lUcont, &ucont); CHKERRQ(ierr);
2608 ierr = DMDAVecRestoreArray(user->fda, user->Bcs.Uch, &uch); CHKERRQ(ierr);
2609 ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
2610 ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
2611 ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
2612 ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
2613
2614 PetscFunctionReturn(0);
2615}
2616
2617#undef __FUNCT__
2618#define __FUNCT__ "Destroy_PeriodicDrivenConstant"
2619/**
2620 * @brief Internal helper implementation: `Destroy_PeriodicDrivenConstant()`.
2621 * @details Local to this translation unit.
2622 */
2623static PetscErrorCode Destroy_PeriodicDrivenConstant(BoundaryCondition *self)
2624{
2625 PetscFunctionBeginUser;
2626
2627 // Check that the handler object and its private data pointer are valid before trying to free.
2628 if (self && self->data) {
2629 // The private data was allocated with PetscNew(), so it must be freed with PetscFree().
2630 PetscFree(self->data);
2631
2632 // It is good practice to nullify the pointer after freeing to prevent
2633 // any accidental use of the dangling pointer (use-after-free).
2634 self->data = NULL;
2635
2636 LOG_ALLOW(LOCAL, LOG_TRACE, "Destroy_PeriodicDrivenConstant: Private data freed successfully.\n");
2637 }
2638
2639 PetscFunctionReturn(0);
2640}
InletParabolicData::cs2_half
PetscReal cs2_half
Half-width (in index space) in cross-stream direction 2.
Definition BC_Handlers.c:710
Initialize_InletProfileFromFile
static PetscErrorCode Initialize_InletProfileFromFile(BoundaryCondition *self, BCContext *ctx)
Initializes a file-prescribed inlet profile handler for one boundary face.
Definition BC_Handlers.c:1344
Destroy_InletProfileFromFile
static PetscErrorCode Destroy_InletProfileFromFile(BoundaryCondition *self)
Releases private storage owned by a file-prescribed inlet profile handler.
Definition BC_Handlers.c:1611
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:2213
Create_InletProfileFromFile
PetscErrorCode Create_InletProfileFromFile(BoundaryCondition *bc)
Implementation of Create_InletProfileFromFile().
Definition BC_Handlers.c:1303
Initialize_InletParabolicProfile
static PetscErrorCode Initialize_InletParabolicProfile(BoundaryCondition *self, BCContext *ctx)
Internal helper implementation: Initialize_InletParabolicProfile().
Definition BC_Handlers.c:750
GetProfileFileExpectedDims
static PetscErrorCode GetProfileFileExpectedDims(UserCtx *user, BCFace face_id, PetscInt *n1, PetscInt *n2)
Computes the expected PICSLICE dimensions for an inlet face.
Definition BC_Handlers.c:1160
PreStep_InletConstantVelocity
static PetscErrorCode PreStep_InletConstantVelocity(BoundaryCondition *self, BCContext *ctx, PetscReal *in, PetscReal *out)
Internal helper implementation: PreStep_InletConstantVelocity().
Definition BC_Handlers.c:402
PreStep_InletProfileFromFile
static PetscErrorCode PreStep_InletProfileFromFile(BoundaryCondition *self, BCContext *ctx, PetscReal *in, PetscReal *out)
Pre-step hook for the static file-prescribed inlet profile handler.
Definition BC_Handlers.c:1389
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:1820
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
InletProfileFileData::n2
PetscInt n2
Definition BC_Handlers.c:1119
DrivenConstantData::targetVolumetricFlux
PetscReal targetVolumetricFlux
Definition BC_Handlers.c:2210
DrivenConstantData::direction
char direction
Definition BC_Handlers.c:2209
InletProfileFileData::profile
PetscReal * profile
Definition BC_Handlers.c:1120
ReadPicSliceProfile
static PetscErrorCode ReadPicSliceProfile(const char *source_file, PetscInt expected_n1, PetscInt expected_n2, InletProfileFileData *data)
Reads and validates a static scalar inlet profile from a canonical PICSLICE file.
Definition BC_Handlers.c:1205
DrivenConstantData::isMasterController
PetscBool isMasterController
Definition BC_Handlers.c:2212
Apply_PeriodicDrivenConstant
static PetscErrorCode Apply_PeriodicDrivenConstant(BoundaryCondition *self, BCContext *ctx)
Internal helper implementation: Apply_PeriodicDrivenConstant().
Definition BC_Handlers.c:2519
Initialize_PeriodicDrivenConstant
static PetscErrorCode Initialize_PeriodicDrivenConstant(BoundaryCondition *self, BCContext *ctx)
Internal helper implementation: Initialize_PeriodicDrivenConstant().
Definition BC_Handlers.c:2269
Create_PeriodicGeometric
PetscErrorCode Create_PeriodicGeometric(BoundaryCondition *bc)
Implementation of Create_PeriodicGeometric().
Definition BC_Handlers.c:2172
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:2623
PostStep_InletProfileFromFile
static PetscErrorCode PostStep_InletProfileFromFile(BoundaryCondition *self, BCContext *ctx, PetscReal *in, PetscReal *out)
Accumulates the applied inlet flux for a file-prescribed profile.
Definition BC_Handlers.c:1535
Create_PeriodicDrivenConstant
PetscErrorCode Create_PeriodicDrivenConstant(BoundaryCondition *bc)
Internal helper implementation: Create_PeriodicDrivenConstant().
Definition BC_Handlers.c:2225
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:2033
InletProfileFileData::max_speed
PetscReal max_speed
Definition BC_Handlers.c:1122
InletConstantData::normal_velocity
PetscReal normal_velocity
Definition BC_Handlers.c:312
ProfileSpeedAt
static PetscReal ProfileSpeedAt(const InletProfileFileData *data, PetscInt a, PetscInt b)
Returns one scalar speed from the flattened PICSLICE profile.
Definition BC_Handlers.c:1290
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
Apply_InletProfileFromFile
static PetscErrorCode Apply_InletProfileFromFile(BoundaryCondition *self, BCContext *ctx)
Applies the loaded PICSLICE scalar profile to Ucont and Ubcs on an inlet face.
Definition BC_Handlers.c:1415
PreStep_PeriodicDrivenConstant
static PetscErrorCode PreStep_PeriodicDrivenConstant(BoundaryCondition *self, BCContext *ctx, PetscReal *in, PetscReal *out)
Internal helper implementation: PreStep_PeriodicDrivenConstant().
Definition BC_Handlers.c:2337
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:1684
DrivenConstantData::boundaryVelocityCorrection
PetscReal boundaryVelocityCorrection
Definition BC_Handlers.c:2211
Create_OutletConservation
PetscErrorCode Create_OutletConservation(BoundaryCondition *bc)
Implementation of Create_OutletConservation().
Definition BC_Handlers.c:1655
InletProfileFileData::n1
PetscInt n1
Definition BC_Handlers.c:1118
GetBCParamStringLocal
static PetscErrorCode GetBCParamStringLocal(BC_Param *params, const char *key, const char **value_out, PetscBool *found)
Looks up a string-valued boundary-condition parameter in a BC_Param list.
Definition BC_Handlers.c:1135
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
InletProfileFileData::min_speed
PetscReal min_speed
Definition BC_Handlers.c:1121
InletProfileFileData::source_file
char * source_file
Definition BC_Handlers.c:1123
DrivenConstantData
Private data structure for the handler.
Definition BC_Handlers.c:2208
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
InletProfileFileData
Definition BC_Handlers.c:1117
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:390
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:411
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:820
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:811
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:848
SimCtx::targetVolumetricFlux
PetscReal targetVolumetricFlux
Definition variables.h:738
UserCtx::lNvert
Vec lNvert
Definition variables.h:856
UserCtx::simCtx
SimCtx * simCtx
Back-pointer to the master simulation context.
Definition variables.h:833
BC_Param_s::next
struct BC_Param_s * next
Definition variables.h:307
UserCtx::KM
PetscInt KM
Definition variables.h:839
UserCtx::lZet
Vec lZet
Definition variables.h:879
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:739
BCContext::face_id
BCFace face_id
Definition variables.h:313
UserCtx::Ucont
Vec Ucont
Definition variables.h:856
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:851
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:879
BoundaryFaceConfig::params
BC_Param * params
Definition variables.h:338
Cmpnts::z
PetscScalar z
Definition variables.h:101
UserCtx::JM
PetscInt JM
Definition variables.h:839
BCContext::global_farfield_inflow_sum
const PetscReal * global_farfield_inflow_sum
Definition variables.h:315
UserCtx::lUcont
Vec lUcont
Definition variables.h:856
SimCtx::AreaOutSum
PetscReal AreaOutSum
Definition variables.h:740
UserCtx::info
DMDALocalInfo info
Definition variables.h:837
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:856
Cmpnts::y
PetscScalar y
Definition variables.h:101
UserCtx::IM
PetscInt IM
Definition variables.h:839
UserCtx::lEta
Vec lEta
Definition variables.h:879
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
BC_Param_s
A node in a linked list for storing key-value parameters from the bcs.dat file.
Definition variables.h:304
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:651
UserCtx
User-defined context containing data specific to a single computational grid level.
Definition variables.h:830
Generated by doxygen 1.9.8