Boundaries.c File Reference

#include "Boundaries.h"
#include <string.h>
#include <ctype.h>

Include dependency graph for Boundaries.c:

Go to the source code of this file.

Macros
#define	__FUNCT__   "CanRankServiceInletFace"
#define	__FUNCT__   "CanRankServiceFace"
#define	__FUNCT__   "GetDeterministicFaceGridLocation"
#define	__FUNCT__   "GetRandomFCellAndLogicOnInletFace"
#define	__FUNCT__   "EnforceRHSBoundaryConditions"
#define	__FUNCT__   "BoundaryCondition_Create"
#define	__FUNCT__   "BoundarySystem_Validate"
#define	__FUNCT__   "BoundarySystem_Initialize"
#define	__FUNCT__   "PropagateBoundaryConfigToCoarserLevels"
#define	__FUNCT__   "BoundarySystem_ExecuteStep"
#define	__FUNCT__   "BoundarySystem_RefreshUbcs"
#define	__FUNCT__   "BoundarySystem_Destroy"
#define	__FUNCT__   "TransferPeriodicFieldByDirection"
#define	__FUNCT__   "SynchronizePeriodicCellFields"
#define	__FUNCT__   "GetPersistentFaceField"
#define	__FUNCT__   "TransferPeriodicFaceFieldByDirection"
#define	__FUNCT__   "SynchronizePeriodicFaceFields"
#define	__FUNCT__   "GetPersistentStaggeredField"
#define	__FUNCT__   "TransferPeriodicStaggeredFieldByDirection"
#define	__FUNCT__   "SynchronizePeriodicStaggeredFields"
#define	__FUNCT__   "PreparePeriodicQuickStencilFields"
#define	__FUNCT__   "SynchronizePeriodicLocalStaggeredField"
#define	__FUNCT__   "TransferPeriodicField"
#define	__FUNCT__   "TransferPeriodicFaceField"
#define	__FUNCT__   "ApplyMetricsPeriodicBCs"
#define	__FUNCT__   "ApplyPeriodicBCs"
#define	__FUNCT__   "UpdateDummyCells"
#define	__FUNCT__   "UpdateCornerNodes"
#define	__FUNCT__   "UpdatePeriodicCornerNodes"
#define	__FUNCT__   "ApplyWallFunction"
#define	__FUNCT__   "FinalizePostProjectionCellFields"
#define	__FUNCT__   "ApplyBoundaryConditions"

Functions
PetscErrorCode	CanRankServiceInletFace (UserCtx *user, const DMDALocalInfo *info, PetscInt IM_nodes_global, PetscInt JM_nodes_global, PetscInt KM_nodes_global, PetscBool *can_service_inlet_out)
	Internal helper implementation: CanRankServiceInletFace().
PetscErrorCode	CanRankServiceFace (const DMDALocalInfo *info, PetscInt IM_nodes_global, PetscInt JM_nodes_global, PetscInt KM_nodes_global, BCFace face_id, PetscBool *can_service_out)
	Implementation of CanRankServiceFace().
PetscErrorCode	GetDeterministicFaceGridLocation (UserCtx *user, const DMDALocalInfo *info, PetscInt xs_gnode_rank, PetscInt ys_gnode_rank, PetscInt zs_gnode_rank, PetscInt IM_cells_global, PetscInt JM_cells_global, PetscInt KM_cells_global, PetscInt64 particle_global_id, PetscInt *ci_metric_lnode_out, PetscInt *cj_metric_lnode_out, PetscInt *ck_metric_lnode_out, PetscReal *xi_metric_logic_out, PetscReal *eta_metric_logic_out, PetscReal *zta_metric_logic_out, PetscBool *placement_successful_out)
	Internal helper implementation: GetDeterministicFaceGridLocation().
PetscErrorCode	GetRandomCellAndLogicalCoordsOnInletFace (UserCtx *user, const DMDALocalInfo *info, PetscInt xs_gnode_rank, PetscInt ys_gnode_rank, PetscInt zs_gnode_rank, PetscInt IM_nodes_global, PetscInt JM_nodes_global, PetscInt KM_nodes_global, PetscRandom *rand_logic_i_ptr, PetscRandom *rand_logic_j_ptr, PetscRandom *rand_logic_k_ptr, PetscInt *ci_metric_lnode_out, PetscInt *cj_metric_lnode_out, PetscInt *ck_metric_lnode_out, PetscReal *xi_metric_logic_out, PetscReal *eta_metric_logic_out, PetscReal *zta_metric_logic_out)
	Internal helper implementation: GetRandomCellAndLogicalCoordsOnInletFace().
PetscErrorCode	EnforceRHSBoundaryConditions (UserCtx *user)
	Internal helper implementation: EnforceRHSBoundaryConditions().
PetscErrorCode	BoundaryCondition_Create (BCHandlerType handler_type, BoundaryCondition **new_bc_ptr)
	Internal helper implementation: BoundaryCondition_Create().
PetscErrorCode	BoundarySystem_Validate (UserCtx *user)
	Internal helper implementation: BoundarySystem_Validate().
PetscErrorCode	BoundarySystem_Initialize (UserCtx *user, const char *bcs_filename)
	Implementation of BoundarySystem_Initialize().
PetscErrorCode	PropagateBoundaryConfigToCoarserLevels (SimCtx *simCtx)
	Internal helper implementation: PropagateBoundaryConfigToCoarserLevels().
PetscErrorCode	BoundarySystem_ExecuteStep (UserCtx *user)
	Implementation of BoundarySystem_ExecuteStep().
PetscErrorCode	BoundarySystem_RefreshUbcs (UserCtx *user)
	Internal helper implementation: BoundarySystem_RefreshUbcs().
PetscErrorCode	BoundarySystem_Destroy (UserCtx *user)
	Implementation of BoundarySystem_Destroy().
PetscErrorCode	TransferPeriodicFieldByDirection (UserCtx *user, const char *field_name, char direction)
	Internal helper implementation: TransferPeriodicFieldByDirection().
PetscErrorCode	SynchronizePeriodicCellFields (UserCtx *user, PetscInt num_fields, const char *field_names[])
	Implementation of SynchronizePeriodicCellFields().
static PetscErrorCode	GetPersistentFaceField (UserCtx *user, const char *field_name, char face_direction, DM *dm, Vec *global_vec, Vec *local_vec, PetscInt *dof)
	Resolves one registered persistent single-face-family field.
static PetscBool	IsFaceCenterCoordinateField (const char *field_name)
	Returns whether a registered face field stores physical coordinates.
static PetscErrorCode	TranslatePeriodicFaceCenterGhosts (UserCtx *user, Vec local_vec)
	Applies geometric translations to wrapped face-center ghost coordinates.
PetscErrorCode	TransferPeriodicFaceFieldByDirection (UserCtx *user, const char *field_name, char face_direction, char periodic_direction)
	Implementation of TransferPeriodicFaceFieldByDirection().
PetscErrorCode	SynchronizePeriodicFaceFields (UserCtx *user, char face_direction, PetscInt num_fields, const char *field_names[])
	Implementation of SynchronizePeriodicFaceFields().
static PetscErrorCode	GetPersistentStaggeredField (UserCtx *user, const char *field_name, DM *dm, Vec *global_vec, Vec *local_vec)
	Resolves one registered persistent component-staggered field.
PetscErrorCode	TransferPeriodicStaggeredFieldByDirection (UserCtx *user, const char *field_name, char periodic_direction)
	Implementation of TransferPeriodicStaggeredFieldByDirection().
PetscErrorCode	SynchronizePeriodicStaggeredFields (UserCtx *user, PetscInt num_fields, const char *field_names[])
	Implementation of SynchronizePeriodicStaggeredFields().
PetscErrorCode	PreparePeriodicQuickStencilFields (UserCtx *user, Vec local_vector_field, Vec local_scalar_field)
	Implementation of PreparePeriodicQuickStencilFields().
PetscErrorCode	SynchronizePeriodicLocalStaggeredField (UserCtx *user, Vec local_field)
	Implementation of SynchronizePeriodicLocalStaggeredField().
PetscErrorCode	TransferPeriodicField (UserCtx *user, const char *field_name)
	Internal helper implementation: TransferPeriodicField().
PetscErrorCode	TransferPeriodicFaceField (UserCtx *user, const char *field_name)
	Internal helper implementation: TransferPeriodicFaceField().
PetscErrorCode	ApplyMetricsPeriodicBCs (UserCtx *user)
	Internal helper implementation: ApplyMetricsPeriodicBCs().
PetscErrorCode	ApplyPeriodicBCs (UserCtx *user)
	Internal helper implementation: ApplyPeriodicBCs().
PetscErrorCode	UpdateDummyCells (UserCtx *user)
	Internal helper implementation: UpdateDummyCells().
PetscErrorCode	UpdateCornerNodes (UserCtx *user)
	Internal helper implementation: UpdateCornerNodes().
PetscErrorCode	UpdatePeriodicCornerNodes (UserCtx *user, PetscInt num_fields, const char *field_names[])
	Internal helper implementation: UpdatePeriodicCornerNodes().
PetscErrorCode	ApplyWallFunction (UserCtx *user)
	Internal helper implementation: ApplyWallFunction().
PetscErrorCode	FinalizePostProjectionCellFields (UserCtx *user)
	Implementation of FinalizePostProjectionCellFields().
PetscErrorCode	ApplyBoundaryConditions (UserCtx *user)
	Implementation of ApplyBoundaryConditions().

Macro Definition Documentation

__FUNCT__ [1/32]

#define __FUNCT__   "CanRankServiceInletFace"

Definition at line 6 of file Boundaries.c.

__FUNCT__ [2/32]

#define __FUNCT__   "CanRankServiceFace"

Definition at line 6 of file Boundaries.c.

__FUNCT__ [3/32]

#define __FUNCT__   "GetDeterministicFaceGridLocation"

Definition at line 6 of file Boundaries.c.

__FUNCT__ [4/32]

#define __FUNCT__   "GetRandomFCellAndLogicOnInletFace"

Definition at line 6 of file Boundaries.c.

__FUNCT__ [5/32]

#define __FUNCT__   "EnforceRHSBoundaryConditions"

Definition at line 6 of file Boundaries.c.

__FUNCT__ [6/32]

#define __FUNCT__   "BoundaryCondition_Create"

Definition at line 6 of file Boundaries.c.

__FUNCT__ [7/32]

#define __FUNCT__   "BoundarySystem_Validate"

Definition at line 6 of file Boundaries.c.

__FUNCT__ [8/32]

#define __FUNCT__   "BoundarySystem_Initialize"

Definition at line 6 of file Boundaries.c.

__FUNCT__ [9/32]

#define __FUNCT__   "PropagateBoundaryConfigToCoarserLevels"

Definition at line 6 of file Boundaries.c.

__FUNCT__ [10/32]

#define __FUNCT__   "BoundarySystem_ExecuteStep"

Definition at line 6 of file Boundaries.c.

__FUNCT__ [11/32]

#define __FUNCT__   "BoundarySystem_RefreshUbcs"

Definition at line 6 of file Boundaries.c.

__FUNCT__ [12/32]

#define __FUNCT__   "BoundarySystem_Destroy"

Definition at line 6 of file Boundaries.c.

__FUNCT__ [13/32]

#define __FUNCT__   "TransferPeriodicFieldByDirection"

Definition at line 6 of file Boundaries.c.

__FUNCT__ [14/32]

#define __FUNCT__   "SynchronizePeriodicCellFields"

Definition at line 6 of file Boundaries.c.

__FUNCT__ [15/32]

#define __FUNCT__   "GetPersistentFaceField"

Definition at line 6 of file Boundaries.c.

__FUNCT__ [16/32]

#define __FUNCT__   "TransferPeriodicFaceFieldByDirection"

Definition at line 6 of file Boundaries.c.

__FUNCT__ [17/32]

#define __FUNCT__   "SynchronizePeriodicFaceFields"

Definition at line 6 of file Boundaries.c.

__FUNCT__ [18/32]

#define __FUNCT__   "GetPersistentStaggeredField"

Definition at line 6 of file Boundaries.c.

__FUNCT__ [19/32]

#define __FUNCT__   "TransferPeriodicStaggeredFieldByDirection"

Definition at line 6 of file Boundaries.c.

__FUNCT__ [20/32]

#define __FUNCT__   "SynchronizePeriodicStaggeredFields"

Definition at line 6 of file Boundaries.c.

__FUNCT__ [21/32]

#define __FUNCT__   "PreparePeriodicQuickStencilFields"

Definition at line 6 of file Boundaries.c.

__FUNCT__ [22/32]

#define __FUNCT__   "SynchronizePeriodicLocalStaggeredField"

Definition at line 6 of file Boundaries.c.

__FUNCT__ [23/32]

#define __FUNCT__   "TransferPeriodicField"

Definition at line 6 of file Boundaries.c.

__FUNCT__ [24/32]

#define __FUNCT__   "TransferPeriodicFaceField"

Definition at line 6 of file Boundaries.c.

__FUNCT__ [25/32]

#define __FUNCT__   "ApplyMetricsPeriodicBCs"

Definition at line 6 of file Boundaries.c.

__FUNCT__ [26/32]

#define __FUNCT__   "ApplyPeriodicBCs"

Definition at line 6 of file Boundaries.c.

__FUNCT__ [27/32]

#define __FUNCT__   "UpdateDummyCells"

Definition at line 6 of file Boundaries.c.

__FUNCT__ [28/32]

#define __FUNCT__   "UpdateCornerNodes"

Definition at line 6 of file Boundaries.c.

__FUNCT__ [29/32]

#define __FUNCT__   "UpdatePeriodicCornerNodes"

Definition at line 6 of file Boundaries.c.

__FUNCT__ [30/32]

#define __FUNCT__   "ApplyWallFunction"

Definition at line 6 of file Boundaries.c.

__FUNCT__ [31/32]

#define __FUNCT__   "FinalizePostProjectionCellFields"

Definition at line 6 of file Boundaries.c.

__FUNCT__ [32/32]

#define __FUNCT__   "ApplyBoundaryConditions"

Definition at line 6 of file Boundaries.c.

Function Documentation

CanRankServiceInletFace()

PetscErrorCode CanRankServiceInletFace	(	UserCtx *	user,
		const DMDALocalInfo *	info,
		PetscInt	IM_nodes_global,
		PetscInt	JM_nodes_global,
		PetscInt	KM_nodes_global,
		PetscBool *	can_service_inlet_out
	)

Internal helper implementation: CanRankServiceInletFace().

Determines if the current MPI rank owns any part of the globally defined inlet face, making it responsible for placing particles on that portion of the surface.

Local to this translation unit.

Definition at line 11 of file Boundaries.c.

{
    PetscErrorCode ierr;
    PetscMPIInt    rank_for_logging; // For detailed debugging logs
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank_for_logging); CHKERRQ(ierr);
 
    *can_service_inlet_out = PETSC_FALSE; // Default to no service
 
    if (!user->inletFaceDefined) {
        LOG_ALLOW(LOCAL, LOG_DEBUG, "[Rank %d]: Inlet face not defined in user context. Cannot service.\n", rank_for_logging);
        PetscFunctionReturn(0);
    }
 
    // Get the range of cells owned by this rank in each dimension
    PetscInt owned_start_cell_i, num_owned_cells_on_rank_i;
    PetscInt owned_start_cell_j, num_owned_cells_on_rank_j;
    PetscInt owned_start_cell_k, num_owned_cells_on_rank_k;
 
    ierr = GetOwnedCellRange(info, 0, &owned_start_cell_i, &num_owned_cells_on_rank_i); CHKERRQ(ierr);
    ierr = GetOwnedCellRange(info, 1, &owned_start_cell_j, &num_owned_cells_on_rank_j); CHKERRQ(ierr);
    ierr = GetOwnedCellRange(info, 2, &owned_start_cell_k, &num_owned_cells_on_rank_k); CHKERRQ(ierr);
 
    // Determine the global index of the last cell (0-indexed) in each direction.
    // Example: If IM_nodes_global = 11 (nodes 0-10), there are 10 cells (0-9). Last cell index is 9.
    // Formula: global_nodes - 1 (num cells) - 1 (0-indexed) = global_nodes - 2.
    PetscInt last_global_cell_idx_i = (IM_nodes_global > 1) ? (IM_nodes_global - 2) : -1; // -1 if 0 or 1 node (i.e., 0 cells)
    PetscInt last_global_cell_idx_j = (JM_nodes_global > 1) ? (JM_nodes_global - 2) : -1;
    PetscInt last_global_cell_idx_k = (KM_nodes_global > 1) ? (KM_nodes_global - 2) : -1;
 
    switch (user->identifiedInletBCFace) {
        case BC_FACE_NEG_X: // Inlet on the global I-minimum face (face of cell C_i=0)
            // Rank services if its first owned node is global node 0 (info->xs == 0),
            // and it owns cells in I, J, and K directions.
            if (info->xs == 0 && num_owned_cells_on_rank_i > 0 &&
                num_owned_cells_on_rank_j > 0 && num_owned_cells_on_rank_k > 0) {
                *can_service_inlet_out = PETSC_TRUE;
            }
            break;
        case BC_FACE_POS_X: // Inlet on the global I-maximum face (face of cell C_i=last_global_cell_idx_i)
            // Rank services if it owns the last cell in I-direction,
            // and has extent in J and K.
            if (last_global_cell_idx_i >= 0 && /* Check for valid global domain */
                (owned_start_cell_i + num_owned_cells_on_rank_i - 1) == last_global_cell_idx_i && /* Rank's last cell is the global last cell */
                num_owned_cells_on_rank_j > 0 && num_owned_cells_on_rank_k > 0) {
                *can_service_inlet_out = PETSC_TRUE;
            }
            break;
        case BC_FACE_NEG_Y:
            if (info->ys == 0 && num_owned_cells_on_rank_j > 0 &&
                num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_k > 0) {
                *can_service_inlet_out = PETSC_TRUE;
            }
            break;
        case BC_FACE_POS_Y:
            if (last_global_cell_idx_j >= 0 &&
                (owned_start_cell_j + num_owned_cells_on_rank_j - 1) == last_global_cell_idx_j &&
                num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_k > 0) {
                *can_service_inlet_out = PETSC_TRUE;
            }
            break;
        case BC_FACE_NEG_Z:
            if (info->zs == 0 && num_owned_cells_on_rank_k > 0 &&
                num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_j > 0) {
                *can_service_inlet_out = PETSC_TRUE;
            }
            break;
        case BC_FACE_POS_Z:
            if (last_global_cell_idx_k >= 0 &&
                (owned_start_cell_k + num_owned_cells_on_rank_k - 1) == last_global_cell_idx_k &&
                num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_j > 0) {
                *can_service_inlet_out = PETSC_TRUE;
            }
            break;
        default:
             LOG_ALLOW(LOCAL, LOG_WARNING, "[Rank %d]: Unknown inlet face %s.\n", rank_for_logging, BCFaceToString((BCFace)user->identifiedInletBCFace));
            break;
    }
 
      LOG_ALLOW(LOCAL, LOG_TRACE,
      "[Rank %d] Check Service for Inlet %s:\n"
      "    - Local Domain: starts at cell (%d,%d,%d), has (%d,%d,%d) cells.\n"
      "    - Global Domain: has (%d,%d,%d) nodes, so last cell is (%d,%d,%d).\n",
      rank_for_logging,
      BCFaceToString((BCFace)user->identifiedInletBCFace),
      owned_start_cell_i, owned_start_cell_j, owned_start_cell_k,
      num_owned_cells_on_rank_i, num_owned_cells_on_rank_j, num_owned_cells_on_rank_k,
      IM_nodes_global, JM_nodes_global, KM_nodes_global,
      last_global_cell_idx_i, last_global_cell_idx_j, last_global_cell_idx_k);
 
      LOG_ALLOW(LOCAL, LOG_INFO,"[Rank %d] Inlet Face %s Service Check Result: %s | Owned Cells (I,J,K): (%d,%d,%d) | Starts at Cell (%d,%d,%d)\n",
                rank_for_logging,
                BCFaceToString((BCFace)user->identifiedInletBCFace),
                (*can_service_inlet_out) ? "CAN SERVICE" : "CANNOT SERVICE",
                num_owned_cells_on_rank_i, num_owned_cells_on_rank_j, num_owned_cells_on_rank_k,
                owned_start_cell_i, owned_start_cell_j, owned_start_cell_k);
 
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

CanRankServiceFace()

PetscErrorCode CanRankServiceFace	(	const DMDALocalInfo *	info,
		PetscInt	IM_nodes_global,
		PetscInt	JM_nodes_global,
		PetscInt	KM_nodes_global,
		BCFace	face_id,
		PetscBool *	can_service_out
	)

Implementation of CanRankServiceFace().

Determines if the current MPI rank owns any part of a specified global face.

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

See also

CanRankServiceFace()

Definition at line 126 of file Boundaries.c.

{
    PetscErrorCode ierr;
    PetscMPIInt    rank_for_logging;
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank_for_logging); CHKERRQ(ierr);
 
    *can_service_out = PETSC_FALSE; // Default to no service
 
    // Get the range of cells owned by this rank
    PetscInt owned_start_cell_i, num_owned_cells_on_rank_i;
    PetscInt owned_start_cell_j, num_owned_cells_on_rank_j;
    PetscInt owned_start_cell_k, num_owned_cells_on_rank_k;
    ierr = GetOwnedCellRange(info, 0, &owned_start_cell_i, &num_owned_cells_on_rank_i); CHKERRQ(ierr);
    ierr = GetOwnedCellRange(info, 1, &owned_start_cell_j, &num_owned_cells_on_rank_j); CHKERRQ(ierr);
    ierr = GetOwnedCellRange(info, 2, &owned_start_cell_k, &num_owned_cells_on_rank_k); CHKERRQ(ierr);
 
    // Determine the global index of the last cell (0-indexed) in each direction.
    PetscInt last_global_cell_idx_i = (IM_nodes_global > 1) ? (IM_nodes_global - 2) : -1;
    PetscInt last_global_cell_idx_j = (JM_nodes_global > 1) ? (JM_nodes_global - 2) : -1;
    PetscInt last_global_cell_idx_k = (KM_nodes_global > 1) ? (KM_nodes_global - 2) : -1;
 
    switch (face_id) {
        case BC_FACE_NEG_X:
            if (info->xs == 0 && num_owned_cells_on_rank_i > 0 &&
                num_owned_cells_on_rank_j > 0 && num_owned_cells_on_rank_k > 0) {
                *can_service_out = PETSC_TRUE;
            }
            break;
        case BC_FACE_POS_X:
            if (last_global_cell_idx_i >= 0 &&
                (owned_start_cell_i + num_owned_cells_on_rank_i - 1) == last_global_cell_idx_i &&
                num_owned_cells_on_rank_j > 0 && num_owned_cells_on_rank_k > 0) {
                *can_service_out = PETSC_TRUE;
            }
            break;
        case BC_FACE_NEG_Y:
            if (info->ys == 0 && num_owned_cells_on_rank_j > 0 &&
                num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_k > 0) {
                *can_service_out = PETSC_TRUE;
            }
            break;
        case BC_FACE_POS_Y:
            if (last_global_cell_idx_j >= 0 &&
                (owned_start_cell_j + num_owned_cells_on_rank_j - 1) == last_global_cell_idx_j &&
                num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_k > 0) {
                *can_service_out = PETSC_TRUE;
            }
            break;
        case BC_FACE_NEG_Z:
            if (info->zs == 0 && num_owned_cells_on_rank_k > 0 &&
                num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_j > 0) {
                *can_service_out = PETSC_TRUE;
            }
            break;
        case BC_FACE_POS_Z:
            if (last_global_cell_idx_k >= 0 &&
                (owned_start_cell_k + num_owned_cells_on_rank_k - 1) == last_global_cell_idx_k &&
                num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_j > 0) {
                *can_service_out = PETSC_TRUE;
            }
            break;
        default:
             LOG_ALLOW(LOCAL, LOG_WARNING, "Rank %d: Unknown face enum %d. \n", rank_for_logging, face_id);
            break;
    }
 
    LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d check for face %s: Result=%s. \n",
        rank_for_logging, BCFaceToString((BCFace)face_id), (*can_service_out ? "TRUE" : "FALSE"));
 
    PROFILE_FUNCTION_END;
        
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

GetDeterministicFaceGridLocation()

PetscErrorCode GetDeterministicFaceGridLocation	(	UserCtx *	user,
		const DMDALocalInfo *	info,
		PetscInt	xs_gnode_rank,
		PetscInt	ys_gnode_rank,
		PetscInt	zs_gnode_rank,
		PetscInt	IM_cells_global,
		PetscInt	JM_cells_global,
		PetscInt	KM_cells_global,
		PetscInt64	particle_global_id,
		PetscInt *	ci_metric_lnode_out,
		PetscInt *	cj_metric_lnode_out,
		PetscInt *	ck_metric_lnode_out,
		PetscReal *	xi_metric_logic_out,
		PetscReal *	eta_metric_logic_out,
		PetscReal *	zta_metric_logic_out,
		PetscBool *	placement_successful_out
	)

Internal helper implementation: GetDeterministicFaceGridLocation().

Places particles in a deterministic grid/raster pattern on a specified domain face.

Local to this translation unit.

Definition at line 212 of file Boundaries.c.

{
    SimCtx *simCtx = user->simCtx;
    PetscReal global_logic_i = 0.0, global_logic_j = 0.0, global_logic_k = 0.0;
    PetscErrorCode ierr;
    PetscMPIInt rank_for_logging;
 
    PetscFunctionBeginUser;
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank_for_logging); CHKERRQ(ierr);
 
    *placement_successful_out = PETSC_FALSE; // Default to failure
 
    // --- Step 1: Configuration and Input Validation ---
 
    // *** Hardcoded number of grid layers. Change this value to alter the pattern. ***
    const PetscInt grid_layers = 2;
 
    LOG_ALLOW(LOCAL, LOG_DEBUG,
        "[Rank %d] Placing particle %lld on face %s with grid_layers=%d in global domain (%d,%d,%d) cells.\n",
        rank_for_logging, (long long)particle_global_id, BCFaceToString(user->identifiedInletBCFace), grid_layers,
        IM_cells_global, JM_cells_global, KM_cells_global);
 
    const char *face_name = BCFaceToString(user->identifiedInletBCFace);
 
    // Fatal Error Checks: Ensure the requested grid is geometrically possible.
    // The total layers from opposite faces (2 * grid_layers) must be less than the domain size.
    switch (user->identifiedInletBCFace) {
        case BC_FACE_NEG_X: case BC_FACE_POS_X:
            if (JM_cells_global <= 1 || KM_cells_global <= 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot place grid on face %s for a 2D/1D domain (J-cells=%d, K-cells=%d).", face_name, JM_cells_global, KM_cells_global);
            if (2 * grid_layers >= JM_cells_global || 2 * grid_layers >= KM_cells_global) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Grid layers (%d) from opposing J/K faces would overlap in this domain (J-cells=%d, K-cells=%d).", grid_layers, JM_cells_global, KM_cells_global);
            break;
        case BC_FACE_NEG_Y: case BC_FACE_POS_Y:
            if (IM_cells_global <= 1 || KM_cells_global <= 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot place grid on face %s for a 2D/1D domain (I-cells=%d, K-cells=%d).", face_name, IM_cells_global, KM_cells_global);
            if (2 * grid_layers >= IM_cells_global || 2 * grid_layers >= KM_cells_global) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Grid layers (%d) from opposing I/K faces would overlap in this domain (I-cells=%d, K-cells=%d).", grid_layers, IM_cells_global, KM_cells_global);
            break;
        case BC_FACE_NEG_Z: case BC_FACE_POS_Z:
            if (IM_cells_global <= 1 || JM_cells_global <= 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot place grid on face %s for a 2D/1D domain (I-cells=%d, J-cells=%d).", face_name, IM_cells_global, JM_cells_global);
            if (2 * grid_layers >= IM_cells_global || 2 * grid_layers >= JM_cells_global) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Grid layers (%d) from opposing I/J faces would overlap in this domain (I-cells=%d, J-cells=%d).", grid_layers, IM_cells_global, JM_cells_global);
            break;
        default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid identifiedInletBCFace specified: %d", user->identifiedInletBCFace);
    }
 
    const PetscInt num_lines_total = 4 * grid_layers;
    if (simCtx->np < num_lines_total) {
        LOG_ALLOW(GLOBAL, LOG_WARNING, "Warning: Total particle count (%lld) is less than the number of grid lines requested (%d). Some lines may be empty.\n", (long long)simCtx->np, num_lines_total);
    }
    if (simCtx->np > 0 && simCtx->np % num_lines_total != 0) {
        LOG_ALLOW(GLOBAL, LOG_WARNING, "Warning: Total particle count (%lld) is not evenly divisible by the number of grid lines (%d). Distribution will be uneven.\n", (long long)simCtx->np, num_lines_total);
    }
 
    // --- Step 2: Map global particle ID to a line and a point on that line ---
    if (simCtx->np == 0) PetscFunctionReturn(0); // Nothing to do
 
    LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Distributing %lld particles over %d lines on face %s.\n",
        rank_for_logging, (long long)simCtx->np, num_lines_total, face_name);
 
    const PetscInt points_per_line = PetscMax(1, simCtx->np / num_lines_total);
    PetscInt line_index = particle_global_id / points_per_line;
    PetscInt point_index_on_line = particle_global_id % points_per_line;
    line_index = PetscMin(line_index, num_lines_total - 1); // Clamp to handle uneven division
 
    // Decode the line_index into an edge group (0-3) and a layer within that group (0 to grid_layers-1)
    const PetscInt edge_group = line_index / grid_layers;
    const PetscInt layer_index = line_index % grid_layers;
 
    // --- Step 3: Calculate placement coordinates based on the decoded indices ---
    const PetscReal layer_spacing_norm_i = (IM_cells_global > 0) ? 1.0 / (PetscReal)IM_cells_global : 0.0;
    const PetscReal layer_spacing_norm_j = (JM_cells_global > 0) ? 1.0 / (PetscReal)JM_cells_global : 0.0;
    const PetscReal layer_spacing_norm_k = (KM_cells_global > 0) ? 1.0 / (PetscReal)KM_cells_global : 0.0;
 
    // Grid-aware epsilon: scale with minimum cell size to keep particles away from rank boundaries
    const PetscReal min_layer_spacing = PetscMin(layer_spacing_norm_i, PetscMin(layer_spacing_norm_j, layer_spacing_norm_k));
    const PetscReal epsilon = 0.5 * min_layer_spacing;  // Keep particles 10% of cell width from boundaries
 
    PetscReal variable_coord; // The coordinate that varies along a line
    if (points_per_line <= 1) {
        variable_coord = 0.5; // Place single point in the middle
    } else {
        variable_coord = ((PetscReal)point_index_on_line + 0.5)/ (PetscReal)(points_per_line);
    }
    variable_coord = PetscMin(1.0 - epsilon, PetscMax(epsilon, variable_coord)); // Clamp within [eps, 1-eps]
 
    // Main logic switch to determine the three global logical coordinates
    switch (user->identifiedInletBCFace) {
        case BC_FACE_NEG_X:
            global_logic_i = 0.5 * layer_spacing_norm_i; // Place near the face, in the middle of the first cell
            if (edge_group == 0)      { global_logic_j = (PetscReal)layer_index * layer_spacing_norm_j + epsilon; global_logic_k = variable_coord; }
            else if (edge_group == 1) { global_logic_j = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_j) - epsilon; global_logic_k = variable_coord; }
            else if (edge_group == 2) { global_logic_k = (PetscReal)layer_index * layer_spacing_norm_k + epsilon; global_logic_j = variable_coord; }
            else /* edge_group == 3 */ { global_logic_k = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_k) - epsilon; global_logic_j = variable_coord; }
            break;
        case BC_FACE_POS_X:
            global_logic_i = 1.0 - (0.5 * layer_spacing_norm_i); // Place near the face, in the middle of the last cell
            if (edge_group == 0)      { global_logic_j = (PetscReal)layer_index * layer_spacing_norm_j + epsilon; global_logic_k = variable_coord; }
            else if (edge_group == 1) { global_logic_j = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_j) - epsilon; global_logic_k = variable_coord; }
            else if (edge_group == 2) { global_logic_k = (PetscReal)layer_index * layer_spacing_norm_k + epsilon; global_logic_j = variable_coord; }
            else /* edge_group == 3 */ { global_logic_k = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_k) - epsilon; global_logic_j = variable_coord; }
            break;
        case BC_FACE_NEG_Y:
            global_logic_j = 0.5 * layer_spacing_norm_j;
            if (edge_group == 0)      { global_logic_i = (PetscReal)layer_index * layer_spacing_norm_i + epsilon; global_logic_k = variable_coord; }
            else if (edge_group == 1) { global_logic_i = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_i) - epsilon; global_logic_k = variable_coord; }
            else if (edge_group == 2) { global_logic_k = (PetscReal)layer_index * layer_spacing_norm_k + epsilon; global_logic_i = variable_coord; }
            else /* edge_group == 3 */ { global_logic_k = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_k) - epsilon; global_logic_i = variable_coord; }
            break;
        case BC_FACE_POS_Y:
            global_logic_j = 1.0 - (0.5 * layer_spacing_norm_j);
            if (edge_group == 0)      { global_logic_i = (PetscReal)layer_index * layer_spacing_norm_i + epsilon; global_logic_k = variable_coord; }
            else if (edge_group == 1) { global_logic_i = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_i) - epsilon; global_logic_k = variable_coord; }
            else if (edge_group == 2) { global_logic_k = (PetscReal)layer_index * layer_spacing_norm_k + epsilon; global_logic_i = variable_coord; }
            else /* edge_group == 3 */ { global_logic_k = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_k) - epsilon; global_logic_i = variable_coord; }
            break;
        case BC_FACE_NEG_Z:
            global_logic_k = 0.5 * layer_spacing_norm_k;
            if (edge_group == 0)      { global_logic_i = (PetscReal)layer_index * layer_spacing_norm_i + epsilon; global_logic_j = variable_coord; }
            else if (edge_group == 1) { global_logic_i = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_i) - epsilon; global_logic_j = variable_coord; }
            else if (edge_group == 2) { global_logic_j = (PetscReal)layer_index * layer_spacing_norm_j + epsilon; global_logic_i = variable_coord; }
            else /* edge_group == 3 */ { global_logic_j = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_j) - epsilon; global_logic_i = variable_coord; }
            break;
        case BC_FACE_POS_Z:
            global_logic_k = 1.0 - (0.5 * layer_spacing_norm_k);
            if (edge_group == 0)      { global_logic_i = (PetscReal)layer_index * layer_spacing_norm_i + epsilon; global_logic_j = variable_coord; }
            else if (edge_group == 1) { global_logic_i = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_i) - epsilon; global_logic_j = variable_coord; }
            else if (edge_group == 2) { global_logic_j = (PetscReal)layer_index * layer_spacing_norm_j + epsilon; global_logic_i = variable_coord; }
            else /* edge_group == 3 */ { global_logic_j = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_j) - epsilon; global_logic_i = variable_coord; }
            break;
    }
 
    LOG_ALLOW(LOCAL, LOG_TRACE,
        "[Rank %d] Particle %lld assigned to line %d (edge group %d, layer %d) with variable_coord=%.4f.\n"
        "    -> Global logical coords: (i,j,k) = (%.6f, %.6f, %.6f)\n",
        rank_for_logging, (long long)particle_global_id, line_index, edge_group, layer_index, variable_coord,
        global_logic_i, global_logic_j, global_logic_k);
 
    // --- Step 4: Convert global logical coordinate to global cell index and intra-cell logicals ---
    PetscReal global_cell_coord_i = global_logic_i * IM_cells_global;
    PetscInt  I_g = (PetscInt)global_cell_coord_i;
    *xi_metric_logic_out = global_cell_coord_i - I_g;
 
    PetscReal global_cell_coord_j = global_logic_j * JM_cells_global;
    PetscInt  J_g = (PetscInt)global_cell_coord_j;
    *eta_metric_logic_out = global_cell_coord_j - J_g;
 
    PetscReal global_cell_coord_k = global_logic_k * KM_cells_global;
    PetscInt  K_g = (PetscInt)global_cell_coord_k;
    *zta_metric_logic_out = global_cell_coord_k - K_g;
 
    // --- Step 5: Check if this rank owns the target cell and finalize outputs ---
    if ((I_g >= info->xs && I_g < info->xs + info->xm) &&
        (J_g >= info->ys && J_g < info->ys + info->ym) &&
        (K_g >= info->zs && K_g < info->zs + info->zm))
    {
        // Convert global cell index to the local node index for this rank's DA patch
        *ci_metric_lnode_out = (I_g - info->xs) + xs_gnode_rank;
        *cj_metric_lnode_out = (J_g - info->ys) + ys_gnode_rank;
        *ck_metric_lnode_out = (K_g - info->zs) + zs_gnode_rank;
        *placement_successful_out = PETSC_TRUE;
    }
 
    LOG_ALLOW(LOCAL, LOG_VERBOSE,
        "[Rank %d] Particle %lld placement %s.\n",
        rank_for_logging, (long long)particle_global_id,
        (*placement_successful_out ? "SUCCESSFUL" : "NOT ON THIS RANK"));
 
    if(*placement_successful_out){    
        LOG_ALLOW(LOCAL,LOG_TRACE,"Local cell origin node: (I,J,K) = (%d,%d,%d), intra-cell logicals: (xi,eta,zta)=(%.6f,%.6f,%.6f)\n",
                    *ci_metric_lnode_out, *cj_metric_lnode_out, *ck_metric_lnode_out,
                    *xi_metric_logic_out, *eta_metric_logic_out, *zta_metric_logic_out);
        }
        
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

GetRandomCellAndLogicalCoordsOnInletFace()

PetscErrorCode GetRandomCellAndLogicalCoordsOnInletFace	(	UserCtx *	user,
		const DMDALocalInfo *	info,
		PetscInt	xs_gnode_rank,
		PetscInt	ys_gnode_rank,
		PetscInt	zs_gnode_rank,
		PetscInt	IM_nodes_global,
		PetscInt	JM_nodes_global,
		PetscInt	KM_nodes_global,
		PetscRandom *	rand_logic_i_ptr,
		PetscRandom *	rand_logic_j_ptr,
		PetscRandom *	rand_logic_k_ptr,
		PetscInt *	ci_metric_lnode_out,
		PetscInt *	cj_metric_lnode_out,
		PetscInt *	ck_metric_lnode_out,
		PetscReal *	xi_metric_logic_out,
		PetscReal *	eta_metric_logic_out,
		PetscReal *	zta_metric_logic_out
	)

Internal helper implementation: GetRandomCellAndLogicalCoordsOnInletFace().

Assuming the current rank services the inlet face, this function selects a random cell (owned by this rank on that face) and random logical coordinates within that cell, suitable for placing a particle on the inlet surface.

Local to this translation unit.

Definition at line 399 of file Boundaries.c.

{
    PetscErrorCode ierr = 0;
    PetscReal r_val_i_sel, r_val_j_sel, r_val_k_sel;
    PetscInt local_cell_idx_on_face_dim1 = 0; // 0-indexed relative to owned cells on face
    PetscInt local_cell_idx_on_face_dim2 = 0;
    PetscMPIInt rank_for_logging; 
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank_for_logging); CHKERRQ(ierr);
 
    // Get number of cells this rank owns in each dimension (tangential to the face mainly)
    PetscInt owned_start_cell_i, num_owned_cells_on_rank_i;
    PetscInt owned_start_cell_j, num_owned_cells_on_rank_j;
    PetscInt owned_start_cell_k, num_owned_cells_on_rank_k;
 
    ierr = GetOwnedCellRange(info, 0, &owned_start_cell_i, &num_owned_cells_on_rank_i); CHKERRQ(ierr);
    ierr = GetOwnedCellRange(info, 1, &owned_start_cell_j, &num_owned_cells_on_rank_j); CHKERRQ(ierr);
    ierr = GetOwnedCellRange(info, 2, &owned_start_cell_k, &num_owned_cells_on_rank_k); CHKERRQ(ierr);
 
    // Defaults for cell origin node (local index for the rank's DA patch, including ghosts)
    *ci_metric_lnode_out = xs_gnode_rank; *cj_metric_lnode_out = ys_gnode_rank; *ck_metric_lnode_out = zs_gnode_rank;
    // Defaults for logical coordinates
    *xi_metric_logic_out = 0.5; *eta_metric_logic_out = 0.5; *zta_metric_logic_out = 0.5;
 
    // Index of the last cell (0-indexed) in each global direction
    PetscInt last_global_cell_idx_i = (IM_nodes_global > 1) ? (IM_nodes_global - 2) : -1;
    PetscInt last_global_cell_idx_j = (JM_nodes_global > 1) ? (JM_nodes_global - 2) : -1;
    PetscInt last_global_cell_idx_k = (KM_nodes_global > 1) ? (KM_nodes_global - 2) : -1;
    
    LOG_ALLOW(LOCAL, LOG_INFO, "PARTICLE_INIT_DEBUG Rank %d: Inlet face %s.\n"
        "  Owned cells (i,j,k): (%d,%d,%d)\n"
        "  Global nodes (I,J,K): (%d,%d,%d)\n"
        "  info->xs,ys,zs (first owned node GLOBAL): (%d,%d,%d)\n"
        "  info->xm,ym,zm (num owned nodes GLOBAL): (%d,%d,%d)\n"
        "  xs_gnode_rank,ys_gnode_rank,zs_gnode_rank (DMDAGetCorners): (%d,%d,%d)\n"
        "  owned_start_cell (i,j,k) GLOBAL: (%d,%d,%d)\n"
        "  last_global_cell_idx (i,j,k): (%d,%d,%d)\n",
        rank_for_logging, BCFaceToString((BCFace)user->identifiedInletBCFace),
        num_owned_cells_on_rank_i,num_owned_cells_on_rank_j,num_owned_cells_on_rank_k,
        IM_nodes_global,JM_nodes_global,KM_nodes_global,
        info->xs, info->ys, info->zs,
        info->xm, info->ym, info->zm,
        xs_gnode_rank,ys_gnode_rank,zs_gnode_rank,
        owned_start_cell_i, owned_start_cell_j, owned_start_cell_k,
        last_global_cell_idx_i, last_global_cell_idx_j, last_global_cell_idx_k);
 
 
    switch (user->identifiedInletBCFace) {
        case BC_FACE_NEG_X: // Particle on -X face of cell C_0 (origin node N_0)
            // Cell origin node is the first owned node in I by this rank (global index info->xs).
            // Its local index within the rank's DA (incl ghosts) is xs_gnode_rank.
            *ci_metric_lnode_out = xs_gnode_rank;
            *xi_metric_logic_out = 1.0e-6;
 
            // Tangential dimensions are J and K. Select an owned cell randomly on this face.
            // num_owned_cells_on_rank_j/k must be > 0 (checked by CanRankServiceInletFace)
            ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, &r_val_j_sel); CHKERRQ(ierr);
            local_cell_idx_on_face_dim1 = (PetscInt)(r_val_j_sel * num_owned_cells_on_rank_j); // Index among owned J-cells
            local_cell_idx_on_face_dim1 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim1), num_owned_cells_on_rank_j - 1);
            *cj_metric_lnode_out = ys_gnode_rank + local_cell_idx_on_face_dim1; // Offset from start of rank's J-nodes
 
            ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, &r_val_k_sel); CHKERRQ(ierr);
            local_cell_idx_on_face_dim2 = (PetscInt)(r_val_k_sel * num_owned_cells_on_rank_k);
            local_cell_idx_on_face_dim2 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim2), num_owned_cells_on_rank_k - 1);
            *ck_metric_lnode_out = zs_gnode_rank + local_cell_idx_on_face_dim2;
 
            ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, eta_metric_logic_out); CHKERRQ(ierr);
            ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, zta_metric_logic_out); CHKERRQ(ierr);
            break;
 
        case BC_FACE_POS_X: // Particle on +X face of cell C_last_I (origin node N_last_I_origin)
            // Origin node of the last I-cell is global_node_idx = last_global_cell_idx_i.
            // Its local index in rank's DA: (last_global_cell_idx_i - info->xs) + xs_gnode_rank
            *ci_metric_lnode_out = xs_gnode_rank + (last_global_cell_idx_i - info->xs);
            *xi_metric_logic_out = 1.0 - 1.0e-6;
 
            ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, &r_val_j_sel); CHKERRQ(ierr);
            local_cell_idx_on_face_dim1 = (PetscInt)(r_val_j_sel * num_owned_cells_on_rank_j);
            local_cell_idx_on_face_dim1 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim1), num_owned_cells_on_rank_j - 1);
            *cj_metric_lnode_out = ys_gnode_rank + local_cell_idx_on_face_dim1;
 
            ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, &r_val_k_sel); CHKERRQ(ierr);
            local_cell_idx_on_face_dim2 = (PetscInt)(r_val_k_sel * num_owned_cells_on_rank_k);
            local_cell_idx_on_face_dim2 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim2), num_owned_cells_on_rank_k - 1);
            *ck_metric_lnode_out = zs_gnode_rank + local_cell_idx_on_face_dim2;
 
            ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, eta_metric_logic_out); CHKERRQ(ierr);
            ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, zta_metric_logic_out); CHKERRQ(ierr);
            break;
        // ... (Cases for Y and Z faces, following the same pattern) ...
        case BC_FACE_NEG_Y:
            *cj_metric_lnode_out = ys_gnode_rank;
            *eta_metric_logic_out = 1.0e-6;
            ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, &r_val_i_sel); CHKERRQ(ierr);
            local_cell_idx_on_face_dim1 = (PetscInt)(r_val_i_sel * num_owned_cells_on_rank_i);
            local_cell_idx_on_face_dim1 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim1), num_owned_cells_on_rank_i - 1);
            *ci_metric_lnode_out = xs_gnode_rank + local_cell_idx_on_face_dim1;
            ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, &r_val_k_sel); CHKERRQ(ierr);
            local_cell_idx_on_face_dim2 = (PetscInt)(r_val_k_sel * num_owned_cells_on_rank_k);
            local_cell_idx_on_face_dim2 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim2), num_owned_cells_on_rank_k - 1);
            *ck_metric_lnode_out = zs_gnode_rank + local_cell_idx_on_face_dim2;
            ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, xi_metric_logic_out); CHKERRQ(ierr);
            ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, zta_metric_logic_out); CHKERRQ(ierr);
            break;
        case BC_FACE_POS_Y:
            *cj_metric_lnode_out = ys_gnode_rank + (last_global_cell_idx_j - info->ys);
            *eta_metric_logic_out = 1.0 - 1.0e-6;
            ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, &r_val_i_sel); CHKERRQ(ierr);
            local_cell_idx_on_face_dim1 = (PetscInt)(r_val_i_sel * num_owned_cells_on_rank_i);
            local_cell_idx_on_face_dim1 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim1), num_owned_cells_on_rank_i - 1);
            *ci_metric_lnode_out = xs_gnode_rank + local_cell_idx_on_face_dim1;
            ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, &r_val_k_sel); CHKERRQ(ierr);
            local_cell_idx_on_face_dim2 = (PetscInt)(r_val_k_sel * num_owned_cells_on_rank_k);
            local_cell_idx_on_face_dim2 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim2), num_owned_cells_on_rank_k - 1);
            *ck_metric_lnode_out = zs_gnode_rank + local_cell_idx_on_face_dim2;
            ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, xi_metric_logic_out); CHKERRQ(ierr);
            ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, zta_metric_logic_out); CHKERRQ(ierr);
            break;
        case BC_FACE_NEG_Z: // Your example case
            *ck_metric_lnode_out = zs_gnode_rank; // Cell origin is the first owned node in K by this rank
            *zta_metric_logic_out = 1.0e-6;      // Place particle slightly inside this cell from its -Z face
            // Tangential dimensions are I and J
            ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, &r_val_i_sel); CHKERRQ(ierr);
            local_cell_idx_on_face_dim1 = (PetscInt)(r_val_i_sel * num_owned_cells_on_rank_i);
            local_cell_idx_on_face_dim1 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim1), num_owned_cells_on_rank_i - 1);
            *ci_metric_lnode_out = xs_gnode_rank + local_cell_idx_on_face_dim1;
 
            ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, &r_val_j_sel); CHKERRQ(ierr);
            local_cell_idx_on_face_dim2 = (PetscInt)(r_val_j_sel * num_owned_cells_on_rank_j);
            local_cell_idx_on_face_dim2 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim2), num_owned_cells_on_rank_j - 1);
            *cj_metric_lnode_out = ys_gnode_rank + local_cell_idx_on_face_dim2;
 
            ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, xi_metric_logic_out); CHKERRQ(ierr); // Intra-cell logical for I
            ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, eta_metric_logic_out); CHKERRQ(ierr); // Intra-cell logical for J
            break;
        case BC_FACE_POS_Z:
            *ck_metric_lnode_out = zs_gnode_rank + (last_global_cell_idx_k - info->zs);
            *zta_metric_logic_out = 1.0 - 1.0e-6;
            ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, &r_val_i_sel); CHKERRQ(ierr);
            local_cell_idx_on_face_dim1 = (PetscInt)(r_val_i_sel * num_owned_cells_on_rank_i);
            local_cell_idx_on_face_dim1 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim1), num_owned_cells_on_rank_i - 1);
            *ci_metric_lnode_out = xs_gnode_rank + local_cell_idx_on_face_dim1;
            ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, &r_val_j_sel); CHKERRQ(ierr);
            local_cell_idx_on_face_dim2 = (PetscInt)(r_val_j_sel * num_owned_cells_on_rank_j);
            local_cell_idx_on_face_dim2 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim2), num_owned_cells_on_rank_j - 1);
            *cj_metric_lnode_out = ys_gnode_rank + local_cell_idx_on_face_dim2;
            ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, xi_metric_logic_out); CHKERRQ(ierr);
            ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, eta_metric_logic_out); CHKERRQ(ierr);
            break;
        default:
             SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "GetRandomCellAndLogicOnInletFace: Invalid user->identifiedInletBCFace %d. \n", user->identifiedInletBCFace);
    }
 
    PetscReal eps = 1.0e-7;
    if (user->identifiedInletBCFace == BC_FACE_NEG_X || user->identifiedInletBCFace == BC_FACE_POS_X) {
        *eta_metric_logic_out = PetscMin(PetscMax(0.0, *eta_metric_logic_out), 1.0 - eps);
        *zta_metric_logic_out = PetscMin(PetscMax(0.0, *zta_metric_logic_out), 1.0 - eps);
    } else if (user->identifiedInletBCFace == BC_FACE_NEG_Y || user->identifiedInletBCFace == BC_FACE_POS_Y) {
        *xi_metric_logic_out  = PetscMin(PetscMax(0.0, *xi_metric_logic_out),  1.0 - eps);
        *zta_metric_logic_out = PetscMin(PetscMax(0.0, *zta_metric_logic_out), 1.0 - eps);
    } else { 
        *xi_metric_logic_out  = PetscMin(PetscMax(0.0, *xi_metric_logic_out),  1.0 - eps);
        *eta_metric_logic_out = PetscMin(PetscMax(0.0, *eta_metric_logic_out), 1.0 - eps);
    }
    
    LOG_ALLOW(LOCAL, LOG_VERBOSE, "Rank %d: Target Cell Node =(%d,%d,%d). (xi,et,zt)=(%.2e,%.2f,%.2f). \n",
        rank_for_logging, *ci_metric_lnode_out, *cj_metric_lnode_out, *ck_metric_lnode_out,
        *xi_metric_logic_out, *eta_metric_logic_out, *zta_metric_logic_out);
 
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

EnforceRHSBoundaryConditions()

PetscErrorCode EnforceRHSBoundaryConditions

(

UserCtx *

user

)

Internal helper implementation: EnforceRHSBoundaryConditions().

Enforces boundary conditions on the momentum equation's Right-Hand-Side (RHS) vector.

Local to this translation unit.

Definition at line 591 of file Boundaries.c.

{
    PetscErrorCode ierr;
    DMDALocalInfo  info = user->info;
    Cmpnts         ***rhs;
 
    // --- Grid extents for this MPI rank and global grid dimensions ---
    const PetscInt xs = info.xs, xe = xs + info.xm;
    const PetscInt ys = info.ys, ye = ys + info.ym;
    const PetscInt zs = info.zs, ze = zs + info.zm;
    const PetscInt mx = info.mx, my = info.my, mz = info.mz;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    // Get a writable pointer to the local data of the global RHS vector.
    ierr = DMDAVecGetArray(user->fda, user->Rhs, &rhs); CHKERRQ(ierr);
 
    // ========================================================================
    // --- I-DIRECTION (X-FACES) ---
    // ========================================================================
 
    // --- Negative X Face (i=0, the first physical face) ---
    if (xs == 0) {
        // This logic applies ONLY to physical (non-periodic) boundaries.
        if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC) {
            const PetscInt i = 0;
            for (PetscInt k = zs; k < ze; k++) {
                for (PetscInt j = ys; j < ye; j++) {
                    rhs[k][j][i].x = 0.0;
                    rhs[k][j][i].y = 0.0;
                    rhs[k][j][i].z = 0.0;
                }
            }
        }
    }
 
    // --- Positive X Face (physical face i=mx-2, dummy location i=mx-1) ---
    if (xe == mx) {
        // Step 1: Enforce strong BC on the LAST PHYSICAL face (i=mx-2) for non-periodic cases.
        if (user->boundary_faces[BC_FACE_POS_X].mathematical_type != PERIODIC) {
            const PetscInt i = mx - 2;
            for (PetscInt k = zs; k < ze; k++) {
                for (PetscInt j = ys; j < ye; j++) {
                    rhs[k][j][i].x = 0.0;
                }
            }
        }
        // Step 2: Unconditionally sanitize the DUMMY location (i=mx-1).
        const PetscInt i = mx - 1;
        for (PetscInt k = zs; k < ze; k++) {
            for (PetscInt j = ys; j < ye; j++) {
                rhs[k][j][i].x = 0.0;
                rhs[k][j][i].y = 0.0;
                rhs[k][j][i].z = 0.0;
            }
        }
    }
 
    // ========================================================================
    // --- J-DIRECTION (Y-FACES) ---
    // ========================================================================
 
    // --- Negative Y Face (j=0, the first physical face) ---
    if (ys == 0) {
        if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC) {
            const PetscInt j = 0;
            for (PetscInt k = zs; k < ze; k++) {
                for (PetscInt i = xs; i < xe; i++) {
                    rhs[k][j][i].x = 0.0;
                    rhs[k][j][i].y = 0.0;
                    rhs[k][j][i].z = 0.0;
                }
            }
        }
    }
 
    // --- Positive Y Face (physical face j=my-2, dummy location j=my-1) ---
    if (ye == my) {
        if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type != PERIODIC) {
            const PetscInt j = my - 2;
            for (PetscInt k = zs; k < ze; k++) {
                for (PetscInt i = xs; i < xe; i++) {
                    rhs[k][j][i].y = 0.0;
                }
            }
        }
        const PetscInt j = my - 1;
        for (PetscInt k = zs; k < ze; k++) {
            for (PetscInt i = xs; i < xe; i++) {
                rhs[k][j][i].x = 0.0;
                rhs[k][j][i].y = 0.0;
                rhs[k][j][i].z = 0.0;
            }
        }
    }
 
    // ========================================================================
    // --- K-DIRECTION (Z-FACES) ---
    // ========================================================================
 
    // --- Negative Z Face (k=0, the first physical face) ---
    if (zs == 0) {
        if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC) {
            const PetscInt k = 0;
            for (PetscInt j = ys; j < ye; j++) {
                for (PetscInt i = xs; i < xe; i++) {
                    rhs[k][j][i].x = 0.0;
                    rhs[k][j][i].y = 0.0;
                    rhs[k][j][i].z = 0.0;
                }
            }
        }
    }
 
    // --- Positive Z Face (physical face k=mz-2, dummy location k=mz-1) ---
    if (ze == mz) {
        if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type != PERIODIC) {
            const PetscInt k = mz - 2;
            for (PetscInt j = ys; j < ye; j++) {
                for (PetscInt i = xs; i < xe; i++) {
                    rhs[k][j][i].z = 0.0;
                }
            }
        }
        const PetscInt k = mz - 1;
        for (PetscInt j = ys; j < ye; j++) {
            for (PetscInt i = xs; i < xe; i++) {
                rhs[k][j][i].x = 0.0;
                rhs[k][j][i].y = 0.0;
                rhs[k][j][i].z = 0.0;
            }
        }
    }
 
    // --- Release the pointer to the local data ---
    ierr = DMDAVecRestoreArray(user->fda, user->Rhs, &rhs); CHKERRQ(ierr);
 
    LOG_ALLOW(LOCAL, LOG_TRACE, "Rank %d, Block %d: Finished enforcing RHS boundary conditions.\n",
              user->simCtx->rank, user->_this);
    
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

BoundaryCondition_Create()

PetscErrorCode BoundaryCondition_Create	(	BCHandlerType	handler_type,
		BoundaryCondition **	new_bc_ptr
	)

Internal helper implementation: BoundaryCondition_Create().

(Private) Creates and configures a specific BoundaryCondition handler object.

Local to this translation unit.

Definition at line 744 of file Boundaries.c.

{
    PetscErrorCode ierr;
    PetscFunctionBeginUser;
 
    const char* handler_name = BCHandlerTypeToString(handler_type);
    LOG_ALLOW(LOCAL, LOG_DEBUG, "Factory called for handler type %s. \n", handler_name);
 
    ierr = PetscMalloc1(1, new_bc_ptr); CHKERRQ(ierr);
    BoundaryCondition *bc = *new_bc_ptr;
 
    bc->type        = handler_type;
    bc->priority    = -1;  // Default priority; can be overridden in specific handlers
    bc->data        = NULL;
    bc->Initialize  = NULL;
    bc->PreStep     = NULL;
    bc->Apply       = NULL;
    bc->PostStep    = NULL;
    bc->UpdateUbcs  = NULL;
    bc->Destroy     = NULL;
 
    LOG_ALLOW(LOCAL, LOG_DEBUG, "Allocated generic handler object at address %p.\n", (void*)bc);
 
    switch (handler_type) {
 
        case BC_HANDLER_OUTLET_CONSERVATION:
            LOG_ALLOW(LOCAL, LOG_DEBUG, "Dispatching to Create_OutletConservation().\n");
            ierr = Create_OutletConservation(bc); CHKERRQ(ierr);
            break;
 
        case BC_HANDLER_WALL_NOSLIP:
            LOG_ALLOW(LOCAL, LOG_DEBUG, "Dispatching to Create_WallNoSlip().\n");
            ierr = Create_WallNoSlip(bc); CHKERRQ(ierr);
            break;
 
        case BC_HANDLER_INLET_CONSTANT_VELOCITY:
            LOG_ALLOW(LOCAL, LOG_DEBUG, "Dispatching to Create_InletConstantVelocity().\n");
              ierr = Create_InletConstantVelocity(bc); CHKERRQ(ierr);
            break;
 
        case BC_HANDLER_PERIODIC_GEOMETRIC:
            LOG_ALLOW(LOCAL,LOG_DEBUG,"Dispatching to Create_PeriodicGeometric().\n");
            ierr = Create_PeriodicGeometric(bc);
            break;
        
        case BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX:
            LOG_ALLOW(LOCAL,LOG_DEBUG,"Dispatching to Create_PeriodicDrivenConstant().\n");
            ierr = Create_PeriodicDrivenConstant(bc);
            break;
        
        //case BC_HANDLER_PERIODIC_DRIVEN_INITIAL_FLUX:
        //    LOG_ALLOW(LOCAL,LOG_DEBUG,"Dispatching to Create_PeriodicDrivenInitial().\n");
        //    ierr = Create_PeriodicDrivenInitial(bc);
        //    break;
                
        case BC_HANDLER_INLET_PARABOLIC:
            LOG_ALLOW(LOCAL, LOG_DEBUG, "Dispatching to Create_InletParabolicProfile().\n");
            ierr = Create_InletParabolicProfile(bc); CHKERRQ(ierr);
            break;
 
        case BC_HANDLER_INLET_PROFILE_FROM_FILE:
            LOG_ALLOW(LOCAL, LOG_DEBUG, "Dispatching to Create_InletProfileFromFile().\n");
            ierr = Create_InletProfileFromFile(bc); CHKERRQ(ierr);
            break;
        //Add cases for other handlers here in future phases 
        
        default:
            LOG_ALLOW(GLOBAL, LOG_ERROR, "Handler type (%s) is not recognized or implemented in the factory.\n", handler_name);
            SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Boundary handler type %d (%s) not recognized in factory.\n", handler_type, handler_name);
    }
 
    if(bc->priority < 0) {
        LOG_ALLOW(GLOBAL, LOG_ERROR, "Handler type %d (%s) did not set a valid priority during creation.\n", handler_type, handler_name);
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Boundary handler type %d (%s) did not set a valid priority during creation.\n", handler_type, handler_name);
    }
    
    LOG_ALLOW(LOCAL, LOG_DEBUG, "Successfully created and configured handler for %s.\n", handler_name);
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

BoundarySystem_Validate()

PetscErrorCode BoundarySystem_Validate

(

UserCtx *

user

)

Internal helper implementation: BoundarySystem_Validate().

(Public) Validates the consistency and compatibility of the parsed boundary condition system.

Local to this translation unit.

Definition at line 830 of file Boundaries.c.

{
    PetscErrorCode ierr;
    const BCFace neg_faces[3] = {BC_FACE_NEG_X, BC_FACE_NEG_Y, BC_FACE_NEG_Z};
    const BCFace pos_faces[3] = {BC_FACE_POS_X, BC_FACE_POS_Y, BC_FACE_POS_Z};
    const char axis_names[3] = {'X', 'Y', 'Z'};
    DMBoundaryType bx, by, bz;
    PetscBool dm_periodic[3];
    PetscFunctionBeginUser;
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Validating parsed boundary condition configuration...\n");
    ierr = DMDAGetInfo(user->da, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
                       &bx, &by, &bz, NULL); CHKERRQ(ierr);
    dm_periodic[0] = (PetscBool)(bx == DM_BOUNDARY_PERIODIC);
    dm_periodic[1] = (PetscBool)(by == DM_BOUNDARY_PERIODIC);
    dm_periodic[2] = (PetscBool)(bz == DM_BOUNDARY_PERIODIC);
 
    // --- Rule Set 1: Geometric periodic faces must be paired and match the DM topology. ---
    for (PetscInt axis = 0; axis < 3; axis++) {
        const PetscBool neg_periodic =
            user->boundary_faces[neg_faces[axis]].mathematical_type == PERIODIC;
        const PetscBool pos_periodic =
            user->boundary_faces[pos_faces[axis]].mathematical_type == PERIODIC;
 
        PetscCheck(neg_periodic == pos_periodic, PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT,
                   "Configuration Error: Periodic boundaries in the %c direction must be paired; "
                   "%s is %s while %s is %s.",
                   axis_names[axis],
                   BCFaceToString(neg_faces[axis]), neg_periodic ? "PERIODIC" : "not periodic",
                   BCFaceToString(pos_faces[axis]), pos_periodic ? "PERIODIC" : "not periodic");
        PetscCheck(dm_periodic[axis] == neg_periodic, PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT,
                   "Configuration Error: The %c-direction DM periodic flag (%d) does not match "
                   "the paired boundary configuration (%s).",
                   axis_names[axis], (int)dm_periodic[axis], neg_periodic ? "PERIODIC" : "not periodic");
    }
 
    // --- Rule Set 2: Driven Flow Handler Consistency ---
    // This specialized validator will check all rules related to driven flow handlers.
    ierr = Validate_DrivenFlowConfiguration(user); CHKERRQ(ierr);
 
    // --- Rule Set 3: (Future Extension) Overset Interface Consistency ---
    // ierr = Validate_OversetConfiguration(user); CHKERRQ(ierr);
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Boundary configuration is valid.\n");
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

BoundarySystem_Initialize()

PetscErrorCode BoundarySystem_Initialize	(	UserCtx *	user,
		const char *	bcs_filename
	)

Implementation of BoundarySystem_Initialize().

Initializes the entire boundary system.

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

See also

BoundarySystem_Initialize()

Definition at line 891 of file Boundaries.c.

{
    PetscErrorCode ierr;
    PetscFunctionBeginUser;
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Starting creation and initialization of all boundary handlers.\n");
 
    // =========================================================================
    // Step 0: Clear any existing boundary handlers (if re-initializing).
    // This ensures no memory leaks if this function is called multiple times.
    // =========================================================================
    for (int i = 0; i < 6; i++) {
        BoundaryFaceConfig *face_cfg = &user->boundary_faces[i];
        if (face_cfg->handler) {
            LOG_ALLOW(LOCAL, LOG_DEBUG, "Destroying existing handler on Face %s before re-initialization.\n", BCFaceToString((BCFace)i));
            if (face_cfg->handler->Destroy) {
                ierr = face_cfg->handler->Destroy(face_cfg->handler); CHKERRQ(ierr);
            }
            ierr = PetscFree(face_cfg->handler); CHKERRQ(ierr);
            face_cfg->handler = NULL;
        }
    }
    // =========================================================================
 
    // Step 0.1: Initiate flux sums to zero
    user->simCtx->FluxInSum = 0.0;
    user->simCtx->FluxOutSum = 0.0;
    user->simCtx->FarFluxInSum = 0.0;
    user->simCtx->FarFluxOutSum = 0.0;
    // =========================================================================
 
    // Step 1: Parse the configuration file to determine user intent.
    // This function, defined in io.c, populates the configuration enums and parameter
    // lists within the user->boundary_faces array on all MPI ranks.
    ierr = ParseAllBoundaryConditions(user, bcs_filename); CHKERRQ(ierr);
    LOG_ALLOW(GLOBAL, LOG_INFO, "Configuration file '%s' parsed successfully.\n", bcs_filename);
 
    // Step 1.1: Validate the parsed configuration to ensure there are no Boundary Condition conflicts
    ierr = BoundarySystem_Validate(user); CHKERRQ(ierr);
 
    // Step 2: Create and Initialize the handler object for each of the 6 faces.
    for (int i = 0; i < 6; i++) {
        BoundaryFaceConfig *face_cfg = &user->boundary_faces[i];
        
        const char *face_name = BCFaceToString(face_cfg->face_id);
        const char *type_name = BCTypeToString(face_cfg->mathematical_type);
        const char *handler_name = BCHandlerTypeToString(face_cfg->handler_type);
 
        LOG_ALLOW(LOCAL, LOG_DEBUG, "Creating handler for Face %s with Type %s and handler '%s'.\n", face_name, type_name,handler_name);
 
        // Use the private factory to construct the correct handler object based on the parsed type.
        // The factory returns a pointer to the new handler object, which we store in the config struct.
        ierr = BoundaryCondition_Create(face_cfg->handler_type, &face_cfg->handler); CHKERRQ(ierr);
 
        // Step 3: Call the specific Initialize() method for the newly created handler.
        // This allows the handler to perform its own setup, like reading parameters from the
        // face_cfg->params list and setting the initial field values on its face.
        if (face_cfg->handler && face_cfg->handler->Initialize) {
            LOG_ALLOW(LOCAL, LOG_DEBUG, "Calling Initialize() method for handler %s(%s) on Face %s.\n",type_name,handler_name,face_name);
            
            // Prepare the context needed by the Initialize() function.
            BCContext ctx = {
                .user = user,
                .face_id = face_cfg->face_id,
                .global_inflow_sum = &user->simCtx->FluxInSum,  // Global flux sums are not relevant during initialization.
                .global_outflow_sum = &user->simCtx->FluxOutSum,
                .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
                .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
            };
            
            ierr = face_cfg->handler->Initialize(face_cfg->handler, &ctx); CHKERRQ(ierr);
        } else {
            LOG_ALLOW(LOCAL, LOG_DEBUG, "Handler %s(%s) for Face %s has no Initialize() method, skipping.\n", type_name,handler_name,face_name);
        }
    }
    // =========================================================================
    // NO SYNCHRONIZATION NEEDED HERE
    // =========================================================================
    // Initialize() only reads parameters and allocates memory.
    // It does NOT modify field values (Ucat, Ucont, Ubcs).
    // Field values are set by:
    //   1. Initial conditions (before this function)
    //   2. Apply() during timestepping (after this function)
    // The first call to ApplyBoundaryConditions() will handle synchronization.
    // =========================================================================
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "All boundary handlers created and initialized successfully.\n");
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

PropagateBoundaryConfigToCoarserLevels()

PetscErrorCode PropagateBoundaryConfigToCoarserLevels

(

SimCtx *

simCtx

)

Internal helper implementation: PropagateBoundaryConfigToCoarserLevels().

Propagates boundary condition configuration from finest to all coarser multigrid levels.

Local to this translation unit.

Definition at line 988 of file Boundaries.c.

{
    PetscErrorCode ierr;
    UserMG *usermg = &simCtx->usermg;
    
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    
    LOG_ALLOW(GLOBAL, LOG_INFO, "Propagating BC configuration from finest to coarser multigrid levels...\n");
    
    // Loop from second-finest down to coarsest
    for (PetscInt level = usermg->mglevels - 2; level >= 0; level--) {
        for (PetscInt bi = 0; bi < simCtx->block_number; bi++) {
            UserCtx *user_coarse = &usermg->mgctx[level].user[bi];
            UserCtx *user_fine   = &usermg->mgctx[level + 1].user[bi];
            
            LOG_ALLOW_SYNC(LOCAL, LOG_DEBUG, "Rank %d: Copying BC config from level %d to level %d, block %d\n",
                          simCtx->rank, level + 1, level, bi);
            
            // Copy the 6 boundary face configurations
            for (int face_i = 0; face_i < 6; face_i++) {
                user_coarse->boundary_faces[face_i].face_id = user_fine->boundary_faces[face_i].face_id;
                user_coarse->boundary_faces[face_i].mathematical_type = user_fine->boundary_faces[face_i].mathematical_type;
                user_coarse->boundary_faces[face_i].handler_type = user_fine->boundary_faces[face_i].handler_type;
                
                // Copy parameter list (deep copy)
                FreeBC_ParamList(user_coarse->boundary_faces[face_i].params); // Clear any existing
                user_coarse->boundary_faces[face_i].params = NULL;
                
                BC_Param **dst_next = &user_coarse->boundary_faces[face_i].params;
                for (BC_Param *src = user_fine->boundary_faces[face_i].params; src; src = src->next) {
                    BC_Param *new_param;
                    ierr = PetscMalloc1(1, &new_param); CHKERRQ(ierr);
                    ierr = PetscStrallocpy(src->key, &new_param->key); CHKERRQ(ierr);
                    ierr = PetscStrallocpy(src->value, &new_param->value); CHKERRQ(ierr);
                    new_param->next = NULL;
                    *dst_next = new_param;
                    dst_next = &new_param->next;
                }
                
                // IMPORTANT: Do NOT create handler objects for coarser levels
                // Handlers are only needed at finest level for timestepping Apply() calls
                user_coarse->boundary_faces[face_i].handler = NULL;
            }
            
            // Propagate the particle inlet lookup fields to coarse levels as well.
            user_coarse->inletFaceDefined = user_fine->inletFaceDefined;
            user_coarse->identifiedInletBCFace = user_fine->identifiedInletBCFace;
        }
    }
    
    LOG_ALLOW(GLOBAL, LOG_INFO, "BC configuration propagation complete.\n");
    
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

BoundarySystem_ExecuteStep()

PetscErrorCode BoundarySystem_ExecuteStep

(

UserCtx *

user

)

Implementation of BoundarySystem_ExecuteStep().

Executes one full boundary condition update cycle for a time step.

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

See also

BoundarySystem_ExecuteStep()

Definition at line 1059 of file Boundaries.c.

{
    PetscErrorCode ierr;
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    
    LOG_ALLOW(LOCAL, LOG_DEBUG, "Starting.\n");
    
    // =========================================================================
    // PRIORITY 0: INLETS
    // =========================================================================
    
        PetscReal local_inflow_pre = 0.0;
        PetscReal local_inflow_post = 0.0;
        PetscReal global_inflow_pre = 0.0;
        PetscReal global_inflow_post = 0.0;
        PetscInt  num_handlers[3] = {0,0,0};
        
        LOG_ALLOW(LOCAL, LOG_TRACE, " (INLETS): Begin.\n");
        
        // Phase 1: PreStep - Preparation (e.g., calculate profiles, read files)
        for (int i = 0; i < 6; i++) {
            BoundaryCondition *handler = user->boundary_faces[i].handler;
            if (!handler || handler->priority != BC_PRIORITY_INLET) continue;
            if (!handler->PreStep) continue;
            
            num_handlers[0]++;
            BCContext ctx = {
                .user = user,
                .face_id = (BCFace)i,
                .global_inflow_sum = NULL,
                .global_outflow_sum = NULL,
                .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
                .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
            };
            
            LOG_ALLOW(LOCAL, LOG_TRACE, "    PreStep: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
            ierr = handler->PreStep(handler, &ctx, &local_inflow_pre, NULL); CHKERRQ(ierr);
        }
        
        // Optional: Global communication for PreStep (for debugging)
        if (local_inflow_pre != 0.0) {
            ierr = MPI_Allreduce(&local_inflow_pre, &global_inflow_pre, 1, MPIU_REAL,
                                MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
            LOG_ALLOW(GLOBAL, LOG_TRACE, "    PreStep predicted flux: %.6e\n", global_inflow_pre);
        }
        
        // Phase 2: Apply - Set boundary conditions
        for (int i = 0; i < 6; i++) {
            BoundaryCondition *handler = user->boundary_faces[i].handler;
            if (!handler || handler->priority != BC_PRIORITY_INLET) continue;
            if(!handler->Apply) continue; // For example Periodic BCs  
            
            num_handlers[1]++;
 
            BCContext ctx = {
                .user = user,
                .face_id = (BCFace)i,
                .global_inflow_sum = NULL,
                .global_outflow_sum = NULL,
                .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
                .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
            };
            
            LOG_ALLOW(LOCAL, LOG_TRACE, "    Apply: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
            ierr = handler->Apply(handler, &ctx); CHKERRQ(ierr);
        }
        
        // Phase 3: PostStep - Measure actual flux
        for (int i = 0; i < 6; i++) {
            BoundaryCondition *handler = user->boundary_faces[i].handler;
            if (!handler || handler->priority != BC_PRIORITY_INLET) continue;
            if (!handler->PostStep) continue;
            
            num_handlers[2]++;
 
            BCContext ctx = {
                .user = user,
                .face_id = (BCFace)i,
                .global_inflow_sum = NULL,
                .global_outflow_sum = NULL,
                .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
                .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
            };
            
            LOG_ALLOW(LOCAL, LOG_TRACE, "    PostStep: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
            ierr = handler->PostStep(handler, &ctx, &local_inflow_post, NULL); CHKERRQ(ierr);
        }
        
        // Phase 4: Global communication - Sum flux for other priorities to use
        ierr = MPI_Allreduce(&local_inflow_post, &global_inflow_post, 1, MPIU_REAL,
                            MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
        
        // Store for next priority levels
        user->simCtx->FluxInSum = global_inflow_post;
        
        LOG_ALLOW(GLOBAL, LOG_INFO, 
                  "  (INLETS): %d Prestep(s), %d Application(s), %d Poststep(s), FluxInSum = %.6e\n",
                  num_handlers[0],num_handlers[1],num_handlers[2], global_inflow_post);
    
    // =========================================================================
    // PRIORITY 1: FARFIELD
    // =========================================================================
    
        PetscReal local_farfield_in_pre = 0.0;
        PetscReal local_farfield_out_pre = 0.0;
        PetscReal local_farfield_in_post = 0.0;
        PetscReal local_farfield_out_post = 0.0;
        PetscReal global_farfield_in_pre = 0.0;
        PetscReal global_farfield_out_pre = 0.0;
        PetscReal global_farfield_in_post = 0.0;
        PetscReal global_farfield_out_post = 0.0;
        memset(num_handlers,0,sizeof(num_handlers));
        
        LOG_ALLOW(LOCAL, LOG_TRACE, "  (FARFIELD): Begin.\n");
        
        // Phase 1: PreStep - Analyze flow direction, measure initial flux
        for (int i = 0; i < 6; i++) {
            BoundaryCondition *handler = user->boundary_faces[i].handler;
            if (!handler || handler->priority != BC_PRIORITY_FARFIELD) continue;
            if (!handler->PreStep) continue;
            
            num_handlers[0]++;
            BCContext ctx = {
                .user = user,
                .face_id = (BCFace)i,
                .global_inflow_sum = &user->simCtx->FluxInSum,  // Available from Priority 0
                .global_outflow_sum = NULL,
                .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
                .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
            };
            
            LOG_ALLOW(LOCAL, LOG_TRACE, "    PreStep: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
            ierr = handler->PreStep(handler, &ctx, &local_farfield_in_pre, &local_farfield_out_pre);
            CHKERRQ(ierr);
        }
        
        // Phase 2: Global communication (optional, for debugging)
        if (local_farfield_in_pre != 0.0 || local_farfield_out_pre != 0.0) {
            ierr = MPI_Allreduce(&local_farfield_in_pre, &global_farfield_in_pre, 1, MPIU_REAL,
                                MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
            ierr = MPI_Allreduce(&local_farfield_out_pre, &global_farfield_out_pre, 1, MPIU_REAL,
                                MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
            
            LOG_ALLOW(GLOBAL, LOG_DEBUG, 
                      "    Farfield pre-analysis: In=%.6e, Out=%.6e\n",
                      global_farfield_in_pre, global_farfield_out_pre);
        }
        
        // Phase 3: Apply - Set farfield boundary conditions
        for (int i = 0; i < 6; i++) {
            BoundaryCondition *handler = user->boundary_faces[i].handler;
            if (!handler || handler->priority != BC_PRIORITY_FARFIELD) continue;
            if(!handler->Apply) continue; // For example Periodic BCs            
            
            num_handlers[1]++;
 
            BCContext ctx = {
                .user = user,
                .face_id = (BCFace)i,
                .global_inflow_sum = &user->simCtx->FluxInSum,
                .global_outflow_sum = NULL,
                .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
                .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
            };
            
            LOG_ALLOW(LOCAL, LOG_TRACE, "    Apply: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
            ierr = handler->Apply(handler, &ctx); CHKERRQ(ierr);
        }
        
        // Phase 4: PostStep - Measure actual farfield fluxes
        for (int i = 0; i < 6; i++) {
            BoundaryCondition *handler = user->boundary_faces[i].handler;
            if (!handler || handler->priority != BC_PRIORITY_FARFIELD) continue;
            if (!handler->PostStep) continue;
            
            num_handlers[2]++;
 
            BCContext ctx = {
                .user = user,
                .face_id = (BCFace)i,
                .global_inflow_sum = &user->simCtx->FluxInSum,
                .global_outflow_sum = NULL,
                .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
                .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
            };
            
            LOG_ALLOW(LOCAL, LOG_TRACE, "    PostStep: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
            ierr = handler->PostStep(handler, &ctx, &local_farfield_in_post, &local_farfield_out_post);
            CHKERRQ(ierr);
        }
        
        // Phase 5: Global communication - Store for outlet priority
        if (num_handlers > 0) {
            ierr = MPI_Allreduce(&local_farfield_in_post, &global_farfield_in_post, 1, MPIU_REAL,
                                MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
            ierr = MPI_Allreduce(&local_farfield_out_post, &global_farfield_out_post, 1, MPIU_REAL,
                                MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
            
            // Store for outlet handlers to use
            user->simCtx->FarFluxInSum = global_farfield_in_post;
            user->simCtx->FarFluxOutSum = global_farfield_out_post;
            
            LOG_ALLOW(GLOBAL, LOG_INFO, 
                      "  (FARFIELD): %d Prestep(s), %d Application(s), %d Poststep(s) , InFlux=%.6e, OutFlux=%.6e\n",
                      num_handlers[0],num_handlers[1],num_handlers[2], global_farfield_in_post, global_farfield_out_post);
        } else {
            // No farfield handlers - zero out the fluxes
            user->simCtx->FarFluxInSum = 0.0;
            user->simCtx->FarFluxOutSum = 0.0;
        }
    
    
    // =========================================================================
    // PRIORITY 2: WALLS
    // =========================================================================
    
        memset(num_handlers,0,sizeof(num_handlers));
        
        LOG_ALLOW(LOCAL, LOG_TRACE, "  (WALLS): Begin.\n");
        
        // Phase 1: PreStep - Preparation (usually no-op for walls)
        for (int i = 0; i < 6; i++) {
            BoundaryCondition *handler = user->boundary_faces[i].handler;
            if (!handler || handler->priority != BC_PRIORITY_WALL) continue;
            if (!handler->PreStep) continue;
            
            num_handlers[0]++;
            BCContext ctx = {
                .user = user,
                .face_id = (BCFace)i,
                .global_inflow_sum = &user->simCtx->FluxInSum,
                .global_outflow_sum = NULL,
                .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
                .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
            };
            
            LOG_ALLOW(LOCAL, LOG_TRACE, "    PreStep: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
            ierr = handler->PreStep(handler, &ctx, NULL, NULL); CHKERRQ(ierr);
        }
        
        // No global communication needed for walls
        
        // Phase 2: Apply - Set boundary conditions
        for (int i = 0; i < 6; i++) {
            BoundaryCondition *handler = user->boundary_faces[i].handler;
            if (!handler || handler->priority != BC_PRIORITY_WALL) continue;
            if(!handler->Apply) continue; // For example Periodic BCs  
            
            num_handlers[1]++;
 
            BCContext ctx = {
                .user = user,
                .face_id = (BCFace)i,
                .global_inflow_sum = &user->simCtx->FluxInSum,
                .global_outflow_sum = NULL,
                .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
                .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
            };
            
            LOG_ALLOW(LOCAL, LOG_TRACE, "    Apply: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
            ierr = handler->Apply(handler, &ctx); CHKERRQ(ierr);
        }
        
        // Phase 3: PostStep - Post-application processing (usually no-op for walls)
        for (int i = 0; i < 6; i++) {
            BoundaryCondition *handler = user->boundary_faces[i].handler;
            if (!handler || handler->priority != BC_PRIORITY_WALL) continue;
            if (!handler->PostStep) continue;
            
            num_handlers[2]++;
 
            BCContext ctx = {
                .user = user,
                .face_id = (BCFace)i,
                .global_inflow_sum = &user->simCtx->FluxInSum,
                .global_outflow_sum = NULL,
                .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
                .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
            };
            
            LOG_ALLOW(LOCAL, LOG_TRACE, "    PostStep: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
            ierr = handler->PostStep(handler, &ctx, NULL, NULL); CHKERRQ(ierr);
        }
        
        // No global communication needed for walls
        
        LOG_ALLOW(GLOBAL, LOG_INFO, "  (WALLS): %d Prestep(s), %d Application(s), %d Poststep(s) applied.\n",
                  num_handlers[0],num_handlers[1],num_handlers[2]);
    
    
    // =========================================================================
    // PRIORITY 3: OUTLETS
    // =========================================================================
    
        PetscReal local_outflow_pre = 0.0;
        PetscReal local_outflow_post = 0.0;
        PetscReal global_outflow_pre = 0.0;
        PetscReal global_outflow_post = 0.0;
        memset(num_handlers,0,sizeof(num_handlers));
        
        LOG_ALLOW(LOCAL, LOG_TRACE, "  (OUTLETS): Begin.\n");
        
        // Phase 1: PreStep - Measure uncorrected outflow (from ucat)
        for (int i = 0; i < 6; i++) {
            BoundaryCondition *handler = user->boundary_faces[i].handler;
            if (!handler || handler->priority != BC_PRIORITY_OUTLET) continue;
            if (!handler->PreStep) continue;
            
            num_handlers[0]++;
            BCContext ctx = {
                .user = user,
                .face_id = (BCFace)i,
                .global_inflow_sum = &user->simCtx->FluxInSum,      // From Priority 0
                .global_outflow_sum = NULL,
                .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
                .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
            };
            
            LOG_ALLOW(LOCAL, LOG_TRACE, "    PreStep: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
            ierr = handler->PreStep(handler, &ctx, NULL, &local_outflow_pre); CHKERRQ(ierr);
        }
        
        // Phase 2: Global communication - Get uncorrected outflow sum
        ierr = MPI_Allreduce(&local_outflow_pre, &global_outflow_pre, 1, MPIU_REAL,
                            MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
        
        // Calculate total inflow (inlet + farfield inflow)
        PetscReal total_inflow = user->simCtx->FluxInSum + user->simCtx->FarFluxInSum;
        
        LOG_ALLOW(GLOBAL, LOG_DEBUG, 
                  "    Uncorrected outflow: %.6e, Total inflow: %.6e (Inlet: %.6e + Farfield: %.6e)\n",
                  global_outflow_pre, total_inflow, user->simCtx->FluxInSum, 
                  user->simCtx->FarFluxInSum);
        
        // Phase 3: Apply - Set corrected boundary conditions
        for (int i = 0; i < 6; i++) {
            BoundaryCondition *handler = user->boundary_faces[i].handler;
            if (!handler || handler->priority != BC_PRIORITY_OUTLET) continue;
            if(!handler->Apply) continue; // For example Periodic BCs  
            
            num_handlers[1]++;
 
            BCContext ctx = {
                .user = user,
                .face_id = (BCFace)i,
                .global_inflow_sum = &user->simCtx->FluxInSum,      // From Priority 0
                .global_outflow_sum = &global_outflow_pre, // From PreStep above
                .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
                .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum 
            };
            
            LOG_ALLOW(LOCAL, LOG_TRACE, "    Apply: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
            ierr = handler->Apply(handler, &ctx); CHKERRQ(ierr);
        }
        
        // Phase 4: PostStep - Measure corrected outflow (verification)
        for (int i = 0; i < 6; i++) {
            BoundaryCondition *handler = user->boundary_faces[i].handler;
            if (!handler || handler->priority != BC_PRIORITY_OUTLET) continue;
            if (!handler->PostStep) continue;
            
            num_handlers[2]++;
 
            BCContext ctx = {
                .user = user,
                .face_id = (BCFace)i,
                .global_inflow_sum = &user->simCtx->FluxInSum,
                .global_outflow_sum = &global_outflow_pre,
                .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
                .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
            };
            
            LOG_ALLOW(LOCAL, LOG_TRACE, "    PostStep: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
            ierr = handler->PostStep(handler, &ctx, NULL, &local_outflow_post); CHKERRQ(ierr);
        }
        
        // Phase 5: Global communication - Verify conservation
        ierr = MPI_Allreduce(&local_outflow_post, &global_outflow_post, 1, MPIU_REAL,
                            MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
        
        // Store for global reporting.
        user->simCtx->FluxOutSum = global_outflow_post;
        
        // Conservation check (compare total outflow vs total inflow)
        PetscReal total_outflow = global_outflow_post + user->simCtx->FarFluxOutSum;
        PetscReal flux_error = PetscAbsReal(total_outflow - total_inflow);
        PetscReal relative_error = (total_inflow > 1e-16) ? 
                                   flux_error / total_inflow : flux_error;
        
        LOG_ALLOW(GLOBAL, LOG_INFO, 
                  "  (OUTLETS): %d Prestep(s), %d Application(s), %d Poststep(s), FluxOutSum = %.6e\n",
                  num_handlers[0],num_handlers[1],num_handlers[2], global_outflow_post);
        LOG_ALLOW(GLOBAL, LOG_INFO, 
                  "    Conservation: Total In=%.6e, Total Out=%.6e, Error=%.3e (%.2e)%%)\n",
                  total_inflow, total_outflow, flux_error, relative_error * 100.0);
        
        if (relative_error > 1e-6) {
            LOG_ALLOW(GLOBAL, LOG_WARNING, 
                     "    WARNING: Large mass conservation error (%.2e%%)!\n",
                     relative_error * 100.0);
        }
    
    
    LOG_ALLOW(LOCAL, LOG_VERBOSE, "Complete.\n");
    
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

BoundarySystem_RefreshUbcs()

PetscErrorCode BoundarySystem_RefreshUbcs

(

UserCtx *

user

)

Internal helper implementation: BoundarySystem_RefreshUbcs().

(Private) A lightweight execution engine that calls the UpdateUbcs() method on all relevant handlers.

Local to this translation unit.

Definition at line 1481 of file Boundaries.c.

{
    PetscErrorCode ierr;
    PetscFunctionBeginUser;
 
    LOG_ALLOW(GLOBAL, LOG_TRACE, "Refreshing `ubcs` targets for flow-dependent boundaries...\n");
 
    // Loop through all 6 faces of the domain
    for (int i = 0; i < 6; i++) {
        BoundaryCondition *handler = user->boundary_faces[i].handler;
        
        // THE FILTER:
        // This is the core logic. We only act if a handler exists for the face
        // AND that handler has explicitly implemented the `UpdateUbcs` method.
        if (handler && handler->UpdateUbcs) {
            
            const char *face_name = BCFaceToString((BCFace)i);
            LOG_ALLOW(LOCAL, LOG_TRACE, "  Calling UpdateUbcs() for handler on Face %s.\n", face_name);
 
            // Prepare the context. For this refresh step, we don't need to pass flux sums.
            BCContext ctx = {
                .user = user,
                .face_id = (BCFace)i,
                .global_inflow_sum = NULL,
                .global_outflow_sum = NULL,
                .global_farfield_inflow_sum = NULL,
                .global_farfield_outflow_sum = NULL
            };
            
            // Call the handler's specific UpdateUbcs function pointer.
            ierr = handler->UpdateUbcs(handler, &ctx); CHKERRQ(ierr);
        }
    }
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

BoundarySystem_Destroy()

PetscErrorCode BoundarySystem_Destroy

(

UserCtx *

user

)

Implementation of BoundarySystem_Destroy().

Cleans up and destroys all boundary system resources.

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

See also

BoundarySystem_Destroy()

Definition at line 1531 of file Boundaries.c.

{
    PetscErrorCode ierr;
    PetscFunctionBeginUser;
 
    
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Starting destruction of all boundary handlers. \n");
 
    for (int i = 0; i < 6; i++) {
        BoundaryFaceConfig *face_cfg = &user->boundary_faces[i];
        const char *face_name = BCFaceToString(face_cfg->face_id);
 
        // --- Step 1: Free the parameter linked list associated with this face ---
        if (face_cfg->params) {
            LOG_ALLOW(LOCAL, LOG_DEBUG, "  Freeing parameter list for Face %d (%s). \n", i, face_name);
            FreeBC_ParamList(face_cfg->params);
            face_cfg->params = NULL; // Good practice to nullify dangling pointers
        }
 
        // --- Step 2: Destroy the handler object itself ---
        if (face_cfg->handler) {
            const char *handler_name = BCHandlerTypeToString(face_cfg->handler->type);
            LOG_ALLOW(LOCAL, LOG_DEBUG, "  Destroying handler '%s' on Face %d (%s).\n", handler_name, i, face_name);
            
            // Call the handler's specific cleanup function first, if it exists.
            // This will free any memory stored in the handler's private `data` pointer.
            if (face_cfg->handler->Destroy) {
                ierr = face_cfg->handler->Destroy(face_cfg->handler); CHKERRQ(ierr);
            }
 
            // Finally, free the generic BoundaryCondition object itself.
            ierr = PetscFree(face_cfg->handler); CHKERRQ(ierr);
            face_cfg->handler = NULL;
        }
    }
    
    LOG_ALLOW(GLOBAL, LOG_INFO, "Destruction complete.\n");
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TransferPeriodicFieldByDirection()

PetscErrorCode TransferPeriodicFieldByDirection	(	UserCtx *	user,
		const char *	field_name,
		char	direction
	)

Internal helper implementation: TransferPeriodicFieldByDirection().

(Private Worker) Copies periodic data for a SINGLE field in a SINGLE direction.

Local to this translation unit.

Definition at line 1578 of file Boundaries.c.

{
    PetscErrorCode ierr;
    DMDALocalInfo  info = user->info;
    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;
 
    // --- Dispatcher to get DM, Vecs, and DoF for the specified field ---
    DM        dm;
    Vec       global_vec;
    Vec       local_vec;
    PetscInt  dof;
    // (This dispatcher is identical to your TransferPeriodicField function)
    if (strcmp(field_name, "Ucat") == 0) {
        dm = user->fda; global_vec = user->Ucat; local_vec = user->lUcat; dof = 3;
    } else if (strcmp(field_name, "P") == 0) {
        dm = user->da; global_vec = user->P; local_vec = user->lP; dof = 1;
    } else if (strcmp(field_name, "Phi") == 0) {
        dm = user->da; global_vec = user->Phi; local_vec = user->lPhi; dof = 1;
    } else if (strcmp(field_name, "Nvert") == 0) {
        dm = user->da; global_vec = user->Nvert; local_vec = user->lNvert; dof = 1;
    } else if (strcmp(field_name, "Nu_t") == 0 || strcmp(field_name, "Eddy Viscosity") == 0) {
        dm = user->da; global_vec = user->Nu_t; local_vec = user->lNu_t; dof = 1;
    } else if (strcmp(field_name, "CS") == 0 || strcmp(field_name, "Cs") == 0) {
        dm = user->da; global_vec = user->CS; local_vec = user->lCs; dof = 1;
    } else if (strcmp(field_name, "Diffusivity") == 0) {
        dm = user->da; global_vec = user->Diffusivity; local_vec = user->lDiffusivity; dof = 1;
    } else if (strcmp(field_name, "Aj") == 0) {
        dm = user->da; global_vec = user->Aj; local_vec = user->lAj; dof = 1;
    } else {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown field name '%s'", field_name);
    }
 
    PetscFunctionBeginUser;
 
    // --- Execute the copy logic based on DoF and Direction ---
    if (dof == 1) { // --- Handle SCALAR fields (PetscReal) ---
        PetscReal ***g_array, ***l_array;
        ierr = DMDAVecGetArray(dm, global_vec, &g_array); CHKERRQ(ierr);
        ierr = DMDAVecGetArrayRead(dm, local_vec, (void*)&l_array); CHKERRQ(ierr); // Use Read for safety
 
        switch (direction) {
            case 'i':
                if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && xs == 0) for (PetscInt k=zs; k<ze; k++) for (PetscInt j=ys; j<ye; j++) g_array[k][j][xs] = l_array[k][j][xs-2];
                if (user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC && xe == mx) for (PetscInt k=zs; k<ze; k++) for (PetscInt j=ys; j<ye; j++) g_array[k][j][xe-1] = l_array[k][j][xe+1];
                break;
            case 'j':
                if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && ys == 0) for (PetscInt k=zs; k<ze; k++) for (PetscInt i=xs; i<xe; i++) g_array[k][ys][i] = l_array[k][ys-2][i];
                if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC && ye == my) for (PetscInt k=zs; k<ze; k++) for (PetscInt i=xs; i<xe; i++) g_array[k][ye-1][i] = l_array[k][ye+1][i];
                break;
            case 'k':
                if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && zs == 0) for (PetscInt j=ys; j<ye; j++) for (PetscInt i=xs; i<xe; i++) g_array[zs][j][i] = l_array[zs-2][j][i];
                if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC && ze == mz) for (PetscInt j=ys; j<ye; j++) for (PetscInt i=xs; i<xe; i++) g_array[ze-1][j][i] = l_array[ze+1][j][i];
                break;
            default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid direction '%c'", direction);
        }
        ierr = DMDAVecRestoreArray(dm, global_vec, &g_array); CHKERRQ(ierr);
        ierr = DMDAVecRestoreArrayRead(dm, local_vec, (void*)&l_array); CHKERRQ(ierr);
 
    } else if (dof == 3) { // --- Handle VECTOR fields (Cmpnts) ---
        Cmpnts ***g_array, ***l_array;
        ierr = DMDAVecGetArray(dm, global_vec, &g_array); CHKERRQ(ierr);
        ierr = DMDAVecGetArrayRead(dm, local_vec, (void*)&l_array); CHKERRQ(ierr);
 
        switch (direction) {
            case 'i':
                if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && xs == 0) for (PetscInt k=zs; k<ze; k++) for (PetscInt j=ys; j<ye; j++) g_array[k][j][xs] = l_array[k][j][xs-2];
                if (user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC && xe == mx) for (PetscInt k=zs; k<ze; k++) for (PetscInt j=ys; j<ye; j++) g_array[k][j][xe-1] = l_array[k][j][xe+1];
                break;
            case 'j':
                if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && ys == 0) for (PetscInt k=zs; k<ze; k++) for (PetscInt i=xs; i<xe; i++) g_array[k][ys][i] = l_array[k][ys-2][i];
                if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC && ye == my) for (PetscInt k=zs; k<ze; k++) for (PetscInt i=xs; i<xe; i++) g_array[k][ye-1][i] = l_array[k][ye+1][i];
                break;
            case 'k':
                if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && zs == 0) for (PetscInt j=ys; j<ye; j++) for (PetscInt i=xs; i<xe; i++) g_array[zs][j][i] = l_array[zs-2][j][i];
                if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC && ze == mz) for (PetscInt j=ys; j<ye; j++) for (PetscInt i=xs; i<xe; i++) g_array[ze-1][j][i] = l_array[ze+1][j][i];
                break;
            default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid direction '%c'", direction);
        }
        ierr = DMDAVecRestoreArray(dm, global_vec, &g_array); CHKERRQ(ierr);
        ierr = DMDAVecRestoreArrayRead(dm, local_vec, (void*)&l_array); CHKERRQ(ierr);
    }
 
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

SynchronizePeriodicCellFields()

PetscErrorCode SynchronizePeriodicCellFields	(	UserCtx *	user,
		PetscInt	num_fields,
		const char *	field_names[]
	)

Implementation of SynchronizePeriodicCellFields().

Synchronizes periodic endpoint cells for a list of cell-centered fields.

Full API contract is documented with the header declaration in include/Boundaries.h.

Definition at line 1673 of file Boundaries.c.

{
    PetscErrorCode ierr;
    PetscBool      periodic_i;
    PetscBool      periodic_j;
    PetscBool      periodic_k;
 
    PetscFunctionBeginUser;
 
    if (num_fields == 0) PetscFunctionReturn(0);
 
    periodic_i =
        user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC ||
        user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC;
    periodic_j =
        user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC ||
        user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC;
    periodic_k =
        user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC ||
        user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC;
 
    if (!periodic_i && !periodic_j && !periodic_k) PetscFunctionReturn(0);
 
    for (PetscInt field = 0; field < num_fields; field++) {
        ierr = UpdateLocalGhosts(user, field_names[field]); CHKERRQ(ierr);
    }
 
    if (periodic_i) {
        for (PetscInt field = 0; field < num_fields; field++) {
            ierr = TransferPeriodicFieldByDirection(user, field_names[field], 'i'); CHKERRQ(ierr);
        }
        for (PetscInt field = 0; field < num_fields; field++) {
            ierr = UpdateLocalGhosts(user, field_names[field]); CHKERRQ(ierr);
        }
    }
 
    if (periodic_j) {
        for (PetscInt field = 0; field < num_fields; field++) {
            ierr = TransferPeriodicFieldByDirection(user, field_names[field], 'j'); CHKERRQ(ierr);
        }
        for (PetscInt field = 0; field < num_fields; field++) {
            ierr = UpdateLocalGhosts(user, field_names[field]); CHKERRQ(ierr);
        }
    }
 
    if (periodic_k) {
        for (PetscInt field = 0; field < num_fields; field++) {
            ierr = TransferPeriodicFieldByDirection(user, field_names[field], 'k'); CHKERRQ(ierr);
        }
        for (PetscInt field = 0; field < num_fields; field++) {
            ierr = UpdateLocalGhosts(user, field_names[field]); CHKERRQ(ierr);
        }
    }
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

GetPersistentFaceField()

static PetscErrorCode GetPersistentFaceField ( UserCtx *  user, const char *  field_name, char  face_direction, DM *  dm, Vec *  global_vec, Vec *  local_vec, PetscInt *  dof  )

static

Resolves one registered persistent single-face-family field.

Definition at line 1735 of file Boundaries.c.

{
    PetscFunctionBeginUser;
    PetscCheck(face_direction == 'i' || face_direction == 'j' || face_direction == 'k',
               PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
               "Invalid face direction '%c'; expected 'i', 'j', or 'k'.", face_direction);
 
    *dm = NULL;
    *global_vec = NULL;
    *local_vec = NULL;
    *dof = 0;
 
    if (face_direction == 'i') {
        if      (strcmp(field_name, "Centx") == 0) { *dm = user->fda; *global_vec = user->Centx; *local_vec = user->lCentx; *dof = 3; }
        else if (strcmp(field_name, "Csi")  == 0) { *dm = user->fda; *global_vec = user->Csi;  *local_vec = user->lCsi;  *dof = 3; }
        else if (strcmp(field_name, "ICsi") == 0) { *dm = user->fda; *global_vec = user->ICsi; *local_vec = user->lICsi; *dof = 3; }
        else if (strcmp(field_name, "IEta") == 0) { *dm = user->fda; *global_vec = user->IEta; *local_vec = user->lIEta; *dof = 3; }
        else if (strcmp(field_name, "IZet") == 0) { *dm = user->fda; *global_vec = user->IZet; *local_vec = user->lIZet; *dof = 3; }
        else if (strcmp(field_name, "IAj")  == 0) { *dm = user->da;  *global_vec = user->IAj;  *local_vec = user->lIAj;  *dof = 1; }
    } else if (face_direction == 'j') {
        if      (strcmp(field_name, "Centy") == 0) { *dm = user->fda; *global_vec = user->Centy; *local_vec = user->lCenty; *dof = 3; }
        else if (strcmp(field_name, "Eta")  == 0) { *dm = user->fda; *global_vec = user->Eta;  *local_vec = user->lEta;  *dof = 3; }
        else if (strcmp(field_name, "JCsi") == 0) { *dm = user->fda; *global_vec = user->JCsi; *local_vec = user->lJCsi; *dof = 3; }
        else if (strcmp(field_name, "JEta") == 0) { *dm = user->fda; *global_vec = user->JEta; *local_vec = user->lJEta; *dof = 3; }
        else if (strcmp(field_name, "JZet") == 0) { *dm = user->fda; *global_vec = user->JZet; *local_vec = user->lJZet; *dof = 3; }
        else if (strcmp(field_name, "JAj")  == 0) { *dm = user->da;  *global_vec = user->JAj;  *local_vec = user->lJAj;  *dof = 1; }
    } else {
        if      (strcmp(field_name, "Centz") == 0) { *dm = user->fda; *global_vec = user->Centz; *local_vec = user->lCentz; *dof = 3; }
        else if (strcmp(field_name, "Zet")  == 0) { *dm = user->fda; *global_vec = user->Zet;  *local_vec = user->lZet;  *dof = 3; }
        else if (strcmp(field_name, "KCsi") == 0) { *dm = user->fda; *global_vec = user->KCsi; *local_vec = user->lKCsi; *dof = 3; }
        else if (strcmp(field_name, "KEta") == 0) { *dm = user->fda; *global_vec = user->KEta; *local_vec = user->lKEta; *dof = 3; }
        else if (strcmp(field_name, "KZet") == 0) { *dm = user->fda; *global_vec = user->KZet; *local_vec = user->lKZet; *dof = 3; }
        else if (strcmp(field_name, "KAj")  == 0) { *dm = user->da;  *global_vec = user->KAj;  *local_vec = user->lKAj;  *dof = 1; }
    }
 
    PetscCheck(*dm && *global_vec && *local_vec, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE,
               "Field '%s' is not a registered %c-face persistent field.",
               field_name, face_direction);
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

IsFaceCenterCoordinateField()

static PetscBool IsFaceCenterCoordinateField ( const char *  field_name )

static

Returns whether a registered face field stores physical coordinates.

Definition at line 1782 of file Boundaries.c.

{
    return (PetscBool)(strcmp(field_name, "Centx") == 0 ||
                       strcmp(field_name, "Centy") == 0 ||
                       strcmp(field_name, "Centz") == 0);
}

Here is the caller graph for this function:

TranslatePeriodicFaceCenterGhosts()

static PetscErrorCode TranslatePeriodicFaceCenterGhosts ( UserCtx *  user, Vec  local_vec  )

static

Applies geometric translations to wrapped face-center ghost coordinates.

Definition at line 1792 of file Boundaries.c.

{
    DMDALocalInfo info;
    Cmpnts ***array;
    const BCFace negative_faces[3] = {BC_FACE_NEG_X, BC_FACE_NEG_Y, BC_FACE_NEG_Z};
    const BCFace positive_faces[3] = {BC_FACE_POS_X, BC_FACE_POS_Y, BC_FACE_POS_Z};
 
    PetscFunctionBeginUser;
    PetscCall(DMDAGetLocalInfo(user->fda, &info));
    PetscCall(DMDAVecGetArray(user->fda, local_vec, &array));
 
    for (PetscInt axis = 0; axis < 3; axis++) {
        const PetscBool active =
            user->boundary_faces[negative_faces[axis]].mathematical_type == PERIODIC ||
            user->boundary_faces[positive_faces[axis]].mathematical_type == PERIODIC;
        Cmpnts translation;
 
        if (!active) continue;
        PetscCheck(user->periodic_translation_valid[axis], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE,
                   "Periodic face-center synchronization requires validated %c-direction geometry.",
                   "XYZ"[axis]);
        translation = user->periodic_translation[axis];
 
        for (PetscInt k = info.gzs; k < info.gzs + info.gzm; k++) {
            for (PetscInt j = info.gys; j < info.gys + info.gym; j++) {
                for (PetscInt i = info.gxs; i < info.gxs + info.gxm; i++) {
                    PetscReal scale = 0.0;
                    const PetscInt index = axis == 0 ? i : (axis == 1 ? j : k);
                    const PetscInt size = axis == 0 ? info.mx : (axis == 1 ? info.my : info.mz);
                    if (index < 0) scale = -1.0;
                    else if (index >= size) scale = 1.0;
                    if (scale == 0.0) continue;
                    array[k][j][i].x += scale * translation.x;
                    array[k][j][i].y += scale * translation.y;
                    array[k][j][i].z += scale * translation.z;
                }
            }
        }
    }
 
    PetscCall(DMDAVecRestoreArray(user->fda, local_vec, &array));
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

TransferPeriodicFaceFieldByDirection()

PetscErrorCode TransferPeriodicFaceFieldByDirection	(	UserCtx *	user,
		const char *	field_name,
		char	face_direction,
		char	periodic_direction
	)

Implementation of TransferPeriodicFaceFieldByDirection().

Transfers one persistent single-face-family field in one periodic direction.

Definition at line 1841 of file Boundaries.c.

{
    PetscErrorCode ierr;
    DMDALocalInfo info = user->info;
    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;
    DM dm;
    Vec global_vec, local_vec;
    PetscInt dof;
 
    PetscFunctionBeginUser;
    PetscCheck(periodic_direction == 'i' || periodic_direction == 'j' || periodic_direction == 'k',
               PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
               "Invalid periodic direction '%c'; expected 'i', 'j', or 'k'.", periodic_direction);
    PetscCall(GetPersistentFaceField(user, field_name, face_direction, &dm, &global_vec, &local_vec, &dof));
 
    if (dof == 1) {
        PetscReal ***global_array, ***local_array;
        ierr = DMDAVecGetArray(dm, global_vec, &global_array); CHKERRQ(ierr);
        ierr = DMDAVecGetArrayRead(dm, local_vec, &local_array); CHKERRQ(ierr);
 
        if (periodic_direction == 'i') {
            if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && xs == 0)
                for (PetscInt k=zs; k<ze; k++) for (PetscInt j=ys; j<ye; j++) global_array[k][j][0] = local_array[k][j][-2];
            if (user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC && xe == mx)
                for (PetscInt k=zs; k<ze; k++) for (PetscInt j=ys; j<ye; j++) global_array[k][j][mx-1] = local_array[k][j][mx+1];
        } else if (periodic_direction == 'j') {
            if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && ys == 0)
                for (PetscInt k=zs; k<ze; k++) for (PetscInt i=xs; i<xe; i++) global_array[k][0][i] = local_array[k][-2][i];
            if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC && ye == my)
                for (PetscInt k=zs; k<ze; k++) for (PetscInt i=xs; i<xe; i++) global_array[k][my-1][i] = local_array[k][my+1][i];
        } else {
            if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && zs == 0)
                for (PetscInt j=ys; j<ye; j++) for (PetscInt i=xs; i<xe; i++) global_array[0][j][i] = local_array[-2][j][i];
            if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC && ze == mz)
                for (PetscInt j=ys; j<ye; j++) for (PetscInt i=xs; i<xe; i++) global_array[mz-1][j][i] = local_array[mz+1][j][i];
        }
        ierr = DMDAVecRestoreArrayRead(dm, local_vec, &local_array); CHKERRQ(ierr);
        ierr = DMDAVecRestoreArray(dm, global_vec, &global_array); CHKERRQ(ierr);
    } else {
        Cmpnts ***global_array, ***local_array;
        ierr = DMDAVecGetArray(dm, global_vec, &global_array); CHKERRQ(ierr);
        ierr = DMDAVecGetArrayRead(dm, local_vec, &local_array); CHKERRQ(ierr);
 
        if (periodic_direction == 'i') {
            if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && xs == 0)
                for (PetscInt k=zs; k<ze; k++) for (PetscInt j=ys; j<ye; j++) global_array[k][j][0] = local_array[k][j][-2];
            if (user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC && xe == mx)
                for (PetscInt k=zs; k<ze; k++) for (PetscInt j=ys; j<ye; j++) global_array[k][j][mx-1] = local_array[k][j][mx+1];
        } else if (periodic_direction == 'j') {
            if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && ys == 0)
                for (PetscInt k=zs; k<ze; k++) for (PetscInt i=xs; i<xe; i++) global_array[k][0][i] = local_array[k][-2][i];
            if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC && ye == my)
                for (PetscInt k=zs; k<ze; k++) for (PetscInt i=xs; i<xe; i++) global_array[k][my-1][i] = local_array[k][my+1][i];
        } else {
            if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && zs == 0)
                for (PetscInt j=ys; j<ye; j++) for (PetscInt i=xs; i<xe; i++) global_array[0][j][i] = local_array[-2][j][i];
            if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC && ze == mz)
                for (PetscInt j=ys; j<ye; j++) for (PetscInt i=xs; i<xe; i++) global_array[mz-1][j][i] = local_array[mz+1][j][i];
        }
        ierr = DMDAVecRestoreArrayRead(dm, local_vec, &local_array); CHKERRQ(ierr);
        ierr = DMDAVecRestoreArray(dm, global_vec, &global_array); CHKERRQ(ierr);
    }
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

SynchronizePeriodicFaceFields()

PetscErrorCode SynchronizePeriodicFaceFields	(	UserCtx *	user,
		char	face_direction,
		PetscInt	num_fields,
		const char *	field_names[]
	)

Implementation of SynchronizePeriodicFaceFields().

Synchronizes persistent fields belonging to one face family.

Definition at line 1916 of file Boundaries.c.

{
    PetscErrorCode ierr;
    const char periodic_directions[3] = {'i', 'j', 'k'};
    const BCFace negative_faces[3] = {BC_FACE_NEG_X, BC_FACE_NEG_Y, BC_FACE_NEG_Z};
    const BCFace positive_faces[3] = {BC_FACE_POS_X, BC_FACE_POS_Y, BC_FACE_POS_Z};
 
    PetscFunctionBeginUser;
    PetscCheck(num_fields >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE,
               "Number of face fields cannot be negative.");
    if (num_fields == 0) PetscFunctionReturn(0);
    PetscCheck(field_names != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
               "Face field-name array cannot be NULL.");
 
    for (PetscInt field = 0; field < num_fields; field++) {
        DM dm;
        Vec global_vec, local_vec;
        PetscInt dof;
        PetscCall(GetPersistentFaceField(user, field_names[field], face_direction,
                                         &dm, &global_vec, &local_vec, &dof));
        ierr = UpdateLocalGhosts(user, field_names[field]); CHKERRQ(ierr);
        if (IsFaceCenterCoordinateField(field_names[field])) {
            ierr = TranslatePeriodicFaceCenterGhosts(user, local_vec); CHKERRQ(ierr);
        }
    }
 
    for (PetscInt direction = 0; direction < 3; direction++) {
        const PetscBool active =
            user->boundary_faces[negative_faces[direction]].mathematical_type == PERIODIC ||
            user->boundary_faces[positive_faces[direction]].mathematical_type == PERIODIC;
        if (!active) continue;
 
        for (PetscInt field = 0; field < num_fields; field++) {
            ierr = TransferPeriodicFaceFieldByDirection(user, field_names[field], face_direction,
                                                        periodic_directions[direction]); CHKERRQ(ierr);
        }
        for (PetscInt field = 0; field < num_fields; field++) {
            DM dm;
            Vec global_vec, local_vec;
            PetscInt dof;
            ierr = UpdateLocalGhosts(user, field_names[field]); CHKERRQ(ierr);
            if (IsFaceCenterCoordinateField(field_names[field])) {
                PetscCall(GetPersistentFaceField(user, field_names[field], face_direction,
                                                 &dm, &global_vec, &local_vec, &dof));
                ierr = TranslatePeriodicFaceCenterGhosts(user, local_vec); CHKERRQ(ierr);
            }
        }
    }
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

GetPersistentStaggeredField()

static PetscErrorCode GetPersistentStaggeredField ( UserCtx *  user, const char *  field_name, DM *  dm, Vec *  global_vec, Vec *  local_vec  )

static

Resolves one registered persistent component-staggered field.

Definition at line 1974 of file Boundaries.c.

{
    PetscFunctionBeginUser;
    *dm = NULL;
    *global_vec = NULL;
    *local_vec = NULL;
 
    if (strcmp(field_name, "Ucont") == 0) {
        *dm = user->fda;
        *global_vec = user->Ucont;
        *local_vec = user->lUcont;
    }
 
    PetscCheck(*dm && *global_vec && *local_vec, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE,
               "Field '%s' is not a registered persistent component-staggered field.",
               field_name);
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

TransferPeriodicStaggeredFieldByDirection()

PetscErrorCode TransferPeriodicStaggeredFieldByDirection	(	UserCtx *	user,
		const char *	field_name,
		char	periodic_direction
	)

Implementation of TransferPeriodicStaggeredFieldByDirection().

Transfers one persistent component-staggered field in one periodic direction.

Definition at line 1999 of file Boundaries.c.

{
    PetscErrorCode ierr;
    DMDALocalInfo info = user->info;
    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;
    DM dm;
    Vec global_vec, local_vec;
    Cmpnts ***global_array, ***local_array;
 
    PetscFunctionBeginUser;
    PetscCheck(periodic_direction == 'i' || periodic_direction == 'j' || periodic_direction == 'k',
               PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
               "Invalid periodic direction '%c'; expected 'i', 'j', or 'k'.", periodic_direction);
    PetscCall(GetPersistentStaggeredField(user, field_name, &dm, &global_vec, &local_vec));
 
    ierr = DMDAVecGetArray(dm, global_vec, &global_array); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(dm, local_vec, &local_array); CHKERRQ(ierr);
 
    if (periodic_direction == 'i') {
        if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && xs == 0)
            for (PetscInt k=zs; k<ze; k++) for (PetscInt j=ys; j<ye; j++) global_array[k][j][0] = local_array[k][j][-2];
        if (user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC && xe == mx)
            for (PetscInt k=zs; k<ze; k++) for (PetscInt j=ys; j<ye; j++) global_array[k][j][mx-1] = local_array[k][j][mx+1];
    } else if (periodic_direction == 'j') {
        if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && ys == 0)
            for (PetscInt k=zs; k<ze; k++) for (PetscInt i=xs; i<xe; i++) global_array[k][0][i] = local_array[k][-2][i];
        if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC && ye == my)
            for (PetscInt k=zs; k<ze; k++) for (PetscInt i=xs; i<xe; i++) global_array[k][my-1][i] = local_array[k][my+1][i];
    } else {
        if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && zs == 0)
            for (PetscInt j=ys; j<ye; j++) for (PetscInt i=xs; i<xe; i++) global_array[0][j][i] = local_array[-2][j][i];
        if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC && ze == mz)
            for (PetscInt j=ys; j<ye; j++) for (PetscInt i=xs; i<xe; i++) global_array[mz-1][j][i] = local_array[mz+1][j][i];
    }
 
    ierr = DMDAVecRestoreArrayRead(dm, local_vec, &local_array); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(dm, global_vec, &global_array); CHKERRQ(ierr);
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

SynchronizePeriodicStaggeredFields()

PetscErrorCode SynchronizePeriodicStaggeredFields	(	UserCtx *	user,
		PetscInt	num_fields,
		const char *	field_names[]
	)

Implementation of SynchronizePeriodicStaggeredFields().

Synchronizes persistent component-staggered vector fields.

Definition at line 2048 of file Boundaries.c.

{
    PetscErrorCode ierr;
    const char periodic_directions[3] = {'i', 'j', 'k'};
    const BCFace negative_faces[3] = {BC_FACE_NEG_X, BC_FACE_NEG_Y, BC_FACE_NEG_Z};
    const BCFace positive_faces[3] = {BC_FACE_POS_X, BC_FACE_POS_Y, BC_FACE_POS_Z};
 
    PetscFunctionBeginUser;
    PetscCheck(num_fields >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE,
               "Number of staggered fields cannot be negative.");
    if (num_fields == 0) PetscFunctionReturn(0);
    PetscCheck(field_names != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
               "Staggered field-name array cannot be NULL.");
 
    for (PetscInt field = 0; field < num_fields; field++) {
        DM dm;
        Vec global_vec, local_vec;
        PetscCall(GetPersistentStaggeredField(user, field_names[field], &dm, &global_vec, &local_vec));
        ierr = UpdateLocalGhosts(user, field_names[field]); CHKERRQ(ierr);
    }
 
    for (PetscInt direction = 0; direction < 3; direction++) {
        const PetscBool active =
            user->boundary_faces[negative_faces[direction]].mathematical_type == PERIODIC ||
            user->boundary_faces[positive_faces[direction]].mathematical_type == PERIODIC;
        if (!active) continue;
 
        for (PetscInt field = 0; field < num_fields; field++) {
            ierr = TransferPeriodicStaggeredFieldByDirection(user, field_names[field],
                                                             periodic_directions[direction]); CHKERRQ(ierr);
        }
        for (PetscInt field = 0; field < num_fields; field++) {
            ierr = UpdateLocalGhosts(user, field_names[field]); CHKERRQ(ierr);
        }
    }
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

PreparePeriodicQuickStencilFields()

PetscErrorCode PreparePeriodicQuickStencilFields	(	UserCtx *	user,
		Vec	local_vector_field,
		Vec	local_scalar_field
	)

Implementation of PreparePeriodicQuickStencilFields().

Repairs the outer adjacent periodic ghosts used by QUICK cell stencils.

Definition at line 2093 of file Boundaries.c.

{
    DMDALocalInfo info = user->info;
    Cmpnts ***vector_array;
    PetscReal ***scalar_array;
    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 gxs = info.gxs, gxe = info.gxs + info.gxm;
    const PetscInt gys = info.gys, gye = info.gys + info.gym;
    const PetscInt gzs = info.gzs, gze = info.gzs + info.gzm;
 
    PetscFunctionBeginUser;
    PetscCheck(local_vector_field && local_scalar_field, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
               "QUICK stencil repair requires both local vector and scalar fields.");
    PetscCall(DMDAVecGetArray(user->fda, local_vector_field, &vector_array));
    PetscCall(DMDAVecGetArray(user->da, local_scalar_field, &scalar_array));
 
    if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && xs == 0) {
        for (PetscInt k = gzs; k < gze; k++) for (PetscInt j = gys; j < gye; j++) {
            vector_array[k][j][-1] = vector_array[k][j][-3];
            scalar_array[k][j][-1] = scalar_array[k][j][-3];
        }
    }
    if (user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC && xe == info.mx) {
        for (PetscInt k = gzs; k < gze; k++) for (PetscInt j = gys; j < gye; j++) {
            vector_array[k][j][info.mx] = vector_array[k][j][info.mx + 2];
            scalar_array[k][j][info.mx] = scalar_array[k][j][info.mx + 2];
        }
    }
    if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && ys == 0) {
        for (PetscInt k = gzs; k < gze; k++) for (PetscInt i = gxs; i < gxe; i++) {
            vector_array[k][-1][i] = vector_array[k][-3][i];
            scalar_array[k][-1][i] = scalar_array[k][-3][i];
        }
    }
    if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC && ye == info.my) {
        for (PetscInt k = gzs; k < gze; k++) for (PetscInt i = gxs; i < gxe; i++) {
            vector_array[k][info.my][i] = vector_array[k][info.my + 2][i];
            scalar_array[k][info.my][i] = scalar_array[k][info.my + 2][i];
        }
    }
    if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && zs == 0) {
        for (PetscInt j = gys; j < gye; j++) for (PetscInt i = gxs; i < gxe; i++) {
            vector_array[-1][j][i] = vector_array[-3][j][i];
            scalar_array[-1][j][i] = scalar_array[-3][j][i];
        }
    }
    if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC && ze == info.mz) {
        for (PetscInt j = gys; j < gye; j++) for (PetscInt i = gxs; i < gxe; i++) {
            vector_array[info.mz][j][i] = vector_array[info.mz + 2][j][i];
            scalar_array[info.mz][j][i] = scalar_array[info.mz + 2][j][i];
        }
    }
 
    PetscCall(DMDAVecRestoreArray(user->da, local_scalar_field, &scalar_array));
    PetscCall(DMDAVecRestoreArray(user->fda, local_vector_field, &vector_array));
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

SynchronizePeriodicLocalStaggeredField()

PetscErrorCode SynchronizePeriodicLocalStaggeredField	(	UserCtx *	user,
		Vec	local_field
	)

Implementation of SynchronizePeriodicLocalStaggeredField().

Synchronizes one local-only component-staggered periodic work field.

Definition at line 2159 of file Boundaries.c.

{
    DMDALocalInfo info = user->info;
    Cmpnts ***array;
    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;
 
    PetscFunctionBeginUser;
    PetscCheck(local_field, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
               "Local staggered field cannot be NULL.");
    PetscCall(DMLocalToLocalBegin(user->fda, local_field, INSERT_VALUES, local_field));
    PetscCall(DMLocalToLocalEnd(user->fda, local_field, INSERT_VALUES, local_field));
    PetscCall(DMDAVecGetArray(user->fda, local_field, &array));
 
    if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && xs == 0)
        for (PetscInt k = zs; k < ze; k++) for (PetscInt j = ys; j < ye; j++) array[k][j][0].x = array[k][j][-2].x;
    if (user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC && xe == info.mx)
        for (PetscInt k = zs; k < ze; k++) for (PetscInt j = ys; j < ye; j++) array[k][j][info.mx - 1].x = array[k][j][info.mx + 1].x;
    if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && ys == 0)
        for (PetscInt k = zs; k < ze; k++) for (PetscInt i = xs; i < xe; i++) array[k][0][i].y = array[k][-2][i].y;
    if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC && ye == info.my)
        for (PetscInt k = zs; k < ze; k++) for (PetscInt i = xs; i < xe; i++) array[k][info.my - 1][i].y = array[k][info.my + 1][i].y;
    if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && zs == 0)
        for (PetscInt j = ys; j < ye; j++) for (PetscInt i = xs; i < xe; i++) array[0][j][i].z = array[-2][j][i].z;
    if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC && ze == info.mz)
        for (PetscInt j = ys; j < ye; j++) for (PetscInt i = xs; i < xe; i++) array[info.mz - 1][j][i].z = array[info.mz + 1][j][i].z;
 
    PetscCall(DMDAVecRestoreArray(user->fda, local_field, &array));
    PetscCall(DMLocalToLocalBegin(user->fda, local_field, INSERT_VALUES, local_field));
    PetscCall(DMLocalToLocalEnd(user->fda, local_field, INSERT_VALUES, local_field));
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

TransferPeriodicField()

PetscErrorCode TransferPeriodicField	(	UserCtx *	user,
		const char *	field_name
	)

Internal helper implementation: TransferPeriodicField().

Legacy monolithic periodic endpoint transfer for one cell-centered field.

Local to this translation unit.

Definition at line 2199 of file Boundaries.c.

{
    PetscErrorCode ierr;
    DMDALocalInfo  info = user->info;
    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;
 
    // --- Local variables to hold the specific details of the chosen field ---
    DM        dm;
    Vec       global_vec;
    Vec       local_vec;
    PetscInt  dof;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    // --- STEP 1: Dispatcher - Set the specific DM, Vecs, and dof based on field_name ---
    // Compatibility implementation retained for legacy callers. New periodic
    // cell fields should be added to the directional dispatcher used by
    // SynchronizePeriodicCellFields().
    if (strcmp(field_name, "Ucat") == 0) {
        dm         = user->fda;
        global_vec = user->Ucat;
        local_vec  = user->lUcat;
        dof        = 3;
    } else if (strcmp(field_name, "P") == 0) {
        dm         = user->da;
        global_vec = user->P;
        local_vec  = user->lP;
        dof        = 1;
    } else if (strcmp(field_name, "Nvert") == 0) {
        dm         = user->da;
        global_vec = user->Nvert;
        local_vec  = user->lNvert;
        dof        = 1;
    } else if (strcmp(field_name, "Eddy Viscosity") == 0) {
        dm         = user->da;
        global_vec = user->Nu_t;
        local_vec  = user->lNu_t;
        dof        = 1;
    }
    /*
    // Example for future extension:
    else if (strcmp(field_name, "Temperature") == 0) {
        dm         = user->da; // Assuming Temperature is scalar
        global_vec = user->T;
        local_vec  = user->lT;
        dof        = 1;
    }
    */
    else {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown field name '%s' in TransferPeriodicFieldByName.", field_name);
    }
 
    LOG_ALLOW(GLOBAL,LOG_TRACE,"Periodic Transform being performed for field: %s with %d DoF.\n",field_name,dof);
    // --- STEP 2: Execute the copy logic using the dispatched variables ---
    if (dof == 1) { // --- Handle SCALAR fields (PetscReal) ---
        PetscReal ***g_array, ***l_array;
        ierr = DMDAVecGetArray(dm, global_vec, &g_array); CHKERRQ(ierr);
        ierr = DMDAVecGetArray(dm, local_vec, &l_array); CHKERRQ(ierr);
 
        if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && xs == 0) for (PetscInt k=zs; k<ze; k++) for (PetscInt j=ys; j<ye; j++) g_array[k][j][xs] = l_array[k][j][xs-2];
        if (user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC && xe == mx) for (PetscInt k=zs; k<ze; k++) for (PetscInt j=ys; j<ye; j++) g_array[k][j][xe-1] = l_array[k][j][xe+1];
        if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && ys == 0) for (PetscInt k=zs; k<ze; k++) for (PetscInt i=xs; i<xe; i++) g_array[k][ys][i] = l_array[k][ys-2][i];
        if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC && ye == my) for (PetscInt k=zs; k<ze; k++) for (PetscInt i=xs; i<xe; i++) g_array[k][ye-1][i] = l_array[k][ye+1][i];
        if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && zs == 0) for (PetscInt j=ys; j<ye; j++) for (PetscInt i=xs; i<xe; i++) g_array[zs][j][i] = l_array[zs-2][j][i];
        if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC && ze == mz) for (PetscInt j=ys; j<ye; j++) for (PetscInt i=xs; i<xe; i++) g_array[ze-1][j][i] = l_array[ze+1][j][i];
        
        ierr = DMDAVecRestoreArray(dm, global_vec, &g_array); CHKERRQ(ierr);
        ierr = DMDAVecRestoreArray(dm, local_vec, &l_array); CHKERRQ(ierr);
 
    } else if (dof == 3) { // --- Handle VECTOR fields (Cmpnts) ---
        Cmpnts ***g_array, ***l_array;
        ierr = DMDAVecGetArray(dm, global_vec, &g_array); CHKERRQ(ierr);
        ierr = DMDAVecGetArray(dm, local_vec, &l_array); CHKERRQ(ierr);
 
        LOG_ALLOW(GLOBAL,LOG_VERBOSE,"Array %s read successfully (Global and Local).\n",field_name);
 
        if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && xs == 0) for (PetscInt k=zs; k<ze; k++) for (PetscInt j=ys; j<ye; j++) g_array[k][j][xs] = l_array[k][j][xs-2];
        if (user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC && xe == mx) for (PetscInt k=zs; k<ze; k++) for (PetscInt j=ys; j<ye; j++) g_array[k][j][xe-1] = l_array[k][j][xe+1];
        if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && ys == 0) for (PetscInt k=zs; k<ze; k++) for (PetscInt i=xs; i<xe; i++) g_array[k][ys][i] = l_array[k][ys-2][i];
        if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC && ye == my) for (PetscInt k=zs; k<ze; k++) for (PetscInt i=xs; i<xe; i++) g_array[k][ye-1][i] = l_array[k][ye+1][i];
        if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && zs == 0) for (PetscInt j=ys; j<ye; j++) for (PetscInt i=xs; i<xe; i++) g_array[zs][j][i] = l_array[zs-2][j][i];
        if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC && ze == mz) for (PetscInt j=ys; j<ye; j++) for (PetscInt i=xs; i<xe; i++) g_array[ze-1][j][i] = l_array[ze+1][j][i];
        
        ierr = DMDAVecRestoreArray(dm, global_vec, &g_array); CHKERRQ(ierr);
        ierr = DMDAVecRestoreArray(dm, local_vec, &l_array); CHKERRQ(ierr);
    }
    else{
         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "This function only accepts Fields with 1 or 3 DoF.");
    }
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

TransferPeriodicFaceField()

PetscErrorCode TransferPeriodicFaceField	(	UserCtx *	user,
		const char *	field_name
	)

Internal helper implementation: TransferPeriodicFaceField().

(Primitive) Copies periodic data from the interior to the local ghost cell region for a single field.

Local to this translation unit.

Definition at line 2303 of file Boundaries.c.

{
    PetscErrorCode ierr;
    DMDALocalInfo  info = user->info;
    PetscInt       gxs = info.gxs, gxe = info.gxs + info.gxm;
    PetscInt       gys = info.gys, gye = info.gys + info.gym;
    PetscInt       gzs = info.gzs, gze = info.gzs + info.gzm;
    PetscInt       mx = info.mx, my = info.my, mz = info.mz;
 
    // --- Dispatcher to get the correct DM, Vec, and DoF for the specified field ---
    // Field-extension note: add one case here when a new field needs periodic
    // ghost-face transfer (typically for stencil-heavy operators).
    DM          dm;
    Vec         local_vec;
    PetscInt    dof;
    // (This dispatcher contains all 17 potential fields)
    if      (strcmp(field_name, "Ucont") == 0) { dm = user->fda; local_vec = user->lUcont; dof = 3; }
    else if (strcmp(field_name, "Csi")   == 0) { dm = user->fda; local_vec = user->lCsi;   dof = 3; }
    else if (strcmp(field_name, "Eta")   == 0) { dm = user->fda; local_vec = user->lEta;   dof = 3; }
    else if (strcmp(field_name, "Zet")   == 0) { dm = user->fda; local_vec = user->lZet;   dof = 3; }
    else if (strcmp(field_name, "ICsi")   == 0) { dm = user->fda; local_vec = user->lICsi;   dof = 3; }
    else if (strcmp(field_name, "IEta")   == 0) { dm = user->fda; local_vec = user->lIEta;   dof = 3; }
    else if (strcmp(field_name, "IZet")   == 0) { dm = user->fda; local_vec = user->lIZet;   dof = 3; }
    else if (strcmp(field_name, "JCsi")   == 0) { dm = user->fda; local_vec = user->lJCsi;   dof = 3; }
    else if (strcmp(field_name, "JEta")   == 0) { dm = user->fda; local_vec = user->lJEta;   dof = 3; }
    else if (strcmp(field_name, "JZet")   == 0) { dm = user->fda; local_vec = user->lJZet;   dof = 3; }
    else if (strcmp(field_name, "KCsi")   == 0) { dm = user->fda; local_vec = user->lKCsi;   dof = 3; }
    else if (strcmp(field_name, "KEta")   == 0) { dm = user->fda; local_vec = user->lKEta;   dof = 3; }
    else if (strcmp(field_name, "KZet")   == 0) { dm = user->fda; local_vec = user->lKZet;   dof = 3; }
    else if (strcmp(field_name, "Aj")     == 0) { dm = user->da;  local_vec = user->lAj;   dof = 1; }
    else if (strcmp(field_name, "IAj")   == 0) { dm = user->da;  local_vec = user->lIAj;   dof = 1; }
    else if (strcmp(field_name, "JAj")   == 0) { dm = user->da;  local_vec = user->lJAj;   dof = 1; }
    else if (strcmp(field_name, "KAj")   == 0) { dm = user->da;  local_vec = user->lKAj;   dof = 1; }
    else {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown field name '%s' in TransferPeriodicFaceField.", field_name);
    }
 
    PetscFunctionBeginUser;
 
    void *l_array_ptr;
    ierr = DMDAVecGetArray(dm, local_vec, &l_array_ptr); CHKERRQ(ierr);
 
    // --- I-DIRECTION ---
    if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC) {
        for (PetscInt k=gzs; k<gze; k++) for (PetscInt j=gys; j<gye; j++) {
            if (dof == 1) {
                PetscReal ***arr = (PetscReal***)l_array_ptr;
                arr[k][j][0]   = arr[k][j][mx-2];
                arr[k][j][-1]  = arr[k][j][mx-3];
            } else {
                Cmpnts ***arr = (Cmpnts***)l_array_ptr;
                arr[k][j][0]   = arr[k][j][mx-2];
                arr[k][j][-1]  = arr[k][j][mx-3];
            }
        }
    }
    if (user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC) {
        for (PetscInt k=gzs; k<gze; k++) for (PetscInt j=gys; j<gye; j++) {
             if (dof == 1) {
                PetscReal ***arr = (PetscReal***)l_array_ptr;
                arr[k][j][mx-1] = arr[k][j][1];
                arr[k][j][mx]   = arr[k][j][2];
            } else {
                Cmpnts ***arr = (Cmpnts***)l_array_ptr;
                arr[k][j][mx-1] = arr[k][j][1];
                arr[k][j][mx]   = arr[k][j][2];
            }
        }
    }
 
    // --- J-DIRECTION ---
    if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC) {
        for (PetscInt k=gzs; k<gze; k++) for (PetscInt i=gxs; i<gxe; i++) {
             if (dof == 1) {
                PetscReal ***arr = (PetscReal***)l_array_ptr;
                arr[k][0][i]   = arr[k][my-2][i];
                arr[k][-1][i]  = arr[k][my-3][i];
            } else {
                Cmpnts ***arr = (Cmpnts***)l_array_ptr;
                arr[k][0][i]   = arr[k][my-2][i];
                arr[k][-1][i]  = arr[k][my-3][i];
            }
        }
    }
    if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC) {
        for (PetscInt k=gzs; k<gze; k++) for (PetscInt i=gxs; i<gxe; i++) {
             if (dof == 1) {
                PetscReal ***arr = (PetscReal***)l_array_ptr;
                arr[k][my-1][i] = arr[k][1][i];
                arr[k][my][i]   = arr[k][2][i];
            } else {
                Cmpnts ***arr = (Cmpnts***)l_array_ptr;
                arr[k][my-1][i] = arr[k][1][i];
                arr[k][my][i]   = arr[k][2][i];
            }
        }
    }
 
    // --- K-DIRECTION ---
    if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC) {
        for (PetscInt j=gys; j<gye; j++) for (PetscInt i=gxs; i<gxe; i++) {
             if (dof == 1) {
                PetscReal ***arr = (PetscReal***)l_array_ptr;
                arr[0][j][i]   = arr[mz-2][j][i];
                arr[-1][j][i]  = arr[mz-3][j][i];
            } else {
                Cmpnts ***arr = (Cmpnts***)l_array_ptr;
                arr[0][j][i]   = arr[mz-2][j][i];
                arr[-1][j][i]  = arr[mz-3][j][i];
            }
        }
    }
    if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC) {
        for (PetscInt j=gys; j<gye; j++) for (PetscInt i=gxs; i<gxe; i++) {
             if (dof == 1) {
                PetscReal ***arr = (PetscReal***)l_array_ptr;
                arr[mz-1][j][i] = arr[1][j][i];
                arr[mz][j][i]   = arr[2][j][i];
            } else {
                Cmpnts ***arr = (Cmpnts***)l_array_ptr;
                arr[mz-1][j][i] = arr[1][j][i];
                arr[mz][j][i]   = arr[2][j][i];
            }
        }
    }
 
    ierr = DMDAVecRestoreArray(dm, local_vec, &l_array_ptr); CHKERRQ(ierr);
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

ApplyMetricsPeriodicBCs()

PetscErrorCode ApplyMetricsPeriodicBCs

(

UserCtx *

user

)

Internal helper implementation: ApplyMetricsPeriodicBCs().

(Orchestrator) Updates all metric-related fields in the local ghost cell regions for periodic boundaries.

Local to this translation unit.

Definition at line 2439 of file Boundaries.c.

{
    PetscErrorCode ierr;
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    const char *cell_fields[] = {"Aj"};
    const char *i_face_fields[] = {"Centx", "Csi", "ICsi", "IEta", "IZet", "IAj"};
    const char *j_face_fields[] = {"Centy", "Eta", "JCsi", "JEta", "JZet", "JAj"};
    const char *k_face_fields[] = {"Centz", "Zet", "KCsi", "KEta", "KZet", "KAj"};
 
    ierr = SynchronizePeriodicCellFields(user, 1, cell_fields); CHKERRQ(ierr);
    ierr = SynchronizePeriodicFaceFields(user, 'i', 6, i_face_fields); CHKERRQ(ierr);
    ierr = SynchronizePeriodicFaceFields(user, 'j', 6, j_face_fields); CHKERRQ(ierr);
    ierr = SynchronizePeriodicFaceFields(user, 'k', 6, k_face_fields); CHKERRQ(ierr);
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

ApplyPeriodicBCs()

PetscErrorCode ApplyPeriodicBCs

(

UserCtx *

user

)

Internal helper implementation: ApplyPeriodicBCs().

Applies periodic boundary conditions by copying data across domain boundaries for all relevant fields.

Local to this translation unit.

Definition at line 2465 of file Boundaries.c.

{
    PetscErrorCode ierr;
    PetscBool is_any_periodic = PETSC_FALSE;
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    for (int i = 0; i < 6; i++) {
        if (user->boundary_faces[i].mathematical_type == PERIODIC) {
            is_any_periodic = PETSC_TRUE;
            break;
        }
    }
 
    if (!is_any_periodic) {
        LOG_ALLOW(GLOBAL,LOG_TRACE, "No periodic boundaries defined; skipping ApplyPeriodicBCs.\n");
        PROFILE_FUNCTION_END;
        PetscFunctionReturn(0);
    }
 
    LOG_ALLOW(GLOBAL, LOG_TRACE, "Applying periodic boundary conditions for all fields.\n");
 
    // STEP 1: Synchronize periodic cell-centered fields in deterministic direction order.
    const char *cell_fields[] = {"Ucat", "P", "Nvert"};
    ierr = SynchronizePeriodicCellFields(user, 3, cell_fields); CHKERRQ(ierr);
 
    /* if (user->solve_temperature) { ierr = UpdateLocalGhosts(user, "Temperature"); CHKERRQ(ierr); } */
 
    // STEP 2: Synchronize persistent staggered endpoints and repair local
    // component-normal ghosts through UpdateLocalGhosts().
    const char *staggered_fields[] = {"Ucont"};
    ierr = SynchronizePeriodicStaggeredFields(user, 1, staggered_fields); CHKERRQ(ierr);
 
    // FUTURE EXTENSION: Add new cell fields through SynchronizePeriodicCellFields().
    /*
    if (user->solve_temperature) {
        const char *temperature_field[] = {"Temperature"};
        ierr = SynchronizePeriodicCellFields(user, 1, temperature_field); CHKERRQ(ierr);
    }
    */
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

UpdateDummyCells()

PetscErrorCode UpdateDummyCells

(

UserCtx *

user

)

Internal helper implementation: UpdateDummyCells().

Updates the dummy cells (ghost nodes) on the faces of the local domain for NON-PERIODIC boundaries.

Local to this translation unit.

Definition at line 2518 of file Boundaries.c.

{
    PetscErrorCode ierr;
    DM            fda = user->fda;
    DMDALocalInfo info = user->info;
    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;
 
    // --- Calculate shrunken loop ranges to avoid edges and corners ---
    PetscInt lxs = xs, lxe = xe;
    PetscInt lys = ys, lye = ye;
    PetscInt lzs = zs, lze = ze;
 
    if (xs == 0) lxs = xs + 1;
    if (ys == 0) lys = ys + 1;
    if (zs == 0) lzs = zs + 1;
 
    if (xe == mx) lxe = xe - 1;
    if (ye == my) lye = ye - 1;
    if (ze == mz) lze = ze - 1;
 
    Cmpnts        ***ucat, ***ubcs;
    PetscFunctionBeginUser;
 
    ierr = DMDAVecGetArray(fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
    ierr = DMDAVecGetArray(fda, user->Ucat, &ucat); CHKERRQ(ierr);
 
    // -X Face
    if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC && xs == 0) {
        for (PetscInt k = lzs; k < lze; k++) for (PetscInt j = lys; j < lye; j++) {
            ucat[k][j][xs].x = 2.0 * ubcs[k][j][xs].x - ucat[k][j][xs + 1].x;
            ucat[k][j][xs].y = 2.0 * ubcs[k][j][xs].y - ucat[k][j][xs + 1].y;
            ucat[k][j][xs].z = 2.0 * ubcs[k][j][xs].z - ucat[k][j][xs + 1].z;
        }
    }
    // +X Face
    if (user->boundary_faces[BC_FACE_POS_X].mathematical_type != PERIODIC && xe == mx) {
        for (PetscInt k = lzs; k < lze; k++) for (PetscInt j = lys; j < lye; j++) {
            ucat[k][j][xe-1].x = 2.0 * ubcs[k][j][xe-1].x - ucat[k][j][xe - 2].x;
            ucat[k][j][xe-1].y = 2.0 * ubcs[k][j][xe-1].y - ucat[k][j][xe - 2].y;
            ucat[k][j][xe-1].z = 2.0 * ubcs[k][j][xe-1].z - ucat[k][j][xe - 2].z;
        }
    }
 
    // -Y Face
    if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC && ys == 0) {
        for (PetscInt k = lzs; k < lze; k++) for (PetscInt i = lxs; i < lxe; i++) {
            ucat[k][ys][i].x = 2.0 * ubcs[k][ys][i].x - ucat[k][ys + 1][i].x;
            ucat[k][ys][i].y = 2.0 * ubcs[k][ys][i].y - ucat[k][ys + 1][i].y;
            ucat[k][ys][i].z = 2.0 * ubcs[k][ys][i].z - ucat[k][ys + 1][i].z;
        }
    }
    // +Y Face
    if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type != PERIODIC && ye == my) {
        for (PetscInt k = lzs; k < lze; k++) for (PetscInt i = lxs; i < lxe; i++) {
            ucat[k][ye-1][i].x = 2.0 * ubcs[k][ye-1][i].x - ucat[k][ye-2][i].x;
            ucat[k][ye-1][i].y = 2.0 * ubcs[k][ye-1][i].y - ucat[k][ye-2][i].y;
            ucat[k][ye-1][i].z = 2.0 * ubcs[k][ye-1][i].z - ucat[k][ye-2][i].z;
        }
    }
 
    // -Z Face
    if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC && zs == 0) {
        for (PetscInt j = lys; j < lye; j++) for (PetscInt i = lxs; i < lxe; i++) {
            ucat[zs][j][i].x = 2.0 * ubcs[zs][j][i].x - ucat[zs + 1][j][i].x;
            ucat[zs][j][i].y = 2.0 * ubcs[zs][j][i].y - ucat[zs + 1][j][i].y;
            ucat[zs][j][i].z = 2.0 * ubcs[zs][j][i].z - ucat[zs + 1][j][i].z;
        }
    }
    // +Z Face
    if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type != PERIODIC && ze == mz) {
        for (PetscInt j = lys; j < lye; j++) for (PetscInt i = lxs; i < lxe; i++) {
            ucat[ze-1][j][i].x = 2.0 * ubcs[ze-1][j][i].x - ucat[ze-2][j][i].x;
            ucat[ze-1][j][i].y = 2.0 * ubcs[ze-1][j][i].y - ucat[ze-2][j][i].y;
            ucat[ze-1][j][i].z = 2.0 * ubcs[ze-1][j][i].z - ucat[ze-2][j][i].z;
        }
    }
 
    ierr = DMDAVecRestoreArray(fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(fda, user->Ucat, &ucat); CHKERRQ(ierr);
 
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

UpdateCornerNodes()

PetscErrorCode UpdateCornerNodes

(

UserCtx *

user

)

Internal helper implementation: UpdateCornerNodes().

Updates the corner and edge ghost nodes of the local domain by averaging.

Local to this translation unit.

Definition at line 2610 of file Boundaries.c.

{
    PetscErrorCode ierr;
    DM            da = user->da, fda = user->fda;
    DMDALocalInfo info = user->info;
    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;
 
    Cmpnts        ***ucat;
    PetscReal     ***p;
 
    PetscFunctionBeginUser;
 
    ierr = DMDAVecGetArray(fda, user->Ucat, &ucat); CHKERRQ(ierr);
    ierr = DMDAVecGetArray(da, user->P, &p); CHKERRQ(ierr);
 
    // --- Update Edges and Corners by Averaging ---
    // The order of these blocks ensures that corners (where 3 faces meet) are
    // computed using data from edges (where 2 faces meet), which are computed first.
// Edges connected to the -Z face (k=zs)
  if (zs == 0) {
      if (xs == 0) {
          for (PetscInt j = ys; j < ye; j++) {
              p[zs][j][xs] = 0.5 * (p[zs+1][j][xs] + p[zs][j][xs+1]);
              ucat[zs][j][xs].x = 0.5 * (ucat[zs+1][j][xs].x + ucat[zs][j][xs+1].x);
              ucat[zs][j][xs].y = 0.5 * (ucat[zs+1][j][xs].y + ucat[zs][j][xs+1].y);
              ucat[zs][j][xs].z = 0.5 * (ucat[zs+1][j][xs].z + ucat[zs][j][xs+1].z);
          }
      }
      if (xe == mx) {
          for (PetscInt j = ys; j < ye; j++) {
              p[zs][j][mx-1] = 0.5 * (p[zs+1][j][mx-1] + p[zs][j][mx-2]);
              ucat[zs][j][mx-1].x = 0.5 * (ucat[zs+1][j][mx-1].x + ucat[zs][j][mx-2].x);
              ucat[zs][j][mx-1].y = 0.5 * (ucat[zs+1][j][mx-1].y + ucat[zs][j][mx-2].y);
              ucat[zs][j][mx-1].z = 0.5 * (ucat[zs+1][j][mx-1].z + ucat[zs][j][mx-2].z);
          }
      }
      if (ys == 0) {
          for (PetscInt i = xs; i < xe; i++) {
              p[zs][ys][i] = 0.5 * (p[zs+1][ys][i] + p[zs][ys+1][i]);
              ucat[zs][ys][i].x = 0.5 * (ucat[zs+1][ys][i].x + ucat[zs][ys+1][i].x);
              ucat[zs][ys][i].y = 0.5 * (ucat[zs+1][ys][i].y + ucat[zs][ys+1][i].y);
              ucat[zs][ys][i].z = 0.5 * (ucat[zs+1][ys][i].z + ucat[zs][ys+1][i].z);
          }
      }
      if (ye == my) {
          for (PetscInt i = xs; i < xe; i++) {
              p[zs][my-1][i] = 0.5 * (p[zs+1][my-1][i] + p[zs][my-2][i]);
              ucat[zs][my-1][i].x = 0.5 * (ucat[zs+1][my-1][i].x + ucat[zs][my-2][i].x);
              ucat[zs][my-1][i].y = 0.5 * (ucat[zs+1][my-1][i].y + ucat[zs][my-2][i].y);
              ucat[zs][my-1][i].z = 0.5 * (ucat[zs+1][my-1][i].z + ucat[zs][my-2][i].z);
          }
      }
  }
 
  // Edges connected to the +Z face (k=ze-1)
  if (ze == mz) {
      if (xs == 0) {
          for (PetscInt j = ys; j < ye; j++) {
              p[mz-1][j][xs] = 0.5 * (p[mz-2][j][xs] + p[mz-1][j][xs+1]);
              ucat[mz-1][j][xs].x = 0.5 * (ucat[mz-2][j][xs].x + ucat[mz-1][j][xs+1].x);
              ucat[mz-1][j][xs].y = 0.5 * (ucat[mz-2][j][xs].y + ucat[mz-1][j][xs+1].y);
              ucat[mz-1][j][xs].z = 0.5 * (ucat[mz-2][j][xs].z + ucat[mz-1][j][xs+1].z);
          }
      }
      if (xe == mx) {
          for (PetscInt j = ys; j < ye; j++) {
              p[mz-1][j][mx-1] = 0.5 * (p[mz-2][j][mx-1] + p[mz-1][j][mx-2]);
              ucat[mz-1][j][mx-1].x = 0.5 * (ucat[mz-2][j][mx-1].x + ucat[mz-1][j][mx-2].x);
              ucat[mz-1][j][mx-1].y = 0.5 * (ucat[mz-2][j][mx-1].y + ucat[mz-1][j][mx-2].y);
              ucat[mz-1][j][mx-1].z = 0.5 * (ucat[mz-2][j][mx-1].z + ucat[mz-1][j][mx-2].z);
          }
      }
      if (ys == 0) {
          for (PetscInt i = xs; i < xe; i++) {
              p[mz-1][ys][i] = 0.5 * (p[mz-2][ys][i] + p[mz-1][ys+1][i]);
              ucat[mz-1][ys][i].x = 0.5 * (ucat[mz-2][ys][i].x + ucat[mz-1][ys+1][i].x);
              ucat[mz-1][ys][i].y = 0.5 * (ucat[mz-2][ys][i].y + ucat[mz-1][ys+1][i].y);
              ucat[mz-1][ys][i].z = 0.5 * (ucat[mz-2][ys][i].z + ucat[mz-1][ys+1][i].z);
          }
      }
      if (ye == my) {
          for (PetscInt i = xs; i < xe; i++) {
              p[mz-1][my-1][i] = 0.5 * (p[mz-2][my-1][i] + p[mz-1][my-2][i]);
              ucat[mz-1][my-1][i].x = 0.5 * (ucat[mz-2][my-1][i].x + ucat[mz-1][my-2][i].x);
              ucat[mz-1][my-1][i].y = 0.5 * (ucat[mz-2][my-1][i].y + ucat[mz-1][my-2][i].y);
              ucat[mz-1][my-1][i].z = 0.5 * (ucat[mz-2][my-1][i].z + ucat[mz-1][my-2][i].z);
          }
      }
  }
 
  // Remaining edges on the XY plane (that are not on Z faces)
  if (ys == 0) {
      if (xs == 0) {
          for (PetscInt k = zs; k < ze; k++) {
              p[k][ys][xs] = 0.5 * (p[k][ys+1][xs] + p[k][ys][xs+1]);
              ucat[k][ys][xs].x = 0.5 * (ucat[k][ys+1][xs].x + ucat[k][ys][xs+1].x);
              ucat[k][ys][xs].y = 0.5 * (ucat[k][ys+1][xs].y + ucat[k][ys][xs+1].y);
              ucat[k][ys][xs].z = 0.5 * (ucat[k][ys+1][xs].z + ucat[k][ys][xs+1].z);
          }
      }
      if (xe == mx) {
          for (PetscInt k = zs; k < ze; k++) {
              p[k][ys][mx-1] = 0.5 * (p[k][ys+1][mx-1] + p[k][ys][mx-2]);
              ucat[k][ys][mx-1].x = 0.5 * (ucat[k][ys+1][mx-1].x + ucat[k][ys][mx-2].x);
              ucat[k][ys][mx-1].y = 0.5 * (ucat[k][ys+1][mx-1].y + ucat[k][ys][mx-2].y);
              ucat[k][ys][mx-1].z = 0.5 * (ucat[k][ys+1][mx-1].z + ucat[k][ys][mx-2].z);
          }
      }
  }
 
  if (ye == my) {
      if (xs == 0) {
          for (PetscInt k = zs; k < ze; k++) {
              p[k][my-1][xs] = 0.5 * (p[k][my-2][xs] + p[k][my-1][xs+1]);
              ucat[k][my-1][xs].x = 0.5 * (ucat[k][my-2][xs].x + ucat[k][my-1][xs+1].x);
              ucat[k][my-1][xs].y = 0.5 * (ucat[k][my-2][xs].y + ucat[k][my-1][xs+1].y);
              ucat[k][my-1][xs].z = 0.5 * (ucat[k][my-2][xs].z + ucat[k][my-1][xs+1].z);
          }
      }
      if (xe == mx) {
          for (PetscInt k = zs; k < ze; k++) {
              p[k][my-1][mx-1] = 0.5 * (p[k][my-2][mx-1] + p[k][my-1][mx-2]);
              ucat[k][my-1][mx-1].x = 0.5 * (ucat[k][my-2][mx-1].x + ucat[k][my-1][mx-2].x);
              ucat[k][my-1][mx-1].y = 0.5 * (ucat[k][my-2][mx-1].y + ucat[k][my-1][mx-2].y);
              ucat[k][my-1][mx-1].z = 0.5 * (ucat[k][my-2][mx-1].z + ucat[k][my-1][mx-2].z);
          }
      }
  }
 
    ierr = DMDAVecRestoreArray(fda, user->Ucat, &ucat); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(da, user->P, &p); CHKERRQ(ierr);
 
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

UpdatePeriodicCornerNodes()

PetscErrorCode UpdatePeriodicCornerNodes	(	UserCtx *	user,
		PetscInt	num_fields,
		const char *	field_names[]
	)

Internal helper implementation: UpdatePeriodicCornerNodes().

Legacy sequential periodic-corner update for a list of fields.

Local to this translation unit.

Definition at line 2754 of file Boundaries.c.

{
    PetscErrorCode ierr;
    PetscFunctionBeginUser;
 
    if (num_fields == 0) PetscFunctionReturn(0);
 
    // --- I-DIRECTION ---
    for (PetscInt i = 0; i < num_fields; i++) {
        ierr = TransferPeriodicFieldByDirection(user, field_names[i], 'i'); CHKERRQ(ierr);
    }
    // --- SYNC ---
    for (PetscInt i = 0; i < num_fields; i++) {
        ierr = UpdateLocalGhosts(user, field_names[i]); CHKERRQ(ierr);
    }
 
    // --- J-DIRECTION ---
    for (PetscInt i = 0; i < num_fields; i++) {
        ierr = TransferPeriodicFieldByDirection(user, field_names[i], 'j'); CHKERRQ(ierr);
    }
    // --- SYNC ---
    for (PetscInt i = 0; i < num_fields; i++) {
        ierr = UpdateLocalGhosts(user, field_names[i]); CHKERRQ(ierr);
    }
 
    // --- K-DIRECTION ---
    for (PetscInt i = 0; i < num_fields; i++) {
        ierr = TransferPeriodicFieldByDirection(user, field_names[i], 'k'); CHKERRQ(ierr);
    }
    // --- FINAL SYNC ---
    for (PetscInt i = 0; i < num_fields; i++) {
        ierr = UpdateLocalGhosts(user, field_names[i]); CHKERRQ(ierr);
    }
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

ApplyWallFunction()

PetscErrorCode ApplyWallFunction

(

UserCtx *

user

)

Internal helper implementation: ApplyWallFunction().

Applies wall function modeling to near-wall velocities for all wall-type boundaries.

Local to this translation unit.

Definition at line 2797 of file Boundaries.c.

{
    PetscErrorCode ierr;
    SimCtx *simCtx = user->simCtx;
    DMDALocalInfo *info = &user->info;
    
    PetscFunctionBeginUser;
    
    // =========================================================================
    // STEP 0: Early exit if wall functions are disabled
    // =========================================================================
    if (!simCtx->wallfunction) {
        PetscFunctionReturn(0);
    }
    
    LOG_ALLOW(LOCAL, LOG_DEBUG, "Processing wall function boundaries.\n");
    
    // =========================================================================
    // STEP 1: Get read/write access to all necessary field arrays
    // =========================================================================
    Cmpnts ***velocity_cartesian;           // Cartesian velocity (modified)
    Cmpnts ***velocity_contravariant;       // Contravariant velocity (set to zero at walls)
    Cmpnts ***velocity_boundary;            // Boundary condition velocity (kept at zero)
    Cmpnts ***csi, ***eta, ***zet;         // Metric tensor components (face normals)
    PetscReal ***node_vertex_flag;          // Fluid/solid indicator (0=fluid, 1=solid)
    PetscReal ***cell_jacobian;             // Grid Jacobian (1/volume)
    PetscReal ***friction_velocity;         // u_tau (friction velocity field)
    
    ierr = DMDAVecGetArray(user->fda, user->Ucat, &velocity_cartesian); CHKERRQ(ierr);
    ierr = DMDAVecGetArray(user->fda, user->Ucont, &velocity_contravariant); CHKERRQ(ierr);
    ierr = DMDAVecGetArray(user->fda, user->Bcs.Ubcs, &velocity_boundary); 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***)&node_vertex_flag); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->da, user->lAj, (const PetscReal***)&cell_jacobian); CHKERRQ(ierr);
    ierr = DMDAVecGetArray(user->da, user->lFriction_Velocity, &friction_velocity); CHKERRQ(ierr);
    
    // =========================================================================
    // STEP 2: Define loop bounds (owned portion of the grid for this MPI rank)
    // =========================================================================
    PetscInt grid_start_i = info->xs, grid_end_i = info->xs + info->xm;
    PetscInt grid_start_j = info->ys, grid_end_j = info->ys + info->ym;
    PetscInt grid_start_k = info->zs, grid_end_k = info->zs + info->zm;
    PetscInt grid_size_i = info->mx, grid_size_j = info->my, grid_size_k = info->mz;
    
    // Shrunken loop bounds: exclude domain edges and corners to avoid double-counting
    PetscInt loop_start_i = grid_start_i, loop_end_i = grid_end_i;
    PetscInt loop_start_j = grid_start_j, loop_end_j = grid_end_j;
    PetscInt loop_start_k = grid_start_k, loop_end_k = grid_end_k;
    
    if (grid_start_i == 0) loop_start_i = grid_start_i + 1;
    if (grid_end_i == grid_size_i) loop_end_i = grid_end_i - 1;
    if (grid_start_j == 0) loop_start_j = grid_start_j + 1;
    if (grid_end_j == grid_size_j) loop_end_j = grid_end_j - 1;
    if (grid_start_k == 0) loop_start_k = grid_start_k + 1;
    if (grid_end_k == grid_size_k) loop_end_k = grid_end_k - 1;
    
    // Wall roughness parameter (smooth wall by default, configurable via -wall_roughness).
    const PetscReal wall_roughness_height = user->simCtx->wall_roughness_height;
    
    // =========================================================================
    // STEP 3: Process each of the 6 domain faces
    // =========================================================================
    for (int face_index = 0; face_index < 6; face_index++) {
        BCFace current_face_id = (BCFace)face_index;
        BoundaryFaceConfig *face_config = &user->boundary_faces[current_face_id];
        
        // Only process faces that are mathematical walls (applies to no-slip, moving, slip, etc.)
        if (face_config->mathematical_type != WALL) {
            continue;
        }
        
        // Check if this MPI rank owns part of this face
        PetscBool rank_owns_this_face;
        ierr = CanRankServiceFace(info, user->IM, user->JM, user->KM, 
                                  current_face_id, &rank_owns_this_face); CHKERRQ(ierr);
        
        if (!rank_owns_this_face) {
            continue;
        }
        
        LOG_ALLOW(LOCAL, LOG_TRACE, "Processing Face %d (%s)\n",
                  current_face_id, BCFaceToString(current_face_id));
        
        // =====================================================================
        // Process each face with appropriate indexing
        // =====================================================================
        switch(current_face_id) {
            
            // =================================================================
            // NEGATIVE X FACE (i = 0, normal points in +X direction)
            // =================================================================
            case BC_FACE_NEG_X: {
                if (grid_start_i == 0) {
                    const PetscInt ghost_cell_index = grid_start_i;
                    const PetscInt first_interior_cell = grid_start_i + 1;
                    const PetscInt second_interior_cell = grid_start_i + 2;
                    
                    for (PetscInt k = loop_start_k; k < loop_end_k; k++) {
                        for (PetscInt j = loop_start_j; j < loop_end_j; j++) {
                            
                            // Skip if this is a solid cell (embedded boundary)
                            if (node_vertex_flag[k][j][first_interior_cell] < 0.1) {
                                
                                // Calculate face area from contravariant metric tensor
                                PetscReal face_area = sqrt(
                                    csi[k][j][ghost_cell_index].x * csi[k][j][ghost_cell_index].x + 
                                    csi[k][j][ghost_cell_index].y * csi[k][j][ghost_cell_index].y + 
                                    csi[k][j][ghost_cell_index].z * csi[k][j][ghost_cell_index].z
                                );
                                
                                // Compute wall-normal distances using cell Jacobians
                                // sb = distance from wall to first interior cell center
                                // sc = distance from wall to second interior cell center
                                PetscReal distance_to_first_cell = 0.5 / cell_jacobian[k][j][first_interior_cell] / face_area;
                                PetscReal distance_to_second_cell = 2.0 * distance_to_first_cell + 
                                                                     0.5 / cell_jacobian[k][j][second_interior_cell] / face_area;
                                
                                // Compute unit normal vector pointing INTO the domain
                                PetscReal wall_normal[3];
                                wall_normal[0] = csi[k][j][ghost_cell_index].x / face_area;
                                wall_normal[1] = csi[k][j][ghost_cell_index].y / face_area;
                                wall_normal[2] = csi[k][j][ghost_cell_index].z / face_area;
                                
                                // Define velocities for wall function calculation
                                Cmpnts wall_velocity;      // Ua = velocity at wall (zero for stationary wall)
                                Cmpnts reference_velocity; // Uc = velocity at second interior cell
                                
                                wall_velocity.x = wall_velocity.y = wall_velocity.z = 0.0;
                                reference_velocity = velocity_cartesian[k][j][second_interior_cell];
                                
                                // Step 1: Linear interpolation (provides initial guess)
                                noslip(user, distance_to_second_cell, distance_to_first_cell,
                                      wall_velocity, reference_velocity,
                                      &velocity_cartesian[k][j][first_interior_cell],
                                      wall_normal[0], wall_normal[1], wall_normal[2]);
                                
                                // Step 2: Apply log-law correction (improves near-wall velocity)
                                wall_function_loglaw(user, wall_roughness_height,
                                                    distance_to_second_cell, distance_to_first_cell,
                                                    wall_velocity, reference_velocity,
                                                    &velocity_cartesian[k][j][first_interior_cell],
                                                    &friction_velocity[k][j][first_interior_cell],
                                                    wall_normal[0], wall_normal[1], wall_normal[2]);
                                
                                // Ensure ghost cell BC remains zero (required for proper extrapolation)
                                velocity_boundary[k][j][ghost_cell_index].x = 0.0;
                                velocity_boundary[k][j][ghost_cell_index].y = 0.0;
                                velocity_boundary[k][j][ghost_cell_index].z = 0.0;
                                velocity_contravariant[k][j][ghost_cell_index].x = 0.0;
                            }
                        }
                    }
                }
            } break;
            
            // =================================================================
            // POSITIVE X FACE (i = mx-1, normal points in -X direction)
            // =================================================================
            case BC_FACE_POS_X: {
                if (grid_end_i == grid_size_i) {
                    const PetscInt ghost_cell_index = grid_end_i - 1;
                    const PetscInt first_interior_cell = grid_end_i - 2;
                    const PetscInt second_interior_cell = grid_end_i - 3;
                    
                    for (PetscInt k = loop_start_k; k < loop_end_k; k++) {
                        for (PetscInt j = loop_start_j; j < loop_end_j; j++) {
                            
                            if (node_vertex_flag[k][j][first_interior_cell] < 0.1) {
                                
                                PetscReal face_area = sqrt(
                                    csi[k][j][first_interior_cell].x * csi[k][j][first_interior_cell].x + 
                                    csi[k][j][first_interior_cell].y * csi[k][j][first_interior_cell].y + 
                                    csi[k][j][first_interior_cell].z * csi[k][j][first_interior_cell].z
                                );
                                
                                PetscReal distance_to_first_cell = 0.5 / cell_jacobian[k][j][first_interior_cell] / face_area;
                                PetscReal distance_to_second_cell = 2.0 * distance_to_first_cell + 
                                                                     0.5 / cell_jacobian[k][j][second_interior_cell] / face_area;
                                
                                // Note: Normal flipped for +X face to point INTO domain
                                PetscReal wall_normal[3];
                                wall_normal[0] = -csi[k][j][first_interior_cell].x / face_area;
                                wall_normal[1] = -csi[k][j][first_interior_cell].y / face_area;
                                wall_normal[2] = -csi[k][j][first_interior_cell].z / face_area;
                                
                                Cmpnts wall_velocity, reference_velocity;
                                wall_velocity.x = wall_velocity.y = wall_velocity.z = 0.0;
                                reference_velocity = velocity_cartesian[k][j][second_interior_cell];
                                
                                noslip(user, distance_to_second_cell, distance_to_first_cell,
                                      wall_velocity, reference_velocity,
                                      &velocity_cartesian[k][j][first_interior_cell],
                                      wall_normal[0], wall_normal[1], wall_normal[2]);
                                
                                wall_function_loglaw(user, wall_roughness_height,
                                                    distance_to_second_cell, distance_to_first_cell,
                                                    wall_velocity, reference_velocity,
                                                    &velocity_cartesian[k][j][first_interior_cell],
                                                    &friction_velocity[k][j][first_interior_cell],
                                                    wall_normal[0], wall_normal[1], wall_normal[2]);
                                
                                velocity_boundary[k][j][ghost_cell_index].x = 0.0;
                                velocity_boundary[k][j][ghost_cell_index].y = 0.0;
                                velocity_boundary[k][j][ghost_cell_index].z = 0.0;
                                velocity_contravariant[k][j][first_interior_cell].x = 0.0;
                            }
                        }
                    }
                }
            } break;
            
            // =================================================================
            // NEGATIVE Y FACE (j = 0, normal points in +Y direction)
            // =================================================================
            case BC_FACE_NEG_Y: {
                if (grid_start_j == 0) {
                    const PetscInt ghost_cell_index = grid_start_j;
                    const PetscInt first_interior_cell = grid_start_j + 1;
                    const PetscInt second_interior_cell = grid_start_j + 2;
                    
                    for (PetscInt k = loop_start_k; k < loop_end_k; k++) {
                        for (PetscInt i = loop_start_i; i < loop_end_i; i++) {
                            
                            if (node_vertex_flag[k][first_interior_cell][i] < 0.1) {
                                
                                PetscReal face_area = sqrt(
                                    eta[k][ghost_cell_index][i].x * eta[k][ghost_cell_index][i].x + 
                                    eta[k][ghost_cell_index][i].y * eta[k][ghost_cell_index][i].y + 
                                    eta[k][ghost_cell_index][i].z * eta[k][ghost_cell_index][i].z
                                );
                                
                                PetscReal distance_to_first_cell = 0.5 / cell_jacobian[k][first_interior_cell][i] / face_area;
                                PetscReal distance_to_second_cell = 2.0 * distance_to_first_cell + 
                                                                     0.5 / cell_jacobian[k][second_interior_cell][i] / face_area;
                                
                                PetscReal wall_normal[3];
                                wall_normal[0] = eta[k][ghost_cell_index][i].x / face_area;
                                wall_normal[1] = eta[k][ghost_cell_index][i].y / face_area;
                                wall_normal[2] = eta[k][ghost_cell_index][i].z / face_area;
                                
                                Cmpnts wall_velocity, reference_velocity;
                                wall_velocity.x = wall_velocity.y = wall_velocity.z = 0.0;
                                reference_velocity = velocity_cartesian[k][second_interior_cell][i];
                                
                                noslip(user, distance_to_second_cell, distance_to_first_cell,
                                      wall_velocity, reference_velocity,
                                      &velocity_cartesian[k][first_interior_cell][i],
                                      wall_normal[0], wall_normal[1], wall_normal[2]);
                                
                                wall_function_loglaw(user, wall_roughness_height,
                                                    distance_to_second_cell, distance_to_first_cell,
                                                    wall_velocity, reference_velocity,
                                                    &velocity_cartesian[k][first_interior_cell][i],
                                                    &friction_velocity[k][first_interior_cell][i],
                                                    wall_normal[0], wall_normal[1], wall_normal[2]);
                                
                                velocity_boundary[k][ghost_cell_index][i].x = 0.0;
                                velocity_boundary[k][ghost_cell_index][i].y = 0.0;
                                velocity_boundary[k][ghost_cell_index][i].z = 0.0;
                                velocity_contravariant[k][ghost_cell_index][i].y = 0.0;
                            }
                        }
                    }
                }
            } break;
            
            // =================================================================
            // POSITIVE Y FACE (j = my-1, normal points in -Y direction)
            // =================================================================
            case BC_FACE_POS_Y: {
                if (grid_end_j == grid_size_j) {
                    const PetscInt ghost_cell_index = grid_end_j - 1;
                    const PetscInt first_interior_cell = grid_end_j - 2;
                    const PetscInt second_interior_cell = grid_end_j - 3;
                    
                    for (PetscInt k = loop_start_k; k < loop_end_k; k++) {
                        for (PetscInt i = loop_start_i; i < loop_end_i; i++) {
                            
                            if (node_vertex_flag[k][first_interior_cell][i] < 0.1) {
                                
                                PetscReal face_area = sqrt(
                                    eta[k][first_interior_cell][i].x * eta[k][first_interior_cell][i].x + 
                                    eta[k][first_interior_cell][i].y * eta[k][first_interior_cell][i].y + 
                                    eta[k][first_interior_cell][i].z * eta[k][first_interior_cell][i].z
                                );
                                
                                PetscReal distance_to_first_cell = 0.5 / cell_jacobian[k][first_interior_cell][i] / face_area;
                                PetscReal distance_to_second_cell = 2.0 * distance_to_first_cell + 
                                                                     0.5 / cell_jacobian[k][second_interior_cell][i] / face_area;
                                
                                PetscReal wall_normal[3];
                                wall_normal[0] = -eta[k][first_interior_cell][i].x / face_area;
                                wall_normal[1] = -eta[k][first_interior_cell][i].y / face_area;
                                wall_normal[2] = -eta[k][first_interior_cell][i].z / face_area;
                                
                                Cmpnts wall_velocity, reference_velocity;
                                wall_velocity.x = wall_velocity.y = wall_velocity.z = 0.0;
                                reference_velocity = velocity_cartesian[k][second_interior_cell][i];
                                
                                noslip(user, distance_to_second_cell, distance_to_first_cell,
                                      wall_velocity, reference_velocity,
                                      &velocity_cartesian[k][first_interior_cell][i],
                                      wall_normal[0], wall_normal[1], wall_normal[2]);
                                
                                wall_function_loglaw(user, wall_roughness_height,
                                                    distance_to_second_cell, distance_to_first_cell,
                                                    wall_velocity, reference_velocity,
                                                    &velocity_cartesian[k][first_interior_cell][i],
                                                    &friction_velocity[k][first_interior_cell][i],
                                                    wall_normal[0], wall_normal[1], wall_normal[2]);
                                
                                velocity_boundary[k][ghost_cell_index][i].x = 0.0;
                                velocity_boundary[k][ghost_cell_index][i].y = 0.0;
                                velocity_boundary[k][ghost_cell_index][i].z = 0.0;
                                velocity_contravariant[k][first_interior_cell][i].y = 0.0;
                            }
                        }
                    }
                }
            } break;
            
            // =================================================================
            // NEGATIVE Z FACE (k = 0, normal points in +Z direction)
            // =================================================================
            case BC_FACE_NEG_Z: {
                if (grid_start_k == 0) {
                    const PetscInt ghost_cell_index = grid_start_k;
                    const PetscInt first_interior_cell = grid_start_k + 1;
                    const PetscInt second_interior_cell = grid_start_k + 2;
                    
                    for (PetscInt j = loop_start_j; j < loop_end_j; j++) {
                        for (PetscInt i = loop_start_i; i < loop_end_i; i++) {
                            
                            if (node_vertex_flag[first_interior_cell][j][i] < 0.1) {
                                
                                PetscReal face_area = sqrt(
                                    zet[ghost_cell_index][j][i].x * zet[ghost_cell_index][j][i].x + 
                                    zet[ghost_cell_index][j][i].y * zet[ghost_cell_index][j][i].y + 
                                    zet[ghost_cell_index][j][i].z * zet[ghost_cell_index][j][i].z
                                );
                                
                                PetscReal distance_to_first_cell = 0.5 / cell_jacobian[first_interior_cell][j][i] / face_area;
                                PetscReal distance_to_second_cell = 2.0 * distance_to_first_cell + 
                                                                     0.5 / cell_jacobian[second_interior_cell][j][i] / face_area;
                                
                                PetscReal wall_normal[3];
                                wall_normal[0] = zet[ghost_cell_index][j][i].x / face_area;
                                wall_normal[1] = zet[ghost_cell_index][j][i].y / face_area;
                                wall_normal[2] = zet[ghost_cell_index][j][i].z / face_area;
                                
                                Cmpnts wall_velocity, reference_velocity;
                                wall_velocity.x = wall_velocity.y = wall_velocity.z = 0.0;
                                reference_velocity = velocity_cartesian[second_interior_cell][j][i];
                                
                                noslip(user, distance_to_second_cell, distance_to_first_cell,
                                      wall_velocity, reference_velocity,
                                      &velocity_cartesian[first_interior_cell][j][i],
                                      wall_normal[0], wall_normal[1], wall_normal[2]);
                                
                                wall_function_loglaw(user, wall_roughness_height,
                                                    distance_to_second_cell, distance_to_first_cell,
                                                    wall_velocity, reference_velocity,
                                                    &velocity_cartesian[first_interior_cell][j][i],
                                                    &friction_velocity[first_interior_cell][j][i],
                                                    wall_normal[0], wall_normal[1], wall_normal[2]);
                                
                                velocity_boundary[ghost_cell_index][j][i].x = 0.0;
                                velocity_boundary[ghost_cell_index][j][i].y = 0.0;
                                velocity_boundary[ghost_cell_index][j][i].z = 0.0;
                                velocity_contravariant[ghost_cell_index][j][i].z = 0.0;
                            }
                        }
                    }
                }
            } break;
            
            // =================================================================
            // POSITIVE Z FACE (k = mz-1, normal points in -Z direction)
            // =================================================================
            case BC_FACE_POS_Z: {
                if (grid_end_k == grid_size_k) {
                    const PetscInt ghost_cell_index = grid_end_k - 1;
                    const PetscInt first_interior_cell = grid_end_k - 2;
                    const PetscInt second_interior_cell = grid_end_k - 3;
                    
                    for (PetscInt j = loop_start_j; j < loop_end_j; j++) {
                        for (PetscInt i = loop_start_i; i < loop_end_i; i++) {
                            
                            if (node_vertex_flag[first_interior_cell][j][i] < 0.1) {
                                
                                PetscReal face_area = sqrt(
                                    zet[first_interior_cell][j][i].x * zet[first_interior_cell][j][i].x + 
                                    zet[first_interior_cell][j][i].y * zet[first_interior_cell][j][i].y + 
                                    zet[first_interior_cell][j][i].z * zet[first_interior_cell][j][i].z
                                );
                                
                                PetscReal distance_to_first_cell = 0.5 / cell_jacobian[first_interior_cell][j][i] / face_area;
                                PetscReal distance_to_second_cell = 2.0 * distance_to_first_cell + 
                                                                     0.5 / cell_jacobian[second_interior_cell][j][i] / face_area;
                                
                                PetscReal wall_normal[3];
                                wall_normal[0] = -zet[first_interior_cell][j][i].x / face_area;
                                wall_normal[1] = -zet[first_interior_cell][j][i].y / face_area;
                                wall_normal[2] = -zet[first_interior_cell][j][i].z / face_area;
                                
                                Cmpnts wall_velocity, reference_velocity;
                                wall_velocity.x = wall_velocity.y = wall_velocity.z = 0.0;
                                reference_velocity = velocity_cartesian[second_interior_cell][j][i];
                                
                                noslip(user, distance_to_second_cell, distance_to_first_cell,
                                      wall_velocity, reference_velocity,
                                      &velocity_cartesian[first_interior_cell][j][i],
                                      wall_normal[0], wall_normal[1], wall_normal[2]);
                                
                                wall_function_loglaw(user, wall_roughness_height,
                                                    distance_to_second_cell, distance_to_first_cell,
                                                    wall_velocity, reference_velocity,
                                                    &velocity_cartesian[first_interior_cell][j][i],
                                                    &friction_velocity[first_interior_cell][j][i],
                                                    wall_normal[0], wall_normal[1], wall_normal[2]);
                                
                                velocity_boundary[ghost_cell_index][j][i].x = 0.0;
                                velocity_boundary[ghost_cell_index][j][i].y = 0.0;
                                velocity_boundary[ghost_cell_index][j][i].z = 0.0;
                                velocity_contravariant[first_interior_cell][j][i].z = 0.0;
                            }
                        }
                    }
                }
            } break;
        }
    }
    
    // =========================================================================
    // STEP 4: Restore all arrays and release memory
    // =========================================================================
    ierr = DMDAVecRestoreArray(user->fda, user->Ucat, &velocity_cartesian); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(user->fda, user->Ucont, &velocity_contravariant); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(user->fda, user->Bcs.Ubcs, &velocity_boundary); 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***)&node_vertex_flag); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->da, user->lAj, (const PetscReal***)&cell_jacobian); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(user->da, user->lFriction_Velocity, &friction_velocity); CHKERRQ(ierr);
    
    LOG_ALLOW(LOCAL, LOG_DEBUG, "Complete.\n");
    
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

FinalizePostProjectionCellFields()

PetscErrorCode FinalizePostProjectionCellFields

(

UserCtx *

user

)

Implementation of FinalizePostProjectionCellFields().

Finalizes cell-centered fields after the projection step.

Full API contract is documented with the header declaration in include/Boundaries.h.

Definition at line 3257 of file Boundaries.c.

{
    PetscErrorCode ierr;
    const char    *cell_fields[] = {"Ucat", "P"};
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Finalizing post-projection cell-centered fields.\n");
 
    // Ensure flow-dependent Ubcs handlers see the newly reconstructed Ucat.
    ierr = UpdateLocalGhosts(user, "Ucat"); CHKERRQ(ierr);
    ierr = BoundarySystem_RefreshUbcs(user); CHKERRQ(ierr);
 
    // Establish flat non-periodic faces and periodic endpoints before corners.
    ierr = UpdateDummyCells(user); CHKERRQ(ierr);
    ierr = SynchronizePeriodicCellFields(user, 2, cell_fields); CHKERRQ(ierr);
 
    // Corner averaging can overwrite periodic endpoints, so restore them after.
    ierr = UpdateCornerNodes(user); CHKERRQ(ierr);
    ierr = SynchronizePeriodicCellFields(user, 2, cell_fields); CHKERRQ(ierr);
 
    // Synchronize explicitly because the periodic helper is a no-op when every
    // direction is non-periodic.
    ierr = UpdateLocalGhosts(user, "Ucat"); CHKERRQ(ierr);
    ierr = UpdateLocalGhosts(user, "P"); CHKERRQ(ierr);
 
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Post-projection cell-centered fields finalized.\n");
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

ApplyBoundaryConditions()

PetscErrorCode ApplyBoundaryConditions

(

UserCtx *

user

)

Implementation of ApplyBoundaryConditions().

Main boundary-condition orchestrator executed during solver timestepping.

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

See also

ApplyBoundaryConditions()

Definition at line 3297 of file Boundaries.c.

{
    PetscErrorCode ierr;
    const char *staggered_fields[] = {"Ucont"};
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    LOG_ALLOW(GLOBAL,LOG_TRACE,"Boundary Condition Application begins.\n");
 
    // STEP 1: Main iteration loop for applying and converging non-periodic BCs.
    // The number of iterations (e.g., 3) allows information to propagate
    // between coupled boundaries, like an inlet and a conserving outlet.
    for (PetscInt iter = 0; iter < 3; iter++) {
        // (a) Execute the boundary system. This phase calculates fluxes across
        //     the domain and then applies the physical logic for each non-periodic
        //     handler, setting the `ubcs` (boundary value) array.
        ierr = BoundarySystem_ExecuteStep(user); CHKERRQ(ierr);
 
        LOG_ALLOW(GLOBAL,LOG_VERBOSE,"Boundary Condition Setup Executed.\n");
 
        // (b) Synchronize the updated ghost cells across all processors to ensure
        //     all ucont values are current before updating the dummy cells.
        ierr = SynchronizePeriodicStaggeredFields(user, 1, staggered_fields); CHKERRQ(ierr);
 
        // (c) Convert updated Contravariant velocities to Cartesian velocities.
        ierr = Contra2Cart(user); CHKERRQ(ierr);
 
        // (d) Synchronize the updated Cartesian velocities across all processors
        //     to ensure all ucat values are current before updating the dummy cells.
        ierr = UpdateLocalGhosts(user, "Ucat"); CHKERRQ(ierr);
 
        // (e) If Wall functions are enabled, apply them now to adjust near-wall velocities.
        if(user->simCtx->wallfunction){
          // Apply wall function adjustments to the boundary velocities.
          ierr = ApplyWallFunction(user); CHKERRQ(ierr);
 
          // Synchronize the updated Cartesian velocities after wall function adjustments.
          ierr = UpdateLocalGhosts(user, "Ucat"); CHKERRQ(ierr);
 
          LOG_ALLOW(GLOBAL,LOG_VERBOSE,"Wall Function Applied at Walls.\n");
        }
 
        // (f) Update the first layer of ghost cells for non-periodic faces using
        //     the newly computed `ubcs` values.
        ierr = UpdateDummyCells(user); CHKERRQ(ierr);
 
        LOG_ALLOW(GLOBAL,LOG_VERBOSE,"Dummy Cells/Ghost Cells Updated.\n");
 
        // (g) Handle all periodic boundaries. This is a parallel direct copy
        // that sets the absolute constraints for the rest of the solve.
        // There is a Ghost update happening inside this function.
        ierr = ApplyPeriodicBCs(user); CHKERRQ(ierr);
 
        // (h) Update the corner and edge ghost nodes. This routine calculates
        // values for corners/edges by averaging their neighbors, which have been
        // finalized in the steps above (both periodic and non-periodic).
        ierr = UpdateCornerNodes(user); CHKERRQ(ierr);
        
        // (i) Synchronize the updated edge and corner cells across all processors to ensure
        //     consistency before the next iteration or finalization.
        ierr = UpdateLocalGhosts(user, "P"); CHKERRQ(ierr);
        ierr = UpdateLocalGhosts(user, "Ucat"); CHKERRQ(ierr);
        ierr = SynchronizePeriodicStaggeredFields(user, 1, staggered_fields); CHKERRQ(ierr);
 
        // (j) Ensure All the corners are synchronized with  a well defined protocol  in case of Periodic boundary conditions
        // To avoid race conditions.
        const char* all_fields[] = {"Ucat", "P", "Nvert"};
        ierr = SynchronizePeriodicCellFields(user, 3, all_fields); CHKERRQ(ierr);
 
    }
 
    // STEP 3: Final ghost node synchronization. This ensures all changes made
    // to the global vectors are reflected in the local ghost regions of all
    // processors, making the state fully consistent before the next solver stage.
    ierr = UpdateLocalGhosts(user, "P"); CHKERRQ(ierr);
    ierr = UpdateLocalGhosts(user, "Ucat"); CHKERRQ(ierr);
    ierr = SynchronizePeriodicStaggeredFields(user, 1, staggered_fields); CHKERRQ(ierr);
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function: