PICurv 0.1.0
A Parallel Particle-In-Cell Solver for Curvilinear LES
Loading...
Searching...
No Matches
Macros | Functions
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__   "TransferPeriodicField"
 
#define __FUNCT__   "TransferPeriodicFaceField"
 
#define __FUNCT__   "ApplyMetricsPeriodicBCs"
 
#define __FUNCT__   "ApplyPeriodicBCs"
 
#define __FUNCT__   "ApplyUcontPeriodicBCs"
 
#define __FUNCT__   "EnforceUcontPeriodicity"
 
#define __FUNCT__   "UpdateDummyCells"
 
#define __FUNCT__   "UpdateCornerNodes"
 
#define __FUNCT__   "UpdatePeriodicCornerNodes"
 
#define __FUNCT__   "ApplyWallFunction"
 
#define __FUNCT__   "RefreshBoundaryGhostCells"
 
#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 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 ApplyUcontPeriodicBCs (UserCtx *user)
 Implementation of ApplyUcontPeriodicBCs().
 
PetscErrorCode EnforceUcontPeriodicity (UserCtx *user)
 Internal helper implementation: EnforceUcontPeriodicity().
 
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 RefreshBoundaryGhostCells (UserCtx *user)
 Internal helper implementation: RefreshBoundaryGhostCells().
 
PetscErrorCode ApplyBoundaryConditions (UserCtx *user)
 Implementation of ApplyBoundaryConditions().
 

Macro Definition Documentation

◆ __FUNCT__ [1/25]

#define __FUNCT__   "CanRankServiceInletFace"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [2/25]

#define __FUNCT__   "CanRankServiceFace"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [3/25]

#define __FUNCT__   "GetDeterministicFaceGridLocation"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [4/25]

#define __FUNCT__   "GetRandomFCellAndLogicOnInletFace"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [5/25]

#define __FUNCT__   "EnforceRHSBoundaryConditions"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [6/25]

#define __FUNCT__   "BoundaryCondition_Create"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [7/25]

#define __FUNCT__   "BoundarySystem_Validate"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [8/25]

#define __FUNCT__   "BoundarySystem_Initialize"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [9/25]

#define __FUNCT__   "PropagateBoundaryConfigToCoarserLevels"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [10/25]

#define __FUNCT__   "BoundarySystem_ExecuteStep"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [11/25]

#define __FUNCT__   "BoundarySystem_RefreshUbcs"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [12/25]

#define __FUNCT__   "BoundarySystem_Destroy"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [13/25]

#define __FUNCT__   "TransferPeriodicFieldByDirection"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [14/25]

#define __FUNCT__   "TransferPeriodicField"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [15/25]

#define __FUNCT__   "TransferPeriodicFaceField"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [16/25]

#define __FUNCT__   "ApplyMetricsPeriodicBCs"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [17/25]

#define __FUNCT__   "ApplyPeriodicBCs"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [18/25]

#define __FUNCT__   "ApplyUcontPeriodicBCs"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [19/25]

#define __FUNCT__   "EnforceUcontPeriodicity"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [20/25]

#define __FUNCT__   "UpdateDummyCells"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [21/25]

#define __FUNCT__   "UpdateCornerNodes"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [22/25]

#define __FUNCT__   "UpdatePeriodicCornerNodes"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [23/25]

#define __FUNCT__   "ApplyWallFunction"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [24/25]

#define __FUNCT__   "RefreshBoundaryGhostCells"

Definition at line 6 of file Boundaries.c.

◆ __FUNCT__ [25/25]

#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.

14{
15 PetscErrorCode ierr;
16 PetscMPIInt rank_for_logging; // For detailed debugging logs
17 PetscFunctionBeginUser;
18 PROFILE_FUNCTION_BEGIN;
19
20 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank_for_logging); CHKERRQ(ierr);
21
22 *can_service_inlet_out = PETSC_FALSE; // Default to no service
23
24 if (!user->inletFaceDefined) {
25 LOG_ALLOW(LOCAL, LOG_DEBUG, "[Rank %d]: Inlet face not defined in user context. Cannot service.\n", rank_for_logging);
26 PetscFunctionReturn(0);
27 }
28
29 // Get the range of cells owned by this rank in each dimension
30 PetscInt owned_start_cell_i, num_owned_cells_on_rank_i;
31 PetscInt owned_start_cell_j, num_owned_cells_on_rank_j;
32 PetscInt owned_start_cell_k, num_owned_cells_on_rank_k;
33
34 ierr = GetOwnedCellRange(info, 0, &owned_start_cell_i, &num_owned_cells_on_rank_i); CHKERRQ(ierr);
35 ierr = GetOwnedCellRange(info, 1, &owned_start_cell_j, &num_owned_cells_on_rank_j); CHKERRQ(ierr);
36 ierr = GetOwnedCellRange(info, 2, &owned_start_cell_k, &num_owned_cells_on_rank_k); CHKERRQ(ierr);
37
38 // Determine the global index of the last cell (0-indexed) in each direction.
39 // Example: If IM_nodes_global = 11 (nodes 0-10), there are 10 cells (0-9). Last cell index is 9.
40 // Formula: global_nodes - 1 (num cells) - 1 (0-indexed) = global_nodes - 2.
41 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)
42 PetscInt last_global_cell_idx_j = (JM_nodes_global > 1) ? (JM_nodes_global - 2) : -1;
43 PetscInt last_global_cell_idx_k = (KM_nodes_global > 1) ? (KM_nodes_global - 2) : -1;
44
45 switch (user->identifiedInletBCFace) {
46 case BC_FACE_NEG_X: // Inlet on the global I-minimum face (face of cell C_i=0)
47 // Rank services if its first owned node is global node 0 (info->xs == 0),
48 // and it owns cells in I, J, and K directions.
49 if (info->xs == 0 && num_owned_cells_on_rank_i > 0 &&
50 num_owned_cells_on_rank_j > 0 && num_owned_cells_on_rank_k > 0) {
51 *can_service_inlet_out = PETSC_TRUE;
52 }
53 break;
54 case BC_FACE_POS_X: // Inlet on the global I-maximum face (face of cell C_i=last_global_cell_idx_i)
55 // Rank services if it owns the last cell in I-direction,
56 // and has extent in J and K.
57 if (last_global_cell_idx_i >= 0 && /* Check for valid global domain */
58 (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 */
59 num_owned_cells_on_rank_j > 0 && num_owned_cells_on_rank_k > 0) {
60 *can_service_inlet_out = PETSC_TRUE;
61 }
62 break;
63 case BC_FACE_NEG_Y:
64 if (info->ys == 0 && num_owned_cells_on_rank_j > 0 &&
65 num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_k > 0) {
66 *can_service_inlet_out = PETSC_TRUE;
67 }
68 break;
69 case BC_FACE_POS_Y:
70 if (last_global_cell_idx_j >= 0 &&
71 (owned_start_cell_j + num_owned_cells_on_rank_j - 1) == last_global_cell_idx_j &&
72 num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_k > 0) {
73 *can_service_inlet_out = PETSC_TRUE;
74 }
75 break;
76 case BC_FACE_NEG_Z:
77 if (info->zs == 0 && num_owned_cells_on_rank_k > 0 &&
78 num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_j > 0) {
79 *can_service_inlet_out = PETSC_TRUE;
80 }
81 break;
82 case BC_FACE_POS_Z:
83 if (last_global_cell_idx_k >= 0 &&
84 (owned_start_cell_k + num_owned_cells_on_rank_k - 1) == last_global_cell_idx_k &&
85 num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_j > 0) {
86 *can_service_inlet_out = PETSC_TRUE;
87 }
88 break;
89 default:
90 LOG_ALLOW(LOCAL, LOG_WARNING, "[Rank %d]: Unknown inlet face %s.\n", rank_for_logging, BCFaceToString((BCFace)user->identifiedInletBCFace));
91 break;
92 }
93
94 LOG_ALLOW(LOCAL, LOG_TRACE,
95 "[Rank %d] Check Service for Inlet %s:\n"
96 " - Local Domain: starts at cell (%d,%d,%d), has (%d,%d,%d) cells.\n"
97 " - Global Domain: has (%d,%d,%d) nodes, so last cell is (%d,%d,%d).\n",
98 rank_for_logging,
99 BCFaceToString((BCFace)user->identifiedInletBCFace),
100 owned_start_cell_i, owned_start_cell_j, owned_start_cell_k,
101 num_owned_cells_on_rank_i, num_owned_cells_on_rank_j, num_owned_cells_on_rank_k,
102 IM_nodes_global, JM_nodes_global, KM_nodes_global,
103 last_global_cell_idx_i, last_global_cell_idx_j, last_global_cell_idx_k);
104
105 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",
106 rank_for_logging,
107 BCFaceToString((BCFace)user->identifiedInletBCFace),
108 (*can_service_inlet_out) ? "CAN SERVICE" : "CANNOT SERVICE",
109 num_owned_cells_on_rank_i, num_owned_cells_on_rank_j, num_owned_cells_on_rank_k,
110 owned_start_cell_i, owned_start_cell_j, owned_start_cell_k);
111
112 PROFILE_FUNCTION_END;
113
114 PetscFunctionReturn(0);
115}
LOCAL
#define LOCAL
Logging scope definitions for controlling message output.
Definition logging.h:44
BCFaceToString
const char * BCFaceToString(BCFace face)
Helper function to convert BCFace enum to a string representation.
Definition logging.c:669
LOG_ALLOW
#define LOG_ALLOW(scope, level, fmt,...)
Logging macro that checks both the log level and whether the calling function is in the allowed-funct...
Definition logging.h:199
PROFILE_FUNCTION_END
#define PROFILE_FUNCTION_END
Marks the end of a profiled code block.
Definition logging.h:779
LOG_TRACE
@ LOG_TRACE
Very fine-grained tracing information for in-depth debugging.
Definition logging.h:32
LOG_INFO
@ LOG_INFO
Informational messages about program execution.
Definition logging.h:30
LOG_WARNING
@ LOG_WARNING
Non-critical issues that warrant attention.
Definition logging.h:29
LOG_DEBUG
@ LOG_DEBUG
Detailed debugging information.
Definition logging.h:31
PROFILE_FUNCTION_BEGIN
#define PROFILE_FUNCTION_BEGIN
Marks the beginning of a profiled code block (typically a function).
Definition logging.h:770
GetOwnedCellRange
PetscErrorCode GetOwnedCellRange(const DMDALocalInfo *info_nodes, PetscInt dim, PetscInt *xs_cell_global_out, PetscInt *xm_cell_local_out)
Determines the global starting index and number of CELLS owned by the current processor in a specifie...
Definition setup.c:1883
UserCtx::inletFaceDefined
PetscBool inletFaceDefined
Definition variables.h:830
UserCtx::identifiedInletBCFace
BCFace identifiedInletBCFace
Definition variables.h:831
BCFace
BCFace
Identifies the six logical faces of a structured computational block.
Definition variables.h:244
BC_FACE_NEG_X
@ BC_FACE_NEG_X
Definition variables.h:245
BC_FACE_POS_Z
@ BC_FACE_POS_Z
Definition variables.h:247
BC_FACE_POS_Y
@ BC_FACE_POS_Y
Definition variables.h:246
BC_FACE_NEG_Z
@ BC_FACE_NEG_Z
Definition variables.h:247
BC_FACE_POS_X
@ BC_FACE_POS_X
Definition variables.h:245
BC_FACE_NEG_Y
@ BC_FACE_NEG_Y
Definition variables.h:246
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.

128{
129 PetscErrorCode ierr;
130 PetscMPIInt rank_for_logging;
131 PetscFunctionBeginUser;
132
133 PROFILE_FUNCTION_BEGIN;
134
135 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank_for_logging); CHKERRQ(ierr);
136
137 *can_service_out = PETSC_FALSE; // Default to no service
138
139 // Get the range of cells owned by this rank
140 PetscInt owned_start_cell_i, num_owned_cells_on_rank_i;
141 PetscInt owned_start_cell_j, num_owned_cells_on_rank_j;
142 PetscInt owned_start_cell_k, num_owned_cells_on_rank_k;
143 ierr = GetOwnedCellRange(info, 0, &owned_start_cell_i, &num_owned_cells_on_rank_i); CHKERRQ(ierr);
144 ierr = GetOwnedCellRange(info, 1, &owned_start_cell_j, &num_owned_cells_on_rank_j); CHKERRQ(ierr);
145 ierr = GetOwnedCellRange(info, 2, &owned_start_cell_k, &num_owned_cells_on_rank_k); CHKERRQ(ierr);
146
147 // Determine the global index of the last cell (0-indexed) in each direction.
148 PetscInt last_global_cell_idx_i = (IM_nodes_global > 1) ? (IM_nodes_global - 2) : -1;
149 PetscInt last_global_cell_idx_j = (JM_nodes_global > 1) ? (JM_nodes_global - 2) : -1;
150 PetscInt last_global_cell_idx_k = (KM_nodes_global > 1) ? (KM_nodes_global - 2) : -1;
151
152 switch (face_id) {
153 case BC_FACE_NEG_X:
154 if (info->xs == 0 && num_owned_cells_on_rank_i > 0 &&
155 num_owned_cells_on_rank_j > 0 && num_owned_cells_on_rank_k > 0) {
156 *can_service_out = PETSC_TRUE;
157 }
158 break;
159 case BC_FACE_POS_X:
160 if (last_global_cell_idx_i >= 0 &&
161 (owned_start_cell_i + num_owned_cells_on_rank_i - 1) == last_global_cell_idx_i &&
162 num_owned_cells_on_rank_j > 0 && num_owned_cells_on_rank_k > 0) {
163 *can_service_out = PETSC_TRUE;
164 }
165 break;
166 case BC_FACE_NEG_Y:
167 if (info->ys == 0 && num_owned_cells_on_rank_j > 0 &&
168 num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_k > 0) {
169 *can_service_out = PETSC_TRUE;
170 }
171 break;
172 case BC_FACE_POS_Y:
173 if (last_global_cell_idx_j >= 0 &&
174 (owned_start_cell_j + num_owned_cells_on_rank_j - 1) == last_global_cell_idx_j &&
175 num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_k > 0) {
176 *can_service_out = PETSC_TRUE;
177 }
178 break;
179 case BC_FACE_NEG_Z:
180 if (info->zs == 0 && num_owned_cells_on_rank_k > 0 &&
181 num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_j > 0) {
182 *can_service_out = PETSC_TRUE;
183 }
184 break;
185 case BC_FACE_POS_Z:
186 if (last_global_cell_idx_k >= 0 &&
187 (owned_start_cell_k + num_owned_cells_on_rank_k - 1) == last_global_cell_idx_k &&
188 num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_j > 0) {
189 *can_service_out = PETSC_TRUE;
190 }
191 break;
192 default:
193 LOG_ALLOW(LOCAL, LOG_WARNING, "Rank %d: Unknown face enum %d. \n", rank_for_logging, face_id);
194 break;
195 }
196
197 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d check for face %s: Result=%s. \n",
198 rank_for_logging, BCFaceToString((BCFace)face_id), (*can_service_out ? "TRUE" : "FALSE"));
199
200 PROFILE_FUNCTION_END;
201
202 PetscFunctionReturn(0);
203}
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.

220{
221 SimCtx *simCtx = user->simCtx;
222 PetscReal global_logic_i = 0.0, global_logic_j = 0.0, global_logic_k = 0.0;
223 PetscErrorCode ierr;
224 PetscMPIInt rank_for_logging;
225
226 PetscFunctionBeginUser;
227 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank_for_logging); CHKERRQ(ierr);
228
229 *placement_successful_out = PETSC_FALSE; // Default to failure
230
231 // --- Step 1: Configuration and Input Validation ---
232
233 // *** Hardcoded number of grid layers. Change this value to alter the pattern. ***
234 const PetscInt grid_layers = 2;
235
236 LOG_ALLOW(LOCAL, LOG_DEBUG,
237 "[Rank %d] Placing particle %lld on face %s with grid_layers=%d in global domain (%d,%d,%d) cells.\n",
238 rank_for_logging, (long long)particle_global_id, BCFaceToString(user->identifiedInletBCFace), grid_layers,
239 IM_cells_global, JM_cells_global, KM_cells_global);
240
241 const char *face_name = BCFaceToString(user->identifiedInletBCFace);
242
243 // Fatal Error Checks: Ensure the requested grid is geometrically possible.
244 // The total layers from opposite faces (2 * grid_layers) must be less than the domain size.
245 switch (user->identifiedInletBCFace) {
246 case BC_FACE_NEG_X: case BC_FACE_POS_X:
247 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);
248 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);
249 break;
250 case BC_FACE_NEG_Y: case BC_FACE_POS_Y:
251 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);
252 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);
253 break;
254 case BC_FACE_NEG_Z: case BC_FACE_POS_Z:
255 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);
256 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);
257 break;
258 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid identifiedInletBCFace specified: %d", user->identifiedInletBCFace);
259 }
260
261 const PetscInt num_lines_total = 4 * grid_layers;
262 if (simCtx->np < num_lines_total) {
263 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);
264 }
265 if (simCtx->np > 0 && simCtx->np % num_lines_total != 0) {
266 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);
267 }
268
269 // --- Step 2: Map global particle ID to a line and a point on that line ---
270 if (simCtx->np == 0) PetscFunctionReturn(0); // Nothing to do
271
272 LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Distributing %lld particles over %d lines on face %s.\n",
273 rank_for_logging, (long long)simCtx->np, num_lines_total, face_name);
274
275 const PetscInt points_per_line = PetscMax(1, simCtx->np / num_lines_total);
276 PetscInt line_index = particle_global_id / points_per_line;
277 PetscInt point_index_on_line = particle_global_id % points_per_line;
278 line_index = PetscMin(line_index, num_lines_total - 1); // Clamp to handle uneven division
279
280 // Decode the line_index into an edge group (0-3) and a layer within that group (0 to grid_layers-1)
281 const PetscInt edge_group = line_index / grid_layers;
282 const PetscInt layer_index = line_index % grid_layers;
283
284 // --- Step 3: Calculate placement coordinates based on the decoded indices ---
285 const PetscReal layer_spacing_norm_i = (IM_cells_global > 0) ? 1.0 / (PetscReal)IM_cells_global : 0.0;
286 const PetscReal layer_spacing_norm_j = (JM_cells_global > 0) ? 1.0 / (PetscReal)JM_cells_global : 0.0;
287 const PetscReal layer_spacing_norm_k = (KM_cells_global > 0) ? 1.0 / (PetscReal)KM_cells_global : 0.0;
288
289 // Grid-aware epsilon: scale with minimum cell size to keep particles away from rank boundaries
290 const PetscReal min_layer_spacing = PetscMin(layer_spacing_norm_i, PetscMin(layer_spacing_norm_j, layer_spacing_norm_k));
291 const PetscReal epsilon = 0.5 * min_layer_spacing; // Keep particles 10% of cell width from boundaries
292
293 PetscReal variable_coord; // The coordinate that varies along a line
294 if (points_per_line <= 1) {
295 variable_coord = 0.5; // Place single point in the middle
296 } else {
297 variable_coord = ((PetscReal)point_index_on_line + 0.5)/ (PetscReal)(points_per_line);
298 }
299 variable_coord = PetscMin(1.0 - epsilon, PetscMax(epsilon, variable_coord)); // Clamp within [eps, 1-eps]
300
301 // Main logic switch to determine the three global logical coordinates
302 switch (user->identifiedInletBCFace) {
303 case BC_FACE_NEG_X:
304 global_logic_i = 0.5 * layer_spacing_norm_i; // Place near the face, in the middle of the first cell
305 if (edge_group == 0) { global_logic_j = (PetscReal)layer_index * layer_spacing_norm_j + epsilon; global_logic_k = variable_coord; }
306 else if (edge_group == 1) { global_logic_j = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_j) - epsilon; global_logic_k = variable_coord; }
307 else if (edge_group == 2) { global_logic_k = (PetscReal)layer_index * layer_spacing_norm_k + epsilon; global_logic_j = variable_coord; }
308 else /* edge_group == 3 */ { global_logic_k = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_k) - epsilon; global_logic_j = variable_coord; }
309 break;
310 case BC_FACE_POS_X:
311 global_logic_i = 1.0 - (0.5 * layer_spacing_norm_i); // Place near the face, in the middle of the last cell
312 if (edge_group == 0) { global_logic_j = (PetscReal)layer_index * layer_spacing_norm_j + epsilon; global_logic_k = variable_coord; }
313 else if (edge_group == 1) { global_logic_j = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_j) - epsilon; global_logic_k = variable_coord; }
314 else if (edge_group == 2) { global_logic_k = (PetscReal)layer_index * layer_spacing_norm_k + epsilon; global_logic_j = variable_coord; }
315 else /* edge_group == 3 */ { global_logic_k = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_k) - epsilon; global_logic_j = variable_coord; }
316 break;
317 case BC_FACE_NEG_Y:
318 global_logic_j = 0.5 * layer_spacing_norm_j;
319 if (edge_group == 0) { global_logic_i = (PetscReal)layer_index * layer_spacing_norm_i + epsilon; global_logic_k = variable_coord; }
320 else if (edge_group == 1) { global_logic_i = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_i) - epsilon; global_logic_k = variable_coord; }
321 else if (edge_group == 2) { global_logic_k = (PetscReal)layer_index * layer_spacing_norm_k + epsilon; global_logic_i = variable_coord; }
322 else /* edge_group == 3 */ { global_logic_k = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_k) - epsilon; global_logic_i = variable_coord; }
323 break;
324 case BC_FACE_POS_Y:
325 global_logic_j = 1.0 - (0.5 * layer_spacing_norm_j);
326 if (edge_group == 0) { global_logic_i = (PetscReal)layer_index * layer_spacing_norm_i + epsilon; global_logic_k = variable_coord; }
327 else if (edge_group == 1) { global_logic_i = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_i) - epsilon; global_logic_k = variable_coord; }
328 else if (edge_group == 2) { global_logic_k = (PetscReal)layer_index * layer_spacing_norm_k + epsilon; global_logic_i = variable_coord; }
329 else /* edge_group == 3 */ { global_logic_k = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_k) - epsilon; global_logic_i = variable_coord; }
330 break;
331 case BC_FACE_NEG_Z:
332 global_logic_k = 0.5 * layer_spacing_norm_k;
333 if (edge_group == 0) { global_logic_i = (PetscReal)layer_index * layer_spacing_norm_i + epsilon; global_logic_j = variable_coord; }
334 else if (edge_group == 1) { global_logic_i = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_i) - epsilon; global_logic_j = variable_coord; }
335 else if (edge_group == 2) { global_logic_j = (PetscReal)layer_index * layer_spacing_norm_j + epsilon; global_logic_i = variable_coord; }
336 else /* edge_group == 3 */ { global_logic_j = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_j) - epsilon; global_logic_i = variable_coord; }
337 break;
338 case BC_FACE_POS_Z:
339 global_logic_k = 1.0 - (0.5 * layer_spacing_norm_k);
340 if (edge_group == 0) { global_logic_i = (PetscReal)layer_index * layer_spacing_norm_i + epsilon; global_logic_j = variable_coord; }
341 else if (edge_group == 1) { global_logic_i = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_i) - epsilon; global_logic_j = variable_coord; }
342 else if (edge_group == 2) { global_logic_j = (PetscReal)layer_index * layer_spacing_norm_j + epsilon; global_logic_i = variable_coord; }
343 else /* edge_group == 3 */ { global_logic_j = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_j) - epsilon; global_logic_i = variable_coord; }
344 break;
345 }
346
347 LOG_ALLOW(LOCAL, LOG_TRACE,
348 "[Rank %d] Particle %lld assigned to line %d (edge group %d, layer %d) with variable_coord=%.4f.\n"
349 " -> Global logical coords: (i,j,k) = (%.6f, %.6f, %.6f)\n",
350 rank_for_logging, (long long)particle_global_id, line_index, edge_group, layer_index, variable_coord,
351 global_logic_i, global_logic_j, global_logic_k);
352
353 // --- Step 4: Convert global logical coordinate to global cell index and intra-cell logicals ---
354 PetscReal global_cell_coord_i = global_logic_i * IM_cells_global;
355 PetscInt I_g = (PetscInt)global_cell_coord_i;
356 *xi_metric_logic_out = global_cell_coord_i - I_g;
357
358 PetscReal global_cell_coord_j = global_logic_j * JM_cells_global;
359 PetscInt J_g = (PetscInt)global_cell_coord_j;
360 *eta_metric_logic_out = global_cell_coord_j - J_g;
361
362 PetscReal global_cell_coord_k = global_logic_k * KM_cells_global;
363 PetscInt K_g = (PetscInt)global_cell_coord_k;
364 *zta_metric_logic_out = global_cell_coord_k - K_g;
365
366 // --- Step 5: Check if this rank owns the target cell and finalize outputs ---
367 if ((I_g >= info->xs && I_g < info->xs + info->xm) &&
368 (J_g >= info->ys && J_g < info->ys + info->ym) &&
369 (K_g >= info->zs && K_g < info->zs + info->zm))
370 {
371 // Convert global cell index to the local node index for this rank's DA patch
372 *ci_metric_lnode_out = (I_g - info->xs) + xs_gnode_rank;
373 *cj_metric_lnode_out = (J_g - info->ys) + ys_gnode_rank;
374 *ck_metric_lnode_out = (K_g - info->zs) + zs_gnode_rank;
375 *placement_successful_out = PETSC_TRUE;
376 }
377
378 LOG_ALLOW(LOCAL, LOG_VERBOSE,
379 "[Rank %d] Particle %lld placement %s.\n",
380 rank_for_logging, (long long)particle_global_id,
381 (*placement_successful_out ? "SUCCESSFUL" : "NOT ON THIS RANK"));
382
383 if(*placement_successful_out){
384 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",
385 *ci_metric_lnode_out, *cj_metric_lnode_out, *ck_metric_lnode_out,
386 *xi_metric_logic_out, *eta_metric_logic_out, *zta_metric_logic_out);
387 }
388
389 PetscFunctionReturn(0);
390}
GLOBAL
#define GLOBAL
Scope for global logging across all processes.
Definition logging.h:45
LOG_VERBOSE
@ LOG_VERBOSE
Extremely detailed logs, typically for development use only.
Definition logging.h:33
UserCtx::simCtx
SimCtx * simCtx
Back-pointer to the master simulation context.
Definition variables.h:814
SimCtx::np
PetscInt np
Definition variables.h:739
SimCtx
The master context for the entire simulation.
Definition variables.h:643
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.

406{
407 PetscErrorCode ierr = 0;
408 PetscReal r_val_i_sel, r_val_j_sel, r_val_k_sel;
409 PetscInt local_cell_idx_on_face_dim1 = 0; // 0-indexed relative to owned cells on face
410 PetscInt local_cell_idx_on_face_dim2 = 0;
411 PetscMPIInt rank_for_logging;
412
413 PetscFunctionBeginUser;
414
415 PROFILE_FUNCTION_BEGIN;
416
417 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank_for_logging); CHKERRQ(ierr);
418
419 // Get number of cells this rank owns in each dimension (tangential to the face mainly)
420 PetscInt owned_start_cell_i, num_owned_cells_on_rank_i;
421 PetscInt owned_start_cell_j, num_owned_cells_on_rank_j;
422 PetscInt owned_start_cell_k, num_owned_cells_on_rank_k;
423
424 ierr = GetOwnedCellRange(info, 0, &owned_start_cell_i, &num_owned_cells_on_rank_i); CHKERRQ(ierr);
425 ierr = GetOwnedCellRange(info, 1, &owned_start_cell_j, &num_owned_cells_on_rank_j); CHKERRQ(ierr);
426 ierr = GetOwnedCellRange(info, 2, &owned_start_cell_k, &num_owned_cells_on_rank_k); CHKERRQ(ierr);
427
428 // Defaults for cell origin node (local index for the rank's DA patch, including ghosts)
429 *ci_metric_lnode_out = xs_gnode_rank; *cj_metric_lnode_out = ys_gnode_rank; *ck_metric_lnode_out = zs_gnode_rank;
430 // Defaults for logical coordinates
431 *xi_metric_logic_out = 0.5; *eta_metric_logic_out = 0.5; *zta_metric_logic_out = 0.5;
432
433 // Index of the last cell (0-indexed) in each global direction
434 PetscInt last_global_cell_idx_i = (IM_nodes_global > 1) ? (IM_nodes_global - 2) : -1;
435 PetscInt last_global_cell_idx_j = (JM_nodes_global > 1) ? (JM_nodes_global - 2) : -1;
436 PetscInt last_global_cell_idx_k = (KM_nodes_global > 1) ? (KM_nodes_global - 2) : -1;
437
438 LOG_ALLOW(LOCAL, LOG_INFO, "PARTICLE_INIT_DEBUG Rank %d: Inlet face %s.\n"
439 " Owned cells (i,j,k): (%d,%d,%d)\n"
440 " Global nodes (I,J,K): (%d,%d,%d)\n"
441 " info->xs,ys,zs (first owned node GLOBAL): (%d,%d,%d)\n"
442 " info->xm,ym,zm (num owned nodes GLOBAL): (%d,%d,%d)\n"
443 " xs_gnode_rank,ys_gnode_rank,zs_gnode_rank (DMDAGetCorners): (%d,%d,%d)\n"
444 " owned_start_cell (i,j,k) GLOBAL: (%d,%d,%d)\n"
445 " last_global_cell_idx (i,j,k): (%d,%d,%d)\n",
446 rank_for_logging, BCFaceToString((BCFace)user->identifiedInletBCFace),
447 num_owned_cells_on_rank_i,num_owned_cells_on_rank_j,num_owned_cells_on_rank_k,
448 IM_nodes_global,JM_nodes_global,KM_nodes_global,
449 info->xs, info->ys, info->zs,
450 info->xm, info->ym, info->zm,
451 xs_gnode_rank,ys_gnode_rank,zs_gnode_rank,
452 owned_start_cell_i, owned_start_cell_j, owned_start_cell_k,
453 last_global_cell_idx_i, last_global_cell_idx_j, last_global_cell_idx_k);
454
455
456 switch (user->identifiedInletBCFace) {
457 case BC_FACE_NEG_X: // Particle on -X face of cell C_0 (origin node N_0)
458 // Cell origin node is the first owned node in I by this rank (global index info->xs).
459 // Its local index within the rank's DA (incl ghosts) is xs_gnode_rank.
460 *ci_metric_lnode_out = xs_gnode_rank;
461 *xi_metric_logic_out = 1.0e-6;
462
463 // Tangential dimensions are J and K. Select an owned cell randomly on this face.
464 // num_owned_cells_on_rank_j/k must be > 0 (checked by CanRankServiceInletFace)
465 ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, &r_val_j_sel); CHKERRQ(ierr);
466 local_cell_idx_on_face_dim1 = (PetscInt)(r_val_j_sel * num_owned_cells_on_rank_j); // Index among owned J-cells
467 local_cell_idx_on_face_dim1 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim1), num_owned_cells_on_rank_j - 1);
468 *cj_metric_lnode_out = ys_gnode_rank + local_cell_idx_on_face_dim1; // Offset from start of rank's J-nodes
469
470 ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, &r_val_k_sel); CHKERRQ(ierr);
471 local_cell_idx_on_face_dim2 = (PetscInt)(r_val_k_sel * num_owned_cells_on_rank_k);
472 local_cell_idx_on_face_dim2 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim2), num_owned_cells_on_rank_k - 1);
473 *ck_metric_lnode_out = zs_gnode_rank + local_cell_idx_on_face_dim2;
474
475 ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, eta_metric_logic_out); CHKERRQ(ierr);
476 ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, zta_metric_logic_out); CHKERRQ(ierr);
477 break;
478
479 case BC_FACE_POS_X: // Particle on +X face of cell C_last_I (origin node N_last_I_origin)
480 // Origin node of the last I-cell is global_node_idx = last_global_cell_idx_i.
481 // Its local index in rank's DA: (last_global_cell_idx_i - info->xs) + xs_gnode_rank
482 *ci_metric_lnode_out = xs_gnode_rank + (last_global_cell_idx_i - info->xs);
483 *xi_metric_logic_out = 1.0 - 1.0e-6;
484
485 ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, &r_val_j_sel); CHKERRQ(ierr);
486 local_cell_idx_on_face_dim1 = (PetscInt)(r_val_j_sel * num_owned_cells_on_rank_j);
487 local_cell_idx_on_face_dim1 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim1), num_owned_cells_on_rank_j - 1);
488 *cj_metric_lnode_out = ys_gnode_rank + local_cell_idx_on_face_dim1;
489
490 ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, &r_val_k_sel); CHKERRQ(ierr);
491 local_cell_idx_on_face_dim2 = (PetscInt)(r_val_k_sel * num_owned_cells_on_rank_k);
492 local_cell_idx_on_face_dim2 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim2), num_owned_cells_on_rank_k - 1);
493 *ck_metric_lnode_out = zs_gnode_rank + local_cell_idx_on_face_dim2;
494
495 ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, eta_metric_logic_out); CHKERRQ(ierr);
496 ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, zta_metric_logic_out); CHKERRQ(ierr);
497 break;
498 // ... (Cases for Y and Z faces, following the same pattern) ...
499 case BC_FACE_NEG_Y:
500 *cj_metric_lnode_out = ys_gnode_rank;
501 *eta_metric_logic_out = 1.0e-6;
502 ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, &r_val_i_sel); CHKERRQ(ierr);
503 local_cell_idx_on_face_dim1 = (PetscInt)(r_val_i_sel * num_owned_cells_on_rank_i);
504 local_cell_idx_on_face_dim1 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim1), num_owned_cells_on_rank_i - 1);
505 *ci_metric_lnode_out = xs_gnode_rank + local_cell_idx_on_face_dim1;
506 ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, &r_val_k_sel); CHKERRQ(ierr);
507 local_cell_idx_on_face_dim2 = (PetscInt)(r_val_k_sel * num_owned_cells_on_rank_k);
508 local_cell_idx_on_face_dim2 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim2), num_owned_cells_on_rank_k - 1);
509 *ck_metric_lnode_out = zs_gnode_rank + local_cell_idx_on_face_dim2;
510 ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, xi_metric_logic_out); CHKERRQ(ierr);
511 ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, zta_metric_logic_out); CHKERRQ(ierr);
512 break;
513 case BC_FACE_POS_Y:
514 *cj_metric_lnode_out = ys_gnode_rank + (last_global_cell_idx_j - info->ys);
515 *eta_metric_logic_out = 1.0 - 1.0e-6;
516 ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, &r_val_i_sel); CHKERRQ(ierr);
517 local_cell_idx_on_face_dim1 = (PetscInt)(r_val_i_sel * num_owned_cells_on_rank_i);
518 local_cell_idx_on_face_dim1 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim1), num_owned_cells_on_rank_i - 1);
519 *ci_metric_lnode_out = xs_gnode_rank + local_cell_idx_on_face_dim1;
520 ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, &r_val_k_sel); CHKERRQ(ierr);
521 local_cell_idx_on_face_dim2 = (PetscInt)(r_val_k_sel * num_owned_cells_on_rank_k);
522 local_cell_idx_on_face_dim2 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim2), num_owned_cells_on_rank_k - 1);
523 *ck_metric_lnode_out = zs_gnode_rank + local_cell_idx_on_face_dim2;
524 ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, xi_metric_logic_out); CHKERRQ(ierr);
525 ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, zta_metric_logic_out); CHKERRQ(ierr);
526 break;
527 case BC_FACE_NEG_Z: // Your example case
528 *ck_metric_lnode_out = zs_gnode_rank; // Cell origin is the first owned node in K by this rank
529 *zta_metric_logic_out = 1.0e-6; // Place particle slightly inside this cell from its -Z face
530 // Tangential dimensions are I and J
531 ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, &r_val_i_sel); CHKERRQ(ierr);
532 local_cell_idx_on_face_dim1 = (PetscInt)(r_val_i_sel * num_owned_cells_on_rank_i);
533 local_cell_idx_on_face_dim1 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim1), num_owned_cells_on_rank_i - 1);
534 *ci_metric_lnode_out = xs_gnode_rank + local_cell_idx_on_face_dim1;
535
536 ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, &r_val_j_sel); CHKERRQ(ierr);
537 local_cell_idx_on_face_dim2 = (PetscInt)(r_val_j_sel * num_owned_cells_on_rank_j);
538 local_cell_idx_on_face_dim2 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim2), num_owned_cells_on_rank_j - 1);
539 *cj_metric_lnode_out = ys_gnode_rank + local_cell_idx_on_face_dim2;
540
541 ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, xi_metric_logic_out); CHKERRQ(ierr); // Intra-cell logical for I
542 ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, eta_metric_logic_out); CHKERRQ(ierr); // Intra-cell logical for J
543 break;
544 case BC_FACE_POS_Z:
545 *ck_metric_lnode_out = zs_gnode_rank + (last_global_cell_idx_k - info->zs);
546 *zta_metric_logic_out = 1.0 - 1.0e-6;
547 ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, &r_val_i_sel); CHKERRQ(ierr);
548 local_cell_idx_on_face_dim1 = (PetscInt)(r_val_i_sel * num_owned_cells_on_rank_i);
549 local_cell_idx_on_face_dim1 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim1), num_owned_cells_on_rank_i - 1);
550 *ci_metric_lnode_out = xs_gnode_rank + local_cell_idx_on_face_dim1;
551 ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, &r_val_j_sel); CHKERRQ(ierr);
552 local_cell_idx_on_face_dim2 = (PetscInt)(r_val_j_sel * num_owned_cells_on_rank_j);
553 local_cell_idx_on_face_dim2 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim2), num_owned_cells_on_rank_j - 1);
554 *cj_metric_lnode_out = ys_gnode_rank + local_cell_idx_on_face_dim2;
555 ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, xi_metric_logic_out); CHKERRQ(ierr);
556 ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, eta_metric_logic_out); CHKERRQ(ierr);
557 break;
558 default:
559 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "GetRandomCellAndLogicOnInletFace: Invalid user->identifiedInletBCFace %d. \n", user->identifiedInletBCFace);
560 }
561
562 PetscReal eps = 1.0e-7;
563 if (user->identifiedInletBCFace == BC_FACE_NEG_X || user->identifiedInletBCFace == BC_FACE_POS_X) {
564 *eta_metric_logic_out = PetscMin(PetscMax(0.0, *eta_metric_logic_out), 1.0 - eps);
565 *zta_metric_logic_out = PetscMin(PetscMax(0.0, *zta_metric_logic_out), 1.0 - eps);
566 } else if (user->identifiedInletBCFace == BC_FACE_NEG_Y || user->identifiedInletBCFace == BC_FACE_POS_Y) {
567 *xi_metric_logic_out = PetscMin(PetscMax(0.0, *xi_metric_logic_out), 1.0 - eps);
568 *zta_metric_logic_out = PetscMin(PetscMax(0.0, *zta_metric_logic_out), 1.0 - eps);
569 } else {
570 *xi_metric_logic_out = PetscMin(PetscMax(0.0, *xi_metric_logic_out), 1.0 - eps);
571 *eta_metric_logic_out = PetscMin(PetscMax(0.0, *eta_metric_logic_out), 1.0 - eps);
572 }
573
574 LOG_ALLOW(LOCAL, LOG_VERBOSE, "Rank %d: Target Cell Node =(%d,%d,%d). (xi,et,zt)=(%.2e,%.2f,%.2f). \n",
575 rank_for_logging, *ci_metric_lnode_out, *cj_metric_lnode_out, *ck_metric_lnode_out,
576 *xi_metric_logic_out, *eta_metric_logic_out, *zta_metric_logic_out);
577
578 PROFILE_FUNCTION_END;
579
580 PetscFunctionReturn(0);
581}
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.

592{
593 PetscErrorCode ierr;
594 DMDALocalInfo info = user->info;
595 Cmpnts ***rhs;
596
597 // --- Grid extents for this MPI rank and global grid dimensions ---
598 const PetscInt xs = info.xs, xe = xs + info.xm;
599 const PetscInt ys = info.ys, ye = ys + info.ym;
600 const PetscInt zs = info.zs, ze = zs + info.zm;
601 const PetscInt mx = info.mx, my = info.my, mz = info.mz;
602
603 PetscFunctionBeginUser;
604 PROFILE_FUNCTION_BEGIN;
605
606 // Get a writable pointer to the local data of the global RHS vector.
607 ierr = DMDAVecGetArray(user->fda, user->Rhs, &rhs); CHKERRQ(ierr);
608
609 // ========================================================================
610 // --- I-DIRECTION (X-FACES) ---
611 // ========================================================================
612
613 // --- Negative X Face (i=0, the first physical face) ---
614 if (xs == 0) {
615 // This logic applies ONLY to physical (non-periodic) boundaries.
616 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC) {
617 const PetscInt i = 0;
618 for (PetscInt k = zs; k < ze; k++) {
619 for (PetscInt j = ys; j < ye; j++) {
620 rhs[k][j][i].x = 0.0;
621 rhs[k][j][i].y = 0.0;
622 rhs[k][j][i].z = 0.0;
623 }
624 }
625 }
626 }
627
628 // --- Positive X Face (physical face i=mx-2, dummy location i=mx-1) ---
629 if (xe == mx) {
630 // Step 1: Enforce strong BC on the LAST PHYSICAL face (i=mx-2) for non-periodic cases.
631 if (user->boundary_faces[BC_FACE_POS_X].mathematical_type != PERIODIC) {
632 const PetscInt i = mx - 2;
633 for (PetscInt k = zs; k < ze; k++) {
634 for (PetscInt j = ys; j < ye; j++) {
635 rhs[k][j][i].x = 0.0;
636 }
637 }
638 }
639 // Step 2: Unconditionally sanitize the DUMMY location (i=mx-1).
640 const PetscInt i = mx - 1;
641 for (PetscInt k = zs; k < ze; k++) {
642 for (PetscInt j = ys; j < ye; j++) {
643 rhs[k][j][i].x = 0.0;
644 rhs[k][j][i].y = 0.0;
645 rhs[k][j][i].z = 0.0;
646 }
647 }
648 }
649
650 // ========================================================================
651 // --- J-DIRECTION (Y-FACES) ---
652 // ========================================================================
653
654 // --- Negative Y Face (j=0, the first physical face) ---
655 if (ys == 0) {
656 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC) {
657 const PetscInt j = 0;
658 for (PetscInt k = zs; k < ze; k++) {
659 for (PetscInt i = xs; i < xe; i++) {
660 rhs[k][j][i].x = 0.0;
661 rhs[k][j][i].y = 0.0;
662 rhs[k][j][i].z = 0.0;
663 }
664 }
665 }
666 }
667
668 // --- Positive Y Face (physical face j=my-2, dummy location j=my-1) ---
669 if (ye == my) {
670 if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type != PERIODIC) {
671 const PetscInt j = my - 2;
672 for (PetscInt k = zs; k < ze; k++) {
673 for (PetscInt i = xs; i < xe; i++) {
674 rhs[k][j][i].y = 0.0;
675 }
676 }
677 }
678 const PetscInt j = my - 1;
679 for (PetscInt k = zs; k < ze; k++) {
680 for (PetscInt i = xs; i < xe; i++) {
681 rhs[k][j][i].x = 0.0;
682 rhs[k][j][i].y = 0.0;
683 rhs[k][j][i].z = 0.0;
684 }
685 }
686 }
687
688 // ========================================================================
689 // --- K-DIRECTION (Z-FACES) ---
690 // ========================================================================
691
692 // --- Negative Z Face (k=0, the first physical face) ---
693 if (zs == 0) {
694 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC) {
695 const PetscInt k = 0;
696 for (PetscInt j = ys; j < ye; j++) {
697 for (PetscInt i = xs; i < xe; i++) {
698 rhs[k][j][i].x = 0.0;
699 rhs[k][j][i].y = 0.0;
700 rhs[k][j][i].z = 0.0;
701 }
702 }
703 }
704 }
705
706 // --- Positive Z Face (physical face k=mz-2, dummy location k=mz-1) ---
707 if (ze == mz) {
708 if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type != PERIODIC) {
709 const PetscInt k = mz - 2;
710 for (PetscInt j = ys; j < ye; j++) {
711 for (PetscInt i = xs; i < xe; i++) {
712 rhs[k][j][i].z = 0.0;
713 }
714 }
715 }
716 const PetscInt k = mz - 1;
717 for (PetscInt j = ys; j < ye; j++) {
718 for (PetscInt i = xs; i < xe; i++) {
719 rhs[k][j][i].x = 0.0;
720 rhs[k][j][i].y = 0.0;
721 rhs[k][j][i].z = 0.0;
722 }
723 }
724 }
725
726 // --- Release the pointer to the local data ---
727 ierr = DMDAVecRestoreArray(user->fda, user->Rhs, &rhs); CHKERRQ(ierr);
728
729 LOG_ALLOW(LOCAL, LOG_TRACE, "Rank %d, Block %d: Finished enforcing RHS boundary conditions.\n",
730 user->simCtx->rank, user->_this);
731
732 PROFILE_FUNCTION_END;
733
734 PetscFunctionReturn(0);
735}
PERIODIC
@ PERIODIC
Definition variables.h:260
UserCtx::Rhs
Vec Rhs
Definition variables.h:845
SimCtx::rank
PetscMPIInt rank
Definition variables.h:646
UserCtx::boundary_faces
BoundaryFaceConfig boundary_faces[6]
Definition variables.h:829
UserCtx::_this
PetscInt _this
Definition variables.h:824
Cmpnts::x
PetscScalar x
Definition variables.h:101
Cmpnts::z
PetscScalar z
Definition variables.h:101
UserCtx::info
DMDALocalInfo info
Definition variables.h:818
Cmpnts::y
PetscScalar y
Definition variables.h:101
BoundaryFaceConfig::mathematical_type
BCType mathematical_type
Definition variables.h:336
Cmpnts
A 3D point or vector with PetscScalar components.
Definition variables.h:100
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.

745{
746 PetscErrorCode ierr;
747 PetscFunctionBeginUser;
748
749 const char* handler_name = BCHandlerTypeToString(handler_type);
750 LOG_ALLOW(LOCAL, LOG_DEBUG, "Factory called for handler type %s. \n", handler_name);
751
752 ierr = PetscMalloc1(1, new_bc_ptr); CHKERRQ(ierr);
753 BoundaryCondition *bc = *new_bc_ptr;
754
755 bc->type = handler_type;
756 bc->priority = -1; // Default priority; can be overridden in specific handlers
757 bc->data = NULL;
758 bc->Initialize = NULL;
759 bc->PreStep = NULL;
760 bc->Apply = NULL;
761 bc->PostStep = NULL;
762 bc->UpdateUbcs = NULL;
763 bc->Destroy = NULL;
764
765 LOG_ALLOW(LOCAL, LOG_DEBUG, "Allocated generic handler object at address %p.\n", (void*)bc);
766
767 switch (handler_type) {
768
769 case BC_HANDLER_OUTLET_CONSERVATION:
770 LOG_ALLOW(LOCAL, LOG_DEBUG, "Dispatching to Create_OutletConservation().\n");
771 ierr = Create_OutletConservation(bc); CHKERRQ(ierr);
772 break;
773
774 case BC_HANDLER_WALL_NOSLIP:
775 LOG_ALLOW(LOCAL, LOG_DEBUG, "Dispatching to Create_WallNoSlip().\n");
776 ierr = Create_WallNoSlip(bc); CHKERRQ(ierr);
777 break;
778
779 case BC_HANDLER_INLET_CONSTANT_VELOCITY:
780 LOG_ALLOW(LOCAL, LOG_DEBUG, "Dispatching to Create_InletConstantVelocity().\n");
781 ierr = Create_InletConstantVelocity(bc); CHKERRQ(ierr);
782 break;
783
784 case BC_HANDLER_PERIODIC_GEOMETRIC:
785 LOG_ALLOW(LOCAL,LOG_DEBUG,"Dispatching to Create_PeriodicGeometric().\n");
786 ierr = Create_PeriodicGeometric(bc);
787 break;
788
789 case BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX:
790 LOG_ALLOW(LOCAL,LOG_DEBUG,"Dispatching to Create_PeriodicDrivenConstant().\n");
791 ierr = Create_PeriodicDrivenConstant(bc);
792 break;
793
794 //case BC_HANDLER_PERIODIC_DRIVEN_INITIAL_FLUX:
795 // LOG_ALLOW(LOCAL,LOG_DEBUG,"Dispatching to Create_PeriodicDrivenInitial().\n");
796 // ierr = Create_PeriodicDrivenInitial(bc);
797 // break;
798
799 case BC_HANDLER_INLET_PARABOLIC:
800 LOG_ALLOW(LOCAL, LOG_DEBUG, "Dispatching to Create_InletParabolicProfile().\n");
801 ierr = Create_InletParabolicProfile(bc); CHKERRQ(ierr);
802 break;
803 //Add cases for other handlers here in future phases
804
805 default:
806 LOG_ALLOW(GLOBAL, LOG_ERROR, "Handler type (%s) is not recognized or implemented in the factory.\n", handler_name);
807 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Boundary handler type %d (%s) not recognized in factory.\n", handler_type, handler_name);
808 }
809
810 if(bc->priority < 0) {
811 LOG_ALLOW(GLOBAL, LOG_ERROR, "Handler type %d (%s) did not set a valid priority during creation.\n", handler_type, handler_name);
812 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);
813 }
814
815 LOG_ALLOW(LOCAL, LOG_DEBUG, "Successfully created and configured handler for %s.\n", handler_name);
816 PetscFunctionReturn(0);
817}
Create_InletConstantVelocity
PetscErrorCode Create_InletConstantVelocity(BoundaryCondition *bc)
Configures a BoundaryCondition object to behave as a constant velocity inlet.
Definition BC_Handlers.c:323
Create_PeriodicGeometric
PetscErrorCode Create_PeriodicGeometric(BoundaryCondition *bc)
Configures a BoundaryCondition object for geometric periodic coupling.
Definition BC_Handlers.c:1646
Create_InletParabolicProfile
PetscErrorCode Create_InletParabolicProfile(BoundaryCondition *bc)
Configures a BoundaryCondition object for a parabolic inlet profile.
Definition BC_Handlers.c:721
Create_PeriodicDrivenConstant
PetscErrorCode Create_PeriodicDrivenConstant(BoundaryCondition *bc)
Configures a BoundaryCondition object for periodic driven-flow forcing.
Definition BC_Handlers.c:1699
Create_WallNoSlip
PetscErrorCode Create_WallNoSlip(BoundaryCondition *bc)
Configures a BoundaryCondition object to behave as a no-slip, stationary wall.
Definition BC_Handlers.c:114
Create_OutletConservation
PetscErrorCode Create_OutletConservation(BoundaryCondition *bc)
Configures a BoundaryCondition object for conservative outlet treatment.
Definition BC_Handlers.c:1129
BCHandlerTypeToString
const char * BCHandlerTypeToString(BCHandlerType handler_type)
Converts a BCHandlerType enum to its string representation.
Definition logging.c:771
LOG_ERROR
@ LOG_ERROR
Critical errors that may halt the program.
Definition logging.h:28
BoundaryCondition
The "virtual table" struct for a boundary condition handler object.
Definition variables.h:321
BoundaryCondition::PostStep
PetscErrorCode(* PostStep)(BoundaryCondition *self, BCContext *ctx, PetscReal *local_inflow, PetscReal *local_outflow)
Definition variables.h:328
BoundaryCondition::PreStep
PetscErrorCode(* PreStep)(BoundaryCondition *self, BCContext *ctx, PetscReal *local_inflow, PetscReal *local_outflow)
Definition variables.h:326
BoundaryCondition::type
BCHandlerType type
Definition variables.h:322
BoundaryCondition::Destroy
PetscErrorCode(* Destroy)(BoundaryCondition *self)
Definition variables.h:330
BoundaryCondition::Initialize
PetscErrorCode(* Initialize)(BoundaryCondition *self, BCContext *ctx)
Definition variables.h:325
BoundaryCondition::UpdateUbcs
PetscErrorCode(* UpdateUbcs)(BoundaryCondition *self, BCContext *ctx)
Definition variables.h:329
BoundaryCondition::Apply
PetscErrorCode(* Apply)(BoundaryCondition *self, BCContext *ctx)
Definition variables.h:327
BoundaryCondition::priority
BCPriorityType priority
Definition variables.h:323
BC_HANDLER_PERIODIC_GEOMETRIC
@ BC_HANDLER_PERIODIC_GEOMETRIC
Definition variables.h:284
BC_HANDLER_INLET_PARABOLIC
@ BC_HANDLER_INLET_PARABOLIC
Definition variables.h:277
BC_HANDLER_INLET_CONSTANT_VELOCITY
@ BC_HANDLER_INLET_CONSTANT_VELOCITY
Definition variables.h:276
BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX
@ BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX
Definition variables.h:286
BC_HANDLER_WALL_NOSLIP
@ BC_HANDLER_WALL_NOSLIP
Definition variables.h:273
BC_HANDLER_OUTLET_CONSERVATION
@ BC_HANDLER_OUTLET_CONSERVATION
Definition variables.h:282
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 825 of file Boundaries.c.

826{
827 PetscErrorCode ierr;
828 PetscFunctionBeginUser;
829
830 LOG_ALLOW(GLOBAL, LOG_INFO, "Validating parsed boundary condition configuration...\n");
831
832 // --- Rule Set 1: Driven Flow Handler Consistency ---
833 // This specialized validator will check all rules related to driven flow handlers.
834 ierr = Validate_DrivenFlowConfiguration(user); CHKERRQ(ierr);
835
836 // --- Rule Set 2: (Future Extension) Overset Interface Consistency ---
837 // ierr = Validate_OversetConfiguration(user); CHKERRQ(ierr);
838
839 LOG_ALLOW(GLOBAL, LOG_INFO, "Boundary configuration is valid.\n");
840
841 PetscFunctionReturn(0);
842}
Validate_DrivenFlowConfiguration
PetscErrorCode Validate_DrivenFlowConfiguration(UserCtx *user)
(Private) Validates all consistency rules for a driven flow (channel/pipe) setup.
Definition BC_Handlers.c:15
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 857 of file Boundaries.c.

858{
859 PetscErrorCode ierr;
860 PetscFunctionBeginUser;
861
862 LOG_ALLOW(GLOBAL, LOG_INFO, "Starting creation and initialization of all boundary handlers.\n");
863
864 // =========================================================================
865 // Step 0: Clear any existing boundary handlers (if re-initializing).
866 // This ensures no memory leaks if this function is called multiple times.
867 // =========================================================================
868 for (int i = 0; i < 6; i++) {
869 BoundaryFaceConfig *face_cfg = &user->boundary_faces[i];
870 if (face_cfg->handler) {
871 LOG_ALLOW(LOCAL, LOG_DEBUG, "Destroying existing handler on Face %s before re-initialization.\n", BCFaceToString((BCFace)i));
872 if (face_cfg->handler->Destroy) {
873 ierr = face_cfg->handler->Destroy(face_cfg->handler); CHKERRQ(ierr);
874 }
875 ierr = PetscFree(face_cfg->handler); CHKERRQ(ierr);
876 face_cfg->handler = NULL;
877 }
878 }
879 // =========================================================================
880
881 // Step 0.1: Initiate flux sums to zero
882 user->simCtx->FluxInSum = 0.0;
883 user->simCtx->FluxOutSum = 0.0;
884 user->simCtx->FarFluxInSum = 0.0;
885 user->simCtx->FarFluxOutSum = 0.0;
886 // =========================================================================
887
888 // Step 1: Parse the configuration file to determine user intent.
889 // This function, defined in io.c, populates the configuration enums and parameter
890 // lists within the user->boundary_faces array on all MPI ranks.
891 ierr = ParseAllBoundaryConditions(user, bcs_filename); CHKERRQ(ierr);
892 LOG_ALLOW(GLOBAL, LOG_INFO, "Configuration file '%s' parsed successfully.\n", bcs_filename);
893
894 // Step 1.1: Validate the parsed configuration to ensure there are no Boundary Condition conflicts
895 ierr = BoundarySystem_Validate(user); CHKERRQ(ierr);
896
897 // Step 2: Create and Initialize the handler object for each of the 6 faces.
898 for (int i = 0; i < 6; i++) {
899 BoundaryFaceConfig *face_cfg = &user->boundary_faces[i];
900
901 const char *face_name = BCFaceToString(face_cfg->face_id);
902 const char *type_name = BCTypeToString(face_cfg->mathematical_type);
903 const char *handler_name = BCHandlerTypeToString(face_cfg->handler_type);
904
905 LOG_ALLOW(LOCAL, LOG_DEBUG, "Creating handler for Face %s with Type %s and handler '%s'.\n", face_name, type_name,handler_name);
906
907 // Use the private factory to construct the correct handler object based on the parsed type.
908 // The factory returns a pointer to the new handler object, which we store in the config struct.
909 ierr = BoundaryCondition_Create(face_cfg->handler_type, &face_cfg->handler); CHKERRQ(ierr);
910
911 // Step 3: Call the specific Initialize() method for the newly created handler.
912 // This allows the handler to perform its own setup, like reading parameters from the
913 // face_cfg->params list and setting the initial field values on its face.
914 if (face_cfg->handler && face_cfg->handler->Initialize) {
915 LOG_ALLOW(LOCAL, LOG_DEBUG, "Calling Initialize() method for handler %s(%s) on Face %s.\n",type_name,handler_name,face_name);
916
917 // Prepare the context needed by the Initialize() function.
918 BCContext ctx = {
919 .user = user,
920 .face_id = face_cfg->face_id,
921 .global_inflow_sum = &user->simCtx->FluxInSum, // Global flux sums are not relevant during initialization.
922 .global_outflow_sum = &user->simCtx->FluxOutSum,
923 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
924 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
925 };
926
927 ierr = face_cfg->handler->Initialize(face_cfg->handler, &ctx); CHKERRQ(ierr);
928 } else {
929 LOG_ALLOW(LOCAL, LOG_DEBUG, "Handler %s(%s) for Face %s has no Initialize() method, skipping.\n", type_name,handler_name,face_name);
930 }
931 }
932 // =========================================================================
933 // NO SYNCHRONIZATION NEEDED HERE
934 // =========================================================================
935 // Initialize() only reads parameters and allocates memory.
936 // It does NOT modify field values (Ucat, Ucont, Ubcs).
937 // Field values are set by:
938 // 1. Initial conditions (before this function)
939 // 2. Apply() during timestepping (after this function)
940 // The first call to ApplyBoundaryConditions() will handle synchronization.
941 // =========================================================================
942
943 // ====================================================================================
944 // --- NEW: Step 4: Synchronize Vectors After Initialization ---
945 // This is the CRITICAL fix. The Initialize() calls have modified local vector
946 // arrays on some ranks but not others. We must now update the global vector state
947 // and then update all local ghost regions to be consistent.
948 // ====================================================================================
949
950 //LOG_ALLOW(GLOBAL, LOG_DEBUG, "Committing global boundary initializations to local vectors.\n");
951
952 // Commit changes from the global vectors (Ucat, Ucont) to the local vectors (lUcat, lUcont)
953 // NOTE: The Apply functions modified Ucat and Ucont via GetArray, which works on the global
954 // representation.
955 /*
956 ierr = DMGlobalToLocalBegin(user->fda, user->Ucat, INSERT_VALUES, user->lUcat); CHKERRQ(ierr);
957 ierr = DMGlobalToLocalEnd(user->fda, user->Ucat, INSERT_VALUES, user->lUcat); CHKERRQ(ierr);
958
959 ierr = DMGlobalToLocalBegin(user->fda, user->Ucont, INSERT_VALUES, user->lUcont); CHKERRQ(ierr);
960 ierr = DMGlobalToLocalEnd(user->fda, user->Ucont, INSERT_VALUES, user->lUcont); CHKERRQ(ierr);
961
962 // Now, update all local vectors (including ghost cells) from the newly consistent global vectors
963
964 ierr = DMLocalToGlobalBegin(user->fda, user->lUcat, INSERT_VALUES, user->Ucat); CHKERRQ(ierr);
965 ierr = DMLocalToGlobalEnd(user->fda, user->lUcat, INSERT_VALUES, user->Ucat); CHKERRQ(ierr);
966
967 ierr = DMLocalToGlobalBegin(user->fda, user->lUcont, INSERT_VALUES, user->Ucont); CHKERRQ(ierr);
968 ierr = DMLocalToGlobalEnd(user->fda, user->lUcont, INSERT_VALUES, user->Ucont); CHKERRQ(ierr);
969 */
970
971 LOG_ALLOW(GLOBAL, LOG_INFO, "All boundary handlers created and initialized successfully.\n");
972 PetscFunctionReturn(0);
973}
BoundarySystem_Validate
PetscErrorCode BoundarySystem_Validate(UserCtx *user)
Internal helper implementation: BoundarySystem_Validate().
Definition Boundaries.c:825
BoundaryCondition_Create
PetscErrorCode BoundaryCondition_Create(BCHandlerType handler_type, BoundaryCondition **new_bc_ptr)
Internal helper implementation: BoundaryCondition_Create().
Definition Boundaries.c:744
ParseAllBoundaryConditions
PetscErrorCode ParseAllBoundaryConditions(UserCtx *user, const char *bcs_input_filename)
Parses the boundary conditions file to configure the type, handler, and any associated parameters for...
Definition io.c:447
BCTypeToString
const char * BCTypeToString(BCType type)
Helper function to convert BCType enum to a string representation.
Definition logging.c:751
SimCtx::FarFluxInSum
PetscReal FarFluxInSum
Definition variables.h:721
SimCtx::FarFluxOutSum
PetscReal FarFluxOutSum
Definition variables.h:721
SimCtx::FluxOutSum
PetscReal FluxOutSum
Definition variables.h:721
BoundaryFaceConfig::handler_type
BCHandlerType handler_type
Definition variables.h:337
BCContext::user
UserCtx * user
Definition variables.h:312
SimCtx::FluxInSum
PetscReal FluxInSum
Definition variables.h:721
BoundaryFaceConfig::handler
BoundaryCondition * handler
Definition variables.h:339
BCContext
Provides execution context for a boundary condition handler.
Definition variables.h:311
BoundaryFaceConfig
Holds the complete configuration for one of the six boundary faces.
Definition variables.h:334
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 982 of file Boundaries.c.

983{
984 PetscErrorCode ierr;
985 UserMG *usermg = &simCtx->usermg;
986
987 PetscFunctionBeginUser;
988 PROFILE_FUNCTION_BEGIN;
989
990 LOG_ALLOW(GLOBAL, LOG_INFO, "Propagating BC configuration from finest to coarser multigrid levels...\n");
991
992 // Loop from second-finest down to coarsest
993 for (PetscInt level = usermg->mglevels - 2; level >= 0; level--) {
994 for (PetscInt bi = 0; bi < simCtx->block_number; bi++) {
995 UserCtx *user_coarse = &usermg->mgctx[level].user[bi];
996 UserCtx *user_fine = &usermg->mgctx[level + 1].user[bi];
997
998 LOG_ALLOW_SYNC(LOCAL, LOG_DEBUG, "Rank %d: Copying BC config from level %d to level %d, block %d\n",
999 simCtx->rank, level + 1, level, bi);
1000
1001 // Copy the 6 boundary face configurations
1002 for (int face_i = 0; face_i < 6; face_i++) {
1003 user_coarse->boundary_faces[face_i].face_id = user_fine->boundary_faces[face_i].face_id;
1004 user_coarse->boundary_faces[face_i].mathematical_type = user_fine->boundary_faces[face_i].mathematical_type;
1005 user_coarse->boundary_faces[face_i].handler_type = user_fine->boundary_faces[face_i].handler_type;
1006
1007 // Copy parameter list (deep copy)
1008 FreeBC_ParamList(user_coarse->boundary_faces[face_i].params); // Clear any existing
1009 user_coarse->boundary_faces[face_i].params = NULL;
1010
1011 BC_Param **dst_next = &user_coarse->boundary_faces[face_i].params;
1012 for (BC_Param *src = user_fine->boundary_faces[face_i].params; src; src = src->next) {
1013 BC_Param *new_param;
1014 ierr = PetscMalloc1(1, &new_param); CHKERRQ(ierr);
1015 ierr = PetscStrallocpy(src->key, &new_param->key); CHKERRQ(ierr);
1016 ierr = PetscStrallocpy(src->value, &new_param->value); CHKERRQ(ierr);
1017 new_param->next = NULL;
1018 *dst_next = new_param;
1019 dst_next = &new_param->next;
1020 }
1021
1022 // IMPORTANT: Do NOT create handler objects for coarser levels
1023 // Handlers are only needed at finest level for timestepping Apply() calls
1024 user_coarse->boundary_faces[face_i].handler = NULL;
1025 }
1026
1027 // Propagate the particle inlet lookup fields to coarse levels as well.
1028 user_coarse->inletFaceDefined = user_fine->inletFaceDefined;
1029 user_coarse->identifiedInletBCFace = user_fine->identifiedInletBCFace;
1030 }
1031 }
1032
1033 LOG_ALLOW(GLOBAL, LOG_INFO, "BC configuration propagation complete.\n");
1034
1035 PROFILE_FUNCTION_END;
1036 PetscFunctionReturn(0);
1037}
FreeBC_ParamList
void FreeBC_ParamList(BC_Param *head)
Frees an entire linked list of boundary-condition parameters.
Definition io.c:302
LOG_ALLOW_SYNC
#define LOG_ALLOW_SYNC(scope, level, fmt,...)
Synchronized logging macro that checks both the log level and whether the calling function is in the ...
Definition logging.h:252
MGCtx::user
UserCtx * user
Definition variables.h:528
SimCtx::block_number
PetscInt block_number
Definition variables.h:712
BC_Param_s::next
struct BC_Param_s * next
Definition variables.h:307
BC_Param_s::key
char * key
Definition variables.h:305
SimCtx::usermg
UserMG usermg
Definition variables.h:764
BC_Param_s::value
char * value
Definition variables.h:306
BoundaryFaceConfig::params
BC_Param * params
Definition variables.h:338
UserMG::mglevels
PetscInt mglevels
Definition variables.h:535
UserMG::mgctx
MGCtx * mgctx
Definition variables.h:538
BC_Param_s
A node in a linked list for storing key-value parameters from the bcs.dat file.
Definition variables.h:304
UserCtx
User-defined context containing data specific to a single computational grid level.
Definition variables.h:811
UserMG
User-level context for managing the entire multigrid hierarchy.
Definition variables.h:534
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 1053 of file Boundaries.c.

1054{
1055 PetscErrorCode ierr;
1056 PetscFunctionBeginUser;
1057 PROFILE_FUNCTION_BEGIN;
1058
1059 LOG_ALLOW(LOCAL, LOG_DEBUG, "Starting.\n");
1060
1061 // =========================================================================
1062 // PRIORITY 0: INLETS
1063 // =========================================================================
1064
1065 PetscReal local_inflow_pre = 0.0;
1066 PetscReal local_inflow_post = 0.0;
1067 PetscReal global_inflow_pre = 0.0;
1068 PetscReal global_inflow_post = 0.0;
1069 PetscInt num_handlers[3] = {0,0,0};
1070
1071 LOG_ALLOW(LOCAL, LOG_TRACE, " (INLETS): Begin.\n");
1072
1073 // Phase 1: PreStep - Preparation (e.g., calculate profiles, read files)
1074 for (int i = 0; i < 6; i++) {
1075 BoundaryCondition *handler = user->boundary_faces[i].handler;
1076 if (!handler || handler->priority != BC_PRIORITY_INLET) continue;
1077 if (!handler->PreStep) continue;
1078
1079 num_handlers[0]++;
1080 BCContext ctx = {
1081 .user = user,
1082 .face_id = (BCFace)i,
1083 .global_inflow_sum = NULL,
1084 .global_outflow_sum = NULL,
1085 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1086 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1087 };
1088
1089 LOG_ALLOW(LOCAL, LOG_TRACE, " PreStep: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1090 ierr = handler->PreStep(handler, &ctx, &local_inflow_pre, NULL); CHKERRQ(ierr);
1091 }
1092
1093 // Optional: Global communication for PreStep (for debugging)
1094 if (local_inflow_pre != 0.0) {
1095 ierr = MPI_Allreduce(&local_inflow_pre, &global_inflow_pre, 1, MPIU_REAL,
1096 MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
1097 LOG_ALLOW(GLOBAL, LOG_TRACE, " PreStep predicted flux: %.6e\n", global_inflow_pre);
1098 }
1099
1100 // Phase 2: Apply - Set boundary conditions
1101 for (int i = 0; i < 6; i++) {
1102 BoundaryCondition *handler = user->boundary_faces[i].handler;
1103 if (!handler || handler->priority != BC_PRIORITY_INLET) continue;
1104 if(!handler->Apply) continue; // For example Periodic BCs
1105
1106 num_handlers[1]++;
1107
1108 BCContext ctx = {
1109 .user = user,
1110 .face_id = (BCFace)i,
1111 .global_inflow_sum = NULL,
1112 .global_outflow_sum = NULL,
1113 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1114 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1115 };
1116
1117 LOG_ALLOW(LOCAL, LOG_TRACE, " Apply: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1118 ierr = handler->Apply(handler, &ctx); CHKERRQ(ierr);
1119 }
1120
1121 // Phase 3: PostStep - Measure actual flux
1122 for (int i = 0; i < 6; i++) {
1123 BoundaryCondition *handler = user->boundary_faces[i].handler;
1124 if (!handler || handler->priority != BC_PRIORITY_INLET) continue;
1125 if (!handler->PostStep) continue;
1126
1127 num_handlers[2]++;
1128
1129 BCContext ctx = {
1130 .user = user,
1131 .face_id = (BCFace)i,
1132 .global_inflow_sum = NULL,
1133 .global_outflow_sum = NULL,
1134 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1135 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1136 };
1137
1138 LOG_ALLOW(LOCAL, LOG_TRACE, " PostStep: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1139 ierr = handler->PostStep(handler, &ctx, &local_inflow_post, NULL); CHKERRQ(ierr);
1140 }
1141
1142 // Phase 4: Global communication - Sum flux for other priorities to use
1143 ierr = MPI_Allreduce(&local_inflow_post, &global_inflow_post, 1, MPIU_REAL,
1144 MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
1145
1146 // Store for next priority levels
1147 user->simCtx->FluxInSum = global_inflow_post;
1148
1149 LOG_ALLOW(GLOBAL, LOG_INFO,
1150 " (INLETS): %d Prestep(s), %d Application(s), %d Poststep(s), FluxInSum = %.6e\n",
1151 num_handlers[0],num_handlers[1],num_handlers[2], global_inflow_post);
1152
1153 // =========================================================================
1154 // PRIORITY 1: FARFIELD
1155 // =========================================================================
1156
1157 PetscReal local_farfield_in_pre = 0.0;
1158 PetscReal local_farfield_out_pre = 0.0;
1159 PetscReal local_farfield_in_post = 0.0;
1160 PetscReal local_farfield_out_post = 0.0;
1161 PetscReal global_farfield_in_pre = 0.0;
1162 PetscReal global_farfield_out_pre = 0.0;
1163 PetscReal global_farfield_in_post = 0.0;
1164 PetscReal global_farfield_out_post = 0.0;
1165 memset(num_handlers,0,sizeof(num_handlers));
1166
1167 LOG_ALLOW(LOCAL, LOG_TRACE, " (FARFIELD): Begin.\n");
1168
1169 // Phase 1: PreStep - Analyze flow direction, measure initial flux
1170 for (int i = 0; i < 6; i++) {
1171 BoundaryCondition *handler = user->boundary_faces[i].handler;
1172 if (!handler || handler->priority != BC_PRIORITY_FARFIELD) continue;
1173 if (!handler->PreStep) continue;
1174
1175 num_handlers[0]++;
1176 BCContext ctx = {
1177 .user = user,
1178 .face_id = (BCFace)i,
1179 .global_inflow_sum = &user->simCtx->FluxInSum, // Available from Priority 0
1180 .global_outflow_sum = NULL,
1181 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1182 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1183 };
1184
1185 LOG_ALLOW(LOCAL, LOG_TRACE, " PreStep: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1186 ierr = handler->PreStep(handler, &ctx, &local_farfield_in_pre, &local_farfield_out_pre);
1187 CHKERRQ(ierr);
1188 }
1189
1190 // Phase 2: Global communication (optional, for debugging)
1191 if (local_farfield_in_pre != 0.0 || local_farfield_out_pre != 0.0) {
1192 ierr = MPI_Allreduce(&local_farfield_in_pre, &global_farfield_in_pre, 1, MPIU_REAL,
1193 MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
1194 ierr = MPI_Allreduce(&local_farfield_out_pre, &global_farfield_out_pre, 1, MPIU_REAL,
1195 MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
1196
1197 LOG_ALLOW(GLOBAL, LOG_DEBUG,
1198 " Farfield pre-analysis: In=%.6e, Out=%.6e\n",
1199 global_farfield_in_pre, global_farfield_out_pre);
1200 }
1201
1202 // Phase 3: Apply - Set farfield boundary conditions
1203 for (int i = 0; i < 6; i++) {
1204 BoundaryCondition *handler = user->boundary_faces[i].handler;
1205 if (!handler || handler->priority != BC_PRIORITY_FARFIELD) continue;
1206 if(!handler->Apply) continue; // For example Periodic BCs
1207
1208 num_handlers[1]++;
1209
1210 BCContext ctx = {
1211 .user = user,
1212 .face_id = (BCFace)i,
1213 .global_inflow_sum = &user->simCtx->FluxInSum,
1214 .global_outflow_sum = NULL,
1215 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1216 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1217 };
1218
1219 LOG_ALLOW(LOCAL, LOG_TRACE, " Apply: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1220 ierr = handler->Apply(handler, &ctx); CHKERRQ(ierr);
1221 }
1222
1223 // Phase 4: PostStep - Measure actual farfield fluxes
1224 for (int i = 0; i < 6; i++) {
1225 BoundaryCondition *handler = user->boundary_faces[i].handler;
1226 if (!handler || handler->priority != BC_PRIORITY_FARFIELD) continue;
1227 if (!handler->PostStep) continue;
1228
1229 num_handlers[2]++;
1230
1231 BCContext ctx = {
1232 .user = user,
1233 .face_id = (BCFace)i,
1234 .global_inflow_sum = &user->simCtx->FluxInSum,
1235 .global_outflow_sum = NULL,
1236 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1237 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1238 };
1239
1240 LOG_ALLOW(LOCAL, LOG_TRACE, " PostStep: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1241 ierr = handler->PostStep(handler, &ctx, &local_farfield_in_post, &local_farfield_out_post);
1242 CHKERRQ(ierr);
1243 }
1244
1245 // Phase 5: Global communication - Store for outlet priority
1246 if (num_handlers > 0) {
1247 ierr = MPI_Allreduce(&local_farfield_in_post, &global_farfield_in_post, 1, MPIU_REAL,
1248 MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
1249 ierr = MPI_Allreduce(&local_farfield_out_post, &global_farfield_out_post, 1, MPIU_REAL,
1250 MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
1251
1252 // Store for outlet handlers to use
1253 user->simCtx->FarFluxInSum = global_farfield_in_post;
1254 user->simCtx->FarFluxOutSum = global_farfield_out_post;
1255
1256 LOG_ALLOW(GLOBAL, LOG_INFO,
1257 " (FARFIELD): %d Prestep(s), %d Application(s), %d Poststep(s) , InFlux=%.6e, OutFlux=%.6e\n",
1258 num_handlers[0],num_handlers[1],num_handlers[2], global_farfield_in_post, global_farfield_out_post);
1259 } else {
1260 // No farfield handlers - zero out the fluxes
1261 user->simCtx->FarFluxInSum = 0.0;
1262 user->simCtx->FarFluxOutSum = 0.0;
1263 }
1264
1265
1266 // =========================================================================
1267 // PRIORITY 2: WALLS
1268 // =========================================================================
1269
1270 memset(num_handlers,0,sizeof(num_handlers));
1271
1272 LOG_ALLOW(LOCAL, LOG_TRACE, " (WALLS): Begin.\n");
1273
1274 // Phase 1: PreStep - Preparation (usually no-op for walls)
1275 for (int i = 0; i < 6; i++) {
1276 BoundaryCondition *handler = user->boundary_faces[i].handler;
1277 if (!handler || handler->priority != BC_PRIORITY_WALL) continue;
1278 if (!handler->PreStep) continue;
1279
1280 num_handlers[0]++;
1281 BCContext ctx = {
1282 .user = user,
1283 .face_id = (BCFace)i,
1284 .global_inflow_sum = &user->simCtx->FluxInSum,
1285 .global_outflow_sum = NULL,
1286 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1287 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1288 };
1289
1290 LOG_ALLOW(LOCAL, LOG_TRACE, " PreStep: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1291 ierr = handler->PreStep(handler, &ctx, NULL, NULL); CHKERRQ(ierr);
1292 }
1293
1294 // No global communication needed for walls
1295
1296 // Phase 2: Apply - Set boundary conditions
1297 for (int i = 0; i < 6; i++) {
1298 BoundaryCondition *handler = user->boundary_faces[i].handler;
1299 if (!handler || handler->priority != BC_PRIORITY_WALL) continue;
1300 if(!handler->Apply) continue; // For example Periodic BCs
1301
1302 num_handlers[1]++;
1303
1304 BCContext ctx = {
1305 .user = user,
1306 .face_id = (BCFace)i,
1307 .global_inflow_sum = &user->simCtx->FluxInSum,
1308 .global_outflow_sum = NULL,
1309 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1310 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1311 };
1312
1313 LOG_ALLOW(LOCAL, LOG_TRACE, " Apply: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1314 ierr = handler->Apply(handler, &ctx); CHKERRQ(ierr);
1315 }
1316
1317 // Phase 3: PostStep - Post-application processing (usually no-op for walls)
1318 for (int i = 0; i < 6; i++) {
1319 BoundaryCondition *handler = user->boundary_faces[i].handler;
1320 if (!handler || handler->priority != BC_PRIORITY_WALL) continue;
1321 if (!handler->PostStep) continue;
1322
1323 num_handlers[2]++;
1324
1325 BCContext ctx = {
1326 .user = user,
1327 .face_id = (BCFace)i,
1328 .global_inflow_sum = &user->simCtx->FluxInSum,
1329 .global_outflow_sum = NULL,
1330 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1331 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1332 };
1333
1334 LOG_ALLOW(LOCAL, LOG_TRACE, " PostStep: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1335 ierr = handler->PostStep(handler, &ctx, NULL, NULL); CHKERRQ(ierr);
1336 }
1337
1338 // No global communication needed for walls
1339
1340 LOG_ALLOW(GLOBAL, LOG_INFO, " (WALLS): %d Prestep(s), %d Application(s), %d Poststep(s) applied.\n",
1341 num_handlers[0],num_handlers[1],num_handlers[2]);
1342
1343
1344 // =========================================================================
1345 // PRIORITY 3: OUTLETS
1346 // =========================================================================
1347
1348 PetscReal local_outflow_pre = 0.0;
1349 PetscReal local_outflow_post = 0.0;
1350 PetscReal global_outflow_pre = 0.0;
1351 PetscReal global_outflow_post = 0.0;
1352 memset(num_handlers,0,sizeof(num_handlers));
1353
1354 LOG_ALLOW(LOCAL, LOG_TRACE, " (OUTLETS): Begin.\n");
1355
1356 // Phase 1: PreStep - Measure uncorrected outflow (from ucat)
1357 for (int i = 0; i < 6; i++) {
1358 BoundaryCondition *handler = user->boundary_faces[i].handler;
1359 if (!handler || handler->priority != BC_PRIORITY_OUTLET) continue;
1360 if (!handler->PreStep) continue;
1361
1362 num_handlers[0]++;
1363 BCContext ctx = {
1364 .user = user,
1365 .face_id = (BCFace)i,
1366 .global_inflow_sum = &user->simCtx->FluxInSum, // From Priority 0
1367 .global_outflow_sum = NULL,
1368 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1369 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1370 };
1371
1372 LOG_ALLOW(LOCAL, LOG_TRACE, " PreStep: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1373 ierr = handler->PreStep(handler, &ctx, NULL, &local_outflow_pre); CHKERRQ(ierr);
1374 }
1375
1376 // Phase 2: Global communication - Get uncorrected outflow sum
1377 ierr = MPI_Allreduce(&local_outflow_pre, &global_outflow_pre, 1, MPIU_REAL,
1378 MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
1379
1380 // Calculate total inflow (inlet + farfield inflow)
1381 PetscReal total_inflow = user->simCtx->FluxInSum + user->simCtx->FarFluxInSum;
1382
1383 LOG_ALLOW(GLOBAL, LOG_DEBUG,
1384 " Uncorrected outflow: %.6e, Total inflow: %.6e (Inlet: %.6e + Farfield: %.6e)\n",
1385 global_outflow_pre, total_inflow, user->simCtx->FluxInSum,
1386 user->simCtx->FarFluxInSum);
1387
1388 // Phase 3: Apply - Set corrected boundary conditions
1389 for (int i = 0; i < 6; i++) {
1390 BoundaryCondition *handler = user->boundary_faces[i].handler;
1391 if (!handler || handler->priority != BC_PRIORITY_OUTLET) continue;
1392 if(!handler->Apply) continue; // For example Periodic BCs
1393
1394 num_handlers[1]++;
1395
1396 BCContext ctx = {
1397 .user = user,
1398 .face_id = (BCFace)i,
1399 .global_inflow_sum = &user->simCtx->FluxInSum, // From Priority 0
1400 .global_outflow_sum = &global_outflow_pre, // From PreStep above
1401 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1402 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1403 };
1404
1405 LOG_ALLOW(LOCAL, LOG_TRACE, " Apply: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1406 ierr = handler->Apply(handler, &ctx); CHKERRQ(ierr);
1407 }
1408
1409 // Phase 4: PostStep - Measure corrected outflow (verification)
1410 for (int i = 0; i < 6; i++) {
1411 BoundaryCondition *handler = user->boundary_faces[i].handler;
1412 if (!handler || handler->priority != BC_PRIORITY_OUTLET) continue;
1413 if (!handler->PostStep) continue;
1414
1415 num_handlers[2]++;
1416
1417 BCContext ctx = {
1418 .user = user,
1419 .face_id = (BCFace)i,
1420 .global_inflow_sum = &user->simCtx->FluxInSum,
1421 .global_outflow_sum = &global_outflow_pre,
1422 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1423 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1424 };
1425
1426 LOG_ALLOW(LOCAL, LOG_TRACE, " PostStep: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1427 ierr = handler->PostStep(handler, &ctx, NULL, &local_outflow_post); CHKERRQ(ierr);
1428 }
1429
1430 // Phase 5: Global communication - Verify conservation
1431 ierr = MPI_Allreduce(&local_outflow_post, &global_outflow_post, 1, MPIU_REAL,
1432 MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
1433
1434 // Store for global reporting.
1435 user->simCtx->FluxOutSum = global_outflow_post;
1436
1437 // Conservation check (compare total outflow vs total inflow)
1438 PetscReal total_outflow = global_outflow_post + user->simCtx->FarFluxOutSum;
1439 PetscReal flux_error = PetscAbsReal(total_outflow - total_inflow);
1440 PetscReal relative_error = (total_inflow > 1e-16) ?
1441 flux_error / total_inflow : flux_error;
1442
1443 LOG_ALLOW(GLOBAL, LOG_INFO,
1444 " (OUTLETS): %d Prestep(s), %d Application(s), %d Poststep(s), FluxOutSum = %.6e\n",
1445 num_handlers[0],num_handlers[1],num_handlers[2], global_outflow_post);
1446 LOG_ALLOW(GLOBAL, LOG_INFO,
1447 " Conservation: Total In=%.6e, Total Out=%.6e, Error=%.3e (%.2e)%%)\n",
1448 total_inflow, total_outflow, flux_error, relative_error * 100.0);
1449
1450 if (relative_error > 1e-6) {
1451 LOG_ALLOW(GLOBAL, LOG_WARNING,
1452 " WARNING: Large mass conservation error (%.2e%%)!\n",
1453 relative_error * 100.0);
1454 }
1455
1456
1457 LOG_ALLOW(LOCAL, LOG_VERBOSE, "Complete.\n");
1458
1459 PROFILE_FUNCTION_END;
1460 PetscFunctionReturn(0);
1461}
BC_PRIORITY_OUTLET
@ BC_PRIORITY_OUTLET
Definition variables.h:296
BC_PRIORITY_FARFIELD
@ BC_PRIORITY_FARFIELD
Definition variables.h:294
BC_PRIORITY_WALL
@ BC_PRIORITY_WALL
Definition variables.h:295
BC_PRIORITY_INLET
@ BC_PRIORITY_INLET
Definition variables.h:293
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 1475 of file Boundaries.c.

1476{
1477 PetscErrorCode ierr;
1478 PetscFunctionBeginUser;
1479
1480 LOG_ALLOW(GLOBAL, LOG_TRACE, "Refreshing `ubcs` targets for flow-dependent boundaries...\n");
1481
1482 // Loop through all 6 faces of the domain
1483 for (int i = 0; i < 6; i++) {
1484 BoundaryCondition *handler = user->boundary_faces[i].handler;
1485
1486 // THE FILTER:
1487 // This is the core logic. We only act if a handler exists for the face
1488 // AND that handler has explicitly implemented the `UpdateUbcs` method.
1489 if (handler && handler->UpdateUbcs) {
1490
1491 const char *face_name = BCFaceToString((BCFace)i);
1492 LOG_ALLOW(LOCAL, LOG_TRACE, " Calling UpdateUbcs() for handler on Face %s.\n", face_name);
1493
1494 // Prepare the context. For this refresh step, we don't need to pass flux sums.
1495 BCContext ctx = {
1496 .user = user,
1497 .face_id = (BCFace)i,
1498 .global_inflow_sum = NULL,
1499 .global_outflow_sum = NULL,
1500 .global_farfield_inflow_sum = NULL,
1501 .global_farfield_outflow_sum = NULL
1502 };
1503
1504 // Call the handler's specific UpdateUbcs function pointer.
1505 ierr = handler->UpdateUbcs(handler, &ctx); CHKERRQ(ierr);
1506 }
1507 }
1508
1509 PetscFunctionReturn(0);
1510}
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 1525 of file Boundaries.c.

1526{
1527 PetscErrorCode ierr;
1528 PetscFunctionBeginUser;
1529
1530
1531
1532 LOG_ALLOW(GLOBAL, LOG_INFO, "Starting destruction of all boundary handlers. \n");
1533
1534 for (int i = 0; i < 6; i++) {
1535 BoundaryFaceConfig *face_cfg = &user->boundary_faces[i];
1536 const char *face_name = BCFaceToString(face_cfg->face_id);
1537
1538 // --- Step 1: Free the parameter linked list associated with this face ---
1539 if (face_cfg->params) {
1540 LOG_ALLOW(LOCAL, LOG_DEBUG, " Freeing parameter list for Face %d (%s). \n", i, face_name);
1541 FreeBC_ParamList(face_cfg->params);
1542 face_cfg->params = NULL; // Good practice to nullify dangling pointers
1543 }
1544
1545 // --- Step 2: Destroy the handler object itself ---
1546 if (face_cfg->handler) {
1547 const char *handler_name = BCHandlerTypeToString(face_cfg->handler->type);
1548 LOG_ALLOW(LOCAL, LOG_DEBUG, " Destroying handler '%s' on Face %d (%s).\n", handler_name, i, face_name);
1549
1550 // Call the handler's specific cleanup function first, if it exists.
1551 // This will free any memory stored in the handler's private `data` pointer.
1552 if (face_cfg->handler->Destroy) {
1553 ierr = face_cfg->handler->Destroy(face_cfg->handler); CHKERRQ(ierr);
1554 }
1555
1556 // Finally, free the generic BoundaryCondition object itself.
1557 ierr = PetscFree(face_cfg->handler); CHKERRQ(ierr);
1558 face_cfg->handler = NULL;
1559 }
1560 }
1561
1562 LOG_ALLOW(GLOBAL, LOG_INFO, "Destruction complete.\n");
1563 PetscFunctionReturn(0);
1564}
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 1572 of file Boundaries.c.

1573{
1574 PetscErrorCode ierr;
1575 DMDALocalInfo info = user->info;
1576 PetscInt xs = info.xs, xe = info.xs + info.xm;
1577 PetscInt ys = info.ys, ye = info.ys + info.ym;
1578 PetscInt zs = info.zs, ze = info.zs + info.zm;
1579 PetscInt mx = info.mx, my = info.my, mz = info.mz;
1580
1581 // --- Dispatcher to get DM, Vecs, and DoF for the specified field ---
1582 DM dm;
1583 Vec global_vec;
1584 Vec local_vec;
1585 PetscInt dof;
1586 // (This dispatcher is identical to your TransferPeriodicField function)
1587 if (strcmp(field_name, "Ucat") == 0) {
1588 dm = user->fda; global_vec = user->Ucat; local_vec = user->lUcat; dof = 3;
1589 } else if (strcmp(field_name, "P") == 0) {
1590 dm = user->da; global_vec = user->P; local_vec = user->lP; dof = 1;
1591 } else if (strcmp(field_name, "Nvert") == 0) {
1592 dm = user->da; global_vec = user->Nvert; local_vec = user->lNvert; dof = 1;
1593 } else {
1594 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown field name '%s'", field_name);
1595 }
1596
1597 PetscFunctionBeginUser;
1598
1599 // --- Execute the copy logic based on DoF and Direction ---
1600 if (dof == 1) { // --- Handle SCALAR fields (PetscReal) ---
1601 PetscReal ***g_array, ***l_array;
1602 ierr = DMDAVecGetArray(dm, global_vec, &g_array); CHKERRQ(ierr);
1603 ierr = DMDAVecGetArrayRead(dm, local_vec, (void*)&l_array); CHKERRQ(ierr); // Use Read for safety
1604
1605 switch (direction) {
1606 case 'i':
1607 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];
1608 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];
1609 break;
1610 case 'j':
1611 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];
1612 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];
1613 break;
1614 case 'k':
1615 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];
1616 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];
1617 break;
1618 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid direction '%c'", direction);
1619 }
1620 ierr = DMDAVecRestoreArray(dm, global_vec, &g_array); CHKERRQ(ierr);
1621 ierr = DMDAVecRestoreArrayRead(dm, local_vec, (void*)&l_array); CHKERRQ(ierr);
1622
1623 } else if (dof == 3) { // --- Handle VECTOR fields (Cmpnts) ---
1624 Cmpnts ***g_array, ***l_array;
1625 ierr = DMDAVecGetArray(dm, global_vec, &g_array); CHKERRQ(ierr);
1626 ierr = DMDAVecGetArrayRead(dm, local_vec, (void*)&l_array); CHKERRQ(ierr);
1627
1628 switch (direction) {
1629 case 'i':
1630 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];
1631 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];
1632 break;
1633 case 'j':
1634 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];
1635 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];
1636 break;
1637 case 'k':
1638 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];
1639 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];
1640 break;
1641 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid direction '%c'", direction);
1642 }
1643 ierr = DMDAVecRestoreArray(dm, global_vec, &g_array); CHKERRQ(ierr);
1644 ierr = DMDAVecRestoreArrayRead(dm, local_vec, (void*)&l_array); CHKERRQ(ierr);
1645 }
1646
1647 PetscFunctionReturn(0);
1648}
UserCtx::lNvert
Vec lNvert
Definition variables.h:837
UserCtx::Ucat
Vec Ucat
Definition variables.h:837
UserCtx::lUcat
Vec lUcat
Definition variables.h:837
UserCtx::Nvert
Vec Nvert
Definition variables.h:837
Here is the caller graph for this function:

◆ TransferPeriodicField()

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

Internal helper implementation: TransferPeriodicField().

(Orchestrator) Applies periodic transfer for one field across all i/j/k directions.

Local to this translation unit.

Definition at line 1656 of file Boundaries.c.

1657{
1658 PetscErrorCode ierr;
1659 DMDALocalInfo info = user->info;
1660 PetscInt xs = info.xs, xe = info.xs + info.xm;
1661 PetscInt ys = info.ys, ye = info.ys + info.ym;
1662 PetscInt zs = info.zs, ze = info.zs + info.zm;
1663 PetscInt mx = info.mx, my = info.my, mz = info.mz;
1664
1665 // --- Local variables to hold the specific details of the chosen field ---
1666 DM dm;
1667 Vec global_vec;
1668 Vec local_vec;
1669 PetscInt dof;
1670
1671 PetscFunctionBeginUser;
1672 PROFILE_FUNCTION_BEGIN;
1673
1674 // --- STEP 1: Dispatcher - Set the specific DM, Vecs, and dof based on field_name ---
1675 // Field-extension note: to add a new periodic field, add one case here and then
1676 // add one call in ApplyPeriodicBCs()/UpdatePeriodicCornerNodes() where appropriate.
1677 if (strcmp(field_name, "Ucat") == 0) {
1678 dm = user->fda;
1679 global_vec = user->Ucat;
1680 local_vec = user->lUcat;
1681 dof = 3;
1682 } else if (strcmp(field_name, "P") == 0) {
1683 dm = user->da;
1684 global_vec = user->P;
1685 local_vec = user->lP;
1686 dof = 1;
1687 } else if (strcmp(field_name, "Nvert") == 0) {
1688 dm = user->da;
1689 global_vec = user->Nvert;
1690 local_vec = user->lNvert;
1691 dof = 1;
1692 } else if (strcmp(field_name, "Eddy Viscosity") == 0) {
1693 dm = user->da;
1694 global_vec = user->Nu_t;
1695 local_vec = user->lNu_t;
1696 dof = 1;
1697 }
1698 /*
1699 // Example for future extension:
1700 else if (strcmp(field_name, "Temperature") == 0) {
1701 dm = user->da; // Assuming Temperature is scalar
1702 global_vec = user->T;
1703 local_vec = user->lT;
1704 dof = 1;
1705 }
1706 */
1707 else {
1708 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown field name '%s' in TransferPeriodicFieldByName.", field_name);
1709 }
1710
1711 LOG_ALLOW(GLOBAL,LOG_TRACE,"Periodic Transform being performed for field: %s with %d DoF.\n",field_name,dof);
1712 // --- STEP 2: Execute the copy logic using the dispatched variables ---
1713 if (dof == 1) { // --- Handle SCALAR fields (PetscReal) ---
1714 PetscReal ***g_array, ***l_array;
1715 ierr = DMDAVecGetArray(dm, global_vec, &g_array); CHKERRQ(ierr);
1716 ierr = DMDAVecGetArray(dm, local_vec, &l_array); CHKERRQ(ierr);
1717
1718 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];
1719 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];
1720 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];
1721 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];
1722 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];
1723 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];
1724
1725 ierr = DMDAVecRestoreArray(dm, global_vec, &g_array); CHKERRQ(ierr);
1726 ierr = DMDAVecRestoreArray(dm, local_vec, &l_array); CHKERRQ(ierr);
1727
1728 } else if (dof == 3) { // --- Handle VECTOR fields (Cmpnts) ---
1729 Cmpnts ***g_array, ***l_array;
1730 ierr = DMDAVecGetArray(dm, global_vec, &g_array); CHKERRQ(ierr);
1731 ierr = DMDAVecGetArray(dm, local_vec, &l_array); CHKERRQ(ierr);
1732
1733 LOG_ALLOW(GLOBAL,LOG_VERBOSE,"Array %s read successfully (Global and Local).\n",field_name);
1734
1735 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];
1736 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];
1737 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];
1738 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];
1739 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];
1740 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];
1741
1742 ierr = DMDAVecRestoreArray(dm, global_vec, &g_array); CHKERRQ(ierr);
1743 ierr = DMDAVecRestoreArray(dm, local_vec, &l_array); CHKERRQ(ierr);
1744 }
1745 else{
1746 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "This function only accepts Fields with 1 or 3 DoF.");
1747 }
1748
1749 PROFILE_FUNCTION_END;
1750 PetscFunctionReturn(0);
1751}
UserCtx::lNu_t
Vec lNu_t
Definition variables.h:865
UserCtx::Nu_t
Vec Nu_t
Definition variables.h:865
Here is the caller graph for this function:

◆ 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 1759 of file Boundaries.c.

1760{
1761 PetscErrorCode ierr;
1762 DMDALocalInfo info = user->info;
1763 PetscInt gxs = info.gxs, gxe = info.gxs + info.gxm;
1764 PetscInt gys = info.gys, gye = info.gys + info.gym;
1765 PetscInt gzs = info.gzs, gze = info.gzs + info.gzm;
1766 PetscInt mx = info.mx, my = info.my, mz = info.mz;
1767
1768 // --- Dispatcher to get the correct DM, Vec, and DoF for the specified field ---
1769 // Field-extension note: add one case here when a new field needs periodic
1770 // ghost-face transfer (typically for stencil-heavy operators).
1771 DM dm;
1772 Vec local_vec;
1773 PetscInt dof;
1774 // (This dispatcher contains all 17 potential fields)
1775 if (strcmp(field_name, "Ucont") == 0) { dm = user->fda; local_vec = user->lUcont; dof = 3; }
1776 else if (strcmp(field_name, "Csi") == 0) { dm = user->fda; local_vec = user->lCsi; dof = 3; }
1777 else if (strcmp(field_name, "Eta") == 0) { dm = user->fda; local_vec = user->lEta; dof = 3; }
1778 else if (strcmp(field_name, "Zet") == 0) { dm = user->fda; local_vec = user->lZet; dof = 3; }
1779 else if (strcmp(field_name, "ICsi") == 0) { dm = user->fda; local_vec = user->lICsi; dof = 3; }
1780 else if (strcmp(field_name, "IEta") == 0) { dm = user->fda; local_vec = user->lIEta; dof = 3; }
1781 else if (strcmp(field_name, "IZet") == 0) { dm = user->fda; local_vec = user->lIZet; dof = 3; }
1782 else if (strcmp(field_name, "JCsi") == 0) { dm = user->fda; local_vec = user->lJCsi; dof = 3; }
1783 else if (strcmp(field_name, "JEta") == 0) { dm = user->fda; local_vec = user->lJEta; dof = 3; }
1784 else if (strcmp(field_name, "JZet") == 0) { dm = user->fda; local_vec = user->lJZet; dof = 3; }
1785 else if (strcmp(field_name, "KCsi") == 0) { dm = user->fda; local_vec = user->lKCsi; dof = 3; }
1786 else if (strcmp(field_name, "KEta") == 0) { dm = user->fda; local_vec = user->lKEta; dof = 3; }
1787 else if (strcmp(field_name, "KZet") == 0) { dm = user->fda; local_vec = user->lKZet; dof = 3; }
1788 else if (strcmp(field_name, "Aj") == 0) { dm = user->da; local_vec = user->lAj; dof = 1; }
1789 else if (strcmp(field_name, "IAj") == 0) { dm = user->da; local_vec = user->lIAj; dof = 1; }
1790 else if (strcmp(field_name, "JAj") == 0) { dm = user->da; local_vec = user->lJAj; dof = 1; }
1791 else if (strcmp(field_name, "KAj") == 0) { dm = user->da; local_vec = user->lKAj; dof = 1; }
1792 else {
1793 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown field name '%s' in TransferPeriodicFaceField.", field_name);
1794 }
1795
1796 PetscFunctionBeginUser;
1797
1798 void *l_array_ptr;
1799 ierr = DMDAVecGetArray(dm, local_vec, &l_array_ptr); CHKERRQ(ierr);
1800
1801 // --- I-DIRECTION ---
1802 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC) {
1803 for (PetscInt k=gzs; k<gze; k++) for (PetscInt j=gys; j<gye; j++) {
1804 if (dof == 1) {
1805 PetscReal ***arr = (PetscReal***)l_array_ptr;
1806 arr[k][j][0] = arr[k][j][mx-2];
1807 arr[k][j][-1] = arr[k][j][mx-3];
1808 } else {
1809 Cmpnts ***arr = (Cmpnts***)l_array_ptr;
1810 arr[k][j][0] = arr[k][j][mx-2];
1811 arr[k][j][-1] = arr[k][j][mx-3];
1812 }
1813 }
1814 }
1815 if (user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC) {
1816 for (PetscInt k=gzs; k<gze; k++) for (PetscInt j=gys; j<gye; j++) {
1817 if (dof == 1) {
1818 PetscReal ***arr = (PetscReal***)l_array_ptr;
1819 arr[k][j][mx-1] = arr[k][j][1];
1820 arr[k][j][mx] = arr[k][j][2];
1821 } else {
1822 Cmpnts ***arr = (Cmpnts***)l_array_ptr;
1823 arr[k][j][mx-1] = arr[k][j][1];
1824 arr[k][j][mx] = arr[k][j][2];
1825 }
1826 }
1827 }
1828
1829 // --- J-DIRECTION ---
1830 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC) {
1831 for (PetscInt k=gzs; k<gze; k++) for (PetscInt i=gxs; i<gxe; i++) {
1832 if (dof == 1) {
1833 PetscReal ***arr = (PetscReal***)l_array_ptr;
1834 arr[k][0][i] = arr[k][my-2][i];
1835 arr[k][-1][i] = arr[k][my-3][i];
1836 } else {
1837 Cmpnts ***arr = (Cmpnts***)l_array_ptr;
1838 arr[k][0][i] = arr[k][my-2][i];
1839 arr[k][-1][i] = arr[k][my-3][i];
1840 }
1841 }
1842 }
1843 if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC) {
1844 for (PetscInt k=gzs; k<gze; k++) for (PetscInt i=gxs; i<gxe; i++) {
1845 if (dof == 1) {
1846 PetscReal ***arr = (PetscReal***)l_array_ptr;
1847 arr[k][my-1][i] = arr[k][1][i];
1848 arr[k][my][i] = arr[k][2][i];
1849 } else {
1850 Cmpnts ***arr = (Cmpnts***)l_array_ptr;
1851 arr[k][my-1][i] = arr[k][1][i];
1852 arr[k][my][i] = arr[k][2][i];
1853 }
1854 }
1855 }
1856
1857 // --- K-DIRECTION ---
1858 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC) {
1859 for (PetscInt j=gys; j<gye; j++) for (PetscInt i=gxs; i<gxe; i++) {
1860 if (dof == 1) {
1861 PetscReal ***arr = (PetscReal***)l_array_ptr;
1862 arr[0][j][i] = arr[mz-2][j][i];
1863 arr[-1][j][i] = arr[mz-3][j][i];
1864 } else {
1865 Cmpnts ***arr = (Cmpnts***)l_array_ptr;
1866 arr[0][j][i] = arr[mz-2][j][i];
1867 arr[-1][j][i] = arr[mz-3][j][i];
1868 }
1869 }
1870 }
1871 if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC) {
1872 for (PetscInt j=gys; j<gye; j++) for (PetscInt i=gxs; i<gxe; i++) {
1873 if (dof == 1) {
1874 PetscReal ***arr = (PetscReal***)l_array_ptr;
1875 arr[mz-1][j][i] = arr[1][j][i];
1876 arr[mz][j][i] = arr[2][j][i];
1877 } else {
1878 Cmpnts ***arr = (Cmpnts***)l_array_ptr;
1879 arr[mz-1][j][i] = arr[1][j][i];
1880 arr[mz][j][i] = arr[2][j][i];
1881 }
1882 }
1883 }
1884
1885 ierr = DMDAVecRestoreArray(dm, local_vec, &l_array_ptr); CHKERRQ(ierr);
1886 PetscFunctionReturn(0);
1887}
UserCtx::lIEta
Vec lIEta
Definition variables.h:860
UserCtx::lIZet
Vec lIZet
Definition variables.h:860
UserCtx::lZet
Vec lZet
Definition variables.h:858
UserCtx::lIAj
Vec lIAj
Definition variables.h:860
UserCtx::lKEta
Vec lKEta
Definition variables.h:862
UserCtx::lJCsi
Vec lJCsi
Definition variables.h:861
UserCtx::lKZet
Vec lKZet
Definition variables.h:862
UserCtx::lJEta
Vec lJEta
Definition variables.h:861
UserCtx::lCsi
Vec lCsi
Definition variables.h:858
UserCtx::lKCsi
Vec lKCsi
Definition variables.h:862
UserCtx::lJZet
Vec lJZet
Definition variables.h:861
UserCtx::lUcont
Vec lUcont
Definition variables.h:837
UserCtx::lAj
Vec lAj
Definition variables.h:858
UserCtx::lICsi
Vec lICsi
Definition variables.h:860
UserCtx::lEta
Vec lEta
Definition variables.h:858
UserCtx::lJAj
Vec lJAj
Definition variables.h:861
UserCtx::lKAj
Vec lKAj
Definition variables.h:862
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 1895 of file Boundaries.c.

1896{
1897 PetscErrorCode ierr;
1898 PetscFunctionBeginUser;
1899 PROFILE_FUNCTION_BEGIN;
1900
1901 const char* metric_fields[] = {
1902 "Csi", "Eta", "Zet", "ICsi", "JCsi", "KCsi", "IEta", "JEta", "KEta",
1903 "IZet", "JZet", "KZet", "Aj", "IAj", "JAj", "KAj"
1904 };
1905 PetscInt num_fields = sizeof(metric_fields) / sizeof(metric_fields[0]);
1906
1907 for (PetscInt i = 0; i < num_fields; i++) {
1908 ierr = TransferPeriodicFaceField(user, metric_fields[i]); CHKERRQ(ierr);
1909 LOG_ALLOW(GLOBAL,LOG_TRACE,"Periodic Transfer complete for %s.\n",metric_fields[i]);
1910 }
1911
1912 PROFILE_FUNCTION_END;
1913 PetscFunctionReturn(0);
1914}
TransferPeriodicFaceField
PetscErrorCode TransferPeriodicFaceField(UserCtx *user, const char *field_name)
Internal helper implementation: TransferPeriodicFaceField().
Definition Boundaries.c:1759
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 1922 of file Boundaries.c.

1923{
1924 PetscErrorCode ierr;
1925 PetscBool is_any_periodic = PETSC_FALSE;
1926
1927 PetscFunctionBeginUser;
1928
1929 PROFILE_FUNCTION_BEGIN;
1930
1931 for (int i = 0; i < 6; i++) {
1932 if (user->boundary_faces[i].mathematical_type == PERIODIC) {
1933 is_any_periodic = PETSC_TRUE;
1934 break;
1935 }
1936 }
1937
1938 if (!is_any_periodic) {
1939 LOG_ALLOW(GLOBAL,LOG_TRACE, "No periodic boundaries defined; skipping ApplyPeriodicBCs.\n");
1940 PROFILE_FUNCTION_END;
1941 PetscFunctionReturn(0);
1942 }
1943
1944 LOG_ALLOW(GLOBAL, LOG_TRACE, "Applying periodic boundary conditions for all fields.\n");
1945
1946 // STEP 1: Perform the collective communication for periodic core fields.
1947 ierr = UpdateLocalGhosts(user, "Ucat"); CHKERRQ(ierr);
1948 ierr = UpdateLocalGhosts(user, "P"); CHKERRQ(ierr);
1949 ierr = UpdateLocalGhosts(user, "Nvert"); CHKERRQ(ierr);
1950 ierr = UpdateLocalGhosts(user, "Ucont"); CHKERRQ(ierr);
1951 /* if (user->solve_temperature) { ierr = UpdateLocalGhosts(user, "Temperature"); CHKERRQ(ierr); } */
1952
1953 // STEP 2: Call the generic copy routine for each face-centered/cell-centered field.
1954 ierr = TransferPeriodicField(user, "Ucat"); CHKERRQ(ierr);
1955 ierr = TransferPeriodicField(user, "P"); CHKERRQ(ierr);
1956 ierr = TransferPeriodicField(user, "Nvert"); CHKERRQ(ierr);
1957
1958 // STEP 3: Update contravariant flux periodic ghosts and enforce strict seam consistency.
1959 // This keeps the staggered Ucont field robust for QUICK-like stencils and periodic flux closure.
1960 ierr = ApplyUcontPeriodicBCs(user); CHKERRQ(ierr);
1961 ierr = EnforceUcontPeriodicity(user); CHKERRQ(ierr);
1962 ierr = UpdateLocalGhosts(user, "Ucont"); CHKERRQ(ierr);
1963
1964 // FUTURE EXTENSION: Adding a new scalar field like Temperature is now trivial.
1965 /*
1966 if (user->solve_temperature) {
1967 ierr = TransferPeriodicField(user, "Temperature"); CHKERRQ(ierr);
1968 }
1969 */
1970
1971 PROFILE_FUNCTION_END;
1972 PetscFunctionReturn(0);
1973}
ApplyUcontPeriodicBCs
PetscErrorCode ApplyUcontPeriodicBCs(UserCtx *user)
Implementation of ApplyUcontPeriodicBCs().
Definition Boundaries.c:1983
EnforceUcontPeriodicity
PetscErrorCode EnforceUcontPeriodicity(UserCtx *user)
Internal helper implementation: EnforceUcontPeriodicity().
Definition Boundaries.c:2007
TransferPeriodicField
PetscErrorCode TransferPeriodicField(UserCtx *user, const char *field_name)
Internal helper implementation: TransferPeriodicField().
Definition Boundaries.c:1656
UpdateLocalGhosts
PetscErrorCode UpdateLocalGhosts(UserCtx *user, const char *fieldName)
Updates the local vector (including ghost points) from its corresponding global vector.
Definition setup.c:1361
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ApplyUcontPeriodicBCs()

PetscErrorCode ApplyUcontPeriodicBCs ( UserCtx *  user)

Implementation of ApplyUcontPeriodicBCs().

(Orchestrator) Updates the contravariant velocity field in the local ghost cell regions for periodic boundaries.

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

See also
ApplyUcontPeriodicBCs()

Definition at line 1983 of file Boundaries.c.

1984{
1985 PetscErrorCode ierr;
1986 PetscFunctionBeginUser;
1987 PROFILE_FUNCTION_BEGIN;
1988
1989 // Update local periodic ghost faces (two-cells deep) for Ucont.
1990 ierr = TransferPeriodicFaceField(user, "Ucont"); CHKERRQ(ierr);
1991
1992 // Commit periodic-face updates back to global Ucont so follow-up routines
1993 // (including EnforceUcontPeriodicity) can safely refresh from global state.
1994 ierr = DMLocalToGlobalBegin(user->fda, user->lUcont, INSERT_VALUES, user->Ucont); CHKERRQ(ierr);
1995 ierr = DMLocalToGlobalEnd(user->fda, user->lUcont, INSERT_VALUES, user->Ucont); CHKERRQ(ierr);
1996
1997 PROFILE_FUNCTION_END;
1998 PetscFunctionReturn(0);
1999}
UserCtx::Ucont
Vec Ucont
Definition variables.h:837
Here is the call graph for this function:
Here is the caller graph for this function:

◆ EnforceUcontPeriodicity()

PetscErrorCode EnforceUcontPeriodicity ( UserCtx *  user)

Internal helper implementation: EnforceUcontPeriodicity().

Enforces strict periodicity on the interior contravariant velocity field.

Local to this translation unit.

Definition at line 2007 of file Boundaries.c.

2008{
2009 PetscErrorCode ierr;
2010 DMDALocalInfo info = user->info;
2011 PetscInt gxs = info.gxs, gxe = info.gxs + info.gxm;
2012 PetscInt gys = info.gys, gye = info.gys + info.gym;
2013 PetscInt gzs = info.gzs, gze = info.gzs + info.gzm;
2014 PetscInt mx = info.mx, my = info.my, mz = info.mz;
2015 Cmpnts ***ucont;
2016
2017 PetscFunctionBeginUser;
2018 PROFILE_FUNCTION_BEGIN;
2019
2020 // STEP 1: Ensure local ghost cells are up-to-date with current global state.
2021 ierr = DMGlobalToLocalBegin(user->fda, user->Ucont, INSERT_VALUES, user->lUcont); CHKERRQ(ierr);
2022 ierr = DMGlobalToLocalEnd(user->fda, user->Ucont, INSERT_VALUES, user->lUcont); CHKERRQ(ierr);
2023
2024 ierr = DMDAVecGetArray(user->fda, user->lUcont, &ucont); CHKERRQ(ierr);
2025
2026 // STEP 2: Perform the component-wise copy from ghost cells to the last interior faces.
2027 if (user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC) {
2028 for (PetscInt k=gzs; k<gze; k++) for (PetscInt j=gys; j<gye; j++) {
2029 ucont[k][j][mx-2].x = ucont[k][j][mx].x; // Note: ucont[mx] is ghost, gets value from ucont[0]
2030 }
2031 }
2032 if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC) {
2033 for (PetscInt k=gzs; k<gze; k++) for (PetscInt i=gxs; i<gxe; i++) {
2034 ucont[k][my-2][i].y = ucont[k][my][i].y; // Note: ucont[my] is ghost, gets value from ucont[0]
2035 }
2036 }
2037 if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC) {
2038 for (PetscInt j=gys; j<gye; j++) for (PetscInt i=gxs; i<gxe; i++) {
2039 ucont[mz-2][j][i].z = ucont[mz][j][i].z; // Note: ucont[mz] is ghost, gets value from ucont[0]
2040 }
2041 }
2042
2043 ierr = DMDAVecRestoreArray(user->fda, user->lUcont, &ucont); CHKERRQ(ierr);
2044
2045 // STEP 3: Communicate the changes made to the interior back to the global vector.
2046 ierr = DMLocalToGlobalBegin(user->fda, user->lUcont, INSERT_VALUES, user->Ucont); CHKERRQ(ierr);
2047 ierr = DMLocalToGlobalEnd(user->fda, user->lUcont, INSERT_VALUES, user->Ucont); CHKERRQ(ierr);
2048
2049 PROFILE_FUNCTION_END;
2050 PetscFunctionReturn(0);
2051}
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 2059 of file Boundaries.c.

2060{
2061 PetscErrorCode ierr;
2062 DM fda = user->fda;
2063 DMDALocalInfo info = user->info;
2064 PetscInt xs = info.xs, xe = info.xs + info.xm;
2065 PetscInt ys = info.ys, ye = info.ys + info.ym;
2066 PetscInt zs = info.zs, ze = info.zs + info.zm;
2067 PetscInt mx = info.mx, my = info.my, mz = info.mz;
2068
2069 // --- Calculate shrunken loop ranges to avoid edges and corners ---
2070 PetscInt lxs = xs, lxe = xe;
2071 PetscInt lys = ys, lye = ye;
2072 PetscInt lzs = zs, lze = ze;
2073
2074 if (xs == 0) lxs = xs + 1;
2075 if (ys == 0) lys = ys + 1;
2076 if (zs == 0) lzs = zs + 1;
2077
2078 if (xe == mx) lxe = xe - 1;
2079 if (ye == my) lye = ye - 1;
2080 if (ze == mz) lze = ze - 1;
2081
2082 Cmpnts ***ucat, ***ubcs;
2083 PetscFunctionBeginUser;
2084
2085 ierr = DMDAVecGetArray(fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
2086 ierr = DMDAVecGetArray(fda, user->Ucat, &ucat); CHKERRQ(ierr);
2087
2088 // -X Face
2089 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC && xs == 0) {
2090 for (PetscInt k = lzs; k < lze; k++) for (PetscInt j = lys; j < lye; j++) {
2091 ucat[k][j][xs].x = 2.0 * ubcs[k][j][xs].x - ucat[k][j][xs + 1].x;
2092 ucat[k][j][xs].y = 2.0 * ubcs[k][j][xs].y - ucat[k][j][xs + 1].y;
2093 ucat[k][j][xs].z = 2.0 * ubcs[k][j][xs].z - ucat[k][j][xs + 1].z;
2094 }
2095 }
2096 // +X Face
2097 if (user->boundary_faces[BC_FACE_POS_X].mathematical_type != PERIODIC && xe == mx) {
2098 for (PetscInt k = lzs; k < lze; k++) for (PetscInt j = lys; j < lye; j++) {
2099 ucat[k][j][xe-1].x = 2.0 * ubcs[k][j][xe-1].x - ucat[k][j][xe - 2].x;
2100 ucat[k][j][xe-1].y = 2.0 * ubcs[k][j][xe-1].y - ucat[k][j][xe - 2].y;
2101 ucat[k][j][xe-1].z = 2.0 * ubcs[k][j][xe-1].z - ucat[k][j][xe - 2].z;
2102 }
2103 }
2104
2105 // -Y Face
2106 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC && ys == 0) {
2107 for (PetscInt k = lzs; k < lze; k++) for (PetscInt i = lxs; i < lxe; i++) {
2108 ucat[k][ys][i].x = 2.0 * ubcs[k][ys][i].x - ucat[k][ys + 1][i].x;
2109 ucat[k][ys][i].y = 2.0 * ubcs[k][ys][i].y - ucat[k][ys + 1][i].y;
2110 ucat[k][ys][i].z = 2.0 * ubcs[k][ys][i].z - ucat[k][ys + 1][i].z;
2111 }
2112 }
2113 // +Y Face
2114 if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type != PERIODIC && ye == my) {
2115 for (PetscInt k = lzs; k < lze; k++) for (PetscInt i = lxs; i < lxe; i++) {
2116 ucat[k][ye-1][i].x = 2.0 * ubcs[k][ye-1][i].x - ucat[k][ye-2][i].x;
2117 ucat[k][ye-1][i].y = 2.0 * ubcs[k][ye-1][i].y - ucat[k][ye-2][i].y;
2118 ucat[k][ye-1][i].z = 2.0 * ubcs[k][ye-1][i].z - ucat[k][ye-2][i].z;
2119 }
2120 }
2121
2122 // -Z Face
2123 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC && zs == 0) {
2124 for (PetscInt j = lys; j < lye; j++) for (PetscInt i = lxs; i < lxe; i++) {
2125 ucat[zs][j][i].x = 2.0 * ubcs[zs][j][i].x - ucat[zs + 1][j][i].x;
2126 ucat[zs][j][i].y = 2.0 * ubcs[zs][j][i].y - ucat[zs + 1][j][i].y;
2127 ucat[zs][j][i].z = 2.0 * ubcs[zs][j][i].z - ucat[zs + 1][j][i].z;
2128 }
2129 }
2130 // +Z Face
2131 if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type != PERIODIC && ze == mz) {
2132 for (PetscInt j = lys; j < lye; j++) for (PetscInt i = lxs; i < lxe; i++) {
2133 ucat[ze-1][j][i].x = 2.0 * ubcs[ze-1][j][i].x - ucat[ze-2][j][i].x;
2134 ucat[ze-1][j][i].y = 2.0 * ubcs[ze-1][j][i].y - ucat[ze-2][j][i].y;
2135 ucat[ze-1][j][i].z = 2.0 * ubcs[ze-1][j][i].z - ucat[ze-2][j][i].z;
2136 }
2137 }
2138
2139 ierr = DMDAVecRestoreArray(fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
2140 ierr = DMDAVecRestoreArray(fda, user->Ucat, &ucat); CHKERRQ(ierr);
2141
2142 PetscFunctionReturn(0);
2143}
BCS::Ubcs
Vec Ubcs
Physical Cartesian velocity at boundary faces. Full 3D array but only boundary-face entries are meani...
Definition variables.h:121
UserCtx::Bcs
BCS Bcs
Definition variables.h:832
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 2151 of file Boundaries.c.

2152{
2153 PetscErrorCode ierr;
2154 DM da = user->da, fda = user->fda;
2155 DMDALocalInfo info = user->info;
2156 PetscInt xs = info.xs, xe = info.xs + info.xm;
2157 PetscInt ys = info.ys, ye = info.ys + info.ym;
2158 PetscInt zs = info.zs, ze = info.zs + info.zm;
2159 PetscInt mx = info.mx, my = info.my, mz = info.mz;
2160
2161 Cmpnts ***ucat;
2162 PetscReal ***p;
2163
2164 PetscFunctionBeginUser;
2165
2166 ierr = DMDAVecGetArray(fda, user->Ucat, &ucat); CHKERRQ(ierr);
2167 ierr = DMDAVecGetArray(da, user->P, &p); CHKERRQ(ierr);
2168
2169 // --- Update Edges and Corners by Averaging ---
2170 // The order of these blocks ensures that corners (where 3 faces meet) are
2171 // computed using data from edges (where 2 faces meet), which are computed first.
2172// Edges connected to the -Z face (k=zs)
2173 if (zs == 0) {
2174 if (xs == 0) {
2175 for (PetscInt j = ys; j < ye; j++) {
2176 p[zs][j][xs] = 0.5 * (p[zs+1][j][xs] + p[zs][j][xs+1]);
2177 ucat[zs][j][xs].x = 0.5 * (ucat[zs+1][j][xs].x + ucat[zs][j][xs+1].x);
2178 ucat[zs][j][xs].y = 0.5 * (ucat[zs+1][j][xs].y + ucat[zs][j][xs+1].y);
2179 ucat[zs][j][xs].z = 0.5 * (ucat[zs+1][j][xs].z + ucat[zs][j][xs+1].z);
2180 }
2181 }
2182 if (xe == mx) {
2183 for (PetscInt j = ys; j < ye; j++) {
2184 p[zs][j][mx-1] = 0.5 * (p[zs+1][j][mx-1] + p[zs][j][mx-2]);
2185 ucat[zs][j][mx-1].x = 0.5 * (ucat[zs+1][j][mx-1].x + ucat[zs][j][mx-2].x);
2186 ucat[zs][j][mx-1].y = 0.5 * (ucat[zs+1][j][mx-1].y + ucat[zs][j][mx-2].y);
2187 ucat[zs][j][mx-1].z = 0.5 * (ucat[zs+1][j][mx-1].z + ucat[zs][j][mx-2].z);
2188 }
2189 }
2190 if (ys == 0) {
2191 for (PetscInt i = xs; i < xe; i++) {
2192 p[zs][ys][i] = 0.5 * (p[zs+1][ys][i] + p[zs][ys+1][i]);
2193 ucat[zs][ys][i].x = 0.5 * (ucat[zs+1][ys][i].x + ucat[zs][ys+1][i].x);
2194 ucat[zs][ys][i].y = 0.5 * (ucat[zs+1][ys][i].y + ucat[zs][ys+1][i].y);
2195 ucat[zs][ys][i].z = 0.5 * (ucat[zs+1][ys][i].z + ucat[zs][ys+1][i].z);
2196 }
2197 }
2198 if (ye == my) {
2199 for (PetscInt i = xs; i < xe; i++) {
2200 p[zs][my-1][i] = 0.5 * (p[zs+1][my-1][i] + p[zs][my-2][i]);
2201 ucat[zs][my-1][i].x = 0.5 * (ucat[zs+1][my-1][i].x + ucat[zs][my-2][i].x);
2202 ucat[zs][my-1][i].y = 0.5 * (ucat[zs+1][my-1][i].y + ucat[zs][my-2][i].y);
2203 ucat[zs][my-1][i].z = 0.5 * (ucat[zs+1][my-1][i].z + ucat[zs][my-2][i].z);
2204 }
2205 }
2206 }
2207
2208 // Edges connected to the +Z face (k=ze-1)
2209 if (ze == mz) {
2210 if (xs == 0) {
2211 for (PetscInt j = ys; j < ye; j++) {
2212 p[mz-1][j][xs] = 0.5 * (p[mz-2][j][xs] + p[mz-1][j][xs+1]);
2213 ucat[mz-1][j][xs].x = 0.5 * (ucat[mz-2][j][xs].x + ucat[mz-1][j][xs+1].x);
2214 ucat[mz-1][j][xs].y = 0.5 * (ucat[mz-2][j][xs].y + ucat[mz-1][j][xs+1].y);
2215 ucat[mz-1][j][xs].z = 0.5 * (ucat[mz-2][j][xs].z + ucat[mz-1][j][xs+1].z);
2216 }
2217 }
2218 if (xe == mx) {
2219 for (PetscInt j = ys; j < ye; j++) {
2220 p[mz-1][j][mx-1] = 0.5 * (p[mz-2][j][mx-1] + p[mz-1][j][mx-2]);
2221 ucat[mz-1][j][mx-1].x = 0.5 * (ucat[mz-2][j][mx-1].x + ucat[mz-1][j][mx-2].x);
2222 ucat[mz-1][j][mx-1].y = 0.5 * (ucat[mz-2][j][mx-1].y + ucat[mz-1][j][mx-2].y);
2223 ucat[mz-1][j][mx-1].z = 0.5 * (ucat[mz-2][j][mx-1].z + ucat[mz-1][j][mx-2].z);
2224 }
2225 }
2226 if (ys == 0) {
2227 for (PetscInt i = xs; i < xe; i++) {
2228 p[mz-1][ys][i] = 0.5 * (p[mz-2][ys][i] + p[mz-1][ys+1][i]);
2229 ucat[mz-1][ys][i].x = 0.5 * (ucat[mz-2][ys][i].x + ucat[mz-1][ys+1][i].x);
2230 ucat[mz-1][ys][i].y = 0.5 * (ucat[mz-2][ys][i].y + ucat[mz-1][ys+1][i].y);
2231 ucat[mz-1][ys][i].z = 0.5 * (ucat[mz-2][ys][i].z + ucat[mz-1][ys+1][i].z);
2232 }
2233 }
2234 if (ye == my) {
2235 for (PetscInt i = xs; i < xe; i++) {
2236 p[mz-1][my-1][i] = 0.5 * (p[mz-2][my-1][i] + p[mz-1][my-2][i]);
2237 ucat[mz-1][my-1][i].x = 0.5 * (ucat[mz-2][my-1][i].x + ucat[mz-1][my-2][i].x);
2238 ucat[mz-1][my-1][i].y = 0.5 * (ucat[mz-2][my-1][i].y + ucat[mz-1][my-2][i].y);
2239 ucat[mz-1][my-1][i].z = 0.5 * (ucat[mz-2][my-1][i].z + ucat[mz-1][my-2][i].z);
2240 }
2241 }
2242 }
2243
2244 // Remaining edges on the XY plane (that are not on Z faces)
2245 if (ys == 0) {
2246 if (xs == 0) {
2247 for (PetscInt k = zs; k < ze; k++) {
2248 p[k][ys][xs] = 0.5 * (p[k][ys+1][xs] + p[k][ys][xs+1]);
2249 ucat[k][ys][xs].x = 0.5 * (ucat[k][ys+1][xs].x + ucat[k][ys][xs+1].x);
2250 ucat[k][ys][xs].y = 0.5 * (ucat[k][ys+1][xs].y + ucat[k][ys][xs+1].y);
2251 ucat[k][ys][xs].z = 0.5 * (ucat[k][ys+1][xs].z + ucat[k][ys][xs+1].z);
2252 }
2253 }
2254 if (xe == mx) {
2255 for (PetscInt k = zs; k < ze; k++) {
2256 p[k][ys][mx-1] = 0.5 * (p[k][ys+1][mx-1] + p[k][ys][mx-2]);
2257 ucat[k][ys][mx-1].x = 0.5 * (ucat[k][ys+1][mx-1].x + ucat[k][ys][mx-2].x);
2258 ucat[k][ys][mx-1].y = 0.5 * (ucat[k][ys+1][mx-1].y + ucat[k][ys][mx-2].y);
2259 ucat[k][ys][mx-1].z = 0.5 * (ucat[k][ys+1][mx-1].z + ucat[k][ys][mx-2].z);
2260 }
2261 }
2262 }
2263
2264 if (ye == my) {
2265 if (xs == 0) {
2266 for (PetscInt k = zs; k < ze; k++) {
2267 p[k][my-1][xs] = 0.5 * (p[k][my-2][xs] + p[k][my-1][xs+1]);
2268 ucat[k][my-1][xs].x = 0.5 * (ucat[k][my-2][xs].x + ucat[k][my-1][xs+1].x);
2269 ucat[k][my-1][xs].y = 0.5 * (ucat[k][my-2][xs].y + ucat[k][my-1][xs+1].y);
2270 ucat[k][my-1][xs].z = 0.5 * (ucat[k][my-2][xs].z + ucat[k][my-1][xs+1].z);
2271 }
2272 }
2273 if (xe == mx) {
2274 for (PetscInt k = zs; k < ze; k++) {
2275 p[k][my-1][mx-1] = 0.5 * (p[k][my-2][mx-1] + p[k][my-1][mx-2]);
2276 ucat[k][my-1][mx-1].x = 0.5 * (ucat[k][my-2][mx-1].x + ucat[k][my-1][mx-2].x);
2277 ucat[k][my-1][mx-1].y = 0.5 * (ucat[k][my-2][mx-1].y + ucat[k][my-1][mx-2].y);
2278 ucat[k][my-1][mx-1].z = 0.5 * (ucat[k][my-2][mx-1].z + ucat[k][my-1][mx-2].z);
2279 }
2280 }
2281 }
2282
2283 ierr = DMDAVecRestoreArray(fda, user->Ucat, &ucat); CHKERRQ(ierr);
2284 ierr = DMDAVecRestoreArray(da, user->P, &p); CHKERRQ(ierr);
2285
2286 PetscFunctionReturn(0);
2287}
Here is the caller graph for this function:

◆ UpdatePeriodicCornerNodes()

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

Internal helper implementation: UpdatePeriodicCornerNodes().

(Orchestrator) Performs a sequential, deterministic periodic update for a list of fields.

Local to this translation unit.

Definition at line 2295 of file Boundaries.c.

2296{
2297 PetscErrorCode ierr;
2298 PetscFunctionBeginUser;
2299
2300 if (num_fields == 0) PetscFunctionReturn(0);
2301
2302 // --- I-DIRECTION ---
2303 for (PetscInt i = 0; i < num_fields; i++) {
2304 ierr = TransferPeriodicFieldByDirection(user, field_names[i], 'i'); CHKERRQ(ierr);
2305 }
2306 // --- SYNC ---
2307 for (PetscInt i = 0; i < num_fields; i++) {
2308 ierr = UpdateLocalGhosts(user, field_names[i]); CHKERRQ(ierr);
2309 }
2310
2311 // --- J-DIRECTION ---
2312 for (PetscInt i = 0; i < num_fields; i++) {
2313 ierr = TransferPeriodicFieldByDirection(user, field_names[i], 'j'); CHKERRQ(ierr);
2314 }
2315 // --- SYNC ---
2316 for (PetscInt i = 0; i < num_fields; i++) {
2317 ierr = UpdateLocalGhosts(user, field_names[i]); CHKERRQ(ierr);
2318 }
2319
2320 // --- K-DIRECTION ---
2321 for (PetscInt i = 0; i < num_fields; i++) {
2322 ierr = TransferPeriodicFieldByDirection(user, field_names[i], 'k'); CHKERRQ(ierr);
2323 }
2324 // --- FINAL SYNC ---
2325 for (PetscInt i = 0; i < num_fields; i++) {
2326 ierr = UpdateLocalGhosts(user, field_names[i]); CHKERRQ(ierr);
2327 }
2328
2329 PetscFunctionReturn(0);
2330}
TransferPeriodicFieldByDirection
PetscErrorCode TransferPeriodicFieldByDirection(UserCtx *user, const char *field_name, char direction)
Internal helper implementation: TransferPeriodicFieldByDirection().
Definition Boundaries.c:1572
Here is the call graph for this function:
Here is the caller 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 2338 of file Boundaries.c.

2339{
2340 PetscErrorCode ierr;
2341 SimCtx *simCtx = user->simCtx;
2342 DMDALocalInfo *info = &user->info;
2343
2344 PetscFunctionBeginUser;
2345
2346 // =========================================================================
2347 // STEP 0: Early exit if wall functions are disabled
2348 // =========================================================================
2349 if (!simCtx->wallfunction) {
2350 PetscFunctionReturn(0);
2351 }
2352
2353 LOG_ALLOW(LOCAL, LOG_DEBUG, "Processing wall function boundaries.\n");
2354
2355 // =========================================================================
2356 // STEP 1: Get read/write access to all necessary field arrays
2357 // =========================================================================
2358 Cmpnts ***velocity_cartesian; // Cartesian velocity (modified)
2359 Cmpnts ***velocity_contravariant; // Contravariant velocity (set to zero at walls)
2360 Cmpnts ***velocity_boundary; // Boundary condition velocity (kept at zero)
2361 Cmpnts ***csi, ***eta, ***zet; // Metric tensor components (face normals)
2362 PetscReal ***node_vertex_flag; // Fluid/solid indicator (0=fluid, 1=solid)
2363 PetscReal ***cell_jacobian; // Grid Jacobian (1/volume)
2364 PetscReal ***friction_velocity; // u_tau (friction velocity field)
2365
2366 ierr = DMDAVecGetArray(user->fda, user->Ucat, &velocity_cartesian); CHKERRQ(ierr);
2367 ierr = DMDAVecGetArray(user->fda, user->Ucont, &velocity_contravariant); CHKERRQ(ierr);
2368 ierr = DMDAVecGetArray(user->fda, user->Bcs.Ubcs, &velocity_boundary); CHKERRQ(ierr);
2369 ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
2370 ierr = DMDAVecGetArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
2371 ierr = DMDAVecGetArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
2372 ierr = DMDAVecGetArrayRead(user->da, user->lNvert, (const PetscReal***)&node_vertex_flag); CHKERRQ(ierr);
2373 ierr = DMDAVecGetArrayRead(user->da, user->lAj, (const PetscReal***)&cell_jacobian); CHKERRQ(ierr);
2374 ierr = DMDAVecGetArray(user->da, user->lFriction_Velocity, &friction_velocity); CHKERRQ(ierr);
2375
2376 // =========================================================================
2377 // STEP 2: Define loop bounds (owned portion of the grid for this MPI rank)
2378 // =========================================================================
2379 PetscInt grid_start_i = info->xs, grid_end_i = info->xs + info->xm;
2380 PetscInt grid_start_j = info->ys, grid_end_j = info->ys + info->ym;
2381 PetscInt grid_start_k = info->zs, grid_end_k = info->zs + info->zm;
2382 PetscInt grid_size_i = info->mx, grid_size_j = info->my, grid_size_k = info->mz;
2383
2384 // Shrunken loop bounds: exclude domain edges and corners to avoid double-counting
2385 PetscInt loop_start_i = grid_start_i, loop_end_i = grid_end_i;
2386 PetscInt loop_start_j = grid_start_j, loop_end_j = grid_end_j;
2387 PetscInt loop_start_k = grid_start_k, loop_end_k = grid_end_k;
2388
2389 if (grid_start_i == 0) loop_start_i = grid_start_i + 1;
2390 if (grid_end_i == grid_size_i) loop_end_i = grid_end_i - 1;
2391 if (grid_start_j == 0) loop_start_j = grid_start_j + 1;
2392 if (grid_end_j == grid_size_j) loop_end_j = grid_end_j - 1;
2393 if (grid_start_k == 0) loop_start_k = grid_start_k + 1;
2394 if (grid_end_k == grid_size_k) loop_end_k = grid_end_k - 1;
2395
2396 // Wall roughness parameter (currently smooth wall)
2397 const PetscReal wall_roughness_height = 1.e-16;
2398
2399 // =========================================================================
2400 // STEP 3: Process each of the 6 domain faces
2401 // =========================================================================
2402 for (int face_index = 0; face_index < 6; face_index++) {
2403 BCFace current_face_id = (BCFace)face_index;
2404 BoundaryFaceConfig *face_config = &user->boundary_faces[current_face_id];
2405
2406 // Only process faces that are mathematical walls (applies to no-slip, moving, slip, etc.)
2407 if (face_config->mathematical_type != WALL) {
2408 continue;
2409 }
2410
2411 // Check if this MPI rank owns part of this face
2412 PetscBool rank_owns_this_face;
2413 ierr = CanRankServiceFace(info, user->IM, user->JM, user->KM,
2414 current_face_id, &rank_owns_this_face); CHKERRQ(ierr);
2415
2416 if (!rank_owns_this_face) {
2417 continue;
2418 }
2419
2420 LOG_ALLOW(LOCAL, LOG_TRACE, "Processing Face %d (%s)\n",
2421 current_face_id, BCFaceToString(current_face_id));
2422
2423 // =====================================================================
2424 // Process each face with appropriate indexing
2425 // =====================================================================
2426 switch(current_face_id) {
2427
2428 // =================================================================
2429 // NEGATIVE X FACE (i = 0, normal points in +X direction)
2430 // =================================================================
2431 case BC_FACE_NEG_X: {
2432 if (grid_start_i == 0) {
2433 const PetscInt ghost_cell_index = grid_start_i;
2434 const PetscInt first_interior_cell = grid_start_i + 1;
2435 const PetscInt second_interior_cell = grid_start_i + 2;
2436
2437 for (PetscInt k = loop_start_k; k < loop_end_k; k++) {
2438 for (PetscInt j = loop_start_j; j < loop_end_j; j++) {
2439
2440 // Skip if this is a solid cell (embedded boundary)
2441 if (node_vertex_flag[k][j][first_interior_cell] < 0.1) {
2442
2443 // Calculate face area from contravariant metric tensor
2444 PetscReal face_area = sqrt(
2445 csi[k][j][ghost_cell_index].x * csi[k][j][ghost_cell_index].x +
2446 csi[k][j][ghost_cell_index].y * csi[k][j][ghost_cell_index].y +
2447 csi[k][j][ghost_cell_index].z * csi[k][j][ghost_cell_index].z
2448 );
2449
2450 // Compute wall-normal distances using cell Jacobians
2451 // sb = distance from wall to first interior cell center
2452 // sc = distance from wall to second interior cell center
2453 PetscReal distance_to_first_cell = 0.5 / cell_jacobian[k][j][first_interior_cell] / face_area;
2454 PetscReal distance_to_second_cell = 2.0 * distance_to_first_cell +
2455 0.5 / cell_jacobian[k][j][second_interior_cell] / face_area;
2456
2457 // Compute unit normal vector pointing INTO the domain
2458 PetscReal wall_normal[3];
2459 wall_normal[0] = csi[k][j][ghost_cell_index].x / face_area;
2460 wall_normal[1] = csi[k][j][ghost_cell_index].y / face_area;
2461 wall_normal[2] = csi[k][j][ghost_cell_index].z / face_area;
2462
2463 // Define velocities for wall function calculation
2464 Cmpnts wall_velocity; // Ua = velocity at wall (zero for stationary wall)
2465 Cmpnts reference_velocity; // Uc = velocity at second interior cell
2466
2467 wall_velocity.x = wall_velocity.y = wall_velocity.z = 0.0;
2468 reference_velocity = velocity_cartesian[k][j][second_interior_cell];
2469
2470 // Step 1: Linear interpolation (provides initial guess)
2471 noslip(user, distance_to_second_cell, distance_to_first_cell,
2472 wall_velocity, reference_velocity,
2473 &velocity_cartesian[k][j][first_interior_cell],
2474 wall_normal[0], wall_normal[1], wall_normal[2]);
2475
2476 // Step 2: Apply log-law correction (improves near-wall velocity)
2477 wall_function_loglaw(user, wall_roughness_height,
2478 distance_to_second_cell, distance_to_first_cell,
2479 wall_velocity, reference_velocity,
2480 &velocity_cartesian[k][j][first_interior_cell],
2481 &friction_velocity[k][j][first_interior_cell],
2482 wall_normal[0], wall_normal[1], wall_normal[2]);
2483
2484 // Ensure ghost cell BC remains zero (required for proper extrapolation)
2485 velocity_boundary[k][j][ghost_cell_index].x = 0.0;
2486 velocity_boundary[k][j][ghost_cell_index].y = 0.0;
2487 velocity_boundary[k][j][ghost_cell_index].z = 0.0;
2488 velocity_contravariant[k][j][ghost_cell_index].x = 0.0;
2489 }
2490 }
2491 }
2492 }
2493 } break;
2494
2495 // =================================================================
2496 // POSITIVE X FACE (i = mx-1, normal points in -X direction)
2497 // =================================================================
2498 case BC_FACE_POS_X: {
2499 if (grid_end_i == grid_size_i) {
2500 const PetscInt ghost_cell_index = grid_end_i - 1;
2501 const PetscInt first_interior_cell = grid_end_i - 2;
2502 const PetscInt second_interior_cell = grid_end_i - 3;
2503
2504 for (PetscInt k = loop_start_k; k < loop_end_k; k++) {
2505 for (PetscInt j = loop_start_j; j < loop_end_j; j++) {
2506
2507 if (node_vertex_flag[k][j][first_interior_cell] < 0.1) {
2508
2509 PetscReal face_area = sqrt(
2510 csi[k][j][first_interior_cell].x * csi[k][j][first_interior_cell].x +
2511 csi[k][j][first_interior_cell].y * csi[k][j][first_interior_cell].y +
2512 csi[k][j][first_interior_cell].z * csi[k][j][first_interior_cell].z
2513 );
2514
2515 PetscReal distance_to_first_cell = 0.5 / cell_jacobian[k][j][first_interior_cell] / face_area;
2516 PetscReal distance_to_second_cell = 2.0 * distance_to_first_cell +
2517 0.5 / cell_jacobian[k][j][second_interior_cell] / face_area;
2518
2519 // Note: Normal flipped for +X face to point INTO domain
2520 PetscReal wall_normal[3];
2521 wall_normal[0] = -csi[k][j][first_interior_cell].x / face_area;
2522 wall_normal[1] = -csi[k][j][first_interior_cell].y / face_area;
2523 wall_normal[2] = -csi[k][j][first_interior_cell].z / face_area;
2524
2525 Cmpnts wall_velocity, reference_velocity;
2526 wall_velocity.x = wall_velocity.y = wall_velocity.z = 0.0;
2527 reference_velocity = velocity_cartesian[k][j][second_interior_cell];
2528
2529 noslip(user, distance_to_second_cell, distance_to_first_cell,
2530 wall_velocity, reference_velocity,
2531 &velocity_cartesian[k][j][first_interior_cell],
2532 wall_normal[0], wall_normal[1], wall_normal[2]);
2533
2534 wall_function_loglaw(user, wall_roughness_height,
2535 distance_to_second_cell, distance_to_first_cell,
2536 wall_velocity, reference_velocity,
2537 &velocity_cartesian[k][j][first_interior_cell],
2538 &friction_velocity[k][j][first_interior_cell],
2539 wall_normal[0], wall_normal[1], wall_normal[2]);
2540
2541 velocity_boundary[k][j][ghost_cell_index].x = 0.0;
2542 velocity_boundary[k][j][ghost_cell_index].y = 0.0;
2543 velocity_boundary[k][j][ghost_cell_index].z = 0.0;
2544 velocity_contravariant[k][j][first_interior_cell].x = 0.0;
2545 }
2546 }
2547 }
2548 }
2549 } break;
2550
2551 // =================================================================
2552 // NEGATIVE Y FACE (j = 0, normal points in +Y direction)
2553 // =================================================================
2554 case BC_FACE_NEG_Y: {
2555 if (grid_start_j == 0) {
2556 const PetscInt ghost_cell_index = grid_start_j;
2557 const PetscInt first_interior_cell = grid_start_j + 1;
2558 const PetscInt second_interior_cell = grid_start_j + 2;
2559
2560 for (PetscInt k = loop_start_k; k < loop_end_k; k++) {
2561 for (PetscInt i = loop_start_i; i < loop_end_i; i++) {
2562
2563 if (node_vertex_flag[k][first_interior_cell][i] < 0.1) {
2564
2565 PetscReal face_area = sqrt(
2566 eta[k][ghost_cell_index][i].x * eta[k][ghost_cell_index][i].x +
2567 eta[k][ghost_cell_index][i].y * eta[k][ghost_cell_index][i].y +
2568 eta[k][ghost_cell_index][i].z * eta[k][ghost_cell_index][i].z
2569 );
2570
2571 PetscReal distance_to_first_cell = 0.5 / cell_jacobian[k][first_interior_cell][i] / face_area;
2572 PetscReal distance_to_second_cell = 2.0 * distance_to_first_cell +
2573 0.5 / cell_jacobian[k][second_interior_cell][i] / face_area;
2574
2575 PetscReal wall_normal[3];
2576 wall_normal[0] = eta[k][ghost_cell_index][i].x / face_area;
2577 wall_normal[1] = eta[k][ghost_cell_index][i].y / face_area;
2578 wall_normal[2] = eta[k][ghost_cell_index][i].z / face_area;
2579
2580 Cmpnts wall_velocity, reference_velocity;
2581 wall_velocity.x = wall_velocity.y = wall_velocity.z = 0.0;
2582 reference_velocity = velocity_cartesian[k][second_interior_cell][i];
2583
2584 noslip(user, distance_to_second_cell, distance_to_first_cell,
2585 wall_velocity, reference_velocity,
2586 &velocity_cartesian[k][first_interior_cell][i],
2587 wall_normal[0], wall_normal[1], wall_normal[2]);
2588
2589 wall_function_loglaw(user, wall_roughness_height,
2590 distance_to_second_cell, distance_to_first_cell,
2591 wall_velocity, reference_velocity,
2592 &velocity_cartesian[k][first_interior_cell][i],
2593 &friction_velocity[k][first_interior_cell][i],
2594 wall_normal[0], wall_normal[1], wall_normal[2]);
2595
2596 velocity_boundary[k][ghost_cell_index][i].x = 0.0;
2597 velocity_boundary[k][ghost_cell_index][i].y = 0.0;
2598 velocity_boundary[k][ghost_cell_index][i].z = 0.0;
2599 velocity_contravariant[k][ghost_cell_index][i].y = 0.0;
2600 }
2601 }
2602 }
2603 }
2604 } break;
2605
2606 // =================================================================
2607 // POSITIVE Y FACE (j = my-1, normal points in -Y direction)
2608 // =================================================================
2609 case BC_FACE_POS_Y: {
2610 if (grid_end_j == grid_size_j) {
2611 const PetscInt ghost_cell_index = grid_end_j - 1;
2612 const PetscInt first_interior_cell = grid_end_j - 2;
2613 const PetscInt second_interior_cell = grid_end_j - 3;
2614
2615 for (PetscInt k = loop_start_k; k < loop_end_k; k++) {
2616 for (PetscInt i = loop_start_i; i < loop_end_i; i++) {
2617
2618 if (node_vertex_flag[k][first_interior_cell][i] < 0.1) {
2619
2620 PetscReal face_area = sqrt(
2621 eta[k][first_interior_cell][i].x * eta[k][first_interior_cell][i].x +
2622 eta[k][first_interior_cell][i].y * eta[k][first_interior_cell][i].y +
2623 eta[k][first_interior_cell][i].z * eta[k][first_interior_cell][i].z
2624 );
2625
2626 PetscReal distance_to_first_cell = 0.5 / cell_jacobian[k][first_interior_cell][i] / face_area;
2627 PetscReal distance_to_second_cell = 2.0 * distance_to_first_cell +
2628 0.5 / cell_jacobian[k][second_interior_cell][i] / face_area;
2629
2630 PetscReal wall_normal[3];
2631 wall_normal[0] = -eta[k][first_interior_cell][i].x / face_area;
2632 wall_normal[1] = -eta[k][first_interior_cell][i].y / face_area;
2633 wall_normal[2] = -eta[k][first_interior_cell][i].z / face_area;
2634
2635 Cmpnts wall_velocity, reference_velocity;
2636 wall_velocity.x = wall_velocity.y = wall_velocity.z = 0.0;
2637 reference_velocity = velocity_cartesian[k][second_interior_cell][i];
2638
2639 noslip(user, distance_to_second_cell, distance_to_first_cell,
2640 wall_velocity, reference_velocity,
2641 &velocity_cartesian[k][first_interior_cell][i],
2642 wall_normal[0], wall_normal[1], wall_normal[2]);
2643
2644 wall_function_loglaw(user, wall_roughness_height,
2645 distance_to_second_cell, distance_to_first_cell,
2646 wall_velocity, reference_velocity,
2647 &velocity_cartesian[k][first_interior_cell][i],
2648 &friction_velocity[k][first_interior_cell][i],
2649 wall_normal[0], wall_normal[1], wall_normal[2]);
2650
2651 velocity_boundary[k][ghost_cell_index][i].x = 0.0;
2652 velocity_boundary[k][ghost_cell_index][i].y = 0.0;
2653 velocity_boundary[k][ghost_cell_index][i].z = 0.0;
2654 velocity_contravariant[k][first_interior_cell][i].y = 0.0;
2655 }
2656 }
2657 }
2658 }
2659 } break;
2660
2661 // =================================================================
2662 // NEGATIVE Z FACE (k = 0, normal points in +Z direction)
2663 // =================================================================
2664 case BC_FACE_NEG_Z: {
2665 if (grid_start_k == 0) {
2666 const PetscInt ghost_cell_index = grid_start_k;
2667 const PetscInt first_interior_cell = grid_start_k + 1;
2668 const PetscInt second_interior_cell = grid_start_k + 2;
2669
2670 for (PetscInt j = loop_start_j; j < loop_end_j; j++) {
2671 for (PetscInt i = loop_start_i; i < loop_end_i; i++) {
2672
2673 if (node_vertex_flag[first_interior_cell][j][i] < 0.1) {
2674
2675 PetscReal face_area = sqrt(
2676 zet[ghost_cell_index][j][i].x * zet[ghost_cell_index][j][i].x +
2677 zet[ghost_cell_index][j][i].y * zet[ghost_cell_index][j][i].y +
2678 zet[ghost_cell_index][j][i].z * zet[ghost_cell_index][j][i].z
2679 );
2680
2681 PetscReal distance_to_first_cell = 0.5 / cell_jacobian[first_interior_cell][j][i] / face_area;
2682 PetscReal distance_to_second_cell = 2.0 * distance_to_first_cell +
2683 0.5 / cell_jacobian[second_interior_cell][j][i] / face_area;
2684
2685 PetscReal wall_normal[3];
2686 wall_normal[0] = zet[ghost_cell_index][j][i].x / face_area;
2687 wall_normal[1] = zet[ghost_cell_index][j][i].y / face_area;
2688 wall_normal[2] = zet[ghost_cell_index][j][i].z / face_area;
2689
2690 Cmpnts wall_velocity, reference_velocity;
2691 wall_velocity.x = wall_velocity.y = wall_velocity.z = 0.0;
2692 reference_velocity = velocity_cartesian[second_interior_cell][j][i];
2693
2694 noslip(user, distance_to_second_cell, distance_to_first_cell,
2695 wall_velocity, reference_velocity,
2696 &velocity_cartesian[first_interior_cell][j][i],
2697 wall_normal[0], wall_normal[1], wall_normal[2]);
2698
2699 wall_function_loglaw(user, wall_roughness_height,
2700 distance_to_second_cell, distance_to_first_cell,
2701 wall_velocity, reference_velocity,
2702 &velocity_cartesian[first_interior_cell][j][i],
2703 &friction_velocity[first_interior_cell][j][i],
2704 wall_normal[0], wall_normal[1], wall_normal[2]);
2705
2706 velocity_boundary[ghost_cell_index][j][i].x = 0.0;
2707 velocity_boundary[ghost_cell_index][j][i].y = 0.0;
2708 velocity_boundary[ghost_cell_index][j][i].z = 0.0;
2709 velocity_contravariant[ghost_cell_index][j][i].z = 0.0;
2710 }
2711 }
2712 }
2713 }
2714 } break;
2715
2716 // =================================================================
2717 // POSITIVE Z FACE (k = mz-1, normal points in -Z direction)
2718 // =================================================================
2719 case BC_FACE_POS_Z: {
2720 if (grid_end_k == grid_size_k) {
2721 const PetscInt ghost_cell_index = grid_end_k - 1;
2722 const PetscInt first_interior_cell = grid_end_k - 2;
2723 const PetscInt second_interior_cell = grid_end_k - 3;
2724
2725 for (PetscInt j = loop_start_j; j < loop_end_j; j++) {
2726 for (PetscInt i = loop_start_i; i < loop_end_i; i++) {
2727
2728 if (node_vertex_flag[first_interior_cell][j][i] < 0.1) {
2729
2730 PetscReal face_area = sqrt(
2731 zet[first_interior_cell][j][i].x * zet[first_interior_cell][j][i].x +
2732 zet[first_interior_cell][j][i].y * zet[first_interior_cell][j][i].y +
2733 zet[first_interior_cell][j][i].z * zet[first_interior_cell][j][i].z
2734 );
2735
2736 PetscReal distance_to_first_cell = 0.5 / cell_jacobian[first_interior_cell][j][i] / face_area;
2737 PetscReal distance_to_second_cell = 2.0 * distance_to_first_cell +
2738 0.5 / cell_jacobian[second_interior_cell][j][i] / face_area;
2739
2740 PetscReal wall_normal[3];
2741 wall_normal[0] = -zet[first_interior_cell][j][i].x / face_area;
2742 wall_normal[1] = -zet[first_interior_cell][j][i].y / face_area;
2743 wall_normal[2] = -zet[first_interior_cell][j][i].z / face_area;
2744
2745 Cmpnts wall_velocity, reference_velocity;
2746 wall_velocity.x = wall_velocity.y = wall_velocity.z = 0.0;
2747 reference_velocity = velocity_cartesian[second_interior_cell][j][i];
2748
2749 noslip(user, distance_to_second_cell, distance_to_first_cell,
2750 wall_velocity, reference_velocity,
2751 &velocity_cartesian[first_interior_cell][j][i],
2752 wall_normal[0], wall_normal[1], wall_normal[2]);
2753
2754 wall_function_loglaw(user, wall_roughness_height,
2755 distance_to_second_cell, distance_to_first_cell,
2756 wall_velocity, reference_velocity,
2757 &velocity_cartesian[first_interior_cell][j][i],
2758 &friction_velocity[first_interior_cell][j][i],
2759 wall_normal[0], wall_normal[1], wall_normal[2]);
2760
2761 velocity_boundary[ghost_cell_index][j][i].x = 0.0;
2762 velocity_boundary[ghost_cell_index][j][i].y = 0.0;
2763 velocity_boundary[ghost_cell_index][j][i].z = 0.0;
2764 velocity_contravariant[first_interior_cell][j][i].z = 0.0;
2765 }
2766 }
2767 }
2768 }
2769 } break;
2770 }
2771 }
2772
2773 // =========================================================================
2774 // STEP 4: Restore all arrays and release memory
2775 // =========================================================================
2776 ierr = DMDAVecRestoreArray(user->fda, user->Ucat, &velocity_cartesian); CHKERRQ(ierr);
2777 ierr = DMDAVecRestoreArray(user->fda, user->Ucont, &velocity_contravariant); CHKERRQ(ierr);
2778 ierr = DMDAVecRestoreArray(user->fda, user->Bcs.Ubcs, &velocity_boundary); CHKERRQ(ierr);
2779 ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
2780 ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
2781 ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
2782 ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, (const PetscReal***)&node_vertex_flag); CHKERRQ(ierr);
2783 ierr = DMDAVecRestoreArrayRead(user->da, user->lAj, (const PetscReal***)&cell_jacobian); CHKERRQ(ierr);
2784 ierr = DMDAVecRestoreArray(user->da, user->lFriction_Velocity, &friction_velocity); CHKERRQ(ierr);
2785
2786 LOG_ALLOW(LOCAL, LOG_DEBUG, "Complete.\n");
2787
2788 PetscFunctionReturn(0);
2789}
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().
Definition Boundaries.c:126
UserCtx::lFriction_Velocity
Vec lFriction_Velocity
Definition variables.h:833
WALL
@ WALL
Definition variables.h:254
UserCtx::KM
PetscInt KM
Definition variables.h:820
UserCtx::JM
PetscInt JM
Definition variables.h:820
SimCtx::wallfunction
PetscInt wallfunction
Definition variables.h:733
UserCtx::IM
PetscInt IM
Definition variables.h:820
wall_function_loglaw
void wall_function_loglaw(UserCtx *user, double roughness_height, double distance_reference, double distance_boundary, Cmpnts velocity_wall, Cmpnts velocity_reference, Cmpnts *velocity_boundary, PetscReal *friction_velocity, double normal_x, double normal_y, double normal_z)
Applies log-law wall function with roughness correction.
Definition wallfunction.c:704
noslip
void noslip(UserCtx *user, double distance_reference, double distance_boundary, Cmpnts velocity_wall, Cmpnts velocity_reference, Cmpnts *velocity_boundary, double normal_x, double normal_y, double normal_z)
Applies no-slip wall boundary condition with linear interpolation.
Definition wallfunction.c:68
Here is the call graph for this function:
Here is the caller graph for this function:

◆ RefreshBoundaryGhostCells()

PetscErrorCode RefreshBoundaryGhostCells ( UserCtx *  user)

Internal helper implementation: RefreshBoundaryGhostCells().

(Public) Orchestrates the "light" refresh of all boundary ghost cells after the projection step.

Local to this translation unit.

Definition at line 2797 of file Boundaries.c.

2798{
2799 PetscErrorCode ierr;
2800 PetscFunctionBeginUser;
2801 PROFILE_FUNCTION_BEGIN;
2802
2803 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Starting post-projection refresh of boundary ghost cells.\n");
2804
2805 // -------------------------------------------------------------------------
2806 // STEP 1: Refresh Flow-Dependent Boundary Value Targets (`ubcs`)
2807 // -------------------------------------------------------------------------
2808 // This is the "physics" part of the refresh. It calls the lightweight engine
2809 // to update `ubcs` for any BCs (like Symmetry) that depend on the now-corrected
2810 // interior velocity field. This step does NOT touch `ucont`.
2811 ierr = BoundarySystem_RefreshUbcs(user); CHKERRQ(ierr);
2812 LOG_ALLOW(GLOBAL, LOG_VERBOSE, " `ubcs` targets refreshed.\n");
2813
2814 ierr = UpdateLocalGhosts(user,"Ucat");
2815
2816 // STEP 1.5: Apply Wall function if applicable.
2817 if(user->simCtx->wallfunction){
2818
2819 ierr = ApplyWallFunction(user); CHKERRQ(ierr);
2820
2821 }
2822 // -------------------------------------------------------------------------
2823 // STEP 2: Apply Geometric Ghost Cell Updates
2824 // -------------------------------------------------------------------------
2825 // With `ubcs` now fully up-to-date, we execute the purely geometric
2826 // operations to fill the entire ghost cell layer. The order is important.
2827
2828 // (a) Update the ghost cells on the faces of non-periodic boundaries.
2829 ierr = UpdateDummyCells(user); CHKERRQ(ierr);
2830 LOG_ALLOW(GLOBAL, LOG_VERBOSE, " Face ghost cells (dummy cells) updated.\n");
2831
2832 // (b) Handle periodic boundaries first. This is a direct data copy.
2833 // Ghost Update is done inside this function.s
2834 ierr = ApplyPeriodicBCs(user); CHKERRQ(ierr);
2835 LOG_ALLOW(GLOBAL, LOG_VERBOSE, " Periodic boundaries synchronized.\n");
2836
2837 // (c) Update the ghost cells at the edges and corners by averaging.
2838 ierr = UpdateCornerNodes(user); CHKERRQ(ierr);
2839 LOG_ALLOW(GLOBAL, LOG_VERBOSE, " Edge and corner ghost cells updated.\n");
2840
2841 // (d) Synchronize Ucat after setting corner nodes.
2842 ierr = UpdateLocalGhosts(user, "Ucat"); CHKERRQ(ierr);
2843
2844 // (e) Update the corners for periodic conditions sequentially
2845 // ensuring no race conditions are raised.
2846 const char* ucat_only[] = {"Ucat"};
2847 ierr = UpdatePeriodicCornerNodes(user,1,ucat_only);CHKERRQ(ierr);
2848
2849 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Boundary ghost cell refresh complete.\n");
2850 PROFILE_FUNCTION_END;
2851 PetscFunctionReturn(0);
2852}
ApplyPeriodicBCs
PetscErrorCode ApplyPeriodicBCs(UserCtx *user)
Internal helper implementation: ApplyPeriodicBCs().
Definition Boundaries.c:1922
ApplyWallFunction
PetscErrorCode ApplyWallFunction(UserCtx *user)
Internal helper implementation: ApplyWallFunction().
Definition Boundaries.c:2338
UpdateDummyCells
PetscErrorCode UpdateDummyCells(UserCtx *user)
Internal helper implementation: UpdateDummyCells().
Definition Boundaries.c:2059
BoundarySystem_RefreshUbcs
PetscErrorCode BoundarySystem_RefreshUbcs(UserCtx *user)
Internal helper implementation: BoundarySystem_RefreshUbcs().
Definition Boundaries.c:1475
UpdatePeriodicCornerNodes
PetscErrorCode UpdatePeriodicCornerNodes(UserCtx *user, PetscInt num_fields, const char *field_names[])
Internal helper implementation: UpdatePeriodicCornerNodes().
Definition Boundaries.c:2295
UpdateCornerNodes
PetscErrorCode UpdateCornerNodes(UserCtx *user)
Internal helper implementation: UpdateCornerNodes().
Definition Boundaries.c:2151
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 2862 of file Boundaries.c.

2863{
2864 PetscErrorCode ierr;
2865 PetscFunctionBeginUser;
2866 PROFILE_FUNCTION_BEGIN;
2867
2868 LOG_ALLOW(GLOBAL,LOG_TRACE,"Boundary Condition Application begins.\n");
2869
2870 // STEP 1: Main iteration loop for applying and converging non-periodic BCs.
2871 // The number of iterations (e.g., 3) allows information to propagate
2872 // between coupled boundaries, like an inlet and a conserving outlet.
2873 for (PetscInt iter = 0; iter < 3; iter++) {
2874 // (a) Execute the boundary system. This phase calculates fluxes across
2875 // the domain and then applies the physical logic for each non-periodic
2876 // handler, setting the `ubcs` (boundary value) array.
2877 ierr = BoundarySystem_ExecuteStep(user); CHKERRQ(ierr);
2878
2879 LOG_ALLOW(GLOBAL,LOG_VERBOSE,"Boundary Condition Setup Executed.\n");
2880
2881 // (b) Synchronize the updated ghost cells across all processors to ensure
2882 // all ucont values are current before updating the dummy cells.
2883 ierr = UpdateLocalGhosts(user, "Ucont"); CHKERRQ(ierr);
2884
2885 // (c) Convert updated Contravariant velocities to Cartesian velocities.
2886 ierr = Contra2Cart(user); CHKERRQ(ierr);
2887
2888 // (d) Synchronize the updated Cartesian velocities across all processors
2889 // to ensure all ucat values are current before updating the dummy cells.
2890 ierr = UpdateLocalGhosts(user, "Ucat"); CHKERRQ(ierr);
2891
2892 // (e) If Wall functions are enabled, apply them now to adjust near-wall velocities.
2893 if(user->simCtx->wallfunction){
2894 // Apply wall function adjustments to the boundary velocities.
2895 ierr = ApplyWallFunction(user); CHKERRQ(ierr);
2896
2897 // Synchronize the updated Cartesian velocities after wall function adjustments.
2898 ierr = UpdateLocalGhosts(user, "Ucat"); CHKERRQ(ierr);
2899
2900 LOG_ALLOW(GLOBAL,LOG_VERBOSE,"Wall Function Applied at Walls.\n");
2901 }
2902
2903 // (f) Update the first layer of ghost cells for non-periodic faces using
2904 // the newly computed `ubcs` values.
2905 ierr = UpdateDummyCells(user); CHKERRQ(ierr);
2906
2907 LOG_ALLOW(GLOBAL,LOG_VERBOSE,"Dummy Cells/Ghost Cells Updated.\n");
2908
2909 // (g) Handle all periodic boundaries. This is a parallel direct copy
2910 // that sets the absolute constraints for the rest of the solve.
2911 // There is a Ghost update happening inside this function.
2912 ierr = ApplyPeriodicBCs(user); CHKERRQ(ierr);
2913
2914 // (h) Update the corner and edge ghost nodes. This routine calculates
2915 // values for corners/edges by averaging their neighbors, which have been
2916 // finalized in the steps above (both periodic and non-periodic).
2917 ierr = UpdateCornerNodes(user); CHKERRQ(ierr);
2918
2919 // (i) Synchronize the updated edge and corner cells across all processors to ensure
2920 // consistency before the next iteration or finalization.
2921 ierr = UpdateLocalGhosts(user, "P"); CHKERRQ(ierr);
2922 ierr = UpdateLocalGhosts(user, "Ucat"); CHKERRQ(ierr);
2923 ierr = UpdateLocalGhosts(user, "Ucont"); CHKERRQ(ierr);
2924
2925 // (j) Ensure All the corners are synchronized with a well defined protocol in case of Periodic boundary conditions
2926 // To avoid race conditions.
2927 const char* all_fields[] = {"Ucat", "P", "Nvert"};
2928 ierr = UpdatePeriodicCornerNodes(user,3,all_fields); CHKERRQ(ierr);
2929
2930 }
2931
2932 // STEP 3: Final ghost node synchronization. This ensures all changes made
2933 // to the global vectors are reflected in the local ghost regions of all
2934 // processors, making the state fully consistent before the next solver stage.
2935 ierr = UpdateLocalGhosts(user, "P"); CHKERRQ(ierr);
2936 ierr = UpdateLocalGhosts(user, "Ucat"); CHKERRQ(ierr);
2937 ierr = UpdateLocalGhosts(user, "Ucont"); CHKERRQ(ierr);
2938
2939 PROFILE_FUNCTION_END;
2940 PetscFunctionReturn(0);
2941}
BoundarySystem_ExecuteStep
PetscErrorCode BoundarySystem_ExecuteStep(UserCtx *user)
Implementation of BoundarySystem_ExecuteStep().
Definition Boundaries.c:1053
Contra2Cart
PetscErrorCode Contra2Cart(UserCtx *user)
Reconstructs Cartesian velocity (Ucat) at cell centers from contravariant velocity (Ucont) defined on...
Definition setup.c:2247
Here is the call graph for this function:
Here is the caller graph for this function:
Generated by doxygen 1.9.8