PICurv 0.1.0
A Parallel Particle-In-Cell Solver for Curvilinear LES
Loading...
Searching...
No Matches
BC_Handlers.c
Go to the documentation of this file.
1#include "BC_Handlers.h" // The header that declares this file's "constructor" functions
2
3
4//================================================================================
5// VALIDATORS
6//================================================================================
7
8
9#undef __FUNCT__
10#define __FUNCT__ "Validate_DrivenFlowConfiguration"
11/**
12 * @brief Internal helper implementation: `Validate_DrivenFlowConfiguration()`.
13 * @details Local to this translation unit.
14 */
15PetscErrorCode Validate_DrivenFlowConfiguration(UserCtx *user)
16{
17 PetscFunctionBeginUser;
18
19 // --- CHECK 1: Detect if a driven flow is active. ---
20 PetscBool is_driven_flow_active = PETSC_FALSE;
21 char driven_direction = ' ';
22 const char* first_driven_face_name = "";
23
24 for (int i = 0; i < 6; i++) {
25 BCHandlerType handler_type = user->boundary_faces[i].handler_type;
26 if (handler_type == BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX ||
27 handler_type == BC_HANDLER_PERIODIC_DRIVEN_INITIAL_FLUX)
28 {
29 is_driven_flow_active = PETSC_TRUE;
30 first_driven_face_name = BCFaceToString((BCFace)i);
31
32 if (i <= 1) driven_direction = 'X';
33 else if (i <= 3) driven_direction = 'Y';
34 else driven_direction = 'Z';
35
36 break; // Exit loop once we've confirmed it's active and found the direction.
37 }
38 }
39
40 // If no driven flow handler is found, validation for this rule set is complete.
41 if (!is_driven_flow_active) {
42 PetscFunctionReturn(0);
43 }
44
45 LOG_ALLOW(GLOBAL, LOG_DEBUG, " - Driven Flow Handler detected on face %s. Applying driven flow validation rules...\n", first_driven_face_name);
46
47 // --- CHECK 2: Ensure no conflicting BCs (Inlet/Outlet/Far-field) are present. ---
48 LOG_ALLOW(GLOBAL, LOG_DEBUG, " - Checking for incompatible Inlet/Outlet/Far-field BCs...\n");
49 for (int i = 0; i < 6; i++) {
50 BCType math_type = user->boundary_faces[i].mathematical_type;
51 if (math_type == INLET || math_type == OUTLET || math_type == FARFIELD) {
52 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT,
53 "Configuration Error: A DRIVEN flow handler is active, which is incompatible with the %s boundary condition found on face %s.",
54 BCTypeToString(math_type), BCFaceToString((BCFace)i));
55 }
56 }
57 LOG_ALLOW(GLOBAL, LOG_DEBUG, " ... No conflicting BC types found. OK.\n");
58
59 // --- CHECK 3: Ensure both ends of the driven direction have identical, valid setups. ---
60 LOG_ALLOW(GLOBAL, LOG_DEBUG, " - Validating symmetry and mathematical types for the '%c' direction...\n", driven_direction);
61
62 PetscInt neg_face_idx = 0, pos_face_idx = 0;
63 if (driven_direction == 'X') {
64 neg_face_idx = BC_FACE_NEG_X; pos_face_idx = BC_FACE_POS_X;
65 } else if (driven_direction == 'Y') {
66 neg_face_idx = BC_FACE_NEG_Y; pos_face_idx = BC_FACE_POS_Y;
67 } else { // 'Z'
68 neg_face_idx = BC_FACE_NEG_Z; pos_face_idx = BC_FACE_POS_Z;
69 }
70
71 BoundaryFaceConfig *neg_face_cfg = &user->boundary_faces[neg_face_idx];
72 BoundaryFaceConfig *pos_face_cfg = &user->boundary_faces[pos_face_idx];
73
74 // Rule 3a: Both faces must be PERIODIC.
75 if (neg_face_cfg->mathematical_type != PERIODIC || pos_face_cfg->mathematical_type != PERIODIC) {
76 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT,
77 "Configuration Error: For a driven flow in the '%c' direction, both the %s and %s faces must be of mathematical_type PERIODIC.",
78 driven_direction, BCFaceToString((BCFace)neg_face_idx), BCFaceToString((BCFace)pos_face_idx));
79 }
80
81 // Rule 3b: Both faces must use the exact same handler type.
82 if (neg_face_cfg->handler_type != pos_face_cfg->handler_type) {
83 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT,
84 "Configuration Error: The DRIVEN handlers on the %s and %s faces of the '%c' direction do not match. Both must be the same type (e.g., both CONSTANT_FLUX).",
85 BCFaceToString((BCFace)neg_face_idx), BCFaceToString((BCFace)pos_face_idx), driven_direction);
86 }
87
88 LOG_ALLOW(GLOBAL, LOG_DEBUG, " ... Symmetry and mathematical types are valid. OK.\n");
89
90 PetscFunctionReturn(0);
91}
92
93//================================================================================
94//
95// HANDLER IMPLEMENTATION: NO-SLIP WALL
96// (Corresponds to BC_HANDLER_WALL_NOSLIP)
97//
98// This handler implements a stationary, impenetrable wall where the fluid
99// velocity is zero (no-slip condition).
100//
101//================================================================================
102
103// --- FORWARD DECLARATIONS ---
104static PetscErrorCode Apply_WallNoSlip(BoundaryCondition *self, BCContext *ctx);
105
106#undef __FUNCT__
107#define __FUNCT__ "Create_WallNoSlip"
108/**
109 * @brief Implementation of \ref Create_WallNoSlip().
110 * @details Full API contract (arguments, ownership, side effects) is documented with
111 * the header declaration in `include/BC_Handlers.h`.
112 * @see Create_WallNoSlip()
113 */
114PetscErrorCode Create_WallNoSlip(BoundaryCondition *bc)
115{
116 PetscFunctionBeginUser;
117
118 if (!bc) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
119 "Input BoundaryCondition object is NULL in Create_WallNoSlip");
120
121 // ✅ Set priority
122 bc->priority = BC_PRIORITY_WALL;
123
124 // Assign function pointers
125 bc->Initialize = NULL;
126 bc->PreStep = NULL;
127 bc->Apply = Apply_WallNoSlip;
128 bc->PostStep = NULL;
129 bc->UpdateUbcs = NULL;
130 bc->Destroy = NULL;
131
132 // No private data needed for this simple handler
133 bc->data = NULL;
134
135 PetscFunctionReturn(0);
136}
137
138#undef __FUNCT__
139#define __FUNCT__ "Apply_WallNoSlip"
140/**
141 * @brief Apply no-slip velocity values to wall-adjacent cells for this boundary condition.
142 */
143static PetscErrorCode Apply_WallNoSlip(BoundaryCondition *self, BCContext *ctx)
144{
145 PetscErrorCode ierr;
146 UserCtx* user = ctx->user;
147 BCFace face_id = ctx->face_id;
148 PetscBool can_service;
149
150 (void)self; // Unused for simple handlers
151
152 PetscFunctionBeginUser;
153 DMDALocalInfo *info = &user->info;
154 Cmpnts ***ubcs, ***ucont;
155 PetscInt IM_nodes_global, JM_nodes_global,KM_nodes_global;
156
157 IM_nodes_global = user->IM;
158 JM_nodes_global = user->JM;
159 KM_nodes_global = user->KM;
160
161 ierr = CanRankServiceFace(info,IM_nodes_global,JM_nodes_global,KM_nodes_global,face_id,&can_service); CHKERRQ(ierr);
162 // Check if this rank owns part of this boundary face
163 if (!can_service) PetscFunctionReturn(0);
164
165 LOG_ALLOW(LOCAL, LOG_DEBUG, "Apply_WallNoSlip: Applying to Face %d (%s).\n",
166 face_id, BCFaceToString(face_id));
167
168 // Get arrays
169
170 ierr = DMDAVecGetArray(user->fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
171 ierr = DMDAVecGetArray(user->fda, user->Ucont, &ucont); CHKERRQ(ierr);
172
173 PetscInt xs = info->xs, xe = info->xs + info->xm;
174 PetscInt ys = info->ys, ye = info->ys + info->ym;
175 PetscInt zs = info->zs, ze = info->zs + info->zm;
176 PetscInt mx = info->mx, my = info->my, mz = info->mz;
177
178 // ✅ Use shrunken loop bounds (avoids edges/corners like inlet handler)
179 PetscInt lxs = xs, lxe = xe, lys = ys, lye = ye, lzs = zs, lze = ze;
180 if (xs == 0) lxs = xs + 1;
181 if (xe == mx) lxe = xe - 1;
182 if (ys == 0) lys = ys + 1;
183 if (ye == my) lye = ye - 1;
184 if (zs == 0) lzs = zs + 1;
185 if (ze == mz) lze = ze - 1;
186
187 switch (face_id) {
188 case BC_FACE_NEG_X: {
189 if (xs == 0){
190 PetscInt i = xs;
191 for (PetscInt k = lzs; k < lze; k++) {
192 for (PetscInt j = lys; j < lye; j++) {
193 // ✅ Set contravariant flux to zero (no penetration)
194 ucont[k][j][i].x = 0.0;
195
196 // ✅ Set boundary velocity to zero (no slip)
197 ubcs[k][j][i].x = 0.0;
198 ubcs[k][j][i].y = 0.0;
199 ubcs[k][j][i].z = 0.0;
200 }
201 }
202 }
203 break;
204 }
205
206 case BC_FACE_POS_X: {
207 if (xe == mx){
208 PetscInt i = xe - 1;
209 for (PetscInt k = lzs; k < lze; k++) {
210 for (PetscInt j = lys; j < lye; j++) {
211 ucont[k][j][i-1].x = 0.0;
212
213 ubcs[k][j][i].x = 0.0;
214 ubcs[k][j][i].y = 0.0;
215 ubcs[k][j][i].z = 0.0;
216 }
217 }
218 }
219 break;
220 }
221 case BC_FACE_NEG_Y: {
222 if (ys == 0){
223 PetscInt j = ys;
224 for (PetscInt k = lzs; k < lze; k++) {
225 for (PetscInt i = lxs; i < lxe; i++) {
226 ucont[k][j][i].y = 0.0;
227
228 ubcs[k][j][i].x = 0.0;
229 ubcs[k][j][i].y = 0.0;
230 ubcs[k][j][i].z = 0.0;
231 }
232 }
233 }
234 } break;
235
236 case BC_FACE_POS_Y: {
237 if (ye == my){
238 PetscInt j = ye - 1;
239 for (PetscInt k = lzs; k < lze; k++) {
240 for (PetscInt i = lxs; i < lxe; i++) {
241 ucont[k][j-1][i].y = 0.0;
242
243 ubcs[k][j][i].x = 0.0;
244 ubcs[k][j][i].y = 0.0;
245 ubcs[k][j][i].z = 0.0;
246 }
247 }
248 }
249 } break;
250
251 case BC_FACE_NEG_Z: {
252 if (zs == 0){
253 PetscInt k = zs;
254 for (PetscInt j = lys; j < lye; j++) {
255 for (PetscInt i = lxs; i < lxe; i++) {
256 ucont[k][j][i].z = 0.0;
257
258 ubcs[k][j][i].x = 0.0;
259 ubcs[k][j][i].y = 0.0;
260 ubcs[k][j][i].z = 0.0;
261 }
262 }
263 }
264 } break;
265
266 case BC_FACE_POS_Z: {
267 if (ze == mz) {
268 PetscInt k = ze - 1;
269 for (PetscInt j = lys; j < lye; j++) {
270 for (PetscInt i = lxs; i < lxe; i++) {
271 ucont[k-1][j][i].z = 0.0;
272
273 ubcs[k][j][i].x = 0.0;
274 ubcs[k][j][i].y = 0.0;
275 ubcs[k][j][i].z = 0.0;
276 }
277 }
278 }
279 } break;
280 }
281
282 // Restore arrays
283 ierr = DMDAVecRestoreArray(user->fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
284 ierr = DMDAVecRestoreArray(user->fda, user->Ucont, &ucont); CHKERRQ(ierr);
285
286 PetscFunctionReturn(0);
287}
288
289////////////////////////////////////////////////////////////////////////
290
291//================================================================================
292//
293// HANDLER IMPLEMENTATION: CONSTANT VELOCITY INLET
294// (Corresponds to BC_HANDLER_INLET_CONSTANT_VELOCITY)
295//
296//================================================================================
297
298// --- FORWARD DECLARATIONS ---
299static PetscErrorCode Initialize_InletConstantVelocity(BoundaryCondition *self, BCContext *ctx);
300static PetscErrorCode PreStep_InletConstantVelocity(BoundaryCondition *self, BCContext *ctx,
301 PetscReal *in, PetscReal *out);
302static PetscErrorCode Apply_InletConstantVelocity(BoundaryCondition *self, BCContext *ctx);
303static PetscErrorCode PostStep_InletConstantVelocity(BoundaryCondition *self, BCContext *ctx,
304 PetscReal *in, PetscReal *out);
305static PetscErrorCode Destroy_InletConstantVelocity(BoundaryCondition *self);
306
307/**
308 * @brief Private data structure for the Constant Velocity Inlet handler.
309 */
310typedef struct{
311 PetscReal normal_velocity; // The desired Cartesian velocity (vx, vy, vz)
312}InletConstantData;
313
314#undef __FUNCT__
315#define __FUNCT__ "Create_InletConstantVelocity"
316/**
317 * @brief Implementation of \ref Create_InletConstantVelocity().
318 * @details Full API contract (arguments, ownership, side effects) is documented with
319 * the header declaration in `include/BC_Handlers.h`.
320 * @see Create_InletConstantVelocity()
321 */
322PetscErrorCode Create_InletConstantVelocity(BoundaryCondition *bc)
323{
324 PetscErrorCode ierr;
325 PetscFunctionBeginUser;
326
327 if (!bc) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "BoundaryCondition is NULL");
328
329 InletConstantData *data = NULL;
330 ierr = PetscMalloc1(1, &data); CHKERRQ(ierr);
331 bc->data = (void*)data;
332
333 bc->priority = BC_PRIORITY_INLET;
334 bc->Initialize = Initialize_InletConstantVelocity;
335 bc->PreStep = PreStep_InletConstantVelocity;
336 bc->Apply = Apply_InletConstantVelocity;
337 bc->PostStep = PostStep_InletConstantVelocity;
338 bc->UpdateUbcs = NULL;
339 bc->Destroy = Destroy_InletConstantVelocity;
340
341 PetscFunctionReturn(0);
342}
343
344
345#undef __FUNCT__
346#define __FUNCT__ "Initialize_InletConstantVelocity"
347/**
348 * @brief Initialize persistent state for a constant-velocity inlet boundary.
349 */
350static PetscErrorCode Initialize_InletConstantVelocity(BoundaryCondition *self, BCContext *ctx)
351{
352 PetscErrorCode ierr;
353 UserCtx* user = ctx->user;
354 BCFace face_id = ctx->face_id;
355 InletConstantData *data = (InletConstantData*)self->data;
356 PetscBool found;
357
358 PetscFunctionBeginUser;
359 LOG_ALLOW(LOCAL, LOG_DEBUG, "Initialize_InletConstantVelocity: Initializing handler for Face %d. \n", face_id);
360 data->normal_velocity = 0.0;
361
362 switch (face_id) {
363 case BC_FACE_NEG_X:
364 case BC_FACE_POS_X:
365 // For X-faces, read "vx" as normal velocity
366 ierr = GetBCParamReal(user->boundary_faces[face_id].params, "vx",
367 &data->normal_velocity, &found); CHKERRQ(ierr);
368 break;
369
370 case BC_FACE_NEG_Y:
371 case BC_FACE_POS_Y:
372 // For Y-faces, read "vy" as normal velocity
373 ierr = GetBCParamReal(user->boundary_faces[face_id].params, "vy",
374 &data->normal_velocity, &found); CHKERRQ(ierr);
375 break;
376
377 case BC_FACE_NEG_Z:
378 case BC_FACE_POS_Z:
379 // For Z-faces, read "vz" as normal velocity
380 ierr = GetBCParamReal(user->boundary_faces[face_id].params, "vz",
381 &data->normal_velocity, &found); CHKERRQ(ierr);
382 break;
383 }
384
385 LOG_ALLOW(LOCAL, LOG_INFO, " Inlet Face %d: normal velocity = %.4f\n",
386 face_id, data->normal_velocity);
387
388 // Set initial boundary state
389 ierr = Apply_InletConstantVelocity(self, ctx); CHKERRQ(ierr);
390
391 PetscFunctionReturn(0);
392}
393
394#undef __FUNCT__
395#define __FUNCT__ "PreStep_InletConstantVelocity"
396/**
397 * @brief Update constant-inlet data required before the next solver step.
398 */
399static PetscErrorCode PreStep_InletConstantVelocity(BoundaryCondition *self, BCContext *ctx,
400 PetscReal *local_inflow_contribution,
401 PetscReal *local_outflow_contribution)
402{
403 // No preparation needed for constant velocity inlet.
404 // The velocity is already stored in self->data from Initialize.
405 // Apply will set ucont, and PostStep will measure the actual flux.
406
407 (void)self;
408 (void)ctx;
409 (void)local_inflow_contribution;
410 (void)local_outflow_contribution;
411
412 PetscFunctionBeginUser;
413 PetscFunctionReturn(0);
414}
415
416#undef __FUNCT__
417#define __FUNCT__ "Apply_InletConstantVelocity"
418/**
419 * @brief Impose the configured constant velocity on inlet boundary cells.
420 */
421static PetscErrorCode Apply_InletConstantVelocity(BoundaryCondition *self, BCContext *ctx)
422{
423 PetscErrorCode ierr;
424 UserCtx* user = ctx->user;
425 BCFace face_id = ctx->face_id;
426 InletConstantData *data = (InletConstantData*)self->data;
427 PetscBool can_service;
428
429 PetscFunctionBeginUser;
430
431 DMDALocalInfo *info = &user->info;
432 Cmpnts ***ubcs, ***ucont, ***csi, ***eta, ***zet;
433 PetscReal ***nvert;
434 PetscInt IM_nodes_global, JM_nodes_global,KM_nodes_global;
435
436 IM_nodes_global = user->IM;
437 JM_nodes_global = user->JM;
438 KM_nodes_global = user->KM;
439
440 ierr = CanRankServiceFace(info,IM_nodes_global,JM_nodes_global,KM_nodes_global,face_id,&can_service); CHKERRQ(ierr);
441
442 if (!can_service) PetscFunctionReturn(0);
443
444 // Get arrays
445
446 ierr = DMDAVecGetArray(user->fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
447 ierr = DMDAVecGetArray(user->fda, user->Ucont, &ucont); CHKERRQ(ierr);
448 ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
449 ierr = DMDAVecGetArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
450 ierr = DMDAVecGetArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
451 ierr = DMDAVecGetArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
452
453 // Get SCALAR velocity (not vector!)
454 PetscReal uin_this_point = data->normal_velocity;
455
456 PetscInt xs = info->xs, xe = info->xs + info->xm;
457 PetscInt ys = info->ys, ye = info->ys + info->ym;
458 PetscInt zs = info->zs, ze = info->zs + info->zm;
459 PetscInt mx = info->mx, my = info->my, mz = info->mz;
460
461 PetscInt lxs = xs, lxe = xe, lys = ys, lye = ye, lzs = zs, lze = ze;
462 if (xs == 0) lxs = xs + 1;
463 if (xe == mx) lxe = xe - 1;
464 if (ys == 0) lys = ys + 1;
465 if (ye == my) lye = ye - 1;
466 if (zs == 0) lzs = zs + 1;
467 if (ze == mz) lze = ze - 1;
468
469 switch (face_id) {
470 case BC_FACE_NEG_X:
471 case BC_FACE_POS_X: {
472 PetscReal sign = (face_id == BC_FACE_NEG_X) ? 1.0 : -1.0;
473 PetscInt i = (face_id == BC_FACE_NEG_X) ? xs : mx - 2;
474
475 for (PetscInt k = lzs; k < lze; k++) {
476 for (PetscInt j = lys; j < lye; j++) {
477 if ((sign > 0 && nvert[k][j][i+1] > 0.1) ||
478 (sign < 0 && nvert[k][j][i] > 0.1)) continue;
479
480 PetscReal CellArea = sqrt(csi[k][j][i].x * csi[k][j][i].x +
481 csi[k][j][i].y * csi[k][j][i].y +
482 csi[k][j][i].z * csi[k][j][i].z);
483
484 ucont[k][j][i].x = sign * uin_this_point * CellArea;
485
486 ubcs[k][j][i + (sign < 0)].x = sign * uin_this_point * csi[k][j][i].x / CellArea;
487 ubcs[k][j][i + (sign < 0)].y = sign * uin_this_point * csi[k][j][i].y / CellArea;
488 ubcs[k][j][i + (sign < 0)].z = sign * uin_this_point * csi[k][j][i].z / CellArea;
489 }
490 }
491 } break;
492
493 case BC_FACE_NEG_Y:
494 case BC_FACE_POS_Y: {
495 PetscReal sign = (face_id == BC_FACE_NEG_Y) ? 1.0 : -1.0;
496 PetscInt j = (face_id == BC_FACE_NEG_Y) ? ys : my - 2;
497
498 for (PetscInt k = lzs; k < lze; k++) {
499 for (PetscInt i = lxs; i < lxe; i++) {
500 if ((sign > 0 && nvert[k][j+1][i] > 0.1) ||
501 (sign < 0 && nvert[k][j][i] > 0.1)) continue;
502
503 PetscReal CellArea = sqrt(eta[k][j][i].x * eta[k][j][i].x +
504 eta[k][j][i].y * eta[k][j][i].y +
505 eta[k][j][i].z * eta[k][j][i].z);
506
507 ucont[k][j][i].y = sign * uin_this_point * CellArea;
508
509 ubcs[k][j + (sign < 0)][i].x = sign * uin_this_point * eta[k][j][i].x / CellArea;
510 ubcs[k][j + (sign < 0)][i].y = sign * uin_this_point * eta[k][j][i].y / CellArea;
511 ubcs[k][j + (sign < 0)][i].z = sign * uin_this_point * eta[k][j][i].z / CellArea;
512 }
513 }
514 } break;
515
516 case BC_FACE_NEG_Z:
517 case BC_FACE_POS_Z: {
518 PetscReal sign = (face_id == BC_FACE_NEG_Z) ? 1.0 : -1.0;
519 PetscInt k = (face_id == BC_FACE_NEG_Z) ? zs : mz - 2;
520
521 for (PetscInt j = lys; j < lye; j++) {
522 for (PetscInt i = lxs; i < lxe; i++) {
523 if ((sign > 0 && nvert[k+1][j][i] > 0.1) ||
524 (sign < 0 && nvert[k][j][i] > 0.1)) continue;
525
526 PetscReal CellArea = sqrt(zet[k][j][i].x * zet[k][j][i].x +
527 zet[k][j][i].y * zet[k][j][i].y +
528 zet[k][j][i].z * zet[k][j][i].z);
529
530 ucont[k][j][i].z = sign * uin_this_point * CellArea;
531
532 ubcs[k + (sign < 0)][j][i].x = sign * uin_this_point * zet[k][j][i].x / CellArea;
533 ubcs[k + (sign < 0)][j][i].y = sign * uin_this_point * zet[k][j][i].y / CellArea;
534 ubcs[k + (sign < 0)][j][i].z = sign * uin_this_point * zet[k][j][i].z / CellArea;
535 }
536 }
537 } break;
538 }
539
540 // Restore arrays
541 ierr = DMDAVecRestoreArray(user->fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
542 ierr = DMDAVecRestoreArray(user->fda, user->Ucont, &ucont); CHKERRQ(ierr);
543 ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
544 ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
545 ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
546 ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
547
548 PetscFunctionReturn(0);
549}
550
551#undef __FUNCT__
552#define __FUNCT__ "PostStep_InletConstantVelocity"
553/**
554 * @brief Perform post-step bookkeeping for a constant-velocity inlet boundary.
555 */
556static PetscErrorCode PostStep_InletConstantVelocity(BoundaryCondition *self, BCContext *ctx,
557 PetscReal *local_inflow_contribution,
558 PetscReal *local_outflow_contribution)
559{
560 PetscErrorCode ierr;
561 UserCtx* user = ctx->user;
562 BCFace face_id = ctx->face_id;
563 PetscBool can_service;
564
565 (void)self;
566 (void)local_outflow_contribution;
567
568 PetscFunctionBeginUser;
569
570 DMDALocalInfo *info = &user->info;
571 Cmpnts ***ucont;
572
573 PetscInt IM_nodes_global, JM_nodes_global,KM_nodes_global;
574
575 IM_nodes_global = user->IM;
576 JM_nodes_global = user->JM;
577 KM_nodes_global = user->KM;
578
579 ierr = CanRankServiceFace(info,IM_nodes_global,JM_nodes_global,KM_nodes_global,face_id,&can_service); CHKERRQ(ierr);
580
581
582 if (!can_service) PetscFunctionReturn(0);
583
584 ierr = DMDAVecGetArrayRead(user->fda, user->Ucont, (const Cmpnts***)&ucont); CHKERRQ(ierr);
585
586 PetscReal local_flux = 0.0;
587
588 PetscInt xs = info->xs, xe = info->xs + info->xm;
589 PetscInt ys = info->ys, ye = info->ys + info->ym;
590 PetscInt zs = info->zs, ze = info->zs + info->zm;
591 PetscInt mx = info->mx, my = info->my, mz = info->mz;
592
593 PetscInt lxs = xs, lxe = xe, lys = ys, lye = ye, lzs = zs, lze = ze;
594 if (xs == 0) lxs = xs + 1;
595 if (xe == mx) lxe = xe - 1;
596 if (ys == 0) lys = ys + 1;
597 if (ye == my) lye = ye - 1;
598 if (zs == 0) lzs = zs + 1;
599 if (ze == mz) lze = ze - 1;
600
601 // Sum ucont components
602 switch (face_id) {
603 case BC_FACE_NEG_X:
604 case BC_FACE_POS_X: {
605 PetscInt i = (face_id == BC_FACE_NEG_X) ? xs : mx - 2;
606 for (PetscInt k = lzs; k < lze; k++) {
607 for (PetscInt j = lys; j < lye; j++) {
608 local_flux += ucont[k][j][i].x;
609 }
610 }
611 } break;
612
613 case BC_FACE_NEG_Y:
614 case BC_FACE_POS_Y: {
615 PetscInt j = (face_id == BC_FACE_NEG_Y) ? ys : my - 2;
616 for (PetscInt k = lzs; k < lze; k++) {
617 for (PetscInt i = lxs; i < lxe; i++) {
618 local_flux += ucont[k][j][i].y;
619 }
620 }
621 } break;
622
623 case BC_FACE_NEG_Z:
624 case BC_FACE_POS_Z: {
625 PetscInt k = (face_id == BC_FACE_NEG_Z) ? zs : mz - 2;
626 for (PetscInt j = lys; j < lye; j++) {
627 for (PetscInt i = lxs; i < lxe; i++) {
628 local_flux += ucont[k][j][i].z;
629 }
630 }
631 } break;
632 }
633
634 ierr = DMDAVecRestoreArrayRead(user->fda, user->Ucont, (const Cmpnts***)&ucont); CHKERRQ(ierr);
635
636 *local_inflow_contribution += local_flux;
637
638 LOG_ALLOW(LOCAL, LOG_DEBUG, "PostStep_InletConstantVelocity: Face %d, flux = %.6e\n",
639 face_id, local_flux);
640
641 PetscFunctionReturn(0);
642}
643
644
645#undef __FUNCT__
646#define __FUNCT__ "Destroy_InletConstantVelocity"
647/**
648 * @brief Release resources owned by a constant-velocity inlet boundary.
649 */
650static PetscErrorCode Destroy_InletConstantVelocity(BoundaryCondition *self)
651{
652 PetscFunctionBeginUser;
653 if (self && self->data) {
654 PetscFree(self->data);
655 self->data = NULL;
656 }
657 PetscFunctionReturn(0);
658}
659
660//================================================================================
661//
662// HANDLER IMPLEMENTATION: PARABOLIC VELOCITY INLET (POISEUILLE PROFILE)
663// (Corresponds to BC_HANDLER_INLET_PARABOLIC)
664//
665// This handler enforces a fully-developed parabolic (Poiseuille) velocity profile
666// on a rectangular/square inlet face. The profile shape is:
667//
668// V(cs1, cs2) = v_max * (1 - cs1_norm^2) * (1 - cs2_norm^2)
669//
670// where cs1 and cs2 are the two cross-stream index directions for the given face,
671// normalized to [-1, +1] across the interior nodes. The profile is zero at the
672// walls and peaks at v_max at the center.
673//
674// Workflow:
675// Constructor -> Allocate private data, wire function pointers.
676// Initialize -> Parse v_max from params, compute cross-stream geometry, call Apply.
677// PreStep -> No-op (profile is static in time).
678// Apply -> Set ucont and ubcs on the inlet face using the parabolic profile.
679// PostStep -> Measure actual volumetric flux through the face.
680// Destroy -> Free private data.
681//
682//================================================================================
683
684// --- FORWARD DECLARATIONS ---
685static PetscErrorCode Initialize_InletParabolicProfile(BoundaryCondition *self, BCContext *ctx);
686static PetscErrorCode PreStep_InletParabolicProfile(BoundaryCondition *self, BCContext *ctx,
687 PetscReal *in, PetscReal *out);
688static PetscErrorCode Apply_InletParabolicProfile(BoundaryCondition *self, BCContext *ctx);
689static PetscErrorCode PostStep_InletParabolicProfile(BoundaryCondition *self, BCContext *ctx,
690 PetscReal *in, PetscReal *out);
691static PetscErrorCode Destroy_InletParabolicProfile(BoundaryCondition *self);
692
693/**
694 * @brief Private data structure for the Parabolic Velocity Inlet handler.
695 *
696 * Stores the peak velocity and pre-computed cross-stream geometry needed
697 * to evaluate the parabolic profile at each boundary node.
698 */
699typedef struct {
700 PetscReal v_max; /**< Peak centerline velocity (from user params). */
701 PetscReal cs1_center; /**< Center index in cross-stream direction 1. */
702 PetscReal cs2_center; /**< Center index in cross-stream direction 2. */
703 PetscReal cs1_half; /**< Half-width (in index space) in cross-stream direction 1. */
704 PetscReal cs2_half; /**< Half-width (in index space) in cross-stream direction 2. */
705} InletParabolicData;
706
707#undef __FUNCT__
708#define __FUNCT__ "Create_InletParabolicProfile"
709/**
710 * @brief Implementation of \ref Create_InletParabolicProfile().
711 * @details Full API contract (arguments, ownership, side effects) is documented with
712 * the header declaration in `include/BC_Handlers.h`.
713 * @see Create_InletParabolicProfile()
714 */
715PetscErrorCode Create_InletParabolicProfile(BoundaryCondition *bc)
716{
717 PetscErrorCode ierr;
718 PetscFunctionBeginUser;
719
720 if (!bc) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "BoundaryCondition is NULL");
721
722 InletParabolicData *data = NULL;
723 ierr = PetscMalloc1(1, &data); CHKERRQ(ierr);
724 bc->data = (void*)data;
725
726 bc->priority = BC_PRIORITY_INLET;
727 bc->Initialize = Initialize_InletParabolicProfile;
728 bc->PreStep = PreStep_InletParabolicProfile;
729 bc->Apply = Apply_InletParabolicProfile;
730 bc->PostStep = PostStep_InletParabolicProfile;
731 bc->UpdateUbcs = NULL;
732 bc->Destroy = Destroy_InletParabolicProfile;
733
734 PetscFunctionReturn(0);
735}
736
737
738#undef __FUNCT__
739#define __FUNCT__ "Initialize_InletParabolicProfile"
740/**
741 * @brief Initialize the geometric data used to evaluate a parabolic inlet profile.
742 */
743static PetscErrorCode Initialize_InletParabolicProfile(BoundaryCondition *self, BCContext *ctx)
744{
745 PetscErrorCode ierr;
746 UserCtx* user = ctx->user;
747 BCFace face_id = ctx->face_id;
748 InletParabolicData *data = (InletParabolicData*)self->data;
749 PetscBool found;
750
751 PetscFunctionBeginUser;
752 LOG_ALLOW(LOCAL, LOG_DEBUG, "Initialize_InletParabolicProfile: Initializing handler for Face %d.\n", face_id);
753
754 // --- Parse v_max from boundary condition parameters ---
755 data->v_max = 0.0;
756 ierr = GetBCParamReal(user->boundary_faces[face_id].params, "v_max",
757 &data->v_max, &found); CHKERRQ(ierr);
758 if (!found) {
759 LOG_ALLOW(GLOBAL, LOG_WARNING, "Initialize_InletParabolicProfile: 'v_max' not found in params for face %d. Defaulting to 0.0.\n", face_id);
760 }
761
762 // --- Determine cross-stream dimensions based on face orientation ---
763 PetscReal cs1_dim, cs2_dim;
764 switch (face_id) {
765 case BC_FACE_NEG_X:
766 case BC_FACE_POS_X:
767 cs1_dim = (PetscReal)user->JM; // j-direction
768 cs2_dim = (PetscReal)user->KM; // k-direction
769 break;
770 case BC_FACE_NEG_Y:
771 case BC_FACE_POS_Y:
772 cs1_dim = (PetscReal)user->IM; // i-direction
773 cs2_dim = (PetscReal)user->KM; // k-direction
774 break;
775 case BC_FACE_NEG_Z:
776 case BC_FACE_POS_Z:
777 default:
778 cs1_dim = (PetscReal)user->IM; // i-direction
779 cs2_dim = (PetscReal)user->JM; // j-direction
780 break;
781 }
782
783 // Interior width = dim - 2 (nodes 1 through dim-2 are interior)
784 PetscReal cs1_width = cs1_dim - 2.0;
785 PetscReal cs2_width = cs2_dim - 2.0;
786
787 data->cs1_center = 1.0 + cs1_width / 2.0;
788 data->cs2_center = 1.0 + cs2_width / 2.0;
789 data->cs1_half = cs1_width / 2.0;
790 data->cs2_half = cs2_width / 2.0;
791
792 LOG_ALLOW(LOCAL, LOG_INFO, " Inlet Face %d (Parabolic): v_max = %.4f\n", face_id, data->v_max);
793 LOG_ALLOW(LOCAL, LOG_DEBUG, " Cross-stream 1: center=%.1f, half=%.1f\n", data->cs1_center, data->cs1_half);
794 LOG_ALLOW(LOCAL, LOG_DEBUG, " Cross-stream 2: center=%.1f, half=%.1f\n", data->cs2_center, data->cs2_half);
795
796 // Set initial boundary state
797 ierr = Apply_InletParabolicProfile(self, ctx); CHKERRQ(ierr);
798
799 PetscFunctionReturn(0);
800}
801
802
803#undef __FUNCT__
804#define __FUNCT__ "PreStep_InletParabolicProfile"
805/**
806 * @brief Refresh parabolic-inlet values required before the solver step.
807 */
808static PetscErrorCode PreStep_InletParabolicProfile(BoundaryCondition *self, BCContext *ctx,
809 PetscReal *local_inflow_contribution,
810 PetscReal *local_outflow_contribution)
811{
812 (void)self;
813 (void)ctx;
814 (void)local_inflow_contribution;
815 (void)local_outflow_contribution;
816
817 PetscFunctionBeginUser;
818 PetscFunctionReturn(0);
819}
820
821
822#undef __FUNCT__
823#define __FUNCT__ "Apply_InletParabolicProfile"
824/**
825 * @brief Impose the evaluated parabolic velocity profile on inlet cells.
826 */
827static PetscErrorCode Apply_InletParabolicProfile(BoundaryCondition *self, BCContext *ctx)
828{
829 PetscErrorCode ierr;
830 UserCtx* user = ctx->user;
831 BCFace face_id = ctx->face_id;
832 InletParabolicData *data = (InletParabolicData*)self->data;
833 PetscBool can_service;
834
835 PetscFunctionBeginUser;
836
837 DMDALocalInfo *info = &user->info;
838 Cmpnts ***ubcs, ***ucont, ***csi, ***eta, ***zet;
839 PetscReal ***nvert;
840 PetscInt IM_nodes_global, JM_nodes_global, KM_nodes_global;
841
842 IM_nodes_global = user->IM;
843 JM_nodes_global = user->JM;
844 KM_nodes_global = user->KM;
845
846 ierr = CanRankServiceFace(info, IM_nodes_global, JM_nodes_global, KM_nodes_global,
847 face_id, &can_service); CHKERRQ(ierr);
848
849 if (!can_service) PetscFunctionReturn(0);
850
851 // Get arrays
852 ierr = DMDAVecGetArray(user->fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
853 ierr = DMDAVecGetArray(user->fda, user->Ucont, &ucont); CHKERRQ(ierr);
854 ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
855 ierr = DMDAVecGetArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
856 ierr = DMDAVecGetArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
857 ierr = DMDAVecGetArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
858
859 PetscInt xs = info->xs, xe = info->xs + info->xm;
860 PetscInt ys = info->ys, ye = info->ys + info->ym;
861 PetscInt zs = info->zs, ze = info->zs + info->zm;
862 PetscInt mx = info->mx, my = info->my, mz = info->mz;
863
864 PetscInt lxs = xs, lxe = xe, lys = ys, lye = ye, lzs = zs, lze = ze;
865 if (xs == 0) lxs = xs + 1;
866 if (xe == mx) lxe = xe - 1;
867 if (ys == 0) lys = ys + 1;
868 if (ye == my) lye = ye - 1;
869 if (zs == 0) lzs = zs + 1;
870 if (ze == mz) lze = ze - 1;
871
872 switch (face_id) {
873 case BC_FACE_NEG_X:
874 case BC_FACE_POS_X: {
875 // X-faces: normal = i, cross-stream = (j, k)
876 PetscReal sign = (face_id == BC_FACE_NEG_X) ? 1.0 : -1.0;
877 PetscInt i = (face_id == BC_FACE_NEG_X) ? xs : mx - 2;
878
879 for (PetscInt k = lzs; k < lze; k++) {
880 for (PetscInt j = lys; j < lye; j++) {
881 if ((sign > 0 && nvert[k][j][i+1] > 0.1) ||
882 (sign < 0 && nvert[k][j][i] > 0.1)) continue;
883
884 // Evaluate parabolic profile: cs1 = j, cs2 = k
885 PetscReal cs1_norm = ((PetscReal)j - data->cs1_center) / data->cs1_half;
886 PetscReal cs2_norm = ((PetscReal)k - data->cs2_center) / data->cs2_half;
887 PetscReal profile = PetscMax(0.0, 1.0 - cs1_norm * cs1_norm)
888 * PetscMax(0.0, 1.0 - cs2_norm * cs2_norm);
889 PetscReal uin_local = data->v_max * profile;
890
891 PetscReal CellArea = sqrt(csi[k][j][i].x * csi[k][j][i].x +
892 csi[k][j][i].y * csi[k][j][i].y +
893 csi[k][j][i].z * csi[k][j][i].z);
894
895 ucont[k][j][i].x = sign * uin_local * CellArea;
896
897 ubcs[k][j][i + (sign < 0)].x = sign * uin_local * csi[k][j][i].x / CellArea;
898 ubcs[k][j][i + (sign < 0)].y = sign * uin_local * csi[k][j][i].y / CellArea;
899 ubcs[k][j][i + (sign < 0)].z = sign * uin_local * csi[k][j][i].z / CellArea;
900 }
901 }
902 } break;
903
904 case BC_FACE_NEG_Y:
905 case BC_FACE_POS_Y: {
906 // Y-faces: normal = j, cross-stream = (i, k)
907 PetscReal sign = (face_id == BC_FACE_NEG_Y) ? 1.0 : -1.0;
908 PetscInt j = (face_id == BC_FACE_NEG_Y) ? ys : my - 2;
909
910 for (PetscInt k = lzs; k < lze; k++) {
911 for (PetscInt i = lxs; i < lxe; i++) {
912 if ((sign > 0 && nvert[k][j+1][i] > 0.1) ||
913 (sign < 0 && nvert[k][j][i] > 0.1)) continue;
914
915 // Evaluate parabolic profile: cs1 = i, cs2 = k
916 PetscReal cs1_norm = ((PetscReal)i - data->cs1_center) / data->cs1_half;
917 PetscReal cs2_norm = ((PetscReal)k - data->cs2_center) / data->cs2_half;
918 PetscReal profile = PetscMax(0.0, 1.0 - cs1_norm * cs1_norm)
919 * PetscMax(0.0, 1.0 - cs2_norm * cs2_norm);
920 PetscReal uin_local = data->v_max * profile;
921
922 PetscReal CellArea = sqrt(eta[k][j][i].x * eta[k][j][i].x +
923 eta[k][j][i].y * eta[k][j][i].y +
924 eta[k][j][i].z * eta[k][j][i].z);
925
926 ucont[k][j][i].y = sign * uin_local * CellArea;
927
928 ubcs[k][j + (sign < 0)][i].x = sign * uin_local * eta[k][j][i].x / CellArea;
929 ubcs[k][j + (sign < 0)][i].y = sign * uin_local * eta[k][j][i].y / CellArea;
930 ubcs[k][j + (sign < 0)][i].z = sign * uin_local * eta[k][j][i].z / CellArea;
931 }
932 }
933 } break;
934
935 case BC_FACE_NEG_Z:
936 case BC_FACE_POS_Z: {
937 // Z-faces: normal = k, cross-stream = (i, j)
938 PetscReal sign = (face_id == BC_FACE_NEG_Z) ? 1.0 : -1.0;
939 PetscInt k = (face_id == BC_FACE_NEG_Z) ? zs : mz - 2;
940
941 for (PetscInt j = lys; j < lye; j++) {
942 for (PetscInt i = lxs; i < lxe; i++) {
943 if ((sign > 0 && nvert[k+1][j][i] > 0.1) ||
944 (sign < 0 && nvert[k][j][i] > 0.1)) continue;
945
946 // Evaluate parabolic profile: cs1 = i, cs2 = j
947 PetscReal cs1_norm = ((PetscReal)i - data->cs1_center) / data->cs1_half;
948 PetscReal cs2_norm = ((PetscReal)j - data->cs2_center) / data->cs2_half;
949 PetscReal profile = PetscMax(0.0, 1.0 - cs1_norm * cs1_norm)
950 * PetscMax(0.0, 1.0 - cs2_norm * cs2_norm);
951 PetscReal uin_local = data->v_max * profile;
952
953 PetscReal CellArea = sqrt(zet[k][j][i].x * zet[k][j][i].x +
954 zet[k][j][i].y * zet[k][j][i].y +
955 zet[k][j][i].z * zet[k][j][i].z);
956
957 ucont[k][j][i].z = sign * uin_local * CellArea;
958
959 ubcs[k + (sign < 0)][j][i].x = sign * uin_local * zet[k][j][i].x / CellArea;
960 ubcs[k + (sign < 0)][j][i].y = sign * uin_local * zet[k][j][i].y / CellArea;
961 ubcs[k + (sign < 0)][j][i].z = sign * uin_local * zet[k][j][i].z / CellArea;
962 }
963 }
964 } break;
965 }
966
967 // Restore arrays
968 ierr = DMDAVecRestoreArray(user->fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
969 ierr = DMDAVecRestoreArray(user->fda, user->Ucont, &ucont); CHKERRQ(ierr);
970 ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
971 ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
972 ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
973 ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
974
975 PetscFunctionReturn(0);
976}
977
978
979#undef __FUNCT__
980#define __FUNCT__ "PostStep_InletParabolicProfile"
981/**
982 * @brief Perform post-step bookkeeping for a parabolic inlet boundary.
983 */
984static PetscErrorCode PostStep_InletParabolicProfile(BoundaryCondition *self, BCContext *ctx,
985 PetscReal *local_inflow_contribution,
986 PetscReal *local_outflow_contribution)
987{
988 PetscErrorCode ierr;
989 UserCtx* user = ctx->user;
990 BCFace face_id = ctx->face_id;
991 PetscBool can_service;
992
993 (void)self;
994 (void)local_outflow_contribution;
995
996 PetscFunctionBeginUser;
997
998 DMDALocalInfo *info = &user->info;
999 Cmpnts ***ucont;
1000
1001 PetscInt IM_nodes_global, JM_nodes_global, KM_nodes_global;
1002
1003 IM_nodes_global = user->IM;
1004 JM_nodes_global = user->JM;
1005 KM_nodes_global = user->KM;
1006
1007 ierr = CanRankServiceFace(info, IM_nodes_global, JM_nodes_global, KM_nodes_global,
1008 face_id, &can_service); CHKERRQ(ierr);
1009
1010 if (!can_service) PetscFunctionReturn(0);
1011
1012 ierr = DMDAVecGetArrayRead(user->fda, user->Ucont, (const Cmpnts***)&ucont); CHKERRQ(ierr);
1013
1014 PetscReal local_flux = 0.0;
1015
1016 PetscInt xs = info->xs, xe = info->xs + info->xm;
1017 PetscInt ys = info->ys, ye = info->ys + info->ym;
1018 PetscInt zs = info->zs, ze = info->zs + info->zm;
1019 PetscInt mx = info->mx, my = info->my, mz = info->mz;
1020
1021 PetscInt lxs = xs, lxe = xe, lys = ys, lye = ye, lzs = zs, lze = ze;
1022 if (xs == 0) lxs = xs + 1;
1023 if (xe == mx) lxe = xe - 1;
1024 if (ys == 0) lys = ys + 1;
1025 if (ye == my) lye = ye - 1;
1026 if (zs == 0) lzs = zs + 1;
1027 if (ze == mz) lze = ze - 1;
1028
1029 switch (face_id) {
1030 case BC_FACE_NEG_X:
1031 case BC_FACE_POS_X: {
1032 PetscInt i = (face_id == BC_FACE_NEG_X) ? xs : mx - 2;
1033 for (PetscInt k = lzs; k < lze; k++) {
1034 for (PetscInt j = lys; j < lye; j++) {
1035 local_flux += ucont[k][j][i].x;
1036 }
1037 }
1038 } break;
1039
1040 case BC_FACE_NEG_Y:
1041 case BC_FACE_POS_Y: {
1042 PetscInt j = (face_id == BC_FACE_NEG_Y) ? ys : my - 2;
1043 for (PetscInt k = lzs; k < lze; k++) {
1044 for (PetscInt i = lxs; i < lxe; i++) {
1045 local_flux += ucont[k][j][i].y;
1046 }
1047 }
1048 } break;
1049
1050 case BC_FACE_NEG_Z:
1051 case BC_FACE_POS_Z: {
1052 PetscInt k = (face_id == BC_FACE_NEG_Z) ? zs : mz - 2;
1053 for (PetscInt j = lys; j < lye; j++) {
1054 for (PetscInt i = lxs; i < lxe; i++) {
1055 local_flux += ucont[k][j][i].z;
1056 }
1057 }
1058 } break;
1059 }
1060
1061 ierr = DMDAVecRestoreArrayRead(user->fda, user->Ucont, (const Cmpnts***)&ucont); CHKERRQ(ierr);
1062
1063 *local_inflow_contribution += local_flux;
1064
1065 LOG_ALLOW(LOCAL, LOG_DEBUG, "PostStep_InletParabolicProfile: Face %d, flux = %.6e\n",
1066 face_id, local_flux);
1067
1068 PetscFunctionReturn(0);
1069}
1070
1071
1072#undef __FUNCT__
1073#define __FUNCT__ "Destroy_InletParabolicProfile"
1074/**
1075 * @brief Release resources owned by a parabolic inlet boundary.
1076 */
1077static PetscErrorCode Destroy_InletParabolicProfile(BoundaryCondition *self)
1078{
1079 PetscFunctionBeginUser;
1080 if (self && self->data) {
1081 PetscFree(self->data);
1082 self->data = NULL;
1083 }
1084 PetscFunctionReturn(0);
1085}
1086
1087//================================================================================
1088//
1089// HANDLER IMPLEMENTATION: PRESCRIBED INLET PROFILE FROM FILE
1090// (Corresponds to BC_HANDLER_INLET_PROFILE_FROM_FILE)
1091//
1092// This handler reads positive scalar normal speeds from a canonical PICSLICE file.
1093// It then applies those speeds through the same face sign and metric conversion
1094// used by the constant and parabolic inlet handlers.
1095//
1096//================================================================================
1097
1098static PetscErrorCode Initialize_InletProfileFromFile(BoundaryCondition *self, BCContext *ctx);
1099static PetscErrorCode PreStep_InletProfileFromFile(BoundaryCondition *self, BCContext *ctx,
1100 PetscReal *in, PetscReal *out);
1101static PetscErrorCode Apply_InletProfileFromFile(BoundaryCondition *self, BCContext *ctx);
1102static PetscErrorCode PostStep_InletProfileFromFile(BoundaryCondition *self, BCContext *ctx,
1103 PetscReal *in, PetscReal *out);
1104static PetscErrorCode Destroy_InletProfileFromFile(BoundaryCondition *self);
1105
1106typedef struct {
1107 PetscInt n1;
1108 PetscInt n2;
1109 PetscReal *profile;
1110 PetscReal min_speed;
1111 PetscReal max_speed;
1112 char *source_file;
1113} InletProfileFileData;
1114
1115/**
1116 * @brief Looks up a string-valued boundary-condition parameter in a BC_Param list.
1117 *
1118 * @param params Head of the boundary-condition parameter linked list.
1119 * @param key Case-insensitive key to search for.
1120 * @param[out] value_out Borrowed pointer to the matching value string, or NULL if absent.
1121 * @param[out] found PETSC_TRUE when the key is present, PETSC_FALSE otherwise.
1122 * @return PetscErrorCode 0 on success.
1123 */
1124static PetscErrorCode GetBCParamStringLocal(BC_Param *params, const char *key,
1125 const char **value_out, PetscBool *found)
1126{
1127 PetscFunctionBeginUser;
1128 *found = PETSC_FALSE;
1129 *value_out = NULL;
1130 for (BC_Param *param = params; param; param = param->next) {
1131 if (strcasecmp(param->key, key) == 0) {
1132 *value_out = param->value;
1133 *found = PETSC_TRUE;
1134 PetscFunctionReturn(0);
1135 }
1136 }
1137 PetscFunctionReturn(0);
1138}
1139
1140/**
1141 * @brief Computes the expected PICSLICE dimensions for an inlet face.
1142 *
1143 * @param user User context containing global grid node counts.
1144 * @param face_id Boundary face whose tangential profile dimensions are requested.
1145 * @param[out] n1 First PICSLICE dimension in handler storage order.
1146 * @param[out] n2 Second PICSLICE dimension in handler storage order.
1147 * @return PetscErrorCode 0 on success, or a PETSc error for unsupported faces or invalid dimensions.
1148 */
1149static PetscErrorCode GetProfileFileExpectedDims(UserCtx *user, BCFace face_id,
1150 PetscInt *n1, PetscInt *n2)
1151{
1152 PetscFunctionBeginUser;
1153 switch (face_id) {
1154 case BC_FACE_NEG_X:
1155 case BC_FACE_POS_X:
1156 *n1 = user->KM - 1;
1157 *n2 = user->JM - 1;
1158 break;
1159 case BC_FACE_NEG_Y:
1160 case BC_FACE_POS_Y:
1161 *n1 = user->KM - 1;
1162 *n2 = user->IM - 1;
1163 break;
1164 case BC_FACE_NEG_Z:
1165 case BC_FACE_POS_Z:
1166 *n1 = user->JM - 1;
1167 *n2 = user->IM - 1;
1168 break;
1169 default:
1170 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE,
1171 "Unsupported face id %d for inlet profile dimensions.", face_id);
1172 }
1173 if (*n1 <= 0 || *n2 <= 0) {
1174 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
1175 "Invalid inlet profile dimensions (%d, %d) for grid (%d, %d, %d).",
1176 *n1, *n2, user->IM, user->JM, user->KM);
1177 }
1178 PetscFunctionReturn(0);
1179}
1180
1181/**
1182 * @brief Reads and validates a static scalar inlet profile from a canonical PICSLICE file.
1183 *
1184 * @details The file must contain magic token `PICSLICE`, frame count 1, the expected
1185 * two-dimensional face shape, and exactly one finite nonnegative scalar speed
1186 * value per interior face slot. Values are stored row-major in `data->profile`.
1187 *
1188 * @param source_file Path to the PICSLICE profile file.
1189 * @param expected_n1 Expected first slice dimension.
1190 * @param expected_n2 Expected second slice dimension.
1191 * @param[in,out] data Handler-private storage that receives dimensions, profile values, and min/max speeds.
1192 * @return PetscErrorCode 0 on success, or a PETSc file/validation error on malformed input.
1193 */
1194static PetscErrorCode ReadPicSliceProfile(const char *source_file, PetscInt expected_n1,
1195 PetscInt expected_n2, InletProfileFileData *data)
1196{
1197 PetscErrorCode ierr;
1198 FILE *fd = NULL;
1199 char magic[32] = {0};
1200 PetscInt frame_count = 0, n1 = 0, n2 = 0;
1201
1202 PetscFunctionBeginUser;
1203 fd = fopen(source_file, "r");
1204 if (!fd) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN,
1205 "Cannot open PICSLICE inlet profile file: %s", source_file);
1206
1207 if (fscanf(fd, "%31s", magic) != 1 || strcmp(magic, "PICSLICE") != 0) {
1208 fclose(fd);
1209 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_READ,
1210 "PICSLICE inlet profile file %s must begin with PICSLICE header.", source_file);
1211 }
1212 if (fscanf(fd, "%d", &frame_count) != 1) {
1213 fclose(fd);
1214 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_READ,
1215 "PICSLICE inlet profile file %s missing frame count.", source_file);
1216 }
1217 if (frame_count != 1) {
1218 fclose(fd);
1219 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED,
1220 "PICSLICE inlet profile file %s has %d frames; static handler requires 1.",
1221 source_file, frame_count);
1222 }
1223 if (fscanf(fd, "%d %d", &n1, &n2) != 2) {
1224 fclose(fd);
1225 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_READ,
1226 "PICSLICE inlet profile file %s missing slice dimensions.", source_file);
1227 }
1228 if (n1 != expected_n1 || n2 != expected_n2) {
1229 fclose(fd);
1230 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED,
1231 "PICSLICE inlet profile dimensions mismatch for %s: expected (%d, %d), found (%d, %d).",
1232 source_file, expected_n1, expected_n2, n1, n2);
1233 }
1234
1235 data->n1 = n1;
1236 data->n2 = n2;
1237 ierr = PetscMalloc1(n1 * n2, &data->profile); CHKERRQ(ierr);
1238 data->min_speed = PETSC_MAX_REAL;
1239 data->max_speed = -PETSC_MAX_REAL;
1240
1241 for (PetscInt idx = 0; idx < n1 * n2; idx++) {
1242 PetscReal value = 0.0;
1243 if (fscanf(fd, "%le", &value) != 1) {
1244 fclose(fd);
1245 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_READ,
1246 "PICSLICE inlet profile file %s ended early: expected %d values.",
1247 source_file, n1 * n2);
1248 }
1249 if (PetscIsInfOrNanReal(value) || value < 0.0) {
1250 fclose(fd);
1251 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED,
1252 "PICSLICE inlet profile file %s contains invalid speed %.6e at flat index %d.",
1253 source_file, (double)value, idx);
1254 }
1255 data->profile[idx] = value;
1256 data->min_speed = PetscMin(data->min_speed, value);
1257 data->max_speed = PetscMax(data->max_speed, value);
1258 }
1259
1260 char extra[64];
1261 if (fscanf(fd, "%63s", extra) == 1) {
1262 fclose(fd);
1263 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED,
1264 "PICSLICE inlet profile file %s has extra token after %d values: %s",
1265 source_file, n1 * n2, extra);
1266 }
1267 fclose(fd);
1268 PetscFunctionReturn(0);
1269}
1270
1271/**
1272 * @brief Returns one scalar speed from the flattened PICSLICE profile.
1273 *
1274 * @param data Handler-private profile storage.
1275 * @param a First profile index in face-specific storage order.
1276 * @param b Second profile index in face-specific storage order.
1277 * @return Scalar normal speed at `(a, b)`.
1278 */
1279static inline PetscReal ProfileSpeedAt(const InletProfileFileData *data, PetscInt a, PetscInt b)
1280{
1281 return data->profile[a * data->n2 + b];
1282}
1283
1284#undef __FUNCT__
1285#define __FUNCT__ "Create_InletProfileFromFile"
1286/**
1287 * @brief Implementation of \ref Create_InletProfileFromFile().
1288 * @details Full API contract (arguments, ownership, side effects) is documented with
1289 * the header declaration in `include/BC_Handlers.h`.
1290 * @see Create_InletProfileFromFile()
1291 */
1292PetscErrorCode Create_InletProfileFromFile(BoundaryCondition *bc)
1293{
1294 PetscErrorCode ierr;
1295 PetscFunctionBeginUser;
1296
1297 if (!bc) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "BoundaryCondition is NULL");
1298
1299 InletProfileFileData *data = NULL;
1300 ierr = PetscMalloc1(1, &data); CHKERRQ(ierr);
1301 data->n1 = 0;
1302 data->n2 = 0;
1303 data->profile = NULL;
1304 data->min_speed = 0.0;
1305 data->max_speed = 0.0;
1306 data->source_file = NULL;
1307 bc->data = (void*)data;
1308
1309 bc->priority = BC_PRIORITY_INLET;
1310 bc->Initialize = Initialize_InletProfileFromFile;
1311 bc->PreStep = PreStep_InletProfileFromFile;
1312 bc->Apply = Apply_InletProfileFromFile;
1313 bc->PostStep = PostStep_InletProfileFromFile;
1314 bc->UpdateUbcs = NULL;
1315 bc->Destroy = Destroy_InletProfileFromFile;
1316
1317 PetscFunctionReturn(0);
1318}
1319
1320#undef __FUNCT__
1321#define __FUNCT__ "Initialize_InletProfileFromFile"
1322/**
1323 * @brief Initializes a file-prescribed inlet profile handler for one boundary face.
1324 *
1325 * @details Reads the `source_file` BC parameter, validates the target face dimensions,
1326 * loads the PICSLICE scalar speed profile, records summary statistics for logging,
1327 * and applies the profile once to initialize boundary state.
1328 *
1329 * @param self BoundaryCondition object configured by Create_InletProfileFromFile().
1330 * @param ctx Runtime boundary context containing the UserCtx and face id.
1331 * @return PetscErrorCode 0 on success, or a PETSc error for missing parameters or malformed files.
1332 */
1333static PetscErrorCode Initialize_InletProfileFromFile(BoundaryCondition *self, BCContext *ctx)
1334{
1335 PetscErrorCode ierr;
1336 UserCtx *user = ctx->user;
1337 BCFace face_id = ctx->face_id;
1338 InletProfileFileData *data = (InletProfileFileData*)self->data;
1339 PetscBool found = PETSC_FALSE;
1340 const char *source_file = NULL;
1341 PetscInt expected_n1 = 0, expected_n2 = 0;
1342
1343 PetscFunctionBeginUser;
1344 ierr = GetBCParamStringLocal(user->boundary_faces[face_id].params, "source_file",
1345 &source_file, &found); CHKERRQ(ierr);
1346 if (!found || !source_file || source_file[0] == '\0') {
1347 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
1348 "InletProfileFromFile requires source_file parameter for face %d.", face_id);
1349 }
1350
1351 ierr = PetscStrallocpy(source_file, &data->source_file); CHKERRQ(ierr);
1352 ierr = GetProfileFileExpectedDims(user, face_id, &expected_n1, &expected_n2); CHKERRQ(ierr);
1353 ierr = ReadPicSliceProfile(source_file, expected_n1, expected_n2, data); CHKERRQ(ierr);
1354
1355 LOG_ALLOW(LOCAL, LOG_INFO,
1356 " Inlet Face %d (Prescribed Flow): source=%s dims=(%d,%d) speed[min,max]=[%.6e, %.6e]\n",
1357 face_id, data->source_file, data->n1, data->n2,
1358 (double)data->min_speed, (double)data->max_speed);
1359
1360 ierr = Apply_InletProfileFromFile(self, ctx); CHKERRQ(ierr);
1361 PetscFunctionReturn(0);
1362}
1363
1364#undef __FUNCT__
1365#define __FUNCT__ "PreStep_InletProfileFromFile"
1366/**
1367 * @brief Pre-step hook for the static file-prescribed inlet profile handler.
1368 *
1369 * @details Static profiles require no per-step preparation. The hook is implemented
1370 * so future time-varying profile support can reuse the same handler lifecycle.
1371 *
1372 * @param self BoundaryCondition object for this inlet handler.
1373 * @param ctx Runtime boundary context.
1374 * @param local_inflow_contribution Inflow accumulator, intentionally unchanged.
1375 * @param local_outflow_contribution Outflow accumulator, intentionally unchanged.
1376 * @return PetscErrorCode 0 on success.
1377 */
1378static PetscErrorCode PreStep_InletProfileFromFile(BoundaryCondition *self, BCContext *ctx,
1379 PetscReal *local_inflow_contribution,
1380 PetscReal *local_outflow_contribution)
1381{
1382 (void)self;
1383 (void)ctx;
1384 (void)local_inflow_contribution;
1385 (void)local_outflow_contribution;
1386 PetscFunctionBeginUser;
1387 PetscFunctionReturn(0);
1388}
1389
1390#undef __FUNCT__
1391#define __FUNCT__ "Apply_InletProfileFromFile"
1392/**
1393 * @brief Applies the loaded PICSLICE scalar profile to Ucont and Ubcs on an inlet face.
1394 *
1395 * @details Each stored scalar is treated as a positive normal speed magnitude. The routine
1396 * uses the same negative/positive face sign convention, immersed-boundary skip
1397 * checks, metric vectors, and CellArea conversion used by the constant and
1398 * parabolic inlet handlers.
1399 *
1400 * @param self BoundaryCondition object with InletProfileFileData storage.
1401 * @param ctx Runtime boundary context containing arrays and face id.
1402 * @return PetscErrorCode 0 on success, or a PETSc error from DMDA array access.
1403 */
1404static PetscErrorCode Apply_InletProfileFromFile(BoundaryCondition *self, BCContext *ctx)
1405{
1406 PetscErrorCode ierr;
1407 UserCtx *user = ctx->user;
1408 BCFace face_id = ctx->face_id;
1409 InletProfileFileData *data = (InletProfileFileData*)self->data;
1410 PetscBool can_service;
1411
1412 PetscFunctionBeginUser;
1413 DMDALocalInfo *info = &user->info;
1414 Cmpnts ***ubcs, ***ucont, ***csi, ***eta, ***zet;
1415 PetscReal ***nvert;
1416
1417 ierr = CanRankServiceFace(info, user->IM, user->JM, user->KM, face_id, &can_service); CHKERRQ(ierr);
1418 if (!can_service) PetscFunctionReturn(0);
1419
1420 ierr = DMDAVecGetArray(user->fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
1421 ierr = DMDAVecGetArray(user->fda, user->Ucont, &ucont); CHKERRQ(ierr);
1422 ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
1423 ierr = DMDAVecGetArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
1424 ierr = DMDAVecGetArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
1425 ierr = DMDAVecGetArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
1426
1427 PetscInt xs = info->xs, xe = info->xs + info->xm;
1428 PetscInt ys = info->ys, ye = info->ys + info->ym;
1429 PetscInt zs = info->zs, ze = info->zs + info->zm;
1430 PetscInt mx = info->mx, my = info->my, mz = info->mz;
1431
1432 PetscInt lxs = xs, lxe = xe, lys = ys, lye = ye, lzs = zs, lze = ze;
1433 if (xs == 0) lxs = xs + 1;
1434 if (xe == mx) lxe = xe - 1;
1435 if (ys == 0) lys = ys + 1;
1436 if (ye == my) lye = ye - 1;
1437 if (zs == 0) lzs = zs + 1;
1438 if (ze == mz) lze = ze - 1;
1439
1440 switch (face_id) {
1441 case BC_FACE_NEG_X:
1442 case BC_FACE_POS_X: {
1443 PetscReal sign = (face_id == BC_FACE_NEG_X) ? 1.0 : -1.0;
1444 PetscInt i = (face_id == BC_FACE_NEG_X) ? xs : mx - 2;
1445 for (PetscInt k = lzs; k < lze; k++) {
1446 for (PetscInt j = lys; j < lye; j++) {
1447 if ((sign > 0 && nvert[k][j][i+1] > 0.1) ||
1448 (sign < 0 && nvert[k][j][i] > 0.1)) continue;
1449 PetscReal uin_local = ProfileSpeedAt(data, k - 1, j - 1);
1450 PetscReal CellArea = sqrt(csi[k][j][i].x * csi[k][j][i].x +
1451 csi[k][j][i].y * csi[k][j][i].y +
1452 csi[k][j][i].z * csi[k][j][i].z);
1453 ucont[k][j][i].x = sign * uin_local * CellArea;
1454 ubcs[k][j][i + (sign < 0)].x = sign * uin_local * csi[k][j][i].x / CellArea;
1455 ubcs[k][j][i + (sign < 0)].y = sign * uin_local * csi[k][j][i].y / CellArea;
1456 ubcs[k][j][i + (sign < 0)].z = sign * uin_local * csi[k][j][i].z / CellArea;
1457 }
1458 }
1459 } break;
1460 case BC_FACE_NEG_Y:
1461 case BC_FACE_POS_Y: {
1462 PetscReal sign = (face_id == BC_FACE_NEG_Y) ? 1.0 : -1.0;
1463 PetscInt j = (face_id == BC_FACE_NEG_Y) ? ys : my - 2;
1464 for (PetscInt k = lzs; k < lze; k++) {
1465 for (PetscInt i = lxs; i < lxe; i++) {
1466 if ((sign > 0 && nvert[k][j+1][i] > 0.1) ||
1467 (sign < 0 && nvert[k][j][i] > 0.1)) continue;
1468 PetscReal uin_local = ProfileSpeedAt(data, k - 1, i - 1);
1469 PetscReal CellArea = sqrt(eta[k][j][i].x * eta[k][j][i].x +
1470 eta[k][j][i].y * eta[k][j][i].y +
1471 eta[k][j][i].z * eta[k][j][i].z);
1472 ucont[k][j][i].y = sign * uin_local * CellArea;
1473 ubcs[k][j + (sign < 0)][i].x = sign * uin_local * eta[k][j][i].x / CellArea;
1474 ubcs[k][j + (sign < 0)][i].y = sign * uin_local * eta[k][j][i].y / CellArea;
1475 ubcs[k][j + (sign < 0)][i].z = sign * uin_local * eta[k][j][i].z / CellArea;
1476 }
1477 }
1478 } break;
1479 case BC_FACE_NEG_Z:
1480 case BC_FACE_POS_Z: {
1481 PetscReal sign = (face_id == BC_FACE_NEG_Z) ? 1.0 : -1.0;
1482 PetscInt k = (face_id == BC_FACE_NEG_Z) ? zs : mz - 2;
1483 for (PetscInt j = lys; j < lye; j++) {
1484 for (PetscInt i = lxs; i < lxe; i++) {
1485 if ((sign > 0 && nvert[k+1][j][i] > 0.1) ||
1486 (sign < 0 && nvert[k][j][i] > 0.1)) continue;
1487 PetscReal uin_local = ProfileSpeedAt(data, j - 1, i - 1);
1488 PetscReal CellArea = sqrt(zet[k][j][i].x * zet[k][j][i].x +
1489 zet[k][j][i].y * zet[k][j][i].y +
1490 zet[k][j][i].z * zet[k][j][i].z);
1491 ucont[k][j][i].z = sign * uin_local * CellArea;
1492 ubcs[k + (sign < 0)][j][i].x = sign * uin_local * zet[k][j][i].x / CellArea;
1493 ubcs[k + (sign < 0)][j][i].y = sign * uin_local * zet[k][j][i].y / CellArea;
1494 ubcs[k + (sign < 0)][j][i].z = sign * uin_local * zet[k][j][i].z / CellArea;
1495 }
1496 }
1497 } break;
1498 }
1499
1500 ierr = DMDAVecRestoreArray(user->fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
1501 ierr = DMDAVecRestoreArray(user->fda, user->Ucont, &ucont); CHKERRQ(ierr);
1502 ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
1503 ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
1504 ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
1505 ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
1506
1507 PetscFunctionReturn(0);
1508}
1509
1510#undef __FUNCT__
1511#define __FUNCT__ "PostStep_InletProfileFromFile"
1512/**
1513 * @brief Accumulates the applied inlet flux for a file-prescribed profile.
1514 *
1515 * @details Sums the face-normal Ucont component over the same interior face slots
1516 * populated by Apply_InletProfileFromFile().
1517 *
1518 * @param self BoundaryCondition object for this inlet handler.
1519 * @param ctx Runtime boundary context containing the UserCtx and face id.
1520 * @param local_inflow_contribution Accumulator incremented by the measured inlet flux.
1521 * @param local_outflow_contribution Outflow accumulator, intentionally unchanged.
1522 * @return PetscErrorCode 0 on success, or a PETSc error from DMDA array access.
1523 */
1524static PetscErrorCode PostStep_InletProfileFromFile(BoundaryCondition *self, BCContext *ctx,
1525 PetscReal *local_inflow_contribution,
1526 PetscReal *local_outflow_contribution)
1527{
1528 PetscErrorCode ierr;
1529 UserCtx *user = ctx->user;
1530 BCFace face_id = ctx->face_id;
1531 PetscBool can_service;
1532
1533 (void)self;
1534 (void)local_outflow_contribution;
1535
1536 PetscFunctionBeginUser;
1537 DMDALocalInfo *info = &user->info;
1538 Cmpnts ***ucont;
1539
1540 ierr = CanRankServiceFace(info, user->IM, user->JM, user->KM, face_id, &can_service); CHKERRQ(ierr);
1541 if (!can_service) PetscFunctionReturn(0);
1542
1543 ierr = DMDAVecGetArrayRead(user->fda, user->Ucont, (const Cmpnts***)&ucont); CHKERRQ(ierr);
1544 PetscReal local_flux = 0.0;
1545
1546 PetscInt xs = info->xs, xe = info->xs + info->xm;
1547 PetscInt ys = info->ys, ye = info->ys + info->ym;
1548 PetscInt zs = info->zs, ze = info->zs + info->zm;
1549 PetscInt mx = info->mx, my = info->my, mz = info->mz;
1550
1551 PetscInt lxs = xs, lxe = xe, lys = ys, lye = ye, lzs = zs, lze = ze;
1552 if (xs == 0) lxs = xs + 1;
1553 if (xe == mx) lxe = xe - 1;
1554 if (ys == 0) lys = ys + 1;
1555 if (ye == my) lye = ye - 1;
1556 if (zs == 0) lzs = zs + 1;
1557 if (ze == mz) lze = ze - 1;
1558
1559 switch (face_id) {
1560 case BC_FACE_NEG_X:
1561 case BC_FACE_POS_X: {
1562 PetscInt i = (face_id == BC_FACE_NEG_X) ? xs : mx - 2;
1563 for (PetscInt k = lzs; k < lze; k++)
1564 for (PetscInt j = lys; j < lye; j++)
1565 local_flux += ucont[k][j][i].x;
1566 } break;
1567 case BC_FACE_NEG_Y:
1568 case BC_FACE_POS_Y: {
1569 PetscInt j = (face_id == BC_FACE_NEG_Y) ? ys : my - 2;
1570 for (PetscInt k = lzs; k < lze; k++)
1571 for (PetscInt i = lxs; i < lxe; i++)
1572 local_flux += ucont[k][j][i].y;
1573 } break;
1574 case BC_FACE_NEG_Z:
1575 case BC_FACE_POS_Z: {
1576 PetscInt k = (face_id == BC_FACE_NEG_Z) ? zs : mz - 2;
1577 for (PetscInt j = lys; j < lye; j++)
1578 for (PetscInt i = lxs; i < lxe; i++)
1579 local_flux += ucont[k][j][i].z;
1580 } break;
1581 }
1582
1583 ierr = DMDAVecRestoreArrayRead(user->fda, user->Ucont, (const Cmpnts***)&ucont); CHKERRQ(ierr);
1584 *local_inflow_contribution += local_flux;
1585
1586 LOG_ALLOW(LOCAL, LOG_DEBUG, "PostStep_InletProfileFromFile: Face %d, flux = %.6e\n",
1587 face_id, local_flux);
1588
1589 PetscFunctionReturn(0);
1590}
1591
1592#undef __FUNCT__
1593#define __FUNCT__ "Destroy_InletProfileFromFile"
1594/**
1595 * @brief Releases private storage owned by a file-prescribed inlet profile handler.
1596 *
1597 * @param self BoundaryCondition object whose `data` field stores InletProfileFileData.
1598 * @return PetscErrorCode 0 on success.
1599 */
1600static PetscErrorCode Destroy_InletProfileFromFile(BoundaryCondition *self)
1601{
1602 PetscFunctionBeginUser;
1603 if (self && self->data) {
1604 InletProfileFileData *data = (InletProfileFileData*)self->data;
1605 PetscFree(data->profile);
1606 PetscFree(data->source_file);
1607 PetscFree(self->data);
1608 self->data = NULL;
1609 }
1610 PetscFunctionReturn(0);
1611}
1612
1613//================================================================================
1614//
1615// HANDLER IMPLEMENTATION: OUTLET WITH MASS CONSERVATION
1616// (Corresponds to BC_HANDLER_OUTLET_CONSERVATION)
1617//
1618// This handler ensures that the total flux leaving through outlet boundaries
1619// balances the total flux entering through inlet and far-field boundaries.
1620//
1621// Workflow:
1622// 1. PreStep: Measures the *uncorrected* flux based on interior velocities.
1623// 2. Apply: Calculates a global correction factor based on the flux imbalance
1624// and applies it to the contravariant velocity (ucont) on the outlet face.
1625// 3. PostStep: Measures the *corrected* flux for verification and logging.
1626//
1627//================================================================================
1628
1629// --- 1. FORWARD DECLARATIONS ---
1630static PetscErrorCode PreStep_OutletConservation(BoundaryCondition *self, BCContext *ctx,
1631 PetscReal *local_inflow_contribution, PetscReal *local_outflow_contribution);
1632static PetscErrorCode Apply_OutletConservation(BoundaryCondition *self, BCContext *ctx);
1633static PetscErrorCode PostStep_OutletConservation(BoundaryCondition *self, BCContext *ctx,
1634 PetscReal *in, PetscReal *out);
1635
1636#undef __FUNCT__
1637#define __FUNCT__ "Create_OutletConservation"
1638/**
1639 * @brief Implementation of \ref Create_OutletConservation().
1640 * @details Full API contract (arguments, ownership, side effects) is documented with
1641 * the header declaration in `include/BC_Handlers.h`.
1642 * @see Create_OutletConservation()
1643 */
1644PetscErrorCode Create_OutletConservation(BoundaryCondition *bc)
1645{
1646 PetscFunctionBeginUser;
1647
1648 if (!bc) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Input BoundaryCondition is NULL");
1649
1650 // This handler has the highest priority to ensure it runs after
1651 // all inflow fluxes have been calculated.
1652 bc->priority = BC_PRIORITY_OUTLET;
1653
1654 // Assign function pointers
1655 bc->Initialize = NULL; // No initialization needed
1656 bc->PreStep = PreStep_OutletConservation;
1657 bc->Apply = Apply_OutletConservation;
1658 bc->PostStep = PostStep_OutletConservation;
1659 bc->UpdateUbcs = NULL;
1660 bc->Destroy = NULL; // No private data to destroy
1661
1662 bc->data = NULL;
1663
1664 PetscFunctionReturn(0);
1665}
1666
1667#undef __FUNCT__
1668#define __FUNCT__ "PreStep_OutletConservation"
1669/**
1670 * @brief Prepare the outlet-conservation correction before advancing the solver.
1671 */
1672static PetscErrorCode PreStep_OutletConservation(BoundaryCondition *self, BCContext *ctx,
1673 PetscReal *local_inflow_contribution, PetscReal *local_outflow_contribution)
1674{
1675 PetscErrorCode ierr;
1676 UserCtx* user = ctx->user;
1677 BCFace face_id = ctx->face_id;
1678 DMDALocalInfo* info = &user->info;
1679 PetscBool can_service;
1680
1681 // Suppress unused parameter warnings for clarity.
1682 (void)self;
1683 (void)local_inflow_contribution;
1684
1685 PetscFunctionBeginUser;
1686
1687 // Step 1: Use the robust utility function to determine if this MPI rank owns a computable
1688 // portion of the specified boundary face. If not, there is no work to do, so we exit immediately.
1689 const PetscInt IM_nodes_global = user->IM;
1690 const PetscInt JM_nodes_global = user->JM;
1691 const PetscInt KM_nodes_global = user->KM;
1692 ierr = CanRankServiceFace(info, IM_nodes_global, JM_nodes_global, KM_nodes_global, face_id, &can_service); CHKERRQ(ierr);
1693
1694 if (!can_service) {
1695 PetscFunctionReturn(0);
1696 }
1697
1698 // Step 2: Get read-only access to the necessary PETSc arrays.
1699 // We use the local versions (`lUcat`, `lNvert`) which include ghost cell data,
1700 // ensuring we have the correct interior values adjacent to the boundary.
1701 Cmpnts ***ucat, ***csi, ***eta, ***zet;
1702 PetscReal ***nvert;
1703 ierr = DMDAVecGetArrayRead(user->fda, user->lUcat, (const Cmpnts***)&ucat); CHKERRQ(ierr);
1704 ierr = DMDAVecGetArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
1705 ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
1706 ierr = DMDAVecGetArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
1707 ierr = DMDAVecGetArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
1708
1709 PetscReal local_flux_out = 0.0;
1710 const PetscInt xs=info->xs, xe=info->xs+info->xm;
1711 const PetscInt ys=info->ys, ye=info->ys+info->ym;
1712 const PetscInt zs=info->zs, ze=info->zs+info->zm;
1713 const PetscInt mx=info->mx, my=info->my, mz=info->mz;
1714
1715 // Step 3: Replicate the legacy shrunk loop bounds to exclude corners and edges.
1716 PetscInt lxs = xs; if (xs == 0) lxs = xs + 1;
1717 PetscInt lxe = xe; if (xe == mx) lxe = xe - 1;
1718 PetscInt lys = ys; if (ys == 0) lys = ys + 1;
1719 PetscInt lye = ye; if (ye == my) lye = ye - 1;
1720 PetscInt lzs = zs; if (zs == 0) lzs = zs + 1;
1721 PetscInt lze = ze; if (ze == mz) lze = ze - 1;
1722
1723 // Step 4: Loop over the specified face using the corrected bounds and indexing to calculate flux.
1724 switch (face_id) {
1725 case BC_FACE_NEG_X: {
1726 const PetscInt i_cell = xs + 1; // Index for first interior cell-centered data
1727 const PetscInt i_face = xs; // Index for the -X face of that cell
1728 for (int k=lzs; k<lze; k++) for (int j=lys; j<lye; j++) {
1729 if (nvert[k][j][i_cell] < 0.1) {
1730 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);
1731 }
1732 }
1733 break;
1734 }
1735 case BC_FACE_POS_X: {
1736 const PetscInt i_cell = xe - 2; // Index for last interior cell-centered data
1737 const PetscInt i_face = xe - 2; // Index for the +X face of that cell
1738 for (int k=lzs; k<lze; k++) for (int j=lys; j<lye; j++) {
1739 if (nvert[k][j][i_cell] < 0.1) {
1740 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);
1741 }
1742 }
1743 break;
1744 }
1745 case BC_FACE_NEG_Y: {
1746 const PetscInt j_cell = ys + 1;
1747 const PetscInt j_face = ys;
1748 for (int k=lzs; k<lze; k++) for (int i=lxs; i<lxe; i++) {
1749 if (nvert[k][j_cell][i] < 0.1) {
1750 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);
1751 }
1752 }
1753 break;
1754 }
1755 case BC_FACE_POS_Y: {
1756 const PetscInt j_cell = ye - 2;
1757 const PetscInt j_face = ye - 2;
1758 for (int k=lzs; k<lze; k++) for (int i=lxs; i<lxe; i++) {
1759 if (nvert[k][j_cell][i] < 0.1) {
1760 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);
1761 }
1762 }
1763 break;
1764 }
1765 case BC_FACE_NEG_Z: {
1766 const PetscInt k_cell = zs + 1;
1767 const PetscInt k_face = zs;
1768 for (int j=lys; j<lye; j++) for (int i=lxs; i<lxe; i++) {
1769 if (nvert[k_cell][j][i] < 0.1) {
1770 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);
1771 }
1772 }
1773 break;
1774 }
1775 case BC_FACE_POS_Z: {
1776 const PetscInt k_cell = ze - 2;
1777 const PetscInt k_face = ze - 2;
1778 for (int j=lys; j<lye; j++) for (int i=lxs; i<lxe; i++) {
1779 if (nvert[k_cell][j][i] < 0.1) {
1780 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);
1781 }
1782 }
1783 break;
1784 }
1785 }
1786
1787 // Step 5: Restore the PETSc arrays.
1788 ierr = DMDAVecRestoreArrayRead(user->fda, user->lUcat, (const Cmpnts***)&ucat); CHKERRQ(ierr);
1789 ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
1790 ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
1791 ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
1792 ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
1793
1794 // Step 6: Add this face's calculated flux to the accumulator for this rank.
1795 *local_outflow_contribution += local_flux_out;
1796
1797 PetscFunctionReturn(0);
1798}
1799
1800#undef __FUNCT__
1801#define __FUNCT__ "Apply_OutletConservation"
1802/**
1803 * @brief (Handler Action) Applies mass conservation correction to the outlet face.
1804 *
1805 * This function calculates a global correction factor based on the total inflow and outflow fluxes
1806 * and applies it to the contravariant velocity (`ucont`) on the outlet face to ensure mass conservation.
1807 */
1808static PetscErrorCode Apply_OutletConservation(BoundaryCondition *self, BCContext *ctx)
1809{
1810 PetscErrorCode ierr;
1811 (void)self;
1812 UserCtx* user = ctx->user;
1813 BCFace face_id = ctx->face_id;
1814 DMDALocalInfo* info = &user->info;
1815 PetscBool can_service;
1816
1817 PetscFunctionBeginUser;
1818 PROFILE_FUNCTION_BEGIN;
1819
1820 const PetscInt IM_nodes_global = user->IM;
1821 const PetscInt JM_nodes_global = user->JM;
1822 const PetscInt KM_nodes_global = user->KM;
1823 ierr = CanRankServiceFace(info, IM_nodes_global, JM_nodes_global, KM_nodes_global, face_id, &can_service); CHKERRQ(ierr);
1824
1825 if (!can_service) {
1826 PROFILE_FUNCTION_END;
1827 PetscFunctionReturn(0);
1828 }
1829
1830 // --- STEP 1: Calculate the correction factor using pre-calculated area ---
1831 PetscReal total_inflow = *ctx->global_inflow_sum + *ctx->global_farfield_inflow_sum;
1832 PetscReal flux_imbalance = total_inflow - *ctx->global_outflow_sum;
1833
1834 // Directly use the pre-calculated area from the simulation context.
1835 PetscReal velocity_correction = (PetscAbsReal(user->simCtx->AreaOutSum) > 1e-12)
1836 ? flux_imbalance / user->simCtx->AreaOutSum
1837 : 0.0;
1838
1839 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Outlet Correction on Face %d: Imbalance=%.4e, Pre-calc Area=%.4e, V_corr=%.4e\n",
1840 face_id, flux_imbalance, user->simCtx->AreaOutSum, velocity_correction);
1841
1842 // --- STEP 2: Apply the correction to ucont on the outlet face ---
1843
1844 // Get read/write access to necessary arrays
1845
1846 Cmpnts ***ubcs, ***ucont, ***csi, ***eta, ***zet, ***ucat;
1847 PetscReal ***nvert;
1848 ierr = DMDAVecGetArray(user->fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
1849 ierr = DMDAVecGetArray(user->fda, user->Ucont, &ucont); CHKERRQ(ierr);
1850 ierr = DMDAVecGetArrayRead(user->fda,user->lUcat, (const Cmpnts***)&ucat); CHKERRQ(ierr);
1851 ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
1852 ierr = DMDAVecGetArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
1853 ierr = DMDAVecGetArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
1854 ierr = DMDAVecGetArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
1855
1856 // Get local grid bounds to exclude corners/edges
1857 PetscInt xs = info->xs, xe = info->xs + info->xm;
1858 PetscInt ys = info->ys, ye = info->ys + info->ym;
1859 PetscInt zs = info->zs, ze = info->zs + info->zm;
1860 PetscInt mx = info->mx, my = info->my, mz = info->mz;
1861 PetscInt lxs = xs, lxe = xe, lys = ys, lye = ye, lzs = zs, lze = ze;
1862
1863 if (xs == 0) lxs = xs + 1;
1864 if (xe == mx) lxe = xe - 1;
1865 if (ys == 0) lys = ys + 1;
1866 if (ye == my) lye = ye - 1;
1867 if (zs == 0) lzs = zs + 1;
1868 if (ze == mz) lze = ze - 1;
1869
1870 // Loop over faces and apply correction
1871 switch(face_id){
1872 case BC_FACE_NEG_X:{
1873 const PetscInt i_cell = xs + 1;
1874 const PetscInt i_face = xs;
1875 const PetscInt i_dummy = xs;
1876 for (PetscInt k = lzs; k < lze; k++) {
1877 for (PetscInt j = lys; j < lye; j++) {
1878 if (nvert[k][j][i_cell] < 0.1) {
1879 // Set ubcs
1880 ubcs[k][j][i_dummy] = ucat[k][j][i_cell];
1881
1882 // Calculate Local uncorrected original flux
1883 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);
1884
1885 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));
1886
1887 PetscReal Correction_flux = velocity_correction*Cell_Area;
1888
1889 ucont[k][j][i_face].x = Uncorrected_local_flux + Correction_flux;
1890 }
1891 }
1892 }
1893 break;
1894 }
1895 case BC_FACE_POS_X:{
1896 const PetscInt i_cell = xe - 2;
1897 const PetscInt i_face = xe - 2;
1898 const PetscInt i_dummy = xe - 1;
1899 for(PetscInt k = lzs; k < lze; k++) for (PetscInt j = lys; j < lye; j++){
1900 if(nvert[k][j][i_cell]<0.1){
1901 // Set ubcs
1902 ubcs[k][j][i_dummy] = ucat[k][j][i_cell];
1903
1904 // Calculate Local uncorrected original flux
1905 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);
1906
1907 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));
1908
1909 PetscReal Correction_flux = velocity_correction*Cell_Area;
1910
1911 ucont[k][j][i_face].x = Uncorrected_local_flux + Correction_flux;
1912 }
1913 }
1914 break;
1915 }
1916 case BC_FACE_NEG_Y:{
1917 const PetscInt j_cell = ys + 1;
1918 const PetscInt j_face = ys;
1919 const PetscInt j_dummy = ys;
1920 for(PetscInt k = lzs; k < lze; k++) for (PetscInt i = lxs; i < lxe; i++){
1921 if(nvert[k][j_cell][i]<0.1){
1922 // Set ubcs
1923 ubcs[k][j_dummy][i] = ucat[k][j_cell][i];
1924
1925 // Calculate Local uncorrected original flux
1926 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);
1927
1928 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));
1929
1930 PetscReal Correction_flux = velocity_correction*Cell_Area;
1931
1932 ucont[k][j_face][i].y = Uncorrected_local_flux + Correction_flux;
1933 }
1934 }
1935 break;
1936 }
1937 case BC_FACE_POS_Y:{
1938 const PetscInt j_cell = ye - 2;
1939 const PetscInt j_face = ye - 2;
1940 const PetscInt j_dummy = ye - 1;
1941 for(PetscInt k = lzs; k < lze; k++) for (PetscInt i = lxs; i < lxe; i++){
1942 if(nvert[k][j_cell][i]<0.1){
1943 // Set ubcs
1944 ubcs[k][j_dummy][i] = ucat[k][j_cell][i];
1945
1946 // Calculate Local uncorrected original flux
1947 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);
1948
1949 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));
1950
1951 PetscReal Correction_flux = velocity_correction*Cell_Area;
1952
1953 ucont[k][j_face][i].y = Uncorrected_local_flux + Correction_flux;
1954 }
1955 }
1956 break;
1957 }
1958 case BC_FACE_NEG_Z:{
1959 const PetscInt k_cell = zs + 1;
1960 const PetscInt k_face = zs;
1961 const PetscInt k_dummy = zs;
1962 for(PetscInt j = lys; j < lye; j++) for (PetscInt i = lxs; i < lxe; i++){
1963 if(nvert[k_cell][j][i]<0.1){
1964 // Set ubcs
1965 ubcs[k_dummy][j][i] = ucat[k_cell][j][i];
1966
1967 // Calculate Local uncorrected original flux
1968 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));
1969
1970 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));
1971
1972 PetscReal Correction_flux = velocity_correction*Cell_Area;
1973
1974 ucont[k_face][j][i].z = Uncorrected_local_flux + Correction_flux;
1975 }
1976 }
1977 break;
1978 }
1979 case BC_FACE_POS_Z:{
1980 const PetscInt k_cell = ze - 2;
1981 const PetscInt k_face = ze - 2;
1982 const PetscInt k_dummy = ze - 1;
1983 for(PetscInt j = lys; j < lye; j++) for (PetscInt i = lxs; i < lxe; i++){
1984 if(nvert[k_cell][j][i]<0.1){
1985 // Set ubcs
1986 ubcs[k_dummy][j][i] = ucat[k_cell][j][i];
1987
1988 // Calculate Local uncorrected original flux
1989 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));
1990
1991 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));
1992
1993 PetscReal Correction_flux = velocity_correction*Cell_Area;
1994
1995 ucont[k_face][j][i].z = Uncorrected_local_flux + Correction_flux;
1996 }
1997 }
1998 break;
1999 }
2000 }
2001
2002 // Restore all arrays
2003 ierr = DMDAVecRestoreArray(user->fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
2004 ierr = DMDAVecRestoreArray(user->fda, user->Ucont, &ucont); CHKERRQ(ierr);
2005 ierr = DMDAVecRestoreArrayRead(user->fda,user->lUcat, (const Cmpnts***)&ucat); CHKERRQ(ierr);
2006 ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
2007 ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
2008 ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
2009 ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
2010
2011 PROFILE_FUNCTION_END;
2012 PetscFunctionReturn(0);
2013}
2014
2015#undef __FUNCT__
2016#define __FUNCT__ "PostStep_OutletConservation"
2017/**
2018 * @brief Update outlet-conservation state after a completed solver step.
2019 */
2020static PetscErrorCode PostStep_OutletConservation(BoundaryCondition *self, BCContext *ctx,
2021 PetscReal *local_inflow_contribution,
2022 PetscReal *local_outflow_contribution)
2023{
2024 PetscErrorCode ierr;
2025 UserCtx* user = ctx->user;
2026 BCFace face_id = ctx->face_id;
2027 DMDALocalInfo* info = &user->info;
2028 PetscBool can_service;
2029
2030 (void)self;
2031 (void)local_inflow_contribution;
2032
2033 PetscFunctionBeginUser;
2034 const PetscInt IM_nodes_global = user->IM;
2035 const PetscInt JM_nodes_global = user->JM;
2036 const PetscInt KM_nodes_global = user->KM;
2037 ierr = CanRankServiceFace(info, IM_nodes_global, JM_nodes_global, KM_nodes_global, face_id, &can_service); CHKERRQ(ierr);
2038
2039 if (!can_service) PetscFunctionReturn(0);
2040
2041 // Get arrays (need both ucont and nvert)
2042 Cmpnts ***ucont;
2043 PetscReal ***nvert; // ✅ ADD nvert
2044
2045 ierr = DMDAVecGetArrayRead(user->fda, user->Ucont, (const Cmpnts***)&ucont); CHKERRQ(ierr);
2046 ierr = DMDAVecGetArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr); // ✅ ADD
2047
2048 PetscReal local_flux = 0.0;
2049
2050 PetscInt xs = info->xs, xe = info->xs + info->xm;
2051 PetscInt ys = info->ys, ye = info->ys + info->ym;
2052 PetscInt zs = info->zs, ze = info->zs + info->zm;
2053 PetscInt mx = info->mx, my = info->my, mz = info->mz;
2054
2055 PetscInt lxs = xs, lxe = xe, lys = ys, lye = ye, lzs = zs, lze = ze;
2056 if (xs == 0) lxs = xs + 1;
2057 if (xe == mx) lxe = xe - 1;
2058 if (ys == 0) lys = ys + 1;
2059 if (ye == my) lye = ye - 1;
2060 if (zs == 0) lzs = zs + 1;
2061 if (ze == mz) lze = ze - 1;
2062
2063 // Sum ucont components, skipping solid cells (same indices as PreStep)
2064 switch (face_id) {
2065 case BC_FACE_NEG_X: {
2066 const PetscInt i_cell = xs + 1; // ✅ Match PreStep
2067 const PetscInt i_face = xs;
2068 for (PetscInt k = lzs; k < lze; k++) {
2069 for (PetscInt j = lys; j < lye; j++) {
2070 if (nvert[k][j][i_cell] < 0.1) { // ✅ Skip solid cells
2071 local_flux += ucont[k][j][i_face].x;
2072 }
2073 }
2074 }
2075 } break;
2076
2077 case BC_FACE_POS_X: {
2078 const PetscInt i_cell = xe - 2; // ✅ Match PreStep
2079 const PetscInt i_face = xe - 2;
2080 for (PetscInt k = lzs; k < lze; k++) {
2081 for (PetscInt j = lys; j < lye; j++) {
2082 if (nvert[k][j][i_cell] < 0.1) { // ✅ Skip solid cells
2083 local_flux += ucont[k][j][i_face].x;
2084 }
2085 }
2086 }
2087 } break;
2088
2089 case BC_FACE_NEG_Y: {
2090 const PetscInt j_cell = ys + 1; // ✅ Match PreStep
2091 const PetscInt j_face = ys;
2092 for (PetscInt k = lzs; k < lze; k++) {
2093 for (PetscInt i = lxs; i < lxe; i++) {
2094 if (nvert[k][j_cell][i] < 0.1) { // ✅ Skip solid cells
2095 local_flux += ucont[k][j_face][i].y;
2096 }
2097 }
2098 }
2099 } break;
2100
2101 case BC_FACE_POS_Y: {
2102 const PetscInt j_cell = ye - 2; // ✅ Match PreStep
2103 const PetscInt j_face = ye - 2;
2104 for (PetscInt k = lzs; k < lze; k++) {
2105 for (PetscInt i = lxs; i < lxe; i++) {
2106 if (nvert[k][j_cell][i] < 0.1) { // ✅ Skip solid cells
2107 local_flux += ucont[k][j_face][i].y;
2108 }
2109 }
2110 }
2111 } break;
2112
2113 case BC_FACE_NEG_Z: {
2114 const PetscInt k_cell = zs + 1; // ✅ Match PreStep
2115 const PetscInt k_face = zs;
2116 for (PetscInt j = lys; j < lye; j++) {
2117 for (PetscInt i = lxs; i < lxe; i++) {
2118 if (nvert[k_cell][j][i] < 0.1) { // ✅ Skip solid cells
2119 local_flux += ucont[k_face][j][i].z;
2120 }
2121 }
2122 }
2123 } break;
2124
2125 case BC_FACE_POS_Z: {
2126 const PetscInt k_cell = ze - 2; // ✅ Match PreStep
2127 const PetscInt k_face = ze - 2;
2128 for (PetscInt j = lys; j < lye; j++) {
2129 for (PetscInt i = lxs; i < lxe; i++) {
2130 if (nvert[k_cell][j][i] < 0.1) { // ✅ Skip solid cells
2131 local_flux += ucont[k_face][j][i].z;
2132 }
2133 }
2134 }
2135 } break;
2136 }
2137
2138 // Restore arrays
2139 ierr = DMDAVecRestoreArrayRead(user->fda, user->Ucont, (const Cmpnts***)&ucont); CHKERRQ(ierr);
2140 ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr); // ✅ ADD
2141
2142 // Add to accumulator
2143 *local_outflow_contribution += local_flux;
2144
2145 LOG_ALLOW(LOCAL, LOG_DEBUG, "PostStep_OutletConservation: Face %d, corrected flux = %.6e\n",
2146 face_id, local_flux);
2147
2148 PetscFunctionReturn(0);
2149}
2150
2151
2152/**
2153 * @brief Implementation of \ref Create_PeriodicGeometric().
2154 * @details Full API contract (arguments, ownership, side effects) is documented with
2155 * the header declaration in `include/BC_Handlers.h`.
2156 * @see Create_PeriodicGeometric()
2157 */
2158
2159PetscErrorCode Create_PeriodicGeometric(BoundaryCondition *bc){
2160 PetscFunctionBeginUser;
2161
2162 if (!bc) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Input BoundaryCondition is NULL");
2163 bc->priority = BC_PRIORITY_WALL;
2164
2165 // Assign function pointers
2166 bc->Initialize = NULL; // No initialization needed
2167 bc->PreStep = NULL;
2168 bc->Apply = NULL;
2169 bc->PostStep = NULL;
2170 bc->UpdateUbcs = NULL;
2171 bc->Destroy = NULL; // No private data to destroy
2172
2173 bc->data = NULL;
2174
2175 PetscFunctionReturn(0);
2176}
2177
2178
2179#undef __FUNCT__
2180#define __FUNCT__ "MeasureDrivenFluxes"
2181/**
2182 * @brief Measure the two volumetric fluxes the driven-flow controller senses.
2183 *
2184 * Both periodic driven handlers steer on the same pair of measurements, so they
2185 * are taken here in a single sweep of `lUcont`:
2186 *
2187 * - `*boundaryFlux` is the flux through the single periodic boundary plane.
2188 * It is fast and responsive but noisy, and drives the boundary trim.
2189 * - `*planarAverageFlux` is the flux averaged over every cross-sectional plane
2190 * in the driven direction. It is stable and inertial, and drives the
2191 * momentum source. It is also the quantity `initial_flux` latches at t=0.
2192 *
2193 * @param[in] user Block context supplying `lUcont`, `lNvert` and `info`.
2194 * @param[in] direction Driven direction, 'X', 'Y' or 'Z'.
2195 * @param[out] boundaryFlux Globally reduced flux through the boundary plane.
2196 * @param[out] planarAverageFlux Globally reduced plane-averaged flux.
2197 * @return PetscErrorCode 0 on success.
2198 */
2199static PetscErrorCode MeasureDrivenFluxes(UserCtx *user, char direction,
2200 PetscReal *boundaryFlux,
2201 PetscReal *planarAverageFlux)
2202{
2203 PetscErrorCode ierr;
2204 DMDALocalInfo info = user->info;
2205 PetscInt i, j, k;
2206
2207 PetscFunctionBeginUser;
2208
2209 // --- Get read-only access to necessary field data ---
2210 Cmpnts ***ucont;
2211 PetscReal ***nvert;
2212 ierr = DMDAVecGetArrayRead(user->fda, user->lUcont, (const Cmpnts***)&ucont); CHKERRQ(ierr);
2213 ierr = DMDAVecGetArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
2214
2215 // --- Define local loop bounds ---
2216 PetscInt lxs = (info.xs == 0) ? 1 : info.xs;
2217 PetscInt lys = (info.ys == 0) ? 1 : info.ys;
2218 PetscInt lzs = (info.zs == 0) ? 1 : info.zs;
2219 PetscInt lxe = (info.xs + info.xm == info.mx) ? info.mx - 1 : info.xs + info.xm;
2220 PetscInt lye = (info.ys + info.ym == info.my) ? info.my - 1 : info.ys + info.ym;
2221 PetscInt lze = (info.zs + info.zm == info.mz) ? info.mz - 1 : info.zs + info.zm;
2222
2223 // --- Initialize local accumulators ---
2224 PetscReal localCurrentBoundaryFlux = 0.0;
2225 PetscReal localAveragePlanarVolumetricFluxTerm = 0.0;
2226
2227 // --- Measure local contributions to the two flux types, generalized by direction ---
2228 switch (direction) {
2229 case 'X':
2230 if (info.xs == 0) { // Only the rank on the negative face contributes to boundary flux
2231 i = 0;
2232 for (k = lzs; k < lze; k++) for (j = lys; j < lye; j++) {
2233 if (nvert[k][j][i + 1] < 0.1) localCurrentBoundaryFlux += ucont[k][j][i].x;
2234 }
2235 }
2236 for (i = info.xs; i < lxe; i++) {
2237 for (k = lzs; k < lze; k++) for (j = lys; j < lye; j++) {
2238 if (nvert[k][j][i + 1] < 0.1) localAveragePlanarVolumetricFluxTerm += ucont[k][j][i].x / (PetscReal)(info.mx - 1);
2239 }
2240 }
2241 break;
2242 case 'Y':
2243 if (info.ys == 0) {
2244 j = 0;
2245 for (k = lzs; k < lze; k++) for (i = lxs; i < lxe; i++) {
2246 if (nvert[k][j + 1][i] < 0.1) localCurrentBoundaryFlux += ucont[k][j][i].y;
2247 }
2248 }
2249 for (j = info.ys; j < lye; j++) {
2250 for (k = lzs; k < lze; k++) for (i = lxs; i < lxe; i++) {
2251 if (nvert[k][j + 1][i] < 0.1) localAveragePlanarVolumetricFluxTerm += ucont[k][j][i].y / (PetscReal)(info.my - 1);
2252 }
2253 }
2254 break;
2255 case 'Z':
2256 if (info.zs == 0) {
2257 k = 0;
2258 for (j = lys; j < lye; j++) for (i = lxs; i < lxe; i++) {
2259 if (nvert[k + 1][j][i] < 0.1) localCurrentBoundaryFlux += ucont[k][j][i].z;
2260 }
2261 }
2262 for (k = info.zs; k < lze; k++) {
2263 for (j = lys; j < lye; j++) for (i = lxs; i < lxe; i++) {
2264 if (nvert[k + 1][j][i] < 0.1) localAveragePlanarVolumetricFluxTerm += ucont[k][j][i].z / (PetscReal)(info.mz - 1);
2265 }
2266 }
2267 break;
2268 default:
2269 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
2270 "MeasureDrivenFluxes received an unknown driven direction '%c'.", direction);
2271 }
2272
2273 // --- Release array access as soon as possible ---
2274 ierr = DMDAVecRestoreArrayRead(user->fda, user->lUcont, (const Cmpnts***)&ucont); CHKERRQ(ierr);
2275 ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
2276
2277 // --- Perform global reductions to get the final flux values ---
2278 ierr = MPI_Allreduce(&localCurrentBoundaryFlux, boundaryFlux, 1, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
2279 ierr = MPI_Allreduce(&localAveragePlanarVolumetricFluxTerm, planarAverageFlux, 1, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
2280
2281 PetscFunctionReturn(0);
2282}
2283
2284// ===============================================================================
2285//
2286// HANDLER IMPLEMENTATION: PERIODIC DRIVEN CONSTANT FLUX
2287// (Corresponds to BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX)
2288//
2289// ===============================================================================
2290
2291// --- 1. FORWARD DECLARATIONS & PRIVATE DATA ---
2292
2293// Forward declarations for the static functions that implement this handler's behavior.
2294static PetscErrorCode Initialize_PeriodicDrivenConstant(BoundaryCondition *self, BCContext *ctx);
2295static PetscErrorCode PreStep_PeriodicDrivenConstant(BoundaryCondition *self, BCContext *ctx, PetscReal *in, PetscReal *out);
2296static PetscErrorCode Apply_PeriodicDrivenConstant(BoundaryCondition *self, BCContext *ctx);
2297static PetscErrorCode Destroy_PeriodicDrivenConstant(BoundaryCondition *self);
2298
2299/**
2300 * @brief Private data structure shared by both periodic driven-flux handlers.
2301 *
2302 * `constant_flux` fills `targetVolumetricFlux` from the bcs file at
2303 * initialization; `initial_flux` latches it from the starting field at the
2304 * first PreStep. Everything downstream of the target is identical, so both
2305 * handlers reuse this struct and the PreStep/Apply/Destroy implementations
2306 * below.
2307 */
2308typedef struct {
2309 char direction; // 'X', 'Y', or 'Z', determined at initialization.
2310 PetscReal targetVolumetricFlux; // The target flux this controller drives to.
2311 PetscBool isMasterController; // Flag: PETSC_TRUE only for the handler on the negative face.
2312 PetscBool enforceSeamFlux; // Flag: PETSC_TRUE to add the seam-flux correction into Ucont.
2313 PetscInt lastBulkCorrectionStep; // Physical step at which the momentum source was last set (-1 = never).
2314} DrivenFluxData;
2315
2316
2317// --- 2. HANDLER CONSTRUCTOR ---
2318
2319#undef __FUNCT__
2320#define __FUNCT__ "Create_PeriodicDrivenConstant"
2321/**
2322 * @brief Internal helper implementation: `Create_PeriodicDrivenConstant()`.
2323 * @details Local to this translation unit.
2324 */
2325PetscErrorCode Create_PeriodicDrivenConstant(BoundaryCondition *bc)
2326{
2327 PetscErrorCode ierr;
2328 PetscFunctionBeginUser;
2329
2330 if (!bc) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Input BoundaryCondition object is NULL in Create_PeriodicDrivenConstantFlux");
2331
2332 // --- Allocate the private data structure ---
2333 DrivenFluxData *data = NULL;
2334 ierr = PetscNew(&data); CHKERRQ(ierr);
2335 // Initialize fields to safe default values
2336 data->direction = ' ';
2337 data->targetVolumetricFlux = 0.0;
2338 data->isMasterController = PETSC_FALSE;
2339 data->enforceSeamFlux = PETSC_FALSE;
2340 data->lastBulkCorrectionStep = -1;
2341
2342 // Attach the private data to the generic handler object
2343 bc->data = (void*)data;
2344
2345 // --- Configure the handler's properties and methods ---
2346
2347 // Set priority: Using BC_PRIORITY_INLET ensures this handler's PreStep runs
2348 // before other handlers (like outlets) that might depend on its calculations.
2349 // It is the caller's responsibility that there are no Inlets called along with driven periodic to avoid clash.
2350 bc->priority = BC_PRIORITY_INLET;
2351
2352 // Assign the function pointers to the implementations in this file.
2353 bc->Initialize = Initialize_PeriodicDrivenConstant;
2354 bc->PreStep = PreStep_PeriodicDrivenConstant;
2355 bc->Apply = Apply_PeriodicDrivenConstant;
2356 bc->PostStep = NULL; // This handler has no action after the main solver step.
2357 bc->UpdateUbcs = NULL; // The boundary value is not flow-dependent (it's periodic).
2358 bc->Destroy = Destroy_PeriodicDrivenConstant;
2359
2360 PetscFunctionReturn(0);
2361}
2362
2363#undef __FUNCT__
2364#define __FUNCT__ "Initialize_PeriodicDrivenConstant"
2365/**
2366 * @brief Initialize constant forcing data for a periodically driven boundary.
2367 */
2368static PetscErrorCode Initialize_PeriodicDrivenConstant(BoundaryCondition *self, BCContext *ctx)
2369{
2370 PetscErrorCode ierr;
2371 DrivenFluxData *data = (DrivenFluxData*)self->data;
2372 BCFace face_id = ctx->face_id;
2373 UserCtx* user = ctx->user;
2374
2375 PetscFunctionBeginUser;
2376
2377 LOG_ALLOW(LOCAL, LOG_DEBUG, "Initializing PERIODIC_DRIVEN_CONSTANT_FLUX handler on Face %s...\n", BCFaceToString(face_id));
2378
2379 // --- 1. Validation: Ensure the mathematical type is PERIODIC ---
2380 if (user->boundary_faces[face_id].mathematical_type != PERIODIC) {
2381 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT,
2382 "Configuration Error: Handler PERIODIC_DRIVEN_CONSTANT_FLUX on Face %s must be applied to a face with mathematical_type PERIODIC.",
2383 BCFaceToString(face_id));
2384 }
2385
2386 // --- 2. Role Assignment: Determine direction and master status ---
2387 data->isMasterController = PETSC_FALSE;
2388 switch (face_id) {
2389 case BC_FACE_NEG_X: data->direction = 'X'; data->isMasterController = PETSC_TRUE; break;
2390 case BC_FACE_POS_X: data->direction = 'X'; break;
2391 case BC_FACE_NEG_Y: data->direction = 'Y'; data->isMasterController = PETSC_TRUE; break;
2392 case BC_FACE_POS_Y: data->direction = 'Y'; break;
2393 case BC_FACE_NEG_Z: data->direction = 'Z'; data->isMasterController = PETSC_TRUE; break;
2394 case BC_FACE_POS_Z: data->direction = 'Z'; break;
2395 }
2396
2397 // --- 3. Parameter Parsing (Master Controller only) ---
2398 if (data->isMasterController) {
2399 PetscBool found;
2400
2401 // Attempt to read the 'target_flux' parameter from the bcs.run file.
2402 ierr = GetBCParamReal(user->boundary_faces[face_id].params, "target_flux",
2403 &data->targetVolumetricFlux, &found); CHKERRQ(ierr);
2404
2405 // If the required parameter is not found, halt with an informative error.
2406 if (!found) {
2407 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT,
2408 "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).",
2409 BCFaceToString(face_id));
2410 }
2411
2412 LOG_ALLOW(GLOBAL, LOG_INFO, "Driven Flow (Dir %c): Constant target volumetric flux set to %le.\n",
2413 data->direction, data->targetVolumetricFlux);
2414
2415 // Store the target flux in the UserCtx. This makes it globally accessible
2416 // to other parts of the solver, such as the `CorrectChannelFluxProfile` enforcer function.
2417 user->simCtx->targetVolumetricFlux = data->targetVolumetricFlux;
2418 // The target is fixed for the run from here on; see the flag's comment in
2419 // variables.h for how initial_flux uses the same latch.
2420 user->simCtx->drivenFluxTargetLatched = PETSC_TRUE;
2421 }
2422
2423 PetscBool trimfound;
2424 // Optional seam-flux enforcement; accepts the deprecated `apply_trim` spelling.
2425 ierr = GetDrivenSeamFluxFlag(user->boundary_faces[face_id].params,
2426 &data->enforceSeamFlux, &trimfound); CHKERRQ(ierr);
2427
2428 if(!trimfound) LOG_ALLOW(GLOBAL,LOG_DEBUG,"Seam-flux enforcement not specified, defaults to %s.\n",data->enforceSeamFlux? "True":"False");
2429
2430 PetscFunctionReturn(0);
2431}
2432
2433#undef __FUNCT__
2434#define __FUNCT__ "PreStep_PeriodicDrivenConstant"
2435/**
2436 * @brief Prepare constant periodic-driving data before the solver step.
2437 */
2438static PetscErrorCode PreStep_PeriodicDrivenConstant(BoundaryCondition *self, BCContext *ctx,
2439 PetscReal *local_inflow_contribution,
2440 PetscReal *local_outflow_contribution)
2441{
2442 PetscErrorCode ierr;
2443 DrivenFluxData *data = (DrivenFluxData*)self->data;
2444 UserCtx* user = ctx->user;
2445 SimCtx* simCtx = user->simCtx;
2446
2447 PetscFunctionBeginUser;
2448
2449 // --- Master Check: Only the handler on the negative face performs calculations ---
2450 if (!data->isMasterController) {
2451 PetscFunctionReturn(0);
2452 }
2453
2454 // The controller senses two fluxes; see MeasureDrivenFluxes() for what each
2455 // one is for and why the controller needs both.
2456 char direction = data->direction;
2457 PetscReal globalCurrentBoundaryFlux, globalAveragePlanarVolumetricFlux;
2458 ierr = MeasureDrivenFluxes(user, direction, &globalCurrentBoundaryFlux,
2459 &globalAveragePlanarVolumetricFlux); CHKERRQ(ierr);
2460
2461 // --- Get cross-sectional area using the dedicated geometry function ---
2462 Cmpnts ignored_center;
2463 PetscReal globalBoundaryArea;
2464 BCFace neg_face_id = (direction == 'X') ? BC_FACE_NEG_X : (direction == 'Y') ? BC_FACE_NEG_Y : BC_FACE_NEG_Z;
2465 ierr = CalculateFaceCenterAndArea(user, neg_face_id, &ignored_center, &globalBoundaryArea); CHKERRQ(ierr);
2466
2467 // --- Calculate the two correction terms ---
2468 //
2469 // These are refreshed on DIFFERENT cadences, and the difference matters.
2470 //
2471 // ApplyBoundaryConditions() runs BoundarySystem_ExecuteStep() -- and hence
2472 // this PreStep -- three times per call, and it is itself called once per
2473 // Jameson RK stage under the Picard solver and once per residual evaluation
2474 // under Newton-Krylov. So "per PreStep" is emphatically not "per timestep".
2475 //
2476 // - bulkVelocityCorrection scales the momentum source in ComputeRHS. Both
2477 // momentum solvers are built on it being FROZEN across a timestep: the
2478 // Picard shadow-Jacobian estimate treats the body force as a constant
2479 // forcing with zero velocity Jacobian, and the Newton solve needs a
2480 // source that does not drift between residual evaluations. It is
2481 // therefore computed once per physical step, from the field at the start
2482 // of that step, and held.
2483 //
2484 // - boundaryVelocityCorrection is the tactical trim applied to the
2485 // boundary fluxes in Apply(). It is deliberately re-measured on every
2486 // pass: Apply() accumulates it into Ucont, and re-measuring is what makes
2487 // that accumulation self-limiting as the seam converges. Freezing it
2488 // would make repeated Apply() calls add the same trim over and over.
2489 if (globalBoundaryArea > 1.0e-12) {
2490 if (data->lastBulkCorrectionStep != simCtx->step) {
2491 simCtx->bulkVelocityCorrection = (data->targetVolumetricFlux - globalAveragePlanarVolumetricFlux) / globalBoundaryArea;
2492 data->lastBulkCorrectionStep = simCtx->step;
2493 }
2494 simCtx->boundaryVelocityCorrection = (data->targetVolumetricFlux - globalCurrentBoundaryFlux) / globalBoundaryArea;
2495 } else {
2496 simCtx->bulkVelocityCorrection = 0.0;
2497 simCtx->boundaryVelocityCorrection = 0.0;
2498 data->lastBulkCorrectionStep = simCtx->step;
2499 }
2500
2501 LOG_ALLOW(GLOBAL, LOG_INFO, "Driven Flow Controller Update (Dir %c):\n", data->direction);
2502 LOG_ALLOW(GLOBAL, LOG_INFO, " - Target Volumetric Flux: %.6e\n", data->targetVolumetricFlux);
2503 LOG_ALLOW(GLOBAL, LOG_INFO, " - Avg Planar Volumetric Flux (Stable): %.6e\n", globalAveragePlanarVolumetricFlux);
2504 LOG_ALLOW(GLOBAL, LOG_INFO, " - Boundary Flux (Fast): %.6e\n", globalCurrentBoundaryFlux);
2505 LOG_ALLOW(GLOBAL, LOG_INFO, " - Bulk Velocity Correction: %.6e (For Momentum Source)\n", simCtx->bulkVelocityCorrection);
2506 LOG_ALLOW(GLOBAL, LOG_INFO, " - Boundary Velocity Correction: %.6e (For Boundary Trim)\n", simCtx->boundaryVelocityCorrection);
2507
2508 // Suppress unused parameter warnings for this handler
2509 (void)local_inflow_contribution;
2510 (void)local_outflow_contribution;
2511
2512 PetscFunctionReturn(0);
2513}
2514
2515#undef __FUNCT__
2516#define __FUNCT__ "Apply_PeriodicDrivenConstant"
2517/**
2518 * @brief Apply the configured constant driving term to the periodic boundary.
2519 */
2520static PetscErrorCode Apply_PeriodicDrivenConstant(BoundaryCondition *self, BCContext *ctx)
2521{
2522 PetscErrorCode ierr;
2523 DrivenFluxData *data = (DrivenFluxData*)self->data;
2524 UserCtx* user = ctx->user;
2525 BCFace face_id = ctx->face_id;
2526 PetscBool can_service;
2527
2528 PetscFunctionBeginUser;
2529
2530 // Check if this rank owns part of this boundary face
2531 ierr = CanRankServiceFace(&user->info, user->IM, user->JM, user->KM, face_id, &can_service); CHKERRQ(ierr);
2532 if (!can_service) {
2533 PetscFunctionReturn(0);
2534 }
2535
2536 // If the correction is negligible, no work is needed.
2537 if (PetscAbsReal(user->simCtx->boundaryVelocityCorrection) < 1e-12) {
2538 PetscFunctionReturn(0);
2539 }
2540
2541 LOG_ALLOW(LOCAL, LOG_TRACE, "Apply_PeriodicDrivenConstant: Applying boundary trim on Face %s...\n", BCFaceToString(face_id));
2542
2543 // --- Get read/write access to necessary arrays ---
2544 DMDALocalInfo info = user->info;
2545 Cmpnts ***ucont, ***uch, ***csi, ***eta, ***zet;
2546 PetscReal ***nvert;
2547
2548 ierr = DMDAVecGetArray(user->fda, user->Ucont, &ucont); CHKERRQ(ierr);
2549 ierr = DMDAVecGetArray(user->fda, user->Bcs.Uch, &uch); CHKERRQ(ierr);
2550 ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
2551 ierr = DMDAVecGetArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
2552 ierr = DMDAVecGetArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
2553 ierr = DMDAVecGetArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
2554
2555 PetscInt lxs = (info.xs == 0) ? 1 : info.xs;
2556 PetscInt lys = (info.ys == 0) ? 1 : info.ys;
2557 PetscInt lzs = (info.zs == 0) ? 1 : info.zs;
2558 PetscInt lxe = (info.xs + info.xm == info.mx) ? info.mx - 1 : info.xs + info.xm;
2559 PetscInt lye = (info.ys + info.ym == info.my) ? info.my - 1 : info.ys + info.ym;
2560 PetscInt lze = (info.zs + info.zm == info.mz) ? info.mz - 1 : info.zs + info.zm;
2561
2562 // --- Apply correction to the appropriate face and velocity component ---
2563 switch (face_id) {
2564 case BC_FACE_NEG_X: case BC_FACE_POS_X: {
2565 PetscInt i_face = (face_id == BC_FACE_NEG_X) ? info.xs : info.mx - 2;
2566 PetscInt i_nvert = (face_id == BC_FACE_NEG_X) ? info.xs + 1 : info.mx - 2;
2567
2568 for (PetscInt k = lzs; k < lze; k++) for (PetscInt j = lys; j < lye; j++) {
2569 if (nvert[k][j][i_nvert] < 0.1) {
2570 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);
2571 PetscReal fluxTrim = user->simCtx->boundaryVelocityCorrection * faceArea;
2572 if(data->enforceSeamFlux) ucont[k][j][i_face].x += fluxTrim;
2573 uch[k][j][i_face].x = fluxTrim; // Store correction for diagnostics
2574 }
2575 }
2576 } break;
2577
2578 case BC_FACE_NEG_Y: case BC_FACE_POS_Y: {
2579 PetscInt j_face = (face_id == BC_FACE_NEG_Y) ? info.ys : info.my - 2;
2580 PetscInt j_nvert = (face_id == BC_FACE_NEG_Y) ? info.ys + 1 : info.my - 2;
2581
2582 for (PetscInt k = lzs; k < lze; k++) for (PetscInt i = lxs; i < lxe; i++) {
2583 if (nvert[k][j_nvert][i] < 0.1) {
2584 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);
2585 PetscReal fluxTrim = user->simCtx->boundaryVelocityCorrection * faceArea;
2586 if(data->enforceSeamFlux) ucont[k][j_face][i].y += fluxTrim;
2587 uch[k][j_face][i].y = fluxTrim;
2588 }
2589 }
2590 } break;
2591
2592 case BC_FACE_NEG_Z: case BC_FACE_POS_Z: {
2593 PetscInt k_face = (face_id == BC_FACE_NEG_Z) ? info.zs : info.mz - 2;
2594 PetscInt k_nvert = (face_id == BC_FACE_NEG_Z) ? info.zs + 1 : info.mz - 2;
2595
2596 for (PetscInt j = lys; j < lye; j++) for (PetscInt i = lxs; i < lxe; i++) {
2597 if (nvert[k_nvert][j][i] < 0.1) {
2598 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);
2599 PetscReal fluxTrim = user->simCtx->boundaryVelocityCorrection * faceArea;
2600 if(data->enforceSeamFlux) ucont[k_face][j][i].z += fluxTrim;
2601 uch[k_face][j][i].z = fluxTrim;
2602 }
2603 }
2604 } break;
2605 }
2606
2607 // --- Restore arrays ---
2608 ierr = DMDAVecRestoreArray(user->fda, user->Ucont, &ucont); CHKERRQ(ierr);
2609 ierr = DMDAVecRestoreArray(user->fda, user->Bcs.Uch, &uch); CHKERRQ(ierr);
2610 ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
2611 ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
2612 ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
2613 ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
2614
2615 PetscFunctionReturn(0);
2616}
2617
2618#undef __FUNCT__
2619#define __FUNCT__ "Destroy_PeriodicDrivenConstant"
2620/**
2621 * @brief Release resources owned by the constant periodic-driving boundary.
2622 */
2623static PetscErrorCode Destroy_PeriodicDrivenConstant(BoundaryCondition *self)
2624{
2625 PetscFunctionBeginUser;
2626
2627 // Check that the handler object and its private data pointer are valid before trying to free.
2628 if (self && self->data) {
2629 // The private data was allocated with PetscNew(), so it must be freed with PetscFree().
2630 PetscFree(self->data);
2631
2632 // It is good practice to nullify the pointer after freeing to prevent
2633 // any accidental use of the dangling pointer (use-after-free).
2634 self->data = NULL;
2635
2636 LOG_ALLOW(LOCAL, LOG_TRACE, "Destroy_PeriodicDrivenConstant: Private data freed successfully.\n");
2637 }
2638
2639 PetscFunctionReturn(0);
2640}
2641
2642
2643// ===============================================================================
2644//
2645// HANDLER IMPLEMENTATION: PERIODIC DRIVEN INITIAL FLUX
2646// (Corresponds to BC_HANDLER_PERIODIC_DRIVEN_INITIAL_FLUX)
2647//
2648// Identical to the CONSTANT_FLUX handler except for where the target comes
2649// from: this one measures the volumetric flux of the field the run starts
2650// with and then holds it, so it takes no `target_flux` parameter. Once the
2651// target is latched the two handlers behave the same, so PreStep delegates
2652// to the constant implementation and Apply/Destroy are shared outright.
2653//
2654// ===============================================================================
2655
2656// --- 1. FORWARD DECLARATIONS ---
2657
2658static PetscErrorCode Initialize_PeriodicDrivenInitial(BoundaryCondition *self, BCContext *ctx);
2659static PetscErrorCode PreStep_PeriodicDrivenInitial(BoundaryCondition *self, BCContext *ctx, PetscReal *in, PetscReal *out);
2660
2661// --- 2. HANDLER CONSTRUCTOR ---
2662
2663#undef __FUNCT__
2664#define __FUNCT__ "Create_PeriodicDrivenInitial"
2665/**
2666 * @brief Internal helper implementation: `Create_PeriodicDrivenInitial()`.
2667 * @details Local to this translation unit.
2668 */
2669PetscErrorCode Create_PeriodicDrivenInitial(BoundaryCondition *bc)
2670{
2671 PetscErrorCode ierr;
2672 PetscFunctionBeginUser;
2673
2674 if (!bc) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Input BoundaryCondition object is NULL in Create_PeriodicDrivenInitial");
2675
2676 // --- Allocate the private data structure ---
2677 DrivenFluxData *data = NULL;
2678 ierr = PetscNew(&data); CHKERRQ(ierr);
2679 // Initialize fields to safe default values
2680 data->direction = ' ';
2681 data->targetVolumetricFlux = 0.0;
2682 data->isMasterController = PETSC_FALSE;
2683 data->enforceSeamFlux = PETSC_FALSE;
2684 data->lastBulkCorrectionStep = -1;
2685
2686 // Attach the private data to the generic handler object
2687 bc->data = (void*)data;
2688
2689 // --- Configure the handler's properties and methods ---
2690
2691 // Same priority reasoning as the constant-flux handler: PreStep must run
2692 // before any handler that depends on the controller's corrections.
2693 bc->priority = BC_PRIORITY_INLET;
2694
2695 bc->Initialize = Initialize_PeriodicDrivenInitial;
2696 bc->PreStep = PreStep_PeriodicDrivenInitial;
2697 bc->Apply = Apply_PeriodicDrivenConstant; // Boundary trim is target-agnostic.
2698 bc->PostStep = NULL; // This handler has no action after the main solver step.
2699 bc->UpdateUbcs = NULL; // The boundary value is not flow-dependent (it's periodic).
2700 bc->Destroy = Destroy_PeriodicDrivenConstant; // Same private data layout.
2701
2702 PetscFunctionReturn(0);
2703}
2704
2705#undef __FUNCT__
2706#define __FUNCT__ "Initialize_PeriodicDrivenInitial"
2707/**
2708 * @brief Initialize forcing data for a periodic boundary driven to its initial flux.
2709 *
2710 * @note The target itself is NOT measured here. Boundary handlers are
2711 * initialized before `InitializeEulerianState()` runs, so at this point
2712 * `Ucont` still holds zeros. The measurement is deferred to the first
2713 * PreStep, by which time either the initial condition has been applied or
2714 * the restart target has been restored from the checkpoint manifest.
2715 */
2716static PetscErrorCode Initialize_PeriodicDrivenInitial(BoundaryCondition *self, BCContext *ctx)
2717{
2718 PetscErrorCode ierr;
2719 DrivenFluxData *data = (DrivenFluxData*)self->data;
2720 BCFace face_id = ctx->face_id;
2721 UserCtx* user = ctx->user;
2722
2723 PetscFunctionBeginUser;
2724
2725 LOG_ALLOW(LOCAL, LOG_DEBUG, "Initializing PERIODIC_DRIVEN_INITIAL_FLUX handler on Face %s...\n", BCFaceToString(face_id));
2726
2727 // --- 1. Validation: Ensure the mathematical type is PERIODIC ---
2728 if (user->boundary_faces[face_id].mathematical_type != PERIODIC) {
2729 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT,
2730 "Configuration Error: Handler PERIODIC_DRIVEN_INITIAL_FLUX on Face %s must be applied to a face with mathematical_type PERIODIC.",
2731 BCFaceToString(face_id));
2732 }
2733
2734 // --- 2. Role Assignment: Determine direction and master status ---
2735 data->isMasterController = PETSC_FALSE;
2736 switch (face_id) {
2737 case BC_FACE_NEG_X: data->direction = 'X'; data->isMasterController = PETSC_TRUE; break;
2738 case BC_FACE_POS_X: data->direction = 'X'; break;
2739 case BC_FACE_NEG_Y: data->direction = 'Y'; data->isMasterController = PETSC_TRUE; break;
2740 case BC_FACE_POS_Y: data->direction = 'Y'; break;
2741 case BC_FACE_NEG_Z: data->direction = 'Z'; data->isMasterController = PETSC_TRUE; break;
2742 case BC_FACE_POS_Z: data->direction = 'Z'; break;
2743 }
2744
2745 // --- 3. Parameter Parsing (Master Controller only) ---
2746 if (data->isMasterController) {
2747 PetscReal unused_flux;
2748 PetscBool found;
2749
2750 // This handler derives its own target, so an explicit one is a config error
2751 // rather than something to silently ignore. Users who want to prescribe the
2752 // flux should select the `constant_flux` handler instead.
2753 ierr = GetBCParamReal(user->boundary_faces[face_id].params, "target_flux",
2754 &unused_flux, &found); CHKERRQ(ierr);
2755 if (found) {
2756 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT,
2757 "Configuration Error: Handler PERIODIC_DRIVEN_INITIAL_FLUX on Face %s takes no 'target_flux' parameter; "
2758 "it measures the flux of the initial condition and holds that. Use handler 'constant_flux' to prescribe a target.",
2759 BCFaceToString(face_id));
2760 }
2761
2762 LOG_ALLOW(GLOBAL, LOG_INFO, "Driven Flow (Dir %c): target volumetric flux will be latched from the initial state.\n",
2763 data->direction);
2764 }
2765
2766 PetscBool trimfound;
2767 // Optional seam-flux enforcement; accepts the deprecated `apply_trim` spelling.
2768 ierr = GetDrivenSeamFluxFlag(user->boundary_faces[face_id].params,
2769 &data->enforceSeamFlux, &trimfound); CHKERRQ(ierr);
2770
2771 if(!trimfound) LOG_ALLOW(GLOBAL,LOG_DEBUG,"Seam-flux enforcement not specified, defaults to %s.\n",data->enforceSeamFlux? "True":"False");
2772
2773 PetscFunctionReturn(0);
2774}
2775
2776#undef __FUNCT__
2777#define __FUNCT__ "PreStep_PeriodicDrivenInitial"
2778/**
2779 * @brief Latch the initial-state flux once, then drive to it like a constant target.
2780 *
2781 * @details The latch is one-shot and guarded by `simCtx->drivenFluxTargetLatched`:
2782 * - Fresh start: the flag is false and `Ucont` now holds the initial
2783 * condition, so the plane-averaged flux is measured and stored.
2784 * - Restart: `ReadSimulationFields()` restored both the target and the
2785 * flag from the checkpoint manifest, so the original target survives
2786 * instead of being re-measured from a drifted field.
2787 */
2788static PetscErrorCode PreStep_PeriodicDrivenInitial(BoundaryCondition *self, BCContext *ctx,
2789 PetscReal *local_inflow_contribution,
2790 PetscReal *local_outflow_contribution)
2791{
2792 PetscErrorCode ierr;
2793 DrivenFluxData *data = (DrivenFluxData*)self->data;
2794 SimCtx* simCtx = ctx->user->simCtx;
2795
2796 PetscFunctionBeginUser;
2797
2798 if (data->isMasterController) {
2799 if (!simCtx->drivenFluxTargetLatched) {
2800 PetscReal boundaryFlux, planarAverageFlux;
2801
2802 /* Boundary handlers are also exercised once during setup, from
2803 * FinalizeBlockState(), and at that point the initial condition has
2804 * been written to Ucont but not yet scattered into the ghosted
2805 * lUcont that MeasureDrivenFluxes() reads. Latching there would
2806 * record a target of zero. Setup runs with step == StartStep, so
2807 * wait for the first PreStep of the first real timestep: by then
2808 * lUcont holds the initial condition (or, on a restart, the state
2809 * the restored target already describes). Do no work at all until
2810 * then, so no bogus correction is derived from a zero target. */
2811 if (simCtx->step <= simCtx->StartStep) {
2812 simCtx->bulkVelocityCorrection = 0.0;
2813 simCtx->boundaryVelocityCorrection = 0.0;
2814 PetscFunctionReturn(0);
2815 }
2816
2817 ierr = MeasureDrivenFluxes(ctx->user, data->direction,
2818 &boundaryFlux, &planarAverageFlux); CHKERRQ(ierr);
2819
2820 simCtx->targetVolumetricFlux = planarAverageFlux;
2821 simCtx->drivenFluxTargetLatched = PETSC_TRUE;
2822
2823 LOG_ALLOW(GLOBAL, LOG_INFO,
2824 "Driven Flow (Dir %c): latched initial volumetric flux %.6e as the target.\n",
2825 data->direction, planarAverageFlux);
2826 }
2827 // Keep the handler's copy in step with the authoritative value in SimCtx,
2828 // which is also where a restart deposits the restored target.
2829 data->targetVolumetricFlux = simCtx->targetVolumetricFlux;
2830 }
2831
2832 ierr = PreStep_PeriodicDrivenConstant(self, ctx,
2833 local_inflow_contribution,
2834 local_outflow_contribution); CHKERRQ(ierr);
2835
2836 PetscFunctionReturn(0);
2837}
InletParabolicData::cs2_half
PetscReal cs2_half
Half-width (in index space) in cross-stream direction 2.
Definition BC_Handlers.c:704
Initialize_InletProfileFromFile
static PetscErrorCode Initialize_InletProfileFromFile(BoundaryCondition *self, BCContext *ctx)
Initializes a file-prescribed inlet profile handler for one boundary face.
Definition BC_Handlers.c:1333
DrivenFluxData::enforceSeamFlux
PetscBool enforceSeamFlux
Definition BC_Handlers.c:2312
Destroy_InletProfileFromFile
static PetscErrorCode Destroy_InletProfileFromFile(BoundaryCondition *self)
Releases private storage owned by a file-prescribed inlet profile handler.
Definition BC_Handlers.c:1600
Create_InletConstantVelocity
PetscErrorCode Create_InletConstantVelocity(BoundaryCondition *bc)
Implementation of Create_InletConstantVelocity().
Definition BC_Handlers.c:322
Apply_WallNoSlip
static PetscErrorCode Apply_WallNoSlip(BoundaryCondition *self, BCContext *ctx)
Apply no-slip velocity values to wall-adjacent cells for this boundary condition.
Definition BC_Handlers.c:143
DrivenFluxData::direction
char direction
Definition BC_Handlers.c:2309
Create_InletProfileFromFile
PetscErrorCode Create_InletProfileFromFile(BoundaryCondition *bc)
Implementation of Create_InletProfileFromFile().
Definition BC_Handlers.c:1292
DrivenFluxData::isMasterController
PetscBool isMasterController
Definition BC_Handlers.c:2311
Initialize_InletParabolicProfile
static PetscErrorCode Initialize_InletParabolicProfile(BoundaryCondition *self, BCContext *ctx)
Initialize the geometric data used to evaluate a parabolic inlet profile.
Definition BC_Handlers.c:743
GetProfileFileExpectedDims
static PetscErrorCode GetProfileFileExpectedDims(UserCtx *user, BCFace face_id, PetscInt *n1, PetscInt *n2)
Computes the expected PICSLICE dimensions for an inlet face.
Definition BC_Handlers.c:1149
PreStep_InletConstantVelocity
static PetscErrorCode PreStep_InletConstantVelocity(BoundaryCondition *self, BCContext *ctx, PetscReal *in, PetscReal *out)
Update constant-inlet data required before the next solver step.
Definition BC_Handlers.c:399
PreStep_InletProfileFromFile
static PetscErrorCode PreStep_InletProfileFromFile(BoundaryCondition *self, BCContext *ctx, PetscReal *in, PetscReal *out)
Pre-step hook for the static file-prescribed inlet profile handler.
Definition BC_Handlers.c:1378
Apply_InletConstantVelocity
static PetscErrorCode Apply_InletConstantVelocity(BoundaryCondition *self, BCContext *ctx)
Impose the configured constant velocity on inlet boundary cells.
Definition BC_Handlers.c:421
DrivenFluxData::lastBulkCorrectionStep
PetscInt lastBulkCorrectionStep
Definition BC_Handlers.c:2313
InletParabolicData::cs2_center
PetscReal cs2_center
Center index in cross-stream direction 2.
Definition BC_Handlers.c:702
Apply_OutletConservation
static PetscErrorCode Apply_OutletConservation(BoundaryCondition *self, BCContext *ctx)
(Handler Action) Applies mass conservation correction to the outlet face.
Definition BC_Handlers.c:1808
Apply_InletParabolicProfile
static PetscErrorCode Apply_InletParabolicProfile(BoundaryCondition *self, BCContext *ctx)
Impose the evaluated parabolic velocity profile on inlet cells.
Definition BC_Handlers.c:827
PostStep_InletParabolicProfile
static PetscErrorCode PostStep_InletParabolicProfile(BoundaryCondition *self, BCContext *ctx, PetscReal *in, PetscReal *out)
Perform post-step bookkeeping for a parabolic inlet boundary.
Definition BC_Handlers.c:984
InletProfileFileData::n2
PetscInt n2
Definition BC_Handlers.c:1108
InletProfileFileData::profile
PetscReal * profile
Definition BC_Handlers.c:1109
ReadPicSliceProfile
static PetscErrorCode ReadPicSliceProfile(const char *source_file, PetscInt expected_n1, PetscInt expected_n2, InletProfileFileData *data)
Reads and validates a static scalar inlet profile from a canonical PICSLICE file.
Definition BC_Handlers.c:1194
Apply_PeriodicDrivenConstant
static PetscErrorCode Apply_PeriodicDrivenConstant(BoundaryCondition *self, BCContext *ctx)
Apply the configured constant driving term to the periodic boundary.
Definition BC_Handlers.c:2520
Initialize_PeriodicDrivenConstant
static PetscErrorCode Initialize_PeriodicDrivenConstant(BoundaryCondition *self, BCContext *ctx)
Initialize constant forcing data for a periodically driven boundary.
Definition BC_Handlers.c:2368
Create_PeriodicGeometric
PetscErrorCode Create_PeriodicGeometric(BoundaryCondition *bc)
Implementation of Create_PeriodicGeometric().
Definition BC_Handlers.c:2159
InletParabolicData::v_max
PetscReal v_max
Peak centerline velocity (from user params).
Definition BC_Handlers.c:700
Validate_DrivenFlowConfiguration
PetscErrorCode Validate_DrivenFlowConfiguration(UserCtx *user)
Internal helper implementation: Validate_DrivenFlowConfiguration().
Definition BC_Handlers.c:15
Create_InletParabolicProfile
PetscErrorCode Create_InletParabolicProfile(BoundaryCondition *bc)
Implementation of Create_InletParabolicProfile().
Definition BC_Handlers.c:715
Initialize_PeriodicDrivenInitial
static PetscErrorCode Initialize_PeriodicDrivenInitial(BoundaryCondition *self, BCContext *ctx)
Initialize forcing data for a periodic boundary driven to its initial flux.
Definition BC_Handlers.c:2716
Destroy_InletConstantVelocity
static PetscErrorCode Destroy_InletConstantVelocity(BoundaryCondition *self)
Release resources owned by a constant-velocity inlet boundary.
Definition BC_Handlers.c:650
Create_PeriodicDrivenInitial
PetscErrorCode Create_PeriodicDrivenInitial(BoundaryCondition *bc)
Internal helper implementation: Create_PeriodicDrivenInitial().
Definition BC_Handlers.c:2669
Destroy_PeriodicDrivenConstant
static PetscErrorCode Destroy_PeriodicDrivenConstant(BoundaryCondition *self)
Release resources owned by the constant periodic-driving boundary.
Definition BC_Handlers.c:2623
PostStep_InletProfileFromFile
static PetscErrorCode PostStep_InletProfileFromFile(BoundaryCondition *self, BCContext *ctx, PetscReal *in, PetscReal *out)
Accumulates the applied inlet flux for a file-prescribed profile.
Definition BC_Handlers.c:1524
Create_PeriodicDrivenConstant
PetscErrorCode Create_PeriodicDrivenConstant(BoundaryCondition *bc)
Internal helper implementation: Create_PeriodicDrivenConstant().
Definition BC_Handlers.c:2325
PreStep_InletParabolicProfile
static PetscErrorCode PreStep_InletParabolicProfile(BoundaryCondition *self, BCContext *ctx, PetscReal *in, PetscReal *out)
Refresh parabolic-inlet values required before the solver step.
Definition BC_Handlers.c:808
MeasureDrivenFluxes
static PetscErrorCode MeasureDrivenFluxes(UserCtx *user, char direction, PetscReal *boundaryFlux, PetscReal *planarAverageFlux)
Measure the two volumetric fluxes the driven-flow controller senses.
Definition BC_Handlers.c:2199
PostStep_OutletConservation
static PetscErrorCode PostStep_OutletConservation(BoundaryCondition *self, BCContext *ctx, PetscReal *in, PetscReal *out)
Update outlet-conservation state after a completed solver step.
Definition BC_Handlers.c:2020
InletProfileFileData::max_speed
PetscReal max_speed
Definition BC_Handlers.c:1111
InletConstantData::normal_velocity
PetscReal normal_velocity
Definition BC_Handlers.c:311
ProfileSpeedAt
static PetscReal ProfileSpeedAt(const InletProfileFileData *data, PetscInt a, PetscInt b)
Returns one scalar speed from the flattened PICSLICE profile.
Definition BC_Handlers.c:1279
Destroy_InletParabolicProfile
static PetscErrorCode Destroy_InletParabolicProfile(BoundaryCondition *self)
Release resources owned by a parabolic inlet boundary.
Definition BC_Handlers.c:1077
Initialize_InletConstantVelocity
static PetscErrorCode Initialize_InletConstantVelocity(BoundaryCondition *self, BCContext *ctx)
Initialize persistent state for a constant-velocity inlet boundary.
Definition BC_Handlers.c:350
Create_WallNoSlip
PetscErrorCode Create_WallNoSlip(BoundaryCondition *bc)
Implementation of Create_WallNoSlip().
Definition BC_Handlers.c:114
Apply_InletProfileFromFile
static PetscErrorCode Apply_InletProfileFromFile(BoundaryCondition *self, BCContext *ctx)
Applies the loaded PICSLICE scalar profile to Ucont and Ubcs on an inlet face.
Definition BC_Handlers.c:1404
PreStep_PeriodicDrivenConstant
static PetscErrorCode PreStep_PeriodicDrivenConstant(BoundaryCondition *self, BCContext *ctx, PetscReal *in, PetscReal *out)
Prepare constant periodic-driving data before the solver step.
Definition BC_Handlers.c:2438
PostStep_InletConstantVelocity
static PetscErrorCode PostStep_InletConstantVelocity(BoundaryCondition *self, BCContext *ctx, PetscReal *in, PetscReal *out)
Perform post-step bookkeeping for a constant-velocity inlet boundary.
Definition BC_Handlers.c:556
PreStep_OutletConservation
static PetscErrorCode PreStep_OutletConservation(BoundaryCondition *self, BCContext *ctx, PetscReal *local_inflow_contribution, PetscReal *local_outflow_contribution)
Prepare the outlet-conservation correction before advancing the solver.
Definition BC_Handlers.c:1672
Create_OutletConservation
PetscErrorCode Create_OutletConservation(BoundaryCondition *bc)
Implementation of Create_OutletConservation().
Definition BC_Handlers.c:1644
InletProfileFileData::n1
PetscInt n1
Definition BC_Handlers.c:1107
PreStep_PeriodicDrivenInitial
static PetscErrorCode PreStep_PeriodicDrivenInitial(BoundaryCondition *self, BCContext *ctx, PetscReal *in, PetscReal *out)
Latch the initial-state flux once, then drive to it like a constant target.
Definition BC_Handlers.c:2788
DrivenFluxData::targetVolumetricFlux
PetscReal targetVolumetricFlux
Definition BC_Handlers.c:2310
GetBCParamStringLocal
static PetscErrorCode GetBCParamStringLocal(BC_Param *params, const char *key, const char **value_out, PetscBool *found)
Looks up a string-valued boundary-condition parameter in a BC_Param list.
Definition BC_Handlers.c:1124
InletParabolicData::cs1_half
PetscReal cs1_half
Half-width (in index space) in cross-stream direction 1.
Definition BC_Handlers.c:703
InletParabolicData::cs1_center
PetscReal cs1_center
Center index in cross-stream direction 1.
Definition BC_Handlers.c:701
InletProfileFileData::min_speed
PetscReal min_speed
Definition BC_Handlers.c:1110
InletProfileFileData::source_file
char * source_file
Definition BC_Handlers.c:1112
DrivenFluxData
Private data structure shared by both periodic driven-flux handlers.
Definition BC_Handlers.c:2308
InletConstantData
Private data structure for the Constant Velocity Inlet handler.
Definition BC_Handlers.c:310
InletParabolicData
Private data structure for the Parabolic Velocity Inlet handler.
Definition BC_Handlers.c:699
InletProfileFileData
Definition BC_Handlers.c:1106
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:127
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:1207
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:752
GetDrivenSeamFluxFlag
PetscErrorCode GetDrivenSeamFluxFlag(BC_Param *params, PetscBool *value_out, PetscBool *found)
Read the driven-flow seam-flux flag, accepting its deprecated apply_trim spelling.
Definition io.c:814
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)
Returns the canonical log token for a boundary-face enum value.
Definition logging.c:671
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:859
BCTypeToString
const char * BCTypeToString(BCType type)
Returns the canonical log token for a boundary mathematical type.
Definition logging.c:773
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:30
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:850
BoundaryCondition
The "virtual table" struct for a boundary condition handler object.
Definition variables.h:353
BoundaryCondition::PostStep
PetscErrorCode(* PostStep)(BoundaryCondition *self, BCContext *ctx, PetscReal *local_inflow, PetscReal *local_outflow)
Definition variables.h:360
BoundaryCondition::PreStep
PetscErrorCode(* PreStep)(BoundaryCondition *self, BCContext *ctx, PetscReal *local_inflow, PetscReal *local_outflow)
Definition variables.h:358
BoundaryCondition::Destroy
PetscErrorCode(* Destroy)(BoundaryCondition *self)
Definition variables.h:362
BoundaryCondition::Initialize
PetscErrorCode(* Initialize)(BoundaryCondition *self, BCContext *ctx)
Definition variables.h:357
BoundaryCondition::UpdateUbcs
PetscErrorCode(* UpdateUbcs)(BoundaryCondition *self, BCContext *ctx)
Definition variables.h:361
BoundaryCondition::Apply
PetscErrorCode(* Apply)(BoundaryCondition *self, BCContext *ctx)
Definition variables.h:359
BoundaryCondition::priority
BCPriorityType priority
Definition variables.h:355
BCContext::global_outflow_sum
const PetscReal * global_outflow_sum
Definition variables.h:349
BCType
BCType
Defines the general mathematical/physical Category of a boundary.
Definition variables.h:283
INLET
@ INLET
Definition variables.h:290
FARFIELD
@ FARFIELD
Definition variables.h:291
OUTLET
@ OUTLET
Definition variables.h:289
PERIODIC
@ PERIODIC
Definition variables.h:292
UserCtx::boundary_faces
BoundaryFaceConfig boundary_faces[6]
Definition variables.h:931
SimCtx::targetVolumetricFlux
PetscReal targetVolumetricFlux
Definition variables.h:807
UserCtx::lNvert
Vec lNvert
Definition variables.h:939
UserCtx::simCtx
SimCtx * simCtx
Back-pointer to the master simulation context.
Definition variables.h:909
BC_Param_s::next
struct BC_Param_s * next
Definition variables.h:339
SimCtx::boundaryVelocityCorrection
PetscReal boundaryVelocityCorrection
Definition variables.h:814
UserCtx::KM
PetscInt KM
Definition variables.h:920
UserCtx::lZet
Vec lZet
Definition variables.h:974
BCHandlerType
BCHandlerType
Defines the specific computational "strategy" for a boundary handler.
Definition variables.h:303
BC_HANDLER_PERIODIC_DRIVEN_INITIAL_FLUX
@ BC_HANDLER_PERIODIC_DRIVEN_INITIAL_FLUX
Definition variables.h:319
BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX
@ BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX
Definition variables.h:318
BoundaryFaceConfig::handler_type
BCHandlerType handler_type
Definition variables.h:369
SimCtx::drivenFluxTargetLatched
PetscBool drivenFluxTargetLatched
Definition variables.h:812
SimCtx::bulkVelocityCorrection
PetscReal bulkVelocityCorrection
Definition variables.h:813
BCContext::face_id
BCFace face_id
Definition variables.h:345
UserCtx::Ucont
Vec Ucont
Definition variables.h:939
SimCtx::StartStep
PetscInt StartStep
Definition variables.h:705
BCS::Ubcs
Vec Ubcs
Physical Cartesian velocity at boundary faces. Full 3D array but only boundary-face entries are meani...
Definition variables.h:123
Cmpnts::x
PetscScalar x
Definition variables.h:103
UserCtx::Bcs
BCS Bcs
Definition variables.h:934
BCContext::user
UserCtx * user
Definition variables.h:344
BCContext::global_inflow_sum
const PetscReal * global_inflow_sum
Definition variables.h:346
UserCtx::lCsi
Vec lCsi
Definition variables.h:974
BoundaryFaceConfig::params
BC_Param * params
Definition variables.h:370
Cmpnts::z
PetscScalar z
Definition variables.h:103
UserCtx::JM
PetscInt JM
Definition variables.h:920
BCContext::global_farfield_inflow_sum
const PetscReal * global_farfield_inflow_sum
Definition variables.h:347
UserCtx::lUcont
Vec lUcont
Definition variables.h:939
SimCtx::step
PetscInt step
Definition variables.h:703
SimCtx::AreaOutSum
PetscReal AreaOutSum
Definition variables.h:815
UserCtx::info
DMDALocalInfo info
Definition variables.h:918
BC_PRIORITY_OUTLET
@ BC_PRIORITY_OUTLET
Definition variables.h:328
BC_PRIORITY_WALL
@ BC_PRIORITY_WALL
Definition variables.h:327
BC_PRIORITY_INLET
@ BC_PRIORITY_INLET
Definition variables.h:325
UserCtx::lUcat
Vec lUcat
Definition variables.h:939
Cmpnts::y
PetscScalar y
Definition variables.h:103
UserCtx::IM
PetscInt IM
Definition variables.h:920
UserCtx::lEta
Vec lEta
Definition variables.h:974
BoundaryFaceConfig::mathematical_type
BCType mathematical_type
Definition variables.h:368
BCS::Uch
Vec Uch
Characteristic velocity for boundary conditions.
Definition variables.h:124
BCFace
BCFace
Identifies the six logical faces of a structured computational block.
Definition variables.h:261
BC_FACE_NEG_X
@ BC_FACE_NEG_X
Definition variables.h:262
BC_FACE_POS_Z
@ BC_FACE_POS_Z
Definition variables.h:264
BC_FACE_POS_Y
@ BC_FACE_POS_Y
Definition variables.h:263
BC_FACE_NEG_Z
@ BC_FACE_NEG_Z
Definition variables.h:264
BC_FACE_POS_X
@ BC_FACE_POS_X
Definition variables.h:262
BC_FACE_NEG_Y
@ BC_FACE_NEG_Y
Definition variables.h:263
BCContext
Provides execution context for a boundary condition handler.
Definition variables.h:343
BC_Param_s
A node in a linked list for storing key-value parameters from the bcs.dat file.
Definition variables.h:336
BoundaryFaceConfig
Holds the complete configuration for one of the six boundary faces.
Definition variables.h:366
Cmpnts
A 3D point or vector with PetscScalar components.
Definition variables.h:102
SimCtx
The master context for the entire simulation.
Definition variables.h:695
UserCtx
User-defined context containing data specific to a single computational grid level.
Definition variables.h:906
Generated by doxygen 1.9.8