BC_Handlers.c File Reference

#include "BC_Handlers.h"

Include dependency graph for BC_Handlers.c:

Go to the source code of this file.

Data Structures
struct	InletConstantData
	Private data structure for the Constant Velocity Inlet handler. More...
struct	InletParabolicData
	Private data structure for the Parabolic Velocity Inlet handler. More...
struct	InletProfileFileData
struct	DrivenFluxData
	Private data structure shared by both periodic driven-flux handlers. More...

Macros
#define	__FUNCT__   "Validate_DrivenFlowConfiguration"
#define	__FUNCT__   "Create_WallNoSlip"
#define	__FUNCT__   "Apply_WallNoSlip"
#define	__FUNCT__   "Create_InletConstantVelocity"
#define	__FUNCT__   "Initialize_InletConstantVelocity"
#define	__FUNCT__   "PreStep_InletConstantVelocity"
#define	__FUNCT__   "Apply_InletConstantVelocity"
#define	__FUNCT__   "PostStep_InletConstantVelocity"
#define	__FUNCT__   "Destroy_InletConstantVelocity"
#define	__FUNCT__   "Create_InletParabolicProfile"
#define	__FUNCT__   "Initialize_InletParabolicProfile"
#define	__FUNCT__   "PreStep_InletParabolicProfile"
#define	__FUNCT__   "Apply_InletParabolicProfile"
#define	__FUNCT__   "PostStep_InletParabolicProfile"
#define	__FUNCT__   "Destroy_InletParabolicProfile"
#define	__FUNCT__   "Create_InletProfileFromFile"
#define	__FUNCT__   "Initialize_InletProfileFromFile"
#define	__FUNCT__   "PreStep_InletProfileFromFile"
#define	__FUNCT__   "Apply_InletProfileFromFile"
#define	__FUNCT__   "PostStep_InletProfileFromFile"
#define	__FUNCT__   "Destroy_InletProfileFromFile"
#define	__FUNCT__   "Create_OutletConservation"
#define	__FUNCT__   "PreStep_OutletConservation"
#define	__FUNCT__   "Apply_OutletConservation"
#define	__FUNCT__   "PostStep_OutletConservation"
#define	__FUNCT__   "MeasureDrivenFluxes"
#define	__FUNCT__   "Create_PeriodicDrivenConstant"
#define	__FUNCT__   "Initialize_PeriodicDrivenConstant"
#define	__FUNCT__   "PreStep_PeriodicDrivenConstant"
#define	__FUNCT__   "Apply_PeriodicDrivenConstant"
#define	__FUNCT__   "Destroy_PeriodicDrivenConstant"
#define	__FUNCT__   "Create_PeriodicDrivenInitial"
#define	__FUNCT__   "Initialize_PeriodicDrivenInitial"
#define	__FUNCT__   "PreStep_PeriodicDrivenInitial"

Functions
PetscErrorCode	Validate_DrivenFlowConfiguration (UserCtx *user)
	Internal helper implementation: Validate_DrivenFlowConfiguration().
static PetscErrorCode	Apply_WallNoSlip (BoundaryCondition *self, BCContext *ctx)
	Apply no-slip velocity values to wall-adjacent cells for this boundary condition.
PetscErrorCode	Create_WallNoSlip (BoundaryCondition *bc)
	Implementation of Create_WallNoSlip().
static PetscErrorCode	Initialize_InletConstantVelocity (BoundaryCondition *self, BCContext *ctx)
	Initialize persistent state for a constant-velocity inlet boundary.
static PetscErrorCode	PreStep_InletConstantVelocity (BoundaryCondition *self, BCContext *ctx, PetscReal *local_inflow_contribution, PetscReal *local_outflow_contribution)
	Update constant-inlet data required before the next solver step.
static PetscErrorCode	Apply_InletConstantVelocity (BoundaryCondition *self, BCContext *ctx)
	Impose the configured constant velocity on inlet boundary cells.
static PetscErrorCode	PostStep_InletConstantVelocity (BoundaryCondition *self, BCContext *ctx, PetscReal *local_inflow_contribution, PetscReal *local_outflow_contribution)
	Perform post-step bookkeeping for a constant-velocity inlet boundary.
static PetscErrorCode	Destroy_InletConstantVelocity (BoundaryCondition *self)
	Release resources owned by a constant-velocity inlet boundary.
PetscErrorCode	Create_InletConstantVelocity (BoundaryCondition *bc)
	Implementation of Create_InletConstantVelocity().
static PetscErrorCode	Initialize_InletParabolicProfile (BoundaryCondition *self, BCContext *ctx)
	Initialize the geometric data used to evaluate a parabolic inlet profile.
static PetscErrorCode	PreStep_InletParabolicProfile (BoundaryCondition *self, BCContext *ctx, PetscReal *local_inflow_contribution, PetscReal *local_outflow_contribution)
	Refresh parabolic-inlet values required before the solver step.
static PetscErrorCode	Apply_InletParabolicProfile (BoundaryCondition *self, BCContext *ctx)
	Impose the evaluated parabolic velocity profile on inlet cells.
static PetscErrorCode	PostStep_InletParabolicProfile (BoundaryCondition *self, BCContext *ctx, PetscReal *local_inflow_contribution, PetscReal *local_outflow_contribution)
	Perform post-step bookkeeping for a parabolic inlet boundary.
static PetscErrorCode	Destroy_InletParabolicProfile (BoundaryCondition *self)
	Release resources owned by a parabolic inlet boundary.
PetscErrorCode	Create_InletParabolicProfile (BoundaryCondition *bc)
	Implementation of Create_InletParabolicProfile().
static PetscErrorCode	Initialize_InletProfileFromFile (BoundaryCondition *self, BCContext *ctx)
	Initializes a file-prescribed inlet profile handler for one boundary face.
static PetscErrorCode	PreStep_InletProfileFromFile (BoundaryCondition *self, BCContext *ctx, PetscReal *local_inflow_contribution, PetscReal *local_outflow_contribution)
	Pre-step hook for the static file-prescribed inlet profile handler.
static PetscErrorCode	Apply_InletProfileFromFile (BoundaryCondition *self, BCContext *ctx)
	Applies the loaded PICSLICE scalar profile to Ucont and Ubcs on an inlet face.
static PetscErrorCode	PostStep_InletProfileFromFile (BoundaryCondition *self, BCContext *ctx, PetscReal *local_inflow_contribution, PetscReal *local_outflow_contribution)
	Accumulates the applied inlet flux for a file-prescribed profile.
static PetscErrorCode	Destroy_InletProfileFromFile (BoundaryCondition *self)
	Releases private storage owned by a file-prescribed inlet profile handler.
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.
static PetscErrorCode	GetProfileFileExpectedDims (UserCtx *user, BCFace face_id, PetscInt *n1, PetscInt *n2)
	Computes the expected PICSLICE dimensions for an inlet face.
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.
static PetscReal	ProfileSpeedAt (const InletProfileFileData *data, PetscInt a, PetscInt b)
	Returns one scalar speed from the flattened PICSLICE profile.
PetscErrorCode	Create_InletProfileFromFile (BoundaryCondition *bc)
	Implementation of Create_InletProfileFromFile().
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.
static PetscErrorCode	Apply_OutletConservation (BoundaryCondition *self, BCContext *ctx)
	(Handler Action) Applies mass conservation correction to the outlet face.
static PetscErrorCode	PostStep_OutletConservation (BoundaryCondition *self, BCContext *ctx, PetscReal *local_inflow_contribution, PetscReal *local_outflow_contribution)
	Update outlet-conservation state after a completed solver step.
PetscErrorCode	Create_OutletConservation (BoundaryCondition *bc)
	Implementation of Create_OutletConservation().
PetscErrorCode	Create_PeriodicGeometric (BoundaryCondition *bc)
	Implementation of Create_PeriodicGeometric().
static PetscErrorCode	MeasureDrivenFluxes (UserCtx *user, char direction, PetscReal *boundaryFlux, PetscReal *planarAverageFlux)
	Measure the two volumetric fluxes the driven-flow controller senses.
static PetscErrorCode	Initialize_PeriodicDrivenConstant (BoundaryCondition *self, BCContext *ctx)
	Initialize constant forcing data for a periodically driven boundary.
static PetscErrorCode	PreStep_PeriodicDrivenConstant (BoundaryCondition *self, BCContext *ctx, PetscReal *local_inflow_contribution, PetscReal *local_outflow_contribution)
	Prepare constant periodic-driving data before the solver step.
static PetscErrorCode	Apply_PeriodicDrivenConstant (BoundaryCondition *self, BCContext *ctx)
	Apply the configured constant driving term to the periodic boundary.
static PetscErrorCode	Destroy_PeriodicDrivenConstant (BoundaryCondition *self)
	Release resources owned by the constant periodic-driving boundary.
PetscErrorCode	Create_PeriodicDrivenConstant (BoundaryCondition *bc)
	Internal helper implementation: Create_PeriodicDrivenConstant().
static PetscErrorCode	Initialize_PeriodicDrivenInitial (BoundaryCondition *self, BCContext *ctx)
	Initialize forcing data for a periodic boundary driven to its initial flux.
static PetscErrorCode	PreStep_PeriodicDrivenInitial (BoundaryCondition *self, BCContext *ctx, PetscReal *local_inflow_contribution, PetscReal *local_outflow_contribution)
	Latch the initial-state flux once, then drive to it like a constant target.
PetscErrorCode	Create_PeriodicDrivenInitial (BoundaryCondition *bc)
	Internal helper implementation: Create_PeriodicDrivenInitial().

---

Data Structure Documentation

InletConstantData

struct InletConstantData

Private data structure for the Constant Velocity Inlet handler.

Definition at line 310 of file BC_Handlers.c.

Data Fields
PetscReal	normal_velocity

InletParabolicData

struct InletParabolicData

Private data structure for the Parabolic Velocity Inlet handler.

Stores the peak velocity and pre-computed cross-stream geometry needed to evaluate the parabolic profile at each boundary node.

Definition at line 699 of file BC_Handlers.c.

Data Fields
PetscReal	v_max	Peak centerline velocity (from user params).
PetscReal	cs1_center	Center index in cross-stream direction 1.
PetscReal	cs2_center	Center index in cross-stream direction 2.
PetscReal	cs1_half	Half-width (in index space) in cross-stream direction 1.
PetscReal	cs2_half	Half-width (in index space) in cross-stream direction 2.

InletProfileFileData

struct InletProfileFileData

Definition at line 1106 of file BC_Handlers.c.

Data Fields
PetscInt	n1
PetscInt	n2
PetscReal *	profile
PetscReal	min_speed
PetscReal	max_speed
char *	source_file

DrivenFluxData

struct DrivenFluxData

Private data structure shared by both periodic driven-flux handlers.

constant_flux fills targetVolumetricFlux from the bcs file at initialization; initial_flux latches it from the starting field at the first PreStep. Everything downstream of the target is identical, so both handlers reuse this struct and the PreStep/Apply/Destroy implementations below.

Definition at line 2308 of file BC_Handlers.c.

Data Fields
char	direction
PetscReal	targetVolumetricFlux
PetscBool	isMasterController
PetscBool	enforceSeamFlux
PetscInt	lastBulkCorrectionStep

Macro Definition Documentation

__FUNCT__ [1/34]

#define __FUNCT__   "Validate_DrivenFlowConfiguration"

Definition at line 10 of file BC_Handlers.c.

__FUNCT__ [2/34]

#define __FUNCT__   "Create_WallNoSlip"

Definition at line 10 of file BC_Handlers.c.

__FUNCT__ [3/34]

#define __FUNCT__   "Apply_WallNoSlip"

Definition at line 10 of file BC_Handlers.c.

__FUNCT__ [4/34]

#define __FUNCT__   "Create_InletConstantVelocity"

Definition at line 10 of file BC_Handlers.c.

__FUNCT__ [5/34]

#define __FUNCT__   "Initialize_InletConstantVelocity"

Definition at line 10 of file BC_Handlers.c.

__FUNCT__ [6/34]

#define __FUNCT__   "PreStep_InletConstantVelocity"

Definition at line 10 of file BC_Handlers.c.

__FUNCT__ [7/34]

#define __FUNCT__   "Apply_InletConstantVelocity"

Definition at line 10 of file BC_Handlers.c.

__FUNCT__ [8/34]

#define __FUNCT__   "PostStep_InletConstantVelocity"

Definition at line 10 of file BC_Handlers.c.

__FUNCT__ [9/34]

#define __FUNCT__   "Destroy_InletConstantVelocity"

Definition at line 10 of file BC_Handlers.c.

__FUNCT__ [10/34]

#define __FUNCT__   "Create_InletParabolicProfile"

Definition at line 10 of file BC_Handlers.c.

__FUNCT__ [11/34]

#define __FUNCT__   "Initialize_InletParabolicProfile"

Definition at line 10 of file BC_Handlers.c.

__FUNCT__ [12/34]

#define __FUNCT__   "PreStep_InletParabolicProfile"

Definition at line 10 of file BC_Handlers.c.

__FUNCT__ [13/34]

#define __FUNCT__   "Apply_InletParabolicProfile"

Definition at line 10 of file BC_Handlers.c.

__FUNCT__ [14/34]

#define __FUNCT__   "PostStep_InletParabolicProfile"

Definition at line 10 of file BC_Handlers.c.

__FUNCT__ [15/34]

#define __FUNCT__   "Destroy_InletParabolicProfile"

Definition at line 10 of file BC_Handlers.c.

__FUNCT__ [16/34]

#define __FUNCT__   "Create_InletProfileFromFile"

Definition at line 10 of file BC_Handlers.c.

__FUNCT__ [17/34]

#define __FUNCT__   "Initialize_InletProfileFromFile"

Definition at line 10 of file BC_Handlers.c.

__FUNCT__ [18/34]

#define __FUNCT__   "PreStep_InletProfileFromFile"

Definition at line 10 of file BC_Handlers.c.

__FUNCT__ [19/34]

#define __FUNCT__   "Apply_InletProfileFromFile"

Definition at line 10 of file BC_Handlers.c.

__FUNCT__ [20/34]

#define __FUNCT__   "PostStep_InletProfileFromFile"

Definition at line 10 of file BC_Handlers.c.

__FUNCT__ [21/34]

#define __FUNCT__   "Destroy_InletProfileFromFile"

Definition at line 10 of file BC_Handlers.c.

__FUNCT__ [22/34]

#define __FUNCT__   "Create_OutletConservation"

Definition at line 10 of file BC_Handlers.c.

__FUNCT__ [23/34]

#define __FUNCT__   "PreStep_OutletConservation"

Definition at line 10 of file BC_Handlers.c.

__FUNCT__ [24/34]

#define __FUNCT__   "Apply_OutletConservation"

Definition at line 10 of file BC_Handlers.c.

__FUNCT__ [25/34]

#define __FUNCT__   "PostStep_OutletConservation"

Definition at line 10 of file BC_Handlers.c.

__FUNCT__ [26/34]

#define __FUNCT__   "MeasureDrivenFluxes"

Definition at line 10 of file BC_Handlers.c.

__FUNCT__ [27/34]

#define __FUNCT__   "Create_PeriodicDrivenConstant"

Definition at line 10 of file BC_Handlers.c.

__FUNCT__ [28/34]

#define __FUNCT__   "Initialize_PeriodicDrivenConstant"

Definition at line 10 of file BC_Handlers.c.

__FUNCT__ [29/34]

#define __FUNCT__   "PreStep_PeriodicDrivenConstant"

Definition at line 10 of file BC_Handlers.c.

__FUNCT__ [30/34]

#define __FUNCT__   "Apply_PeriodicDrivenConstant"

Definition at line 10 of file BC_Handlers.c.

__FUNCT__ [31/34]

#define __FUNCT__   "Destroy_PeriodicDrivenConstant"

Definition at line 10 of file BC_Handlers.c.

__FUNCT__ [32/34]

#define __FUNCT__   "Create_PeriodicDrivenInitial"

Definition at line 10 of file BC_Handlers.c.

__FUNCT__ [33/34]

#define __FUNCT__   "Initialize_PeriodicDrivenInitial"

Definition at line 10 of file BC_Handlers.c.

__FUNCT__ [34/34]

#define __FUNCT__   "PreStep_PeriodicDrivenInitial"

Definition at line 10 of file BC_Handlers.c.

Function Documentation

Validate_DrivenFlowConfiguration()

PetscErrorCode Validate_DrivenFlowConfiguration

(

UserCtx *

user

)

Internal helper implementation: Validate_DrivenFlowConfiguration().

(Private) Validates all consistency rules for a driven flow (channel/pipe) setup.

Local to this translation unit.

Definition at line 15 of file BC_Handlers.c.

{
    PetscFunctionBeginUser;
 
    // --- CHECK 1: Detect if a driven flow is active. ---
    PetscBool is_driven_flow_active = PETSC_FALSE;
    char      driven_direction = ' ';
    const char* first_driven_face_name = "";
 
    for (int i = 0; i < 6; i++) {
        BCHandlerType handler_type = user->boundary_faces[i].handler_type;
        if (handler_type == BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX ||
            handler_type == BC_HANDLER_PERIODIC_DRIVEN_INITIAL_FLUX)
        {
            is_driven_flow_active = PETSC_TRUE;
            first_driven_face_name = BCFaceToString((BCFace)i);
            
            if (i <= 1)      driven_direction = 'X';
            else if (i <= 3) driven_direction = 'Y';
            else             driven_direction = 'Z';
            
            break; // Exit loop once we've confirmed it's active and found the direction.
        }
    }
 
    // If no driven flow handler is found, validation for this rule set is complete.
    if (!is_driven_flow_active) {
        PetscFunctionReturn(0);
    }
    
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "  - Driven Flow Handler detected on face %s. Applying driven flow validation rules...\n", first_driven_face_name);
 
    // --- CHECK 2: Ensure no conflicting BCs (Inlet/Outlet/Far-field) are present. ---
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "  - Checking for incompatible Inlet/Outlet/Far-field BCs...\n");
    for (int i = 0; i < 6; i++) {
        BCType math_type = user->boundary_faces[i].mathematical_type;
        if (math_type == INLET || math_type == OUTLET || math_type == FARFIELD) {
             SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT,
                     "Configuration Error: A DRIVEN flow handler is active, which is incompatible with the %s boundary condition found on face %s.",
                     BCTypeToString(math_type), BCFaceToString((BCFace)i));
        }
    }
    LOG_ALLOW(GLOBAL, LOG_DEBUG, " ... No conflicting BC types found. OK.\n");
 
    // --- CHECK 3: Ensure both ends of the driven direction have identical, valid setups. ---
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "      - Validating symmetry and mathematical types for the '%c' direction...\n", driven_direction);
    
    PetscInt neg_face_idx = 0, pos_face_idx = 0;
    if (driven_direction == 'X') {
        neg_face_idx = BC_FACE_NEG_X; pos_face_idx = BC_FACE_POS_X;
    } else if (driven_direction == 'Y') {
        neg_face_idx = BC_FACE_NEG_Y; pos_face_idx = BC_FACE_POS_Y;
    } else { // 'Z'
        neg_face_idx = BC_FACE_NEG_Z; pos_face_idx = BC_FACE_POS_Z;
    }
 
    BoundaryFaceConfig *neg_face_cfg = &user->boundary_faces[neg_face_idx];
    BoundaryFaceConfig *pos_face_cfg = &user->boundary_faces[pos_face_idx];
 
    // Rule 3a: Both faces must be PERIODIC.
    if (neg_face_cfg->mathematical_type != PERIODIC || pos_face_cfg->mathematical_type != PERIODIC) {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT,
                "Configuration Error: For a driven flow in the '%c' direction, both the %s and %s faces must be of mathematical_type PERIODIC.",
                driven_direction, BCFaceToString((BCFace)neg_face_idx), BCFaceToString((BCFace)pos_face_idx));
    }
 
    // Rule 3b: Both faces must use the exact same handler type.
    if (neg_face_cfg->handler_type != pos_face_cfg->handler_type) {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT,
                "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).",
                BCFaceToString((BCFace)neg_face_idx), BCFaceToString((BCFace)pos_face_idx), driven_direction);
    }
    
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "          ... Symmetry and mathematical types are valid. OK.\n");
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

Apply_WallNoSlip()

static PetscErrorCode Apply_WallNoSlip ( BoundaryCondition *  self, BCContext *  ctx  )

static

Apply no-slip velocity values to wall-adjacent cells for this boundary condition.

Definition at line 143 of file BC_Handlers.c.

{
    PetscErrorCode ierr;
    UserCtx*       user = ctx->user;
    BCFace         face_id = ctx->face_id;
    PetscBool      can_service;
    
    (void)self;  // Unused for simple handlers
    
    PetscFunctionBeginUser;
    DMDALocalInfo *info = &user->info;
    Cmpnts ***ubcs, ***ucont;
    PetscInt IM_nodes_global, JM_nodes_global,KM_nodes_global;
 
    IM_nodes_global = user->IM;
    JM_nodes_global = user->JM;
    KM_nodes_global = user->KM;
    
    ierr = CanRankServiceFace(info,IM_nodes_global,JM_nodes_global,KM_nodes_global,face_id,&can_service); CHKERRQ(ierr);
    // Check if this rank owns part of this boundary face
    if (!can_service) PetscFunctionReturn(0);
 
    LOG_ALLOW(LOCAL, LOG_DEBUG, "Apply_WallNoSlip: Applying to Face %d (%s).\n",
              face_id, BCFaceToString(face_id));
 
    // Get arrays
    
    ierr = DMDAVecGetArray(user->fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
    ierr = DMDAVecGetArray(user->fda, user->Ucont, &ucont); CHKERRQ(ierr);
 
    PetscInt xs = info->xs, xe = info->xs + info->xm;
    PetscInt ys = info->ys, ye = info->ys + info->ym;
    PetscInt zs = info->zs, ze = info->zs + info->zm;
    PetscInt mx = info->mx, my = info->my, mz = info->mz;
    
    // ✅ Use shrunken loop bounds (avoids edges/corners like inlet handler)
    PetscInt lxs = xs, lxe = xe, lys = ys, lye = ye, lzs = zs, lze = ze;
    if (xs == 0) lxs = xs + 1;
    if (xe == mx) lxe = xe - 1;
    if (ys == 0) lys = ys + 1;
    if (ye == my) lye = ye - 1;
    if (zs == 0) lzs = zs + 1;
    if (ze == mz) lze = ze - 1;
 
    switch (face_id) {
        case BC_FACE_NEG_X: {
            if (xs == 0){  
                PetscInt i = xs;
                for (PetscInt k = lzs; k < lze; k++) {
                    for (PetscInt j = lys; j < lye; j++) {
                        // ✅ Set contravariant flux to zero (no penetration)
                        ucont[k][j][i].x = 0.0;
                    
                        // ✅ Set boundary velocity to zero (no slip)
                        ubcs[k][j][i].x = 0.0;
                        ubcs[k][j][i].y = 0.0;
                        ubcs[k][j][i].z = 0.0;
                        }                
                    }
                }
            break;
        }
 
        case BC_FACE_POS_X: {
            if (xe == mx){            
                PetscInt i = xe - 1;
                for (PetscInt k = lzs; k < lze; k++) {
                for (PetscInt j = lys; j < lye; j++) {
                    ucont[k][j][i-1].x = 0.0;
                    
                    ubcs[k][j][i].x = 0.0;
                    ubcs[k][j][i].y = 0.0;
                    ubcs[k][j][i].z = 0.0;
                    }
                }    
            }
        break;
        }    
        case BC_FACE_NEG_Y: {
            if (ys == 0){
            PetscInt j = ys;
            for (PetscInt k = lzs; k < lze; k++) {
                for (PetscInt i = lxs; i < lxe; i++) {
                    ucont[k][j][i].y = 0.0;
                    
                    ubcs[k][j][i].x = 0.0;
                    ubcs[k][j][i].y = 0.0;
                    ubcs[k][j][i].z = 0.0;
                    }
                }
            }
        } break;
 
        case BC_FACE_POS_Y: {
            if (ye == my){
            PetscInt j = ye - 1;
            for (PetscInt k = lzs; k < lze; k++) {
                for (PetscInt i = lxs; i < lxe; i++) {
                    ucont[k][j-1][i].y = 0.0;
                    
                    ubcs[k][j][i].x = 0.0;
                    ubcs[k][j][i].y = 0.0;
                    ubcs[k][j][i].z = 0.0;
                    }
                }
            }
        } break;
            
        case BC_FACE_NEG_Z: {
            if (zs == 0){
            PetscInt k = zs;
            for (PetscInt j = lys; j < lye; j++) {
                for (PetscInt i = lxs; i < lxe; i++) {
                    ucont[k][j][i].z = 0.0;
                    
                    ubcs[k][j][i].x = 0.0;
                    ubcs[k][j][i].y = 0.0;
                    ubcs[k][j][i].z = 0.0;
                    }
                }
            }
        } break;
 
        case BC_FACE_POS_Z: {
            if (ze == mz) {
            PetscInt k = ze - 1;
            for (PetscInt j = lys; j < lye; j++) {
                for (PetscInt i = lxs; i < lxe; i++) {
                    ucont[k-1][j][i].z = 0.0;
                    
                    ubcs[k][j][i].x = 0.0;
                    ubcs[k][j][i].y = 0.0;
                    ubcs[k][j][i].z = 0.0;
                    }
                }
            }    
        } break;
    }
 
    // Restore arrays
    ierr = DMDAVecRestoreArray(user->fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(user->fda, user->Ucont, &ucont); CHKERRQ(ierr);
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

Create_WallNoSlip()

PetscErrorCode Create_WallNoSlip

(

BoundaryCondition *

bc

)

Implementation of Create_WallNoSlip().

Configures a BoundaryCondition object to behave as a no-slip, stationary wall.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/BC_Handlers.h.

See also

Create_WallNoSlip()

Definition at line 114 of file BC_Handlers.c.

{
    PetscFunctionBeginUser;
    
    if (!bc) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, 
                     "Input BoundaryCondition object is NULL in Create_WallNoSlip");
 
    // ✅ Set priority
    bc->priority   = BC_PRIORITY_WALL;
    
    // Assign function pointers
    bc->Initialize = NULL;
    bc->PreStep    = NULL;
    bc->Apply      = Apply_WallNoSlip;
    bc->PostStep   = NULL;
    bc->UpdateUbcs = NULL;
    bc->Destroy    = NULL;
    
    // No private data needed for this simple handler
    bc->data = NULL;
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

Initialize_InletConstantVelocity()

static PetscErrorCode Initialize_InletConstantVelocity ( BoundaryCondition *  self, BCContext *  ctx  )

static

Initialize persistent state for a constant-velocity inlet boundary.

Definition at line 350 of file BC_Handlers.c.

{
    PetscErrorCode ierr;
    UserCtx*       user = ctx->user;
    BCFace         face_id = ctx->face_id;
    InletConstantData *data = (InletConstantData*)self->data;
    PetscBool      found;
    
    PetscFunctionBeginUser;
    LOG_ALLOW(LOCAL, LOG_DEBUG, "Initialize_InletConstantVelocity: Initializing handler for Face %d. \n", face_id);
    data->normal_velocity = 0.0;
    
    switch (face_id) {
        case BC_FACE_NEG_X:
        case BC_FACE_POS_X:
            // For X-faces, read "vx" as normal velocity
            ierr = GetBCParamReal(user->boundary_faces[face_id].params, "vx", 
                                 &data->normal_velocity, &found); CHKERRQ(ierr);
            break;
            
        case BC_FACE_NEG_Y:
        case BC_FACE_POS_Y:
            // For Y-faces, read "vy" as normal velocity
            ierr = GetBCParamReal(user->boundary_faces[face_id].params, "vy", 
                                 &data->normal_velocity, &found); CHKERRQ(ierr);
            break;
            
        case BC_FACE_NEG_Z:
        case BC_FACE_POS_Z:
            // For Z-faces, read "vz" as normal velocity
            ierr = GetBCParamReal(user->boundary_faces[face_id].params, "vz", 
                                 &data->normal_velocity, &found); CHKERRQ(ierr);
            break;
    }
    
    LOG_ALLOW(LOCAL, LOG_INFO, "  Inlet Face %d: normal velocity = %.4f\n",
              face_id, data->normal_velocity);
 
    // Set initial boundary state
    ierr = Apply_InletConstantVelocity(self, ctx); CHKERRQ(ierr);
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

PreStep_InletConstantVelocity()

static PetscErrorCode PreStep_InletConstantVelocity ( BoundaryCondition *  self, BCContext *  ctx, PetscReal *  in, PetscReal *  out  )

static

Update constant-inlet data required before the next solver step.

Definition at line 399 of file BC_Handlers.c.

{
    // No preparation needed for constant velocity inlet.
    // The velocity is already stored in self->data from Initialize.
    // Apply will set ucont, and PostStep will measure the actual flux.
    
    (void)self;
    (void)ctx;
    (void)local_inflow_contribution;
    (void)local_outflow_contribution;
    
    PetscFunctionBeginUser;
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

Apply_InletConstantVelocity()

static PetscErrorCode Apply_InletConstantVelocity ( BoundaryCondition *  self, BCContext *  ctx  )

static

Impose the configured constant velocity on inlet boundary cells.

Definition at line 421 of file BC_Handlers.c.

{
    PetscErrorCode ierr;
    UserCtx*       user = ctx->user;
    BCFace         face_id = ctx->face_id;
    InletConstantData *data = (InletConstantData*)self->data;
    PetscBool      can_service;
    
    PetscFunctionBeginUser;
    
    DMDALocalInfo *info = &user->info;
    Cmpnts ***ubcs, ***ucont, ***csi, ***eta, ***zet;
    PetscReal ***nvert;
    PetscInt IM_nodes_global, JM_nodes_global,KM_nodes_global;
 
    IM_nodes_global = user->IM;
    JM_nodes_global = user->JM;
    KM_nodes_global = user->KM;
    
    ierr = CanRankServiceFace(info,IM_nodes_global,JM_nodes_global,KM_nodes_global,face_id,&can_service); CHKERRQ(ierr);
    
    if (!can_service) PetscFunctionReturn(0);
 
    // Get arrays
 
    ierr = DMDAVecGetArray(user->fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
    ierr = DMDAVecGetArray(user->fda, user->Ucont, &ucont); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
 
    // Get SCALAR velocity (not vector!)
    PetscReal uin_this_point = data->normal_velocity;
    
    PetscInt xs = info->xs, xe = info->xs + info->xm;
    PetscInt ys = info->ys, ye = info->ys + info->ym;
    PetscInt zs = info->zs, ze = info->zs + info->zm;
    PetscInt mx = info->mx, my = info->my, mz = info->mz;
    
    PetscInt lxs = xs, lxe = xe, lys = ys, lye = ye, lzs = zs, lze = ze;
    if (xs == 0) lxs = xs + 1;
    if (xe == mx) lxe = xe - 1;
    if (ys == 0) lys = ys + 1;
    if (ye == my) lye = ye - 1;
    if (zs == 0) lzs = zs + 1;
    if (ze == mz) lze = ze - 1;
 
    switch (face_id) {
        case BC_FACE_NEG_X:
        case BC_FACE_POS_X: {
            PetscReal sign = (face_id == BC_FACE_NEG_X) ? 1.0 : -1.0;
            PetscInt i = (face_id == BC_FACE_NEG_X) ? xs : mx - 2;
            
            for (PetscInt k = lzs; k < lze; k++) {
                for (PetscInt j = lys; j < lye; j++) {
                    if ((sign > 0 && nvert[k][j][i+1] > 0.1) || 
                        (sign < 0 && nvert[k][j][i] > 0.1)) continue;
                    
                    PetscReal CellArea = sqrt(csi[k][j][i].x * csi[k][j][i].x + 
                                             csi[k][j][i].y * csi[k][j][i].y + 
                                             csi[k][j][i].z * csi[k][j][i].z);
                    
                    ucont[k][j][i].x = sign * uin_this_point * CellArea;
                    
                    ubcs[k][j][i + (sign < 0)].x = sign * uin_this_point * csi[k][j][i].x / CellArea;
                    ubcs[k][j][i + (sign < 0)].y = sign * uin_this_point * csi[k][j][i].y / CellArea;
                    ubcs[k][j][i + (sign < 0)].z = sign * uin_this_point * csi[k][j][i].z / CellArea;
                }
            }
        } break;
            
        case BC_FACE_NEG_Y:
        case BC_FACE_POS_Y: {
            PetscReal sign = (face_id == BC_FACE_NEG_Y) ? 1.0 : -1.0;
            PetscInt j = (face_id == BC_FACE_NEG_Y) ? ys : my - 2;
            
            for (PetscInt k = lzs; k < lze; k++) {
                for (PetscInt i = lxs; i < lxe; i++) {
                    if ((sign > 0 && nvert[k][j+1][i] > 0.1) || 
                        (sign < 0 && nvert[k][j][i] > 0.1)) continue;
                    
                    PetscReal CellArea = sqrt(eta[k][j][i].x * eta[k][j][i].x + 
                                             eta[k][j][i].y * eta[k][j][i].y + 
                                             eta[k][j][i].z * eta[k][j][i].z);
                    
                    ucont[k][j][i].y = sign * uin_this_point * CellArea;
                    
                    ubcs[k][j + (sign < 0)][i].x = sign * uin_this_point * eta[k][j][i].x / CellArea;
                    ubcs[k][j + (sign < 0)][i].y = sign * uin_this_point * eta[k][j][i].y / CellArea;
                    ubcs[k][j + (sign < 0)][i].z = sign * uin_this_point * eta[k][j][i].z / CellArea;
                }
            }
        } break;
            
        case BC_FACE_NEG_Z:
        case BC_FACE_POS_Z: {
            PetscReal sign = (face_id == BC_FACE_NEG_Z) ? 1.0 : -1.0;
            PetscInt k = (face_id == BC_FACE_NEG_Z) ? zs : mz - 2;
            
            for (PetscInt j = lys; j < lye; j++) {
                for (PetscInt i = lxs; i < lxe; i++) {
                    if ((sign > 0 && nvert[k+1][j][i] > 0.1) || 
                        (sign < 0 && nvert[k][j][i] > 0.1)) continue;
                    
                    PetscReal CellArea = sqrt(zet[k][j][i].x * zet[k][j][i].x + 
                                             zet[k][j][i].y * zet[k][j][i].y + 
                                             zet[k][j][i].z * zet[k][j][i].z);
                    
                    ucont[k][j][i].z = sign * uin_this_point * CellArea;
                    
                    ubcs[k + (sign < 0)][j][i].x = sign * uin_this_point * zet[k][j][i].x / CellArea;
                    ubcs[k + (sign < 0)][j][i].y = sign * uin_this_point * zet[k][j][i].y / CellArea;
                    ubcs[k + (sign < 0)][j][i].z = sign * uin_this_point * zet[k][j][i].z / CellArea;
                }
            }
        } break;
    }
 
    // Restore arrays
    ierr = DMDAVecRestoreArray(user->fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(user->fda, user->Ucont, &ucont); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

PostStep_InletConstantVelocity()

static PetscErrorCode PostStep_InletConstantVelocity ( BoundaryCondition *  self, BCContext *  ctx, PetscReal *  in, PetscReal *  out  )

static

Perform post-step bookkeeping for a constant-velocity inlet boundary.

Definition at line 556 of file BC_Handlers.c.

{
    PetscErrorCode ierr;
    UserCtx*       user = ctx->user;
    BCFace         face_id = ctx->face_id;
    PetscBool      can_service;
    
    (void)self;
    (void)local_outflow_contribution;
    
    PetscFunctionBeginUser;
 
    DMDALocalInfo *info = &user->info;
    Cmpnts ***ucont;
 
    PetscInt IM_nodes_global, JM_nodes_global,KM_nodes_global;
 
    IM_nodes_global = user->IM;
    JM_nodes_global = user->JM;
    KM_nodes_global = user->KM;
    
    ierr = CanRankServiceFace(info,IM_nodes_global,JM_nodes_global,KM_nodes_global,face_id,&can_service); CHKERRQ(ierr);
    
 
    if (!can_service) PetscFunctionReturn(0);
    
    ierr = DMDAVecGetArrayRead(user->fda, user->Ucont, (const Cmpnts***)&ucont); CHKERRQ(ierr);
 
    PetscReal local_flux = 0.0;
    
    PetscInt xs = info->xs, xe = info->xs + info->xm;
    PetscInt ys = info->ys, ye = info->ys + info->ym;
    PetscInt zs = info->zs, ze = info->zs + info->zm;
    PetscInt mx = info->mx, my = info->my, mz = info->mz;
    
    PetscInt lxs = xs, lxe = xe, lys = ys, lye = ye, lzs = zs, lze = ze;
    if (xs == 0) lxs = xs + 1;
    if (xe == mx) lxe = xe - 1;
    if (ys == 0) lys = ys + 1;
    if (ye == my) lye = ye - 1;
    if (zs == 0) lzs = zs + 1;
    if (ze == mz) lze = ze - 1;
 
    // Sum ucont components
    switch (face_id) {
        case BC_FACE_NEG_X:
        case BC_FACE_POS_X: {
            PetscInt i = (face_id == BC_FACE_NEG_X) ? xs : mx - 2;
            for (PetscInt k = lzs; k < lze; k++) {
                for (PetscInt j = lys; j < lye; j++) {
                    local_flux += ucont[k][j][i].x;
                }
            }
        } break;
        
        case BC_FACE_NEG_Y:
        case BC_FACE_POS_Y: {
            PetscInt j = (face_id == BC_FACE_NEG_Y) ? ys : my - 2;
            for (PetscInt k = lzs; k < lze; k++) {
                for (PetscInt i = lxs; i < lxe; i++) {
                    local_flux += ucont[k][j][i].y;
                }
            }
        } break;
        
        case BC_FACE_NEG_Z:
        case BC_FACE_POS_Z: {
            PetscInt k = (face_id == BC_FACE_NEG_Z) ? zs : mz - 2;
            for (PetscInt j = lys; j < lye; j++) {
                for (PetscInt i = lxs; i < lxe; i++) {
                    local_flux += ucont[k][j][i].z;
                }
            }
        } break;
    }
 
    ierr = DMDAVecRestoreArrayRead(user->fda, user->Ucont, (const Cmpnts***)&ucont); CHKERRQ(ierr);
 
    *local_inflow_contribution += local_flux;
    
    LOG_ALLOW(LOCAL, LOG_DEBUG, "PostStep_InletConstantVelocity: Face %d, flux = %.6e\n",
              face_id, local_flux);
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

Destroy_InletConstantVelocity()

static PetscErrorCode Destroy_InletConstantVelocity ( BoundaryCondition *  self )

static

Release resources owned by a constant-velocity inlet boundary.

Definition at line 650 of file BC_Handlers.c.

{
    PetscFunctionBeginUser;
    if (self && self->data) {
        PetscFree(self->data);
        self->data = NULL;
    }
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

Create_InletConstantVelocity()

PetscErrorCode Create_InletConstantVelocity

(

BoundaryCondition *

bc

)

Implementation of Create_InletConstantVelocity().

Configures a BoundaryCondition object to behave as a constant velocity inlet.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/BC_Handlers.h.

See also

Create_InletConstantVelocity()

Definition at line 322 of file BC_Handlers.c.

{
    PetscErrorCode ierr;
    PetscFunctionBeginUser;
    
    if (!bc) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "BoundaryCondition is NULL");
 
    InletConstantData *data = NULL;
    ierr = PetscMalloc1(1, &data); CHKERRQ(ierr);
    bc->data = (void*)data;
    
    bc->priority   = BC_PRIORITY_INLET;
    bc->Initialize = Initialize_InletConstantVelocity;
    bc->PreStep    = PreStep_InletConstantVelocity;
    bc->Apply      = Apply_InletConstantVelocity;
    bc->PostStep   = PostStep_InletConstantVelocity;
    bc->UpdateUbcs = NULL;
    bc->Destroy    = Destroy_InletConstantVelocity;
    
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

Initialize_InletParabolicProfile()

static PetscErrorCode Initialize_InletParabolicProfile ( BoundaryCondition *  self, BCContext *  ctx  )

static

Initialize the geometric data used to evaluate a parabolic inlet profile.

Definition at line 743 of file BC_Handlers.c.

{
    PetscErrorCode ierr;
    UserCtx*       user = ctx->user;
    BCFace         face_id = ctx->face_id;
    InletParabolicData *data = (InletParabolicData*)self->data;
    PetscBool      found;
 
    PetscFunctionBeginUser;
    LOG_ALLOW(LOCAL, LOG_DEBUG, "Initialize_InletParabolicProfile: Initializing handler for Face %d.\n", face_id);
 
    // --- Parse v_max from boundary condition parameters ---
    data->v_max = 0.0;
    ierr = GetBCParamReal(user->boundary_faces[face_id].params, "v_max",
                          &data->v_max, &found); CHKERRQ(ierr);
    if (!found) {
        LOG_ALLOW(GLOBAL, LOG_WARNING, "Initialize_InletParabolicProfile: 'v_max' not found in params for face %d. Defaulting to 0.0.\n", face_id);
    }
 
    // --- Determine cross-stream dimensions based on face orientation ---
    PetscReal cs1_dim, cs2_dim;
    switch (face_id) {
        case BC_FACE_NEG_X:
        case BC_FACE_POS_X:
            cs1_dim = (PetscReal)user->JM;  // j-direction
            cs2_dim = (PetscReal)user->KM;  // k-direction
            break;
        case BC_FACE_NEG_Y:
        case BC_FACE_POS_Y:
            cs1_dim = (PetscReal)user->IM;  // i-direction
            cs2_dim = (PetscReal)user->KM;  // k-direction
            break;
        case BC_FACE_NEG_Z:
        case BC_FACE_POS_Z:
        default:
            cs1_dim = (PetscReal)user->IM;  // i-direction
            cs2_dim = (PetscReal)user->JM;  // j-direction
            break;
    }
 
    // Interior width = dim - 2 (nodes 1 through dim-2 are interior)
    PetscReal cs1_width = cs1_dim - 2.0;
    PetscReal cs2_width = cs2_dim - 2.0;
 
    data->cs1_center = 1.0 + cs1_width / 2.0;
    data->cs2_center = 1.0 + cs2_width / 2.0;
    data->cs1_half   = cs1_width / 2.0;
    data->cs2_half   = cs2_width / 2.0;
 
    LOG_ALLOW(LOCAL, LOG_INFO, "  Inlet Face %d (Parabolic): v_max = %.4f\n", face_id, data->v_max);
    LOG_ALLOW(LOCAL, LOG_DEBUG, "    Cross-stream 1: center=%.1f, half=%.1f\n", data->cs1_center, data->cs1_half);
    LOG_ALLOW(LOCAL, LOG_DEBUG, "    Cross-stream 2: center=%.1f, half=%.1f\n", data->cs2_center, data->cs2_half);
 
    // Set initial boundary state
    ierr = Apply_InletParabolicProfile(self, ctx); CHKERRQ(ierr);
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

PreStep_InletParabolicProfile()

static PetscErrorCode PreStep_InletParabolicProfile ( BoundaryCondition *  self, BCContext *  ctx, PetscReal *  in, PetscReal *  out  )

static

Refresh parabolic-inlet values required before the solver step.

Definition at line 808 of file BC_Handlers.c.

{
    (void)self;
    (void)ctx;
    (void)local_inflow_contribution;
    (void)local_outflow_contribution;
 
    PetscFunctionBeginUser;
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

Apply_InletParabolicProfile()

static PetscErrorCode Apply_InletParabolicProfile ( BoundaryCondition *  self, BCContext *  ctx  )

static

Impose the evaluated parabolic velocity profile on inlet cells.

Definition at line 827 of file BC_Handlers.c.

{
    PetscErrorCode ierr;
    UserCtx*       user = ctx->user;
    BCFace         face_id = ctx->face_id;
    InletParabolicData *data = (InletParabolicData*)self->data;
    PetscBool      can_service;
 
    PetscFunctionBeginUser;
 
    DMDALocalInfo *info = &user->info;
    Cmpnts ***ubcs, ***ucont, ***csi, ***eta, ***zet;
    PetscReal ***nvert;
    PetscInt IM_nodes_global, JM_nodes_global, KM_nodes_global;
 
    IM_nodes_global = user->IM;
    JM_nodes_global = user->JM;
    KM_nodes_global = user->KM;
 
    ierr = CanRankServiceFace(info, IM_nodes_global, JM_nodes_global, KM_nodes_global,
                              face_id, &can_service); CHKERRQ(ierr);
 
    if (!can_service) PetscFunctionReturn(0);
 
    // Get arrays
    ierr = DMDAVecGetArray(user->fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
    ierr = DMDAVecGetArray(user->fda, user->Ucont, &ucont); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
 
    PetscInt xs = info->xs, xe = info->xs + info->xm;
    PetscInt ys = info->ys, ye = info->ys + info->ym;
    PetscInt zs = info->zs, ze = info->zs + info->zm;
    PetscInt mx = info->mx, my = info->my, mz = info->mz;
 
    PetscInt lxs = xs, lxe = xe, lys = ys, lye = ye, lzs = zs, lze = ze;
    if (xs == 0) lxs = xs + 1;
    if (xe == mx) lxe = xe - 1;
    if (ys == 0) lys = ys + 1;
    if (ye == my) lye = ye - 1;
    if (zs == 0) lzs = zs + 1;
    if (ze == mz) lze = ze - 1;
 
    switch (face_id) {
        case BC_FACE_NEG_X:
        case BC_FACE_POS_X: {
            // X-faces: normal = i, cross-stream = (j, k)
            PetscReal sign = (face_id == BC_FACE_NEG_X) ? 1.0 : -1.0;
            PetscInt i = (face_id == BC_FACE_NEG_X) ? xs : mx - 2;
 
            for (PetscInt k = lzs; k < lze; k++) {
                for (PetscInt j = lys; j < lye; j++) {
                    if ((sign > 0 && nvert[k][j][i+1] > 0.1) ||
                        (sign < 0 && nvert[k][j][i] > 0.1)) continue;
 
                    // Evaluate parabolic profile: cs1 = j, cs2 = k
                    PetscReal cs1_norm = ((PetscReal)j - data->cs1_center) / data->cs1_half;
                    PetscReal cs2_norm = ((PetscReal)k - data->cs2_center) / data->cs2_half;
                    PetscReal profile = PetscMax(0.0, 1.0 - cs1_norm * cs1_norm)
                                      * PetscMax(0.0, 1.0 - cs2_norm * cs2_norm);
                    PetscReal uin_local = data->v_max * profile;
 
                    PetscReal CellArea = sqrt(csi[k][j][i].x * csi[k][j][i].x +
                                             csi[k][j][i].y * csi[k][j][i].y +
                                             csi[k][j][i].z * csi[k][j][i].z);
 
                    ucont[k][j][i].x = sign * uin_local * CellArea;
 
                    ubcs[k][j][i + (sign < 0)].x = sign * uin_local * csi[k][j][i].x / CellArea;
                    ubcs[k][j][i + (sign < 0)].y = sign * uin_local * csi[k][j][i].y / CellArea;
                    ubcs[k][j][i + (sign < 0)].z = sign * uin_local * csi[k][j][i].z / CellArea;
                }
            }
        } break;
 
        case BC_FACE_NEG_Y:
        case BC_FACE_POS_Y: {
            // Y-faces: normal = j, cross-stream = (i, k)
            PetscReal sign = (face_id == BC_FACE_NEG_Y) ? 1.0 : -1.0;
            PetscInt j = (face_id == BC_FACE_NEG_Y) ? ys : my - 2;
 
            for (PetscInt k = lzs; k < lze; k++) {
                for (PetscInt i = lxs; i < lxe; i++) {
                    if ((sign > 0 && nvert[k][j+1][i] > 0.1) ||
                        (sign < 0 && nvert[k][j][i] > 0.1)) continue;
 
                    // Evaluate parabolic profile: cs1 = i, cs2 = k
                    PetscReal cs1_norm = ((PetscReal)i - data->cs1_center) / data->cs1_half;
                    PetscReal cs2_norm = ((PetscReal)k - data->cs2_center) / data->cs2_half;
                    PetscReal profile = PetscMax(0.0, 1.0 - cs1_norm * cs1_norm)
                                      * PetscMax(0.0, 1.0 - cs2_norm * cs2_norm);
                    PetscReal uin_local = data->v_max * profile;
 
                    PetscReal CellArea = sqrt(eta[k][j][i].x * eta[k][j][i].x +
                                             eta[k][j][i].y * eta[k][j][i].y +
                                             eta[k][j][i].z * eta[k][j][i].z);
 
                    ucont[k][j][i].y = sign * uin_local * CellArea;
 
                    ubcs[k][j + (sign < 0)][i].x = sign * uin_local * eta[k][j][i].x / CellArea;
                    ubcs[k][j + (sign < 0)][i].y = sign * uin_local * eta[k][j][i].y / CellArea;
                    ubcs[k][j + (sign < 0)][i].z = sign * uin_local * eta[k][j][i].z / CellArea;
                }
            }
        } break;
 
        case BC_FACE_NEG_Z:
        case BC_FACE_POS_Z: {
            // Z-faces: normal = k, cross-stream = (i, j)
            PetscReal sign = (face_id == BC_FACE_NEG_Z) ? 1.0 : -1.0;
            PetscInt k = (face_id == BC_FACE_NEG_Z) ? zs : mz - 2;
 
            for (PetscInt j = lys; j < lye; j++) {
                for (PetscInt i = lxs; i < lxe; i++) {
                    if ((sign > 0 && nvert[k+1][j][i] > 0.1) ||
                        (sign < 0 && nvert[k][j][i] > 0.1)) continue;
 
                    // Evaluate parabolic profile: cs1 = i, cs2 = j
                    PetscReal cs1_norm = ((PetscReal)i - data->cs1_center) / data->cs1_half;
                    PetscReal cs2_norm = ((PetscReal)j - data->cs2_center) / data->cs2_half;
                    PetscReal profile = PetscMax(0.0, 1.0 - cs1_norm * cs1_norm)
                                      * PetscMax(0.0, 1.0 - cs2_norm * cs2_norm);
                    PetscReal uin_local = data->v_max * profile;
 
                    PetscReal CellArea = sqrt(zet[k][j][i].x * zet[k][j][i].x +
                                             zet[k][j][i].y * zet[k][j][i].y +
                                             zet[k][j][i].z * zet[k][j][i].z);
 
                    ucont[k][j][i].z = sign * uin_local * CellArea;
 
                    ubcs[k + (sign < 0)][j][i].x = sign * uin_local * zet[k][j][i].x / CellArea;
                    ubcs[k + (sign < 0)][j][i].y = sign * uin_local * zet[k][j][i].y / CellArea;
                    ubcs[k + (sign < 0)][j][i].z = sign * uin_local * zet[k][j][i].z / CellArea;
                }
            }
        } break;
    }
 
    // Restore arrays
    ierr = DMDAVecRestoreArray(user->fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(user->fda, user->Ucont, &ucont); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

PostStep_InletParabolicProfile()

static PetscErrorCode PostStep_InletParabolicProfile ( BoundaryCondition *  self, BCContext *  ctx, PetscReal *  in, PetscReal *  out  )

static

Perform post-step bookkeeping for a parabolic inlet boundary.

Definition at line 984 of file BC_Handlers.c.

{
    PetscErrorCode ierr;
    UserCtx*       user = ctx->user;
    BCFace         face_id = ctx->face_id;
    PetscBool      can_service;
 
    (void)self;
    (void)local_outflow_contribution;
 
    PetscFunctionBeginUser;
 
    DMDALocalInfo *info = &user->info;
    Cmpnts ***ucont;
 
    PetscInt IM_nodes_global, JM_nodes_global, KM_nodes_global;
 
    IM_nodes_global = user->IM;
    JM_nodes_global = user->JM;
    KM_nodes_global = user->KM;
 
    ierr = CanRankServiceFace(info, IM_nodes_global, JM_nodes_global, KM_nodes_global,
                              face_id, &can_service); CHKERRQ(ierr);
 
    if (!can_service) PetscFunctionReturn(0);
 
    ierr = DMDAVecGetArrayRead(user->fda, user->Ucont, (const Cmpnts***)&ucont); CHKERRQ(ierr);
 
    PetscReal local_flux = 0.0;
 
    PetscInt xs = info->xs, xe = info->xs + info->xm;
    PetscInt ys = info->ys, ye = info->ys + info->ym;
    PetscInt zs = info->zs, ze = info->zs + info->zm;
    PetscInt mx = info->mx, my = info->my, mz = info->mz;
 
    PetscInt lxs = xs, lxe = xe, lys = ys, lye = ye, lzs = zs, lze = ze;
    if (xs == 0) lxs = xs + 1;
    if (xe == mx) lxe = xe - 1;
    if (ys == 0) lys = ys + 1;
    if (ye == my) lye = ye - 1;
    if (zs == 0) lzs = zs + 1;
    if (ze == mz) lze = ze - 1;
 
    switch (face_id) {
        case BC_FACE_NEG_X:
        case BC_FACE_POS_X: {
            PetscInt i = (face_id == BC_FACE_NEG_X) ? xs : mx - 2;
            for (PetscInt k = lzs; k < lze; k++) {
                for (PetscInt j = lys; j < lye; j++) {
                    local_flux += ucont[k][j][i].x;
                }
            }
        } break;
 
        case BC_FACE_NEG_Y:
        case BC_FACE_POS_Y: {
            PetscInt j = (face_id == BC_FACE_NEG_Y) ? ys : my - 2;
            for (PetscInt k = lzs; k < lze; k++) {
                for (PetscInt i = lxs; i < lxe; i++) {
                    local_flux += ucont[k][j][i].y;
                }
            }
        } break;
 
        case BC_FACE_NEG_Z:
        case BC_FACE_POS_Z: {
            PetscInt k = (face_id == BC_FACE_NEG_Z) ? zs : mz - 2;
            for (PetscInt j = lys; j < lye; j++) {
                for (PetscInt i = lxs; i < lxe; i++) {
                    local_flux += ucont[k][j][i].z;
                }
            }
        } break;
    }
 
    ierr = DMDAVecRestoreArrayRead(user->fda, user->Ucont, (const Cmpnts***)&ucont); CHKERRQ(ierr);
 
    *local_inflow_contribution += local_flux;
 
    LOG_ALLOW(LOCAL, LOG_DEBUG, "PostStep_InletParabolicProfile: Face %d, flux = %.6e\n",
              face_id, local_flux);
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

Destroy_InletParabolicProfile()

static PetscErrorCode Destroy_InletParabolicProfile ( BoundaryCondition *  self )

static

Release resources owned by a parabolic inlet boundary.

Definition at line 1077 of file BC_Handlers.c.

{
    PetscFunctionBeginUser;
    if (self && self->data) {
        PetscFree(self->data);
        self->data = NULL;
    }
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

Create_InletParabolicProfile()

PetscErrorCode Create_InletParabolicProfile

(

BoundaryCondition *

bc

)

Implementation of Create_InletParabolicProfile().

Configures a BoundaryCondition object for a parabolic inlet profile.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/BC_Handlers.h.

See also

Create_InletParabolicProfile()

Definition at line 715 of file BC_Handlers.c.

{
    PetscErrorCode ierr;
    PetscFunctionBeginUser;
 
    if (!bc) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "BoundaryCondition is NULL");
 
    InletParabolicData *data = NULL;
    ierr = PetscMalloc1(1, &data); CHKERRQ(ierr);
    bc->data = (void*)data;
 
    bc->priority   = BC_PRIORITY_INLET;
    bc->Initialize = Initialize_InletParabolicProfile;
    bc->PreStep    = PreStep_InletParabolicProfile;
    bc->Apply      = Apply_InletParabolicProfile;
    bc->PostStep   = PostStep_InletParabolicProfile;
    bc->UpdateUbcs = NULL;
    bc->Destroy    = Destroy_InletParabolicProfile;
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

Initialize_InletProfileFromFile()

static PetscErrorCode Initialize_InletProfileFromFile ( BoundaryCondition *  self, BCContext *  ctx  )

static

Initializes a file-prescribed inlet profile handler for one boundary face.

Reads the source_file BC parameter, validates the target face dimensions, loads the PICSLICE scalar speed profile, records summary statistics for logging, and applies the profile once to initialize boundary state.

Parameters

self	BoundaryCondition object configured by Create_InletProfileFromFile().
ctx	Runtime boundary context containing the UserCtx and face id.

Returns

PetscErrorCode 0 on success, or a PETSc error for missing parameters or malformed files.

Definition at line 1333 of file BC_Handlers.c.

{
    PetscErrorCode ierr;
    UserCtx *user = ctx->user;
    BCFace face_id = ctx->face_id;
    InletProfileFileData *data = (InletProfileFileData*)self->data;
    PetscBool found = PETSC_FALSE;
    const char *source_file = NULL;
    PetscInt expected_n1 = 0, expected_n2 = 0;
 
    PetscFunctionBeginUser;
    ierr = GetBCParamStringLocal(user->boundary_faces[face_id].params, "source_file",
                                 &source_file, &found); CHKERRQ(ierr);
    if (!found || !source_file || source_file[0] == '\0') {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
                "InletProfileFromFile requires source_file parameter for face %d.", face_id);
    }
 
    ierr = PetscStrallocpy(source_file, &data->source_file); CHKERRQ(ierr);
    ierr = GetProfileFileExpectedDims(user, face_id, &expected_n1, &expected_n2); CHKERRQ(ierr);
    ierr = ReadPicSliceProfile(source_file, expected_n1, expected_n2, data); CHKERRQ(ierr);
 
    LOG_ALLOW(LOCAL, LOG_INFO,
              "  Inlet Face %d (Prescribed Flow): source=%s dims=(%d,%d) speed[min,max]=[%.6e, %.6e]\n",
              face_id, data->source_file, data->n1, data->n2,
              (double)data->min_speed, (double)data->max_speed);
 
    ierr = Apply_InletProfileFromFile(self, ctx); CHKERRQ(ierr);
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

PreStep_InletProfileFromFile()

static PetscErrorCode PreStep_InletProfileFromFile ( BoundaryCondition *  self, BCContext *  ctx, PetscReal *  local_inflow_contribution, PetscReal *  local_outflow_contribution  )

static

Pre-step hook for the static file-prescribed inlet profile handler.

Static profiles require no per-step preparation. The hook is implemented so future time-varying profile support can reuse the same handler lifecycle.

Parameters

self	BoundaryCondition object for this inlet handler.
ctx	Runtime boundary context.
local_inflow_contribution	Inflow accumulator, intentionally unchanged.
local_outflow_contribution	Outflow accumulator, intentionally unchanged.

Returns

PetscErrorCode 0 on success.

Definition at line 1378 of file BC_Handlers.c.

{
    (void)self;
    (void)ctx;
    (void)local_inflow_contribution;
    (void)local_outflow_contribution;
    PetscFunctionBeginUser;
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

Apply_InletProfileFromFile()

static PetscErrorCode Apply_InletProfileFromFile ( BoundaryCondition *  self, BCContext *  ctx  )

static

Applies the loaded PICSLICE scalar profile to Ucont and Ubcs on an inlet face.

Each stored scalar is treated as a positive normal speed magnitude. The routine uses the same negative/positive face sign convention, immersed-boundary skip checks, metric vectors, and CellArea conversion used by the constant and parabolic inlet handlers.

Parameters

self	BoundaryCondition object with InletProfileFileData storage.
ctx	Runtime boundary context containing arrays and face id.

Returns

PetscErrorCode 0 on success, or a PETSc error from DMDA array access.

Definition at line 1404 of file BC_Handlers.c.

{
    PetscErrorCode ierr;
    UserCtx *user = ctx->user;
    BCFace face_id = ctx->face_id;
    InletProfileFileData *data = (InletProfileFileData*)self->data;
    PetscBool can_service;
 
    PetscFunctionBeginUser;
    DMDALocalInfo *info = &user->info;
    Cmpnts ***ubcs, ***ucont, ***csi, ***eta, ***zet;
    PetscReal ***nvert;
 
    ierr = CanRankServiceFace(info, user->IM, user->JM, user->KM, face_id, &can_service); CHKERRQ(ierr);
    if (!can_service) PetscFunctionReturn(0);
 
    ierr = DMDAVecGetArray(user->fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
    ierr = DMDAVecGetArray(user->fda, user->Ucont, &ucont); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
 
    PetscInt xs = info->xs, xe = info->xs + info->xm;
    PetscInt ys = info->ys, ye = info->ys + info->ym;
    PetscInt zs = info->zs, ze = info->zs + info->zm;
    PetscInt mx = info->mx, my = info->my, mz = info->mz;
 
    PetscInt lxs = xs, lxe = xe, lys = ys, lye = ye, lzs = zs, lze = ze;
    if (xs == 0) lxs = xs + 1;
    if (xe == mx) lxe = xe - 1;
    if (ys == 0) lys = ys + 1;
    if (ye == my) lye = ye - 1;
    if (zs == 0) lzs = zs + 1;
    if (ze == mz) lze = ze - 1;
 
    switch (face_id) {
        case BC_FACE_NEG_X:
        case BC_FACE_POS_X: {
            PetscReal sign = (face_id == BC_FACE_NEG_X) ? 1.0 : -1.0;
            PetscInt i = (face_id == BC_FACE_NEG_X) ? xs : mx - 2;
            for (PetscInt k = lzs; k < lze; k++) {
                for (PetscInt j = lys; j < lye; j++) {
                    if ((sign > 0 && nvert[k][j][i+1] > 0.1) ||
                        (sign < 0 && nvert[k][j][i] > 0.1)) continue;
                    PetscReal uin_local = ProfileSpeedAt(data, k - 1, j - 1);
                    PetscReal CellArea = sqrt(csi[k][j][i].x * csi[k][j][i].x +
                                             csi[k][j][i].y * csi[k][j][i].y +
                                             csi[k][j][i].z * csi[k][j][i].z);
                    ucont[k][j][i].x = sign * uin_local * CellArea;
                    ubcs[k][j][i + (sign < 0)].x = sign * uin_local * csi[k][j][i].x / CellArea;
                    ubcs[k][j][i + (sign < 0)].y = sign * uin_local * csi[k][j][i].y / CellArea;
                    ubcs[k][j][i + (sign < 0)].z = sign * uin_local * csi[k][j][i].z / CellArea;
                }
            }
        } break;
        case BC_FACE_NEG_Y:
        case BC_FACE_POS_Y: {
            PetscReal sign = (face_id == BC_FACE_NEG_Y) ? 1.0 : -1.0;
            PetscInt j = (face_id == BC_FACE_NEG_Y) ? ys : my - 2;
            for (PetscInt k = lzs; k < lze; k++) {
                for (PetscInt i = lxs; i < lxe; i++) {
                    if ((sign > 0 && nvert[k][j+1][i] > 0.1) ||
                        (sign < 0 && nvert[k][j][i] > 0.1)) continue;
                    PetscReal uin_local = ProfileSpeedAt(data, k - 1, i - 1);
                    PetscReal CellArea = sqrt(eta[k][j][i].x * eta[k][j][i].x +
                                             eta[k][j][i].y * eta[k][j][i].y +
                                             eta[k][j][i].z * eta[k][j][i].z);
                    ucont[k][j][i].y = sign * uin_local * CellArea;
                    ubcs[k][j + (sign < 0)][i].x = sign * uin_local * eta[k][j][i].x / CellArea;
                    ubcs[k][j + (sign < 0)][i].y = sign * uin_local * eta[k][j][i].y / CellArea;
                    ubcs[k][j + (sign < 0)][i].z = sign * uin_local * eta[k][j][i].z / CellArea;
                }
            }
        } break;
        case BC_FACE_NEG_Z:
        case BC_FACE_POS_Z: {
            PetscReal sign = (face_id == BC_FACE_NEG_Z) ? 1.0 : -1.0;
            PetscInt k = (face_id == BC_FACE_NEG_Z) ? zs : mz - 2;
            for (PetscInt j = lys; j < lye; j++) {
                for (PetscInt i = lxs; i < lxe; i++) {
                    if ((sign > 0 && nvert[k+1][j][i] > 0.1) ||
                        (sign < 0 && nvert[k][j][i] > 0.1)) continue;
                    PetscReal uin_local = ProfileSpeedAt(data, j - 1, i - 1);
                    PetscReal CellArea = sqrt(zet[k][j][i].x * zet[k][j][i].x +
                                             zet[k][j][i].y * zet[k][j][i].y +
                                             zet[k][j][i].z * zet[k][j][i].z);
                    ucont[k][j][i].z = sign * uin_local * CellArea;
                    ubcs[k + (sign < 0)][j][i].x = sign * uin_local * zet[k][j][i].x / CellArea;
                    ubcs[k + (sign < 0)][j][i].y = sign * uin_local * zet[k][j][i].y / CellArea;
                    ubcs[k + (sign < 0)][j][i].z = sign * uin_local * zet[k][j][i].z / CellArea;
                }
            }
        } break;
    }
 
    ierr = DMDAVecRestoreArray(user->fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(user->fda, user->Ucont, &ucont); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

PostStep_InletProfileFromFile()

static PetscErrorCode PostStep_InletProfileFromFile ( BoundaryCondition *  self, BCContext *  ctx, PetscReal *  local_inflow_contribution, PetscReal *  local_outflow_contribution  )

static

Accumulates the applied inlet flux for a file-prescribed profile.

Sums the face-normal Ucont component over the same interior face slots populated by Apply_InletProfileFromFile().

Parameters

self	BoundaryCondition object for this inlet handler.
ctx	Runtime boundary context containing the UserCtx and face id.
local_inflow_contribution	Accumulator incremented by the measured inlet flux.
local_outflow_contribution	Outflow accumulator, intentionally unchanged.

Returns

PetscErrorCode 0 on success, or a PETSc error from DMDA array access.

Definition at line 1524 of file BC_Handlers.c.

{
    PetscErrorCode ierr;
    UserCtx *user = ctx->user;
    BCFace face_id = ctx->face_id;
    PetscBool can_service;
 
    (void)self;
    (void)local_outflow_contribution;
 
    PetscFunctionBeginUser;
    DMDALocalInfo *info = &user->info;
    Cmpnts ***ucont;
 
    ierr = CanRankServiceFace(info, user->IM, user->JM, user->KM, face_id, &can_service); CHKERRQ(ierr);
    if (!can_service) PetscFunctionReturn(0);
 
    ierr = DMDAVecGetArrayRead(user->fda, user->Ucont, (const Cmpnts***)&ucont); CHKERRQ(ierr);
    PetscReal local_flux = 0.0;
 
    PetscInt xs = info->xs, xe = info->xs + info->xm;
    PetscInt ys = info->ys, ye = info->ys + info->ym;
    PetscInt zs = info->zs, ze = info->zs + info->zm;
    PetscInt mx = info->mx, my = info->my, mz = info->mz;
 
    PetscInt lxs = xs, lxe = xe, lys = ys, lye = ye, lzs = zs, lze = ze;
    if (xs == 0) lxs = xs + 1;
    if (xe == mx) lxe = xe - 1;
    if (ys == 0) lys = ys + 1;
    if (ye == my) lye = ye - 1;
    if (zs == 0) lzs = zs + 1;
    if (ze == mz) lze = ze - 1;
 
    switch (face_id) {
        case BC_FACE_NEG_X:
        case BC_FACE_POS_X: {
            PetscInt i = (face_id == BC_FACE_NEG_X) ? xs : mx - 2;
            for (PetscInt k = lzs; k < lze; k++)
                for (PetscInt j = lys; j < lye; j++)
                    local_flux += ucont[k][j][i].x;
        } break;
        case BC_FACE_NEG_Y:
        case BC_FACE_POS_Y: {
            PetscInt j = (face_id == BC_FACE_NEG_Y) ? ys : my - 2;
            for (PetscInt k = lzs; k < lze; k++)
                for (PetscInt i = lxs; i < lxe; i++)
                    local_flux += ucont[k][j][i].y;
        } break;
        case BC_FACE_NEG_Z:
        case BC_FACE_POS_Z: {
            PetscInt k = (face_id == BC_FACE_NEG_Z) ? zs : mz - 2;
            for (PetscInt j = lys; j < lye; j++)
                for (PetscInt i = lxs; i < lxe; i++)
                    local_flux += ucont[k][j][i].z;
        } break;
    }
 
    ierr = DMDAVecRestoreArrayRead(user->fda, user->Ucont, (const Cmpnts***)&ucont); CHKERRQ(ierr);
    *local_inflow_contribution += local_flux;
 
    LOG_ALLOW(LOCAL, LOG_DEBUG, "PostStep_InletProfileFromFile: Face %d, flux = %.6e\n",
              face_id, local_flux);
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

Destroy_InletProfileFromFile()

static PetscErrorCode Destroy_InletProfileFromFile ( BoundaryCondition *  self )

static

Releases private storage owned by a file-prescribed inlet profile handler.

Parameters

self	BoundaryCondition object whose data field stores InletProfileFileData.

Returns

PetscErrorCode 0 on success.

Definition at line 1600 of file BC_Handlers.c.

{
    PetscFunctionBeginUser;
    if (self && self->data) {
        InletProfileFileData *data = (InletProfileFileData*)self->data;
        PetscFree(data->profile);
        PetscFree(data->source_file);
        PetscFree(self->data);
        self->data = NULL;
    }
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

GetBCParamStringLocal()

static PetscErrorCode GetBCParamStringLocal ( BC_Param *  params, const char *  key, const char **  value_out, PetscBool *  found  )

static

Looks up a string-valued boundary-condition parameter in a BC_Param list.

Parameters

	params	Head of the boundary-condition parameter linked list.
	key	Case-insensitive key to search for.
[out]	value_out	Borrowed pointer to the matching value string, or NULL if absent.
[out]	found	PETSC_TRUE when the key is present, PETSC_FALSE otherwise.

Returns

PetscErrorCode 0 on success.

Definition at line 1124 of file BC_Handlers.c.

{
    PetscFunctionBeginUser;
    *found = PETSC_FALSE;
    *value_out = NULL;
    for (BC_Param *param = params; param; param = param->next) {
        if (strcasecmp(param->key, key) == 0) {
            *value_out = param->value;
            *found = PETSC_TRUE;
            PetscFunctionReturn(0);
        }
    }
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

GetProfileFileExpectedDims()

static PetscErrorCode GetProfileFileExpectedDims ( UserCtx *  user, BCFace  face_id, PetscInt *  n1, PetscInt *  n2  )

static

Computes the expected PICSLICE dimensions for an inlet face.

Parameters

	user	User context containing global grid node counts.
	face_id	Boundary face whose tangential profile dimensions are requested.
[out]	n1	First PICSLICE dimension in handler storage order.
[out]	n2	Second PICSLICE dimension in handler storage order.

Returns

PetscErrorCode 0 on success, or a PETSc error for unsupported faces or invalid dimensions.

Definition at line 1149 of file BC_Handlers.c.

{
    PetscFunctionBeginUser;
    switch (face_id) {
        case BC_FACE_NEG_X:
        case BC_FACE_POS_X:
            *n1 = user->KM - 1;
            *n2 = user->JM - 1;
            break;
        case BC_FACE_NEG_Y:
        case BC_FACE_POS_Y:
            *n1 = user->KM - 1;
            *n2 = user->IM - 1;
            break;
        case BC_FACE_NEG_Z:
        case BC_FACE_POS_Z:
            *n1 = user->JM - 1;
            *n2 = user->IM - 1;
            break;
        default:
            SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE,
                    "Unsupported face id %d for inlet profile dimensions.", face_id);
    }
    if (*n1 <= 0 || *n2 <= 0) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
                "Invalid inlet profile dimensions (%d, %d) for grid (%d, %d, %d).",
                *n1, *n2, user->IM, user->JM, user->KM);
    }
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

ReadPicSliceProfile()

static PetscErrorCode ReadPicSliceProfile ( const char *  source_file, PetscInt  expected_n1, PetscInt  expected_n2, InletProfileFileData *  data  )

static

Reads and validates a static scalar inlet profile from a canonical PICSLICE file.

The file must contain magic token PICSLICE, frame count 1, the expected two-dimensional face shape, and exactly one finite nonnegative scalar speed value per interior face slot. Values are stored row-major in data->profile.

Parameters

	source_file	Path to the PICSLICE profile file.
	expected_n1	Expected first slice dimension.
	expected_n2	Expected second slice dimension.
[in,out]	data	Handler-private storage that receives dimensions, profile values, and min/max speeds.

Returns

PetscErrorCode 0 on success, or a PETSc file/validation error on malformed input.

Definition at line 1194 of file BC_Handlers.c.

{
    PetscErrorCode ierr;
    FILE *fd = NULL;
    char magic[32] = {0};
    PetscInt frame_count = 0, n1 = 0, n2 = 0;
 
    PetscFunctionBeginUser;
    fd = fopen(source_file, "r");
    if (!fd) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN,
                     "Cannot open PICSLICE inlet profile file: %s", source_file);
 
    if (fscanf(fd, "%31s", magic) != 1 || strcmp(magic, "PICSLICE") != 0) {
        fclose(fd);
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_READ,
                "PICSLICE inlet profile file %s must begin with PICSLICE header.", source_file);
    }
    if (fscanf(fd, "%d", &frame_count) != 1) {
        fclose(fd);
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_READ,
                "PICSLICE inlet profile file %s missing frame count.", source_file);
    }
    if (frame_count != 1) {
        fclose(fd);
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED,
                "PICSLICE inlet profile file %s has %d frames; static handler requires 1.",
                source_file, frame_count);
    }
    if (fscanf(fd, "%d %d", &n1, &n2) != 2) {
        fclose(fd);
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_READ,
                "PICSLICE inlet profile file %s missing slice dimensions.", source_file);
    }
    if (n1 != expected_n1 || n2 != expected_n2) {
        fclose(fd);
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED,
                "PICSLICE inlet profile dimensions mismatch for %s: expected (%d, %d), found (%d, %d).",
                source_file, expected_n1, expected_n2, n1, n2);
    }
 
    data->n1 = n1;
    data->n2 = n2;
    ierr = PetscMalloc1(n1 * n2, &data->profile); CHKERRQ(ierr);
    data->min_speed = PETSC_MAX_REAL;
    data->max_speed = -PETSC_MAX_REAL;
 
    for (PetscInt idx = 0; idx < n1 * n2; idx++) {
        PetscReal value = 0.0;
        if (fscanf(fd, "%le", &value) != 1) {
            fclose(fd);
            SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_READ,
                    "PICSLICE inlet profile file %s ended early: expected %d values.",
                    source_file, n1 * n2);
        }
        if (PetscIsInfOrNanReal(value) || value < 0.0) {
            fclose(fd);
            SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED,
                    "PICSLICE inlet profile file %s contains invalid speed %.6e at flat index %d.",
                    source_file, (double)value, idx);
        }
        data->profile[idx] = value;
        data->min_speed = PetscMin(data->min_speed, value);
        data->max_speed = PetscMax(data->max_speed, value);
    }
 
    char extra[64];
    if (fscanf(fd, "%63s", extra) == 1) {
        fclose(fd);
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED,
                "PICSLICE inlet profile file %s has extra token after %d values: %s",
                source_file, n1 * n2, extra);
    }
    fclose(fd);
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

ProfileSpeedAt()

static PetscReal ProfileSpeedAt ( const InletProfileFileData *  data, PetscInt  a, PetscInt  b  )

inlinestatic

Returns one scalar speed from the flattened PICSLICE profile.

Parameters

data	Handler-private profile storage.
a	First profile index in face-specific storage order.
b	Second profile index in face-specific storage order.

Returns

Scalar normal speed at (a, b).

Definition at line 1279 of file BC_Handlers.c.

{
    return data->profile[a * data->n2 + b];
}

Here is the caller graph for this function:

Create_InletProfileFromFile()

PetscErrorCode Create_InletProfileFromFile

(

BoundaryCondition *

bc

)

Implementation of Create_InletProfileFromFile().

Configures a BoundaryCondition object for a file-prescribed inlet profile.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/BC_Handlers.h.

See also

Create_InletProfileFromFile()

Definition at line 1292 of file BC_Handlers.c.

{
    PetscErrorCode ierr;
    PetscFunctionBeginUser;
 
    if (!bc) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "BoundaryCondition is NULL");
 
    InletProfileFileData *data = NULL;
    ierr = PetscMalloc1(1, &data); CHKERRQ(ierr);
    data->n1 = 0;
    data->n2 = 0;
    data->profile = NULL;
    data->min_speed = 0.0;
    data->max_speed = 0.0;
    data->source_file = NULL;
    bc->data = (void*)data;
 
    bc->priority   = BC_PRIORITY_INLET;
    bc->Initialize = Initialize_InletProfileFromFile;
    bc->PreStep    = PreStep_InletProfileFromFile;
    bc->Apply      = Apply_InletProfileFromFile;
    bc->PostStep   = PostStep_InletProfileFromFile;
    bc->UpdateUbcs = NULL;
    bc->Destroy    = Destroy_InletProfileFromFile;
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

PreStep_OutletConservation()

static PetscErrorCode PreStep_OutletConservation ( BoundaryCondition *  self, BCContext *  ctx, PetscReal *  local_inflow_contribution, PetscReal *  local_outflow_contribution  )

static

Prepare the outlet-conservation correction before advancing the solver.

Definition at line 1672 of file BC_Handlers.c.

{
    PetscErrorCode ierr;
    UserCtx*       user = ctx->user;
    BCFace         face_id = ctx->face_id;
    DMDALocalInfo* info = &user->info;
    PetscBool      can_service;
 
    // Suppress unused parameter warnings for clarity.
    (void)self;
    (void)local_inflow_contribution;
 
    PetscFunctionBeginUser;
 
    // Step 1: Use the robust utility function to determine if this MPI rank owns a computable
    // portion of the specified boundary face. If not, there is no work to do, so we exit immediately.
    const PetscInt IM_nodes_global = user->IM;
    const PetscInt JM_nodes_global = user->JM;
    const PetscInt KM_nodes_global = user->KM;
    ierr = CanRankServiceFace(info, IM_nodes_global, JM_nodes_global, KM_nodes_global, face_id, &can_service); CHKERRQ(ierr);
 
    if (!can_service) {
        PetscFunctionReturn(0);
    }
 
    // Step 2: Get read-only access to the necessary PETSc arrays.
    // We use the local versions (`lUcat`, `lNvert`) which include ghost cell data,
    // ensuring we have the correct interior values adjacent to the boundary.
    Cmpnts ***ucat, ***csi, ***eta, ***zet;
    PetscReal ***nvert;
    ierr = DMDAVecGetArrayRead(user->fda, user->lUcat, (const Cmpnts***)&ucat); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
 
    PetscReal local_flux_out = 0.0;
    const PetscInt xs=info->xs, xe=info->xs+info->xm;
    const PetscInt ys=info->ys, ye=info->ys+info->ym;
    const PetscInt zs=info->zs, ze=info->zs+info->zm;
    const PetscInt mx=info->mx, my=info->my, mz=info->mz;
 
    // Step 3: Replicate the legacy shrunk loop bounds to exclude corners and edges.
    PetscInt lxs = xs; if (xs == 0)    lxs = xs + 1;
    PetscInt lxe = xe; if (xe == mx)   lxe = xe - 1;
    PetscInt lys = ys; if (ys == 0)    lys = ys + 1;
    PetscInt lye = ye; if (ye == my)   lye = ye - 1;
    PetscInt lzs = zs; if (zs == 0)    lzs = zs + 1;
    PetscInt lze = ze; if (ze == mz)   lze = ze - 1;
 
    // Step 4: Loop over the specified face using the corrected bounds and indexing to calculate flux.
    switch (face_id) {
        case BC_FACE_NEG_X: {
            const PetscInt i_cell = xs + 1; // Index for first interior cell-centered data
            const PetscInt i_face = xs;     // Index for the -X face of that cell
            for (int k=lzs; k<lze; k++) for (int j=lys; j<lye; j++) {
                if (nvert[k][j][i_cell] < 0.1) {
                    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);
                }
            }
            break;
        }
        case BC_FACE_POS_X: {
            const PetscInt i_cell = xe - 2; // Index for last interior cell-centered data
            const PetscInt i_face = xe - 2; // Index for the +X face of that cell
            for (int k=lzs; k<lze; k++) for (int j=lys; j<lye; j++) {
                if (nvert[k][j][i_cell] < 0.1) {
                    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);
                }
            }
            break;
        }
        case BC_FACE_NEG_Y: {
            const PetscInt j_cell = ys + 1;
            const PetscInt j_face = ys;
            for (int k=lzs; k<lze; k++) for (int i=lxs; i<lxe; i++) {
                if (nvert[k][j_cell][i] < 0.1) {
                    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);
                }
            }
            break;
        }
        case BC_FACE_POS_Y: {
            const PetscInt j_cell = ye - 2;
            const PetscInt j_face = ye - 2;
            for (int k=lzs; k<lze; k++) for (int i=lxs; i<lxe; i++) {
                if (nvert[k][j_cell][i] < 0.1) {
                    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);
                }
            }
            break;
        }
        case BC_FACE_NEG_Z: {
            const PetscInt k_cell = zs + 1;
            const PetscInt k_face = zs;
            for (int j=lys; j<lye; j++) for (int i=lxs; i<lxe; i++) {
                if (nvert[k_cell][j][i] < 0.1) {
                    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);
                }
            }
            break;
        }
        case BC_FACE_POS_Z: {
            const PetscInt k_cell = ze - 2;
            const PetscInt k_face = ze - 2;
            for (int j=lys; j<lye; j++) for (int i=lxs; i<lxe; i++) {
                if (nvert[k_cell][j][i] < 0.1) {
                    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);
                }
            }
            break;
        }
    }
 
    // Step 5: Restore the PETSc arrays.
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lUcat, (const Cmpnts***)&ucat); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
 
    // Step 6: Add this face's calculated flux to the accumulator for this rank.
    *local_outflow_contribution += local_flux_out;
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

Apply_OutletConservation()

static PetscErrorCode Apply_OutletConservation ( BoundaryCondition *  self, BCContext *  ctx  )

static

(Handler Action) Applies mass conservation correction to the outlet face.

This function calculates a global correction factor based on the total inflow and outflow fluxes and applies it to the contravariant velocity (ucont) on the outlet face to ensure mass conservation.

Definition at line 1808 of file BC_Handlers.c.

{
    PetscErrorCode ierr;
    (void)self;
    UserCtx*       user = ctx->user;
    BCFace         face_id = ctx->face_id;
    DMDALocalInfo* info = &user->info;
    PetscBool      can_service;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    const PetscInt IM_nodes_global = user->IM;
    const PetscInt JM_nodes_global = user->JM;
    const PetscInt KM_nodes_global = user->KM;
    ierr = CanRankServiceFace(info, IM_nodes_global, JM_nodes_global, KM_nodes_global, face_id, &can_service); CHKERRQ(ierr);
    
    if (!can_service) {
        PROFILE_FUNCTION_END;
        PetscFunctionReturn(0);
    }
 
    // --- STEP 1: Calculate the correction factor using pre-calculated area ---
    PetscReal total_inflow = *ctx->global_inflow_sum + *ctx->global_farfield_inflow_sum;
    PetscReal flux_imbalance = total_inflow - *ctx->global_outflow_sum;
    
    // Directly use the pre-calculated area from the simulation context.
    PetscReal velocity_correction = (PetscAbsReal(user->simCtx->AreaOutSum) > 1e-12) 
                                    ? flux_imbalance / user->simCtx->AreaOutSum 
                                    : 0.0;
    
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Outlet Correction on Face %d: Imbalance=%.4e, Pre-calc Area=%.4e, V_corr=%.4e\n",
              face_id, flux_imbalance, user->simCtx->AreaOutSum, velocity_correction);
 
    // --- STEP 2: Apply the correction to ucont on the outlet face ---
 
    // Get read/write access to necessary arrays
 
    Cmpnts ***ubcs, ***ucont, ***csi, ***eta, ***zet, ***ucat;
    PetscReal ***nvert;
    ierr = DMDAVecGetArray(user->fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
    ierr = DMDAVecGetArray(user->fda, user->Ucont, &ucont); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda,user->lUcat, (const Cmpnts***)&ucat); CHKERRQ(ierr);  
    ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
 
    // Get local grid bounds to exclude corners/edges
    PetscInt xs = info->xs, xe = info->xs + info->xm;
    PetscInt ys = info->ys, ye = info->ys + info->ym;
    PetscInt zs = info->zs, ze = info->zs + info->zm;
    PetscInt mx = info->mx, my = info->my, mz = info->mz;
    PetscInt lxs = xs, lxe = xe, lys = ys, lye = ye, lzs = zs, lze = ze;
 
    if (xs == 0) lxs = xs + 1;
    if (xe == mx) lxe = xe - 1;
    if (ys == 0) lys = ys + 1;
    if (ye == my) lye = ye - 1;
    if (zs == 0) lzs = zs + 1;          
    if (ze == mz) lze = ze - 1;
 
    // Loop over faces and apply correction    
    switch(face_id){
        case BC_FACE_NEG_X:{
            const PetscInt i_cell = xs + 1;
            const PetscInt i_face = xs;
            const PetscInt i_dummy = xs;
            for (PetscInt k = lzs; k < lze; k++) {
                for (PetscInt j = lys; j < lye; j++) {
                    if (nvert[k][j][i_cell] < 0.1) {
                        // Set ubcs
                        ubcs[k][j][i_dummy] = ucat[k][j][i_cell];
 
                        // Calculate Local uncorrected original flux
                        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);
                        
                        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));
                                                
                        PetscReal Correction_flux = velocity_correction*Cell_Area;
 
                        ucont[k][j][i_face].x = Uncorrected_local_flux + Correction_flux;
                    }     
                }
            }
            break;        
        }
        case BC_FACE_POS_X:{
            const PetscInt i_cell = xe - 2;
            const PetscInt i_face = xe - 2;
            const PetscInt i_dummy = xe - 1;
            for(PetscInt k = lzs; k < lze; k++) for (PetscInt j = lys; j < lye; j++){
                if(nvert[k][j][i_cell]<0.1){
                    // Set ubcs
                    ubcs[k][j][i_dummy] = ucat[k][j][i_cell];
                    
                    // Calculate Local uncorrected original flux
                    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);
                    
                    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));
                                            
                    PetscReal Correction_flux = velocity_correction*Cell_Area;
 
                    ucont[k][j][i_face].x = Uncorrected_local_flux + Correction_flux;
                }
            }
            break;
        }
        case BC_FACE_NEG_Y:{
            const PetscInt j_cell = ys + 1;
            const PetscInt j_face = ys;
            const PetscInt j_dummy = ys;
            for(PetscInt k = lzs; k < lze; k++) for (PetscInt i = lxs; i < lxe; i++){
                if(nvert[k][j_cell][i]<0.1){
                    // Set ubcs
                    ubcs[k][j_dummy][i] = ucat[k][j_cell][i];
 
                    // Calculate Local uncorrected original flux
                    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);
                    
                    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));
 
                    PetscReal Correction_flux = velocity_correction*Cell_Area;
 
                    ucont[k][j_face][i].y = Uncorrected_local_flux + Correction_flux;
                }
            }
            break;
        }
        case BC_FACE_POS_Y:{
            const PetscInt j_cell = ye - 2;
            const PetscInt j_face = ye - 2;
            const PetscInt j_dummy = ye - 1;
            for(PetscInt k = lzs; k < lze; k++) for (PetscInt i = lxs; i < lxe; i++){
                if(nvert[k][j_cell][i]<0.1){
                    // Set ubcs
                    ubcs[k][j_dummy][i] = ucat[k][j_cell][i];
 
                    // Calculate Local uncorrected original flux
                    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);
                    
                    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));
 
                    PetscReal Correction_flux = velocity_correction*Cell_Area;
 
                    ucont[k][j_face][i].y = Uncorrected_local_flux + Correction_flux;
                }
            }
            break;
        }
        case BC_FACE_NEG_Z:{
            const PetscInt k_cell = zs + 1;
            const PetscInt k_face = zs;
            const PetscInt k_dummy = zs;
            for(PetscInt j = lys; j < lye; j++) for (PetscInt i = lxs; i < lxe; i++){
                if(nvert[k_cell][j][i]<0.1){
                    // Set ubcs
                    ubcs[k_dummy][j][i] = ucat[k_cell][j][i];
 
                    // Calculate Local uncorrected original flux
                    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));
 
                    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));
 
                    PetscReal Correction_flux = velocity_correction*Cell_Area;
 
                    ucont[k_face][j][i].z = Uncorrected_local_flux + Correction_flux;
                }
            }
            break;
        }
        case BC_FACE_POS_Z:{
            const PetscInt k_cell = ze - 2;
            const PetscInt k_face = ze - 2;
            const PetscInt k_dummy = ze - 1;
             for(PetscInt j = lys; j < lye; j++) for (PetscInt i = lxs; i < lxe; i++){
                if(nvert[k_cell][j][i]<0.1){
                    // Set ubcs
                    ubcs[k_dummy][j][i] = ucat[k_cell][j][i];
 
                    // Calculate Local uncorrected original flux
                    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));
 
                    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));
 
                    PetscReal Correction_flux = velocity_correction*Cell_Area;
 
                    ucont[k_face][j][i].z = Uncorrected_local_flux + Correction_flux;
                }
            }
            break;           
        }
    }    
    
    // Restore all arrays
    ierr = DMDAVecRestoreArray(user->fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(user->fda, user->Ucont, &ucont); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda,user->lUcat, (const Cmpnts***)&ucat); CHKERRQ(ierr);  
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

PostStep_OutletConservation()

static PetscErrorCode PostStep_OutletConservation ( BoundaryCondition *  self, BCContext *  ctx, PetscReal *  in, PetscReal *  out  )

static

Update outlet-conservation state after a completed solver step.

Definition at line 2020 of file BC_Handlers.c.

{
    PetscErrorCode ierr;
    UserCtx*       user = ctx->user;
    BCFace         face_id = ctx->face_id;
    DMDALocalInfo* info = &user->info;
    PetscBool      can_service;
    
    (void)self;
    (void)local_inflow_contribution;
    
    PetscFunctionBeginUser;
    const PetscInt IM_nodes_global = user->IM;
    const PetscInt JM_nodes_global = user->JM;
    const PetscInt KM_nodes_global = user->KM;
    ierr = CanRankServiceFace(info, IM_nodes_global, JM_nodes_global, KM_nodes_global, face_id, &can_service); CHKERRQ(ierr);
 
    if (!can_service) PetscFunctionReturn(0);
 
    // Get arrays (need both ucont and nvert)
    Cmpnts ***ucont;
    PetscReal ***nvert;  // ✅ ADD nvert
    
    ierr = DMDAVecGetArrayRead(user->fda, user->Ucont, (const Cmpnts***)&ucont); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);  // ✅ ADD
 
    PetscReal local_flux = 0.0;
    
    PetscInt xs = info->xs, xe = info->xs + info->xm;
    PetscInt ys = info->ys, ye = info->ys + info->ym;
    PetscInt zs = info->zs, ze = info->zs + info->zm;
    PetscInt mx = info->mx, my = info->my, mz = info->mz;
    
    PetscInt lxs = xs, lxe = xe, lys = ys, lye = ye, lzs = zs, lze = ze;
    if (xs == 0) lxs = xs + 1;
    if (xe == mx) lxe = xe - 1;
    if (ys == 0) lys = ys + 1;
    if (ye == my) lye = ye - 1;
    if (zs == 0) lzs = zs + 1;
    if (ze == mz) lze = ze - 1;
 
    // Sum ucont components, skipping solid cells (same indices as PreStep)
    switch (face_id) {
        case BC_FACE_NEG_X: {
            const PetscInt i_cell = xs + 1;  // ✅ Match PreStep
            const PetscInt i_face = xs;
            for (PetscInt k = lzs; k < lze; k++) {
                for (PetscInt j = lys; j < lye; j++) {
                    if (nvert[k][j][i_cell] < 0.1) {  // ✅ Skip solid cells
                        local_flux += ucont[k][j][i_face].x;
                    }
                }
            }
        } break;
        
        case BC_FACE_POS_X: {
            const PetscInt i_cell = xe - 2;  // ✅ Match PreStep
            const PetscInt i_face = xe - 2;
            for (PetscInt k = lzs; k < lze; k++) {
                for (PetscInt j = lys; j < lye; j++) {
                    if (nvert[k][j][i_cell] < 0.1) {  // ✅ Skip solid cells
                        local_flux += ucont[k][j][i_face].x;
                    }
                }
            }
        } break;
        
        case BC_FACE_NEG_Y: {
            const PetscInt j_cell = ys + 1;  // ✅ Match PreStep
            const PetscInt j_face = ys;
            for (PetscInt k = lzs; k < lze; k++) {
                for (PetscInt i = lxs; i < lxe; i++) {
                    if (nvert[k][j_cell][i] < 0.1) {  // ✅ Skip solid cells
                        local_flux += ucont[k][j_face][i].y;
                    }
                }
            }
        } break;
        
        case BC_FACE_POS_Y: {
            const PetscInt j_cell = ye - 2;  // ✅ Match PreStep
            const PetscInt j_face = ye - 2;
            for (PetscInt k = lzs; k < lze; k++) {
                for (PetscInt i = lxs; i < lxe; i++) {
                    if (nvert[k][j_cell][i] < 0.1) {  // ✅ Skip solid cells
                        local_flux += ucont[k][j_face][i].y;
                    }
                }
            }
        } break;
        
        case BC_FACE_NEG_Z: {
            const PetscInt k_cell = zs + 1;  // ✅ Match PreStep
            const PetscInt k_face = zs;
            for (PetscInt j = lys; j < lye; j++) {
                for (PetscInt i = lxs; i < lxe; i++) {
                    if (nvert[k_cell][j][i] < 0.1) {  // ✅ Skip solid cells
                        local_flux += ucont[k_face][j][i].z;
                    }
                }
            }
        } break;
        
        case BC_FACE_POS_Z: {
            const PetscInt k_cell = ze - 2;  // ✅ Match PreStep
            const PetscInt k_face = ze - 2;
            for (PetscInt j = lys; j < lye; j++) {
                for (PetscInt i = lxs; i < lxe; i++) {
                    if (nvert[k_cell][j][i] < 0.1) {  // ✅ Skip solid cells
                        local_flux += ucont[k_face][j][i].z;
                    }
                }
            }
        } break;
    }
 
    // Restore arrays
    ierr = DMDAVecRestoreArrayRead(user->fda, user->Ucont, (const Cmpnts***)&ucont); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);  // ✅ ADD
 
    // Add to accumulator
    *local_outflow_contribution += local_flux;
    
    LOG_ALLOW(LOCAL, LOG_DEBUG, "PostStep_OutletConservation: Face %d, corrected flux = %.6e\n",
              face_id, local_flux);
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

Create_OutletConservation()

PetscErrorCode Create_OutletConservation

(

BoundaryCondition *

bc

)

Implementation of Create_OutletConservation().

Configures a BoundaryCondition object for conservative outlet treatment.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/BC_Handlers.h.

See also

Create_OutletConservation()

Definition at line 1644 of file BC_Handlers.c.

{
    PetscFunctionBeginUser;
    
    if (!bc) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Input BoundaryCondition is NULL");
 
    // This handler has the highest priority to ensure it runs after
    // all inflow fluxes have been calculated.
    bc->priority   = BC_PRIORITY_OUTLET;
    
    // Assign function pointers
    bc->Initialize = NULL; // No initialization needed
    bc->PreStep    = PreStep_OutletConservation;
    bc->Apply      = Apply_OutletConservation;
    bc->PostStep   = PostStep_OutletConservation;
    bc->UpdateUbcs = NULL;
    bc->Destroy    = NULL; // No private data to destroy
    
    bc->data = NULL;
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

Create_PeriodicGeometric()

PetscErrorCode Create_PeriodicGeometric

(

BoundaryCondition *

bc

)

Implementation of Create_PeriodicGeometric().

Configures a BoundaryCondition object for geometric periodic coupling.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/BC_Handlers.h.

See also

Create_PeriodicGeometric()

Definition at line 2159 of file BC_Handlers.c.

                                                              {
    PetscFunctionBeginUser;
    
    if (!bc) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Input BoundaryCondition is NULL");
    bc->priority = BC_PRIORITY_WALL;
 
    // Assign function pointers
    bc->Initialize = NULL; // No initialization needed
    bc->PreStep    = NULL;
    bc->Apply      = NULL;
    bc->PostStep   = NULL;
    bc->UpdateUbcs = NULL;
    bc->Destroy    = NULL; // No private data to destroy
 
    bc->data = NULL;
 
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

MeasureDrivenFluxes()

static PetscErrorCode MeasureDrivenFluxes ( UserCtx *  user, char  direction, PetscReal *  boundaryFlux, PetscReal *  planarAverageFlux  )

static

Measure the two volumetric fluxes the driven-flow controller senses.

Both periodic driven handlers steer on the same pair of measurements, so they are taken here in a single sweep of lUcont:

• *boundaryFlux is the flux through the single periodic boundary plane. It is fast and responsive but noisy, and drives the boundary trim.
• *planarAverageFlux is the flux averaged over every cross-sectional plane in the driven direction. It is stable and inertial, and drives the momentum source. It is also the quantity initial_flux latches at t=0.

Parameters

[in]	user	Block context supplying lUcont, lNvert and info.
[in]	direction	Driven direction, 'X', 'Y' or 'Z'.
[out]	boundaryFlux	Globally reduced flux through the boundary plane.
[out]	planarAverageFlux	Globally reduced plane-averaged flux.

Returns

PetscErrorCode 0 on success.

Definition at line 2199 of file BC_Handlers.c.

{
    PetscErrorCode ierr;
    DMDALocalInfo info = user->info;
    PetscInt i, j, k;
 
    PetscFunctionBeginUser;
 
    // --- Get read-only access to necessary field data ---
    Cmpnts ***ucont;
    PetscReal ***nvert;
    ierr = DMDAVecGetArrayRead(user->fda, user->lUcont, (const Cmpnts***)&ucont); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
 
    // --- Define local loop bounds ---
    PetscInt lxs = (info.xs == 0) ? 1 : info.xs;
    PetscInt lys = (info.ys == 0) ? 1 : info.ys;
    PetscInt lzs = (info.zs == 0) ? 1 : info.zs;
    PetscInt lxe = (info.xs + info.xm == info.mx) ? info.mx - 1 : info.xs + info.xm;
    PetscInt lye = (info.ys + info.ym == info.my) ? info.my - 1 : info.ys + info.ym;
    PetscInt lze = (info.zs + info.zm == info.mz) ? info.mz - 1 : info.zs + info.zm;
 
    // --- Initialize local accumulators ---
    PetscReal localCurrentBoundaryFlux = 0.0;
    PetscReal localAveragePlanarVolumetricFluxTerm = 0.0;
 
    // --- Measure local contributions to the two flux types, generalized by direction ---
    switch (direction) {
        case 'X':
            if (info.xs == 0) { // Only the rank on the negative face contributes to boundary flux
                i = 0;
                for (k = lzs; k < lze; k++) for (j = lys; j < lye; j++) {
                    if (nvert[k][j][i + 1] < 0.1) localCurrentBoundaryFlux += ucont[k][j][i].x;
                }
            }
            for (i = info.xs; i < lxe; i++) {
                for (k = lzs; k < lze; k++) for (j = lys; j < lye; j++) {
                    if (nvert[k][j][i + 1] < 0.1) localAveragePlanarVolumetricFluxTerm += ucont[k][j][i].x / (PetscReal)(info.mx - 1);
                }
            }
            break;
        case 'Y':
            if (info.ys == 0) {
                j = 0;
                for (k = lzs; k < lze; k++) for (i = lxs; i < lxe; i++) {
                    if (nvert[k][j + 1][i] < 0.1) localCurrentBoundaryFlux += ucont[k][j][i].y;
                }
            }
            for (j = info.ys; j < lye; j++) {
                for (k = lzs; k < lze; k++) for (i = lxs; i < lxe; i++) {
                    if (nvert[k][j + 1][i] < 0.1) localAveragePlanarVolumetricFluxTerm += ucont[k][j][i].y / (PetscReal)(info.my - 1);
                }
            }
            break;
        case 'Z':
            if (info.zs == 0) {
                k = 0;
                for (j = lys; j < lye; j++) for (i = lxs; i < lxe; i++) {
                    if (nvert[k + 1][j][i] < 0.1) localCurrentBoundaryFlux += ucont[k][j][i].z;
                }
            }
            for (k = info.zs; k < lze; k++) {
                for (j = lys; j < lye; j++) for (i = lxs; i < lxe; i++) {
                    if (nvert[k + 1][j][i] < 0.1) localAveragePlanarVolumetricFluxTerm += ucont[k][j][i].z / (PetscReal)(info.mz - 1);
                }
            }
            break;
        default:
            SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
                    "MeasureDrivenFluxes received an unknown driven direction '%c'.", direction);
    }
 
    // --- Release array access as soon as possible ---
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lUcont, (const Cmpnts***)&ucont); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
 
    // --- Perform global reductions to get the final flux values ---
    ierr = MPI_Allreduce(&localCurrentBoundaryFlux, boundaryFlux, 1, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
    ierr = MPI_Allreduce(&localAveragePlanarVolumetricFluxTerm, planarAverageFlux, 1, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
 
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

Initialize_PeriodicDrivenConstant()

static PetscErrorCode Initialize_PeriodicDrivenConstant ( BoundaryCondition *  self, BCContext *  ctx  )

static

Initialize constant forcing data for a periodically driven boundary.

Definition at line 2368 of file BC_Handlers.c.

{
    PetscErrorCode ierr;
    DrivenFluxData *data = (DrivenFluxData*)self->data;
    BCFace face_id = ctx->face_id;
    UserCtx* user = ctx->user;
    
    PetscFunctionBeginUser;
 
    LOG_ALLOW(LOCAL, LOG_DEBUG, "Initializing PERIODIC_DRIVEN_CONSTANT_FLUX handler on Face %s...\n", BCFaceToString(face_id));
 
    // --- 1. Validation: Ensure the mathematical type is PERIODIC ---
    if (user->boundary_faces[face_id].mathematical_type != PERIODIC) {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT,
                "Configuration Error: Handler PERIODIC_DRIVEN_CONSTANT_FLUX on Face %s must be applied to a face with mathematical_type PERIODIC.",
                BCFaceToString(face_id));
    }
    
    // --- 2. Role Assignment: Determine direction and master status ---
    data->isMasterController = PETSC_FALSE;
    switch (face_id) {
        case BC_FACE_NEG_X: data->direction = 'X'; data->isMasterController = PETSC_TRUE; break;
        case BC_FACE_POS_X: data->direction = 'X'; break;
        case BC_FACE_NEG_Y: data->direction = 'Y'; data->isMasterController = PETSC_TRUE; break;
        case BC_FACE_POS_Y: data->direction = 'Y'; break;
        case BC_FACE_NEG_Z: data->direction = 'Z'; data->isMasterController = PETSC_TRUE; break;
        case BC_FACE_POS_Z: data->direction = 'Z'; break;
    }
    
    // --- 3. Parameter Parsing (Master Controller only) ---
    if (data->isMasterController) {
        PetscBool found;
 
        // Attempt to read the 'target_flux' parameter from the bcs.run file.
        ierr = GetBCParamReal(user->boundary_faces[face_id].params, "target_flux",
                              &data->targetVolumetricFlux, &found); CHKERRQ(ierr);                      
                
        // If the required parameter is not found, halt with an informative error.
        if (!found) {
            SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT,
                    "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).",
                    BCFaceToString(face_id));
        }
        
        LOG_ALLOW(GLOBAL, LOG_INFO, "Driven Flow (Dir %c): Constant target volumetric flux set to %le.\n",
                  data->direction, data->targetVolumetricFlux);
 
        // Store the target flux in the UserCtx. This makes it globally accessible
        // to other parts of the solver, such as the `CorrectChannelFluxProfile` enforcer function.
        user->simCtx->targetVolumetricFlux = data->targetVolumetricFlux;
        // The target is fixed for the run from here on; see the flag's comment in
        // variables.h for how initial_flux uses the same latch.
        user->simCtx->drivenFluxTargetLatched = PETSC_TRUE;
    }
 
    PetscBool trimfound;
    // Optional seam-flux enforcement; accepts the deprecated `apply_trim` spelling.
    ierr = GetDrivenSeamFluxFlag(user->boundary_faces[face_id].params,
                                 &data->enforceSeamFlux, &trimfound); CHKERRQ(ierr);
 
    if(!trimfound) LOG_ALLOW(GLOBAL,LOG_DEBUG,"Seam-flux enforcement not specified, defaults to %s.\n",data->enforceSeamFlux? "True":"False");
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

PreStep_PeriodicDrivenConstant()

static PetscErrorCode PreStep_PeriodicDrivenConstant ( BoundaryCondition *  self, BCContext *  ctx, PetscReal *  in, PetscReal *  out  )

static

Prepare constant periodic-driving data before the solver step.

Definition at line 2438 of file BC_Handlers.c.

{
    PetscErrorCode ierr;
    DrivenFluxData *data = (DrivenFluxData*)self->data;
    UserCtx* user = ctx->user;
    SimCtx* simCtx = user->simCtx;
 
    PetscFunctionBeginUser;
 
    // --- Master Check: Only the handler on the negative face performs calculations ---
    if (!data->isMasterController) {
        PetscFunctionReturn(0);
    }
    
    // The controller senses two fluxes; see MeasureDrivenFluxes() for what each
    // one is for and why the controller needs both.
    char direction = data->direction;
    PetscReal globalCurrentBoundaryFlux, globalAveragePlanarVolumetricFlux;
    ierr = MeasureDrivenFluxes(user, direction, &globalCurrentBoundaryFlux,
                               &globalAveragePlanarVolumetricFlux); CHKERRQ(ierr);
 
    // --- Get cross-sectional area using the dedicated geometry function ---
    Cmpnts ignored_center;
    PetscReal globalBoundaryArea;
    BCFace neg_face_id = (direction == 'X') ? BC_FACE_NEG_X : (direction == 'Y') ? BC_FACE_NEG_Y : BC_FACE_NEG_Z;
    ierr = CalculateFaceCenterAndArea(user, neg_face_id, &ignored_center, &globalBoundaryArea); CHKERRQ(ierr);
 
    // --- Calculate the two correction terms ---
    //
    // These are refreshed on DIFFERENT cadences, and the difference matters.
    //
    // ApplyBoundaryConditions() runs BoundarySystem_ExecuteStep() -- and hence
    // this PreStep -- three times per call, and it is itself called once per
    // Jameson RK stage under the Picard solver and once per residual evaluation
    // under Newton-Krylov. So "per PreStep" is emphatically not "per timestep".
    //
    //  - bulkVelocityCorrection scales the momentum source in ComputeRHS. Both
    //    momentum solvers are built on it being FROZEN across a timestep: the
    //    Picard shadow-Jacobian estimate treats the body force as a constant
    //    forcing with zero velocity Jacobian, and the Newton solve needs a
    //    source that does not drift between residual evaluations. It is
    //    therefore computed once per physical step, from the field at the start
    //    of that step, and held.
    //
    //  - boundaryVelocityCorrection is the tactical trim applied to the
    //    boundary fluxes in Apply(). It is deliberately re-measured on every
    //    pass: Apply() accumulates it into Ucont, and re-measuring is what makes
    //    that accumulation self-limiting as the seam converges. Freezing it
    //    would make repeated Apply() calls add the same trim over and over.
    if (globalBoundaryArea > 1.0e-12) {
        if (data->lastBulkCorrectionStep != simCtx->step) {
            simCtx->bulkVelocityCorrection = (data->targetVolumetricFlux - globalAveragePlanarVolumetricFlux) / globalBoundaryArea;
            data->lastBulkCorrectionStep = simCtx->step;
        }
        simCtx->boundaryVelocityCorrection = (data->targetVolumetricFlux - globalCurrentBoundaryFlux) / globalBoundaryArea;
    } else {
        simCtx->bulkVelocityCorrection = 0.0;
        simCtx->boundaryVelocityCorrection = 0.0;
        data->lastBulkCorrectionStep = simCtx->step;
    }
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Driven Flow Controller Update (Dir %c):\n", data->direction);
    LOG_ALLOW(GLOBAL, LOG_INFO, "  - Target Volumetric Flux:          %.6e\n", data->targetVolumetricFlux);
    LOG_ALLOW(GLOBAL, LOG_INFO, "  - Avg Planar Volumetric Flux (Stable): %.6e\n", globalAveragePlanarVolumetricFlux);
    LOG_ALLOW(GLOBAL, LOG_INFO, "  - Boundary Flux (Fast):            %.6e\n", globalCurrentBoundaryFlux);
    LOG_ALLOW(GLOBAL, LOG_INFO, "  - Bulk Velocity Correction:     %.6e (For Momentum Source)\n", simCtx->bulkVelocityCorrection);
    LOG_ALLOW(GLOBAL, LOG_INFO, "  - Boundary Velocity Correction: %.6e (For Boundary Trim)\n", simCtx->boundaryVelocityCorrection);
    
    // Suppress unused parameter warnings for this handler
    (void)local_inflow_contribution;
    (void)local_outflow_contribution;
    
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

Apply_PeriodicDrivenConstant()

static PetscErrorCode Apply_PeriodicDrivenConstant ( BoundaryCondition *  self, BCContext *  ctx  )

static

Apply the configured constant driving term to the periodic boundary.

Definition at line 2520 of file BC_Handlers.c.

{
    PetscErrorCode ierr;
    DrivenFluxData *data = (DrivenFluxData*)self->data;
    UserCtx* user = ctx->user;
    BCFace face_id = ctx->face_id;
    PetscBool can_service;
    
    PetscFunctionBeginUser;
    
    // Check if this rank owns part of this boundary face
    ierr = CanRankServiceFace(&user->info, user->IM, user->JM, user->KM, face_id, &can_service); CHKERRQ(ierr);
    if (!can_service) {
        PetscFunctionReturn(0);
    }
    
    // If the correction is negligible, no work is needed.
    if (PetscAbsReal(user->simCtx->boundaryVelocityCorrection) < 1e-12) {
        PetscFunctionReturn(0);
    }
 
    LOG_ALLOW(LOCAL, LOG_TRACE, "Apply_PeriodicDrivenConstant: Applying boundary trim on Face %s...\n", BCFaceToString(face_id));
 
    // --- Get read/write access to necessary arrays ---
    DMDALocalInfo info = user->info;
    Cmpnts ***ucont, ***uch, ***csi, ***eta, ***zet;
    PetscReal ***nvert;
    
    ierr = DMDAVecGetArray(user->fda, user->Ucont, &ucont); CHKERRQ(ierr);
    ierr = DMDAVecGetArray(user->fda, user->Bcs.Uch, &uch); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
    
    PetscInt lxs = (info.xs == 0) ? 1 : info.xs;
    PetscInt lys = (info.ys == 0) ? 1 : info.ys;
    PetscInt lzs = (info.zs == 0) ? 1 : info.zs;
    PetscInt lxe = (info.xs + info.xm == info.mx) ? info.mx - 1 : info.xs + info.xm;
    PetscInt lye = (info.ys + info.ym == info.my) ? info.my - 1 : info.ys + info.ym;
    PetscInt lze = (info.zs + info.zm == info.mz) ? info.mz - 1 : info.zs + info.zm;
 
    // --- Apply correction to the appropriate face and velocity component ---
    switch (face_id) {
        case BC_FACE_NEG_X: case BC_FACE_POS_X: {
            PetscInt i_face = (face_id == BC_FACE_NEG_X) ? info.xs : info.mx - 2;
            PetscInt i_nvert = (face_id == BC_FACE_NEG_X) ? info.xs + 1 : info.mx - 2;
 
            for (PetscInt k = lzs; k < lze; k++) for (PetscInt j = lys; j < lye; j++) {
                if (nvert[k][j][i_nvert] < 0.1) {
                    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);
                    PetscReal fluxTrim = user->simCtx->boundaryVelocityCorrection * faceArea;
                    if(data->enforceSeamFlux) ucont[k][j][i_face].x += fluxTrim;
                    uch[k][j][i_face].x = fluxTrim; // Store correction for diagnostics
                }
            }
        } break;
        
        case BC_FACE_NEG_Y: case BC_FACE_POS_Y: {
            PetscInt j_face = (face_id == BC_FACE_NEG_Y) ? info.ys : info.my - 2;
            PetscInt j_nvert = (face_id == BC_FACE_NEG_Y) ? info.ys + 1 : info.my - 2;
            
            for (PetscInt k = lzs; k < lze; k++) for (PetscInt i = lxs; i < lxe; i++) {
                if (nvert[k][j_nvert][i] < 0.1) {
                    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);
                    PetscReal fluxTrim = user->simCtx->boundaryVelocityCorrection * faceArea;
                    if(data->enforceSeamFlux) ucont[k][j_face][i].y += fluxTrim;
                    uch[k][j_face][i].y = fluxTrim;
                }
            }
        } break;
            
        case BC_FACE_NEG_Z: case BC_FACE_POS_Z: {
            PetscInt k_face = (face_id == BC_FACE_NEG_Z) ? info.zs : info.mz - 2;
            PetscInt k_nvert = (face_id == BC_FACE_NEG_Z) ? info.zs + 1 : info.mz - 2;
            
            for (PetscInt j = lys; j < lye; j++) for (PetscInt i = lxs; i < lxe; i++) {
                if (nvert[k_nvert][j][i] < 0.1) {
                    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);
                    PetscReal fluxTrim = user->simCtx->boundaryVelocityCorrection * faceArea;
                    if(data->enforceSeamFlux) ucont[k_face][j][i].z += fluxTrim;
                    uch[k_face][j][i].z = fluxTrim;
                }
            }
        } break;
    }
 
    // --- Restore arrays ---
    ierr = DMDAVecRestoreArray(user->fda, user->Ucont, &ucont); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(user->fda, user->Bcs.Uch, &uch); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
    
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

Destroy_PeriodicDrivenConstant()

static PetscErrorCode Destroy_PeriodicDrivenConstant ( BoundaryCondition *  self )

static

Release resources owned by the constant periodic-driving boundary.

Definition at line 2623 of file BC_Handlers.c.

{
    PetscFunctionBeginUser;
 
    // Check that the handler object and its private data pointer are valid before trying to free.
    if (self && self->data) {
        // The private data was allocated with PetscNew(), so it must be freed with PetscFree().
        PetscFree(self->data);
        
        // It is good practice to nullify the pointer after freeing to prevent
        // any accidental use of the dangling pointer (use-after-free).
        self->data = NULL;
        
        LOG_ALLOW(LOCAL, LOG_TRACE, "Destroy_PeriodicDrivenConstant: Private data freed successfully.\n");
    }
 
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

Create_PeriodicDrivenConstant()

PetscErrorCode Create_PeriodicDrivenConstant

(

BoundaryCondition *

bc

)

Internal helper implementation: Create_PeriodicDrivenConstant().

Configures a BoundaryCondition object for periodic driven-flow forcing.

Local to this translation unit.

Definition at line 2325 of file BC_Handlers.c.

{
    PetscErrorCode ierr;
    PetscFunctionBeginUser;
    
    if (!bc) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Input BoundaryCondition object is NULL in Create_PeriodicDrivenConstantFlux");
 
    // --- Allocate the private data structure ---
    DrivenFluxData *data = NULL;
    ierr = PetscNew(&data); CHKERRQ(ierr);
    // Initialize fields to safe default values
    data->direction = ' ';
    data->targetVolumetricFlux = 0.0;
    data->isMasterController = PETSC_FALSE;
    data->enforceSeamFlux = PETSC_FALSE;
    data->lastBulkCorrectionStep = -1;
    
    // Attach the private data to the generic handler object
    bc->data = (void*)data;
    
    // --- Configure the handler's properties and methods ---
    
    // Set priority: Using BC_PRIORITY_INLET ensures this handler's PreStep runs
    // before other handlers (like outlets) that might depend on its calculations.
    // It is the caller's responsibility that there are no Inlets called along with driven periodic to avoid clash.
    bc->priority   = BC_PRIORITY_INLET;
    
    // Assign the function pointers to the implementations in this file.
    bc->Initialize = Initialize_PeriodicDrivenConstant;
    bc->PreStep    = PreStep_PeriodicDrivenConstant;
    bc->Apply      = Apply_PeriodicDrivenConstant;
    bc->PostStep   = NULL; // This handler has no action after the main solver step.
    bc->UpdateUbcs = NULL; // The boundary value is not flow-dependent (it's periodic).
    bc->Destroy    = Destroy_PeriodicDrivenConstant;
    
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

Initialize_PeriodicDrivenInitial()

static PetscErrorCode Initialize_PeriodicDrivenInitial ( BoundaryCondition *  self, BCContext *  ctx  )

static

Initialize forcing data for a periodic boundary driven to its initial flux.

Note

The target itself is NOT measured here. Boundary handlers are initialized before InitializeEulerianState() runs, so at this point Ucont still holds zeros. The measurement is deferred to the first PreStep, by which time either the initial condition has been applied or the restart target has been restored from the checkpoint manifest.

Definition at line 2716 of file BC_Handlers.c.

{
    PetscErrorCode ierr;
    DrivenFluxData *data = (DrivenFluxData*)self->data;
    BCFace face_id = ctx->face_id;
    UserCtx* user = ctx->user;
 
    PetscFunctionBeginUser;
 
    LOG_ALLOW(LOCAL, LOG_DEBUG, "Initializing PERIODIC_DRIVEN_INITIAL_FLUX handler on Face %s...\n", BCFaceToString(face_id));
 
    // --- 1. Validation: Ensure the mathematical type is PERIODIC ---
    if (user->boundary_faces[face_id].mathematical_type != PERIODIC) {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT,
                "Configuration Error: Handler PERIODIC_DRIVEN_INITIAL_FLUX on Face %s must be applied to a face with mathematical_type PERIODIC.",
                BCFaceToString(face_id));
    }
 
    // --- 2. Role Assignment: Determine direction and master status ---
    data->isMasterController = PETSC_FALSE;
    switch (face_id) {
        case BC_FACE_NEG_X: data->direction = 'X'; data->isMasterController = PETSC_TRUE; break;
        case BC_FACE_POS_X: data->direction = 'X'; break;
        case BC_FACE_NEG_Y: data->direction = 'Y'; data->isMasterController = PETSC_TRUE; break;
        case BC_FACE_POS_Y: data->direction = 'Y'; break;
        case BC_FACE_NEG_Z: data->direction = 'Z'; data->isMasterController = PETSC_TRUE; break;
        case BC_FACE_POS_Z: data->direction = 'Z'; break;
    }
 
    // --- 3. Parameter Parsing (Master Controller only) ---
    if (data->isMasterController) {
        PetscReal unused_flux;
        PetscBool found;
 
        // This handler derives its own target, so an explicit one is a config error
        // rather than something to silently ignore. Users who want to prescribe the
        // flux should select the `constant_flux` handler instead.
        ierr = GetBCParamReal(user->boundary_faces[face_id].params, "target_flux",
                              &unused_flux, &found); CHKERRQ(ierr);
        if (found) {
            SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT,
                    "Configuration Error: Handler PERIODIC_DRIVEN_INITIAL_FLUX on Face %s takes no 'target_flux' parameter; "
                    "it measures the flux of the initial condition and holds that. Use handler 'constant_flux' to prescribe a target.",
                    BCFaceToString(face_id));
        }
 
        LOG_ALLOW(GLOBAL, LOG_INFO, "Driven Flow (Dir %c): target volumetric flux will be latched from the initial state.\n",
                  data->direction);
    }
 
    PetscBool trimfound;
    // Optional seam-flux enforcement; accepts the deprecated `apply_trim` spelling.
    ierr = GetDrivenSeamFluxFlag(user->boundary_faces[face_id].params,
                                 &data->enforceSeamFlux, &trimfound); CHKERRQ(ierr);
 
    if(!trimfound) LOG_ALLOW(GLOBAL,LOG_DEBUG,"Seam-flux enforcement not specified, defaults to %s.\n",data->enforceSeamFlux? "True":"False");
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

PreStep_PeriodicDrivenInitial()

static PetscErrorCode PreStep_PeriodicDrivenInitial ( BoundaryCondition *  self, BCContext *  ctx, PetscReal *  local_inflow_contribution, PetscReal *  local_outflow_contribution  )

static

Latch the initial-state flux once, then drive to it like a constant target.

The latch is one-shot and guarded by simCtx->drivenFluxTargetLatched:

• Fresh start: the flag is false and Ucont now holds the initial condition, so the plane-averaged flux is measured and stored.
• Restart: ReadSimulationFields() restored both the target and the flag from the checkpoint manifest, so the original target survives instead of being re-measured from a drifted field.

Definition at line 2788 of file BC_Handlers.c.

{
    PetscErrorCode ierr;
    DrivenFluxData *data = (DrivenFluxData*)self->data;
    SimCtx* simCtx = ctx->user->simCtx;
 
    PetscFunctionBeginUser;
 
    if (data->isMasterController) {
        if (!simCtx->drivenFluxTargetLatched) {
            PetscReal boundaryFlux, planarAverageFlux;
 
            /* Boundary handlers are also exercised once during setup, from
             * FinalizeBlockState(), and at that point the initial condition has
             * been written to Ucont but not yet scattered into the ghosted
             * lUcont that MeasureDrivenFluxes() reads. Latching there would
             * record a target of zero. Setup runs with step == StartStep, so
             * wait for the first PreStep of the first real timestep: by then
             * lUcont holds the initial condition (or, on a restart, the state
             * the restored target already describes). Do no work at all until
             * then, so no bogus correction is derived from a zero target. */
            if (simCtx->step <= simCtx->StartStep) {
                simCtx->bulkVelocityCorrection = 0.0;
                simCtx->boundaryVelocityCorrection = 0.0;
                PetscFunctionReturn(0);
            }
 
            ierr = MeasureDrivenFluxes(ctx->user, data->direction,
                                       &boundaryFlux, &planarAverageFlux); CHKERRQ(ierr);
 
            simCtx->targetVolumetricFlux    = planarAverageFlux;
            simCtx->drivenFluxTargetLatched = PETSC_TRUE;
 
            LOG_ALLOW(GLOBAL, LOG_INFO,
                      "Driven Flow (Dir %c): latched initial volumetric flux %.6e as the target.\n",
                      data->direction, planarAverageFlux);
        }
        // Keep the handler's copy in step with the authoritative value in SimCtx,
        // which is also where a restart deposits the restored target.
        data->targetVolumetricFlux = simCtx->targetVolumetricFlux;
    }
 
    ierr = PreStep_PeriodicDrivenConstant(self, ctx,
                                          local_inflow_contribution,
                                          local_outflow_contribution); CHKERRQ(ierr);
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

Create_PeriodicDrivenInitial()

PetscErrorCode Create_PeriodicDrivenInitial

(

BoundaryCondition *

bc

)

Internal helper implementation: Create_PeriodicDrivenInitial().

Configures a BoundaryCondition object for initial-flux periodic driving.

Local to this translation unit.

Definition at line 2669 of file BC_Handlers.c.

{
    PetscErrorCode ierr;
    PetscFunctionBeginUser;
 
    if (!bc) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Input BoundaryCondition object is NULL in Create_PeriodicDrivenInitial");
 
    // --- Allocate the private data structure ---
    DrivenFluxData *data = NULL;
    ierr = PetscNew(&data); CHKERRQ(ierr);
    // Initialize fields to safe default values
    data->direction = ' ';
    data->targetVolumetricFlux = 0.0;
    data->isMasterController = PETSC_FALSE;
    data->enforceSeamFlux = PETSC_FALSE;
    data->lastBulkCorrectionStep = -1;
 
    // Attach the private data to the generic handler object
    bc->data = (void*)data;
 
    // --- Configure the handler's properties and methods ---
 
    // Same priority reasoning as the constant-flux handler: PreStep must run
    // before any handler that depends on the controller's corrections.
    bc->priority   = BC_PRIORITY_INLET;
 
    bc->Initialize = Initialize_PeriodicDrivenInitial;
    bc->PreStep    = PreStep_PeriodicDrivenInitial;
    bc->Apply      = Apply_PeriodicDrivenConstant;  // Boundary trim is target-agnostic.
    bc->PostStep   = NULL; // This handler has no action after the main solver step.
    bc->UpdateUbcs = NULL; // The boundary value is not flow-dependent (it's periodic).
    bc->Destroy    = Destroy_PeriodicDrivenConstant; // Same private data layout.
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function: