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