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