PICurv 0.1.0
A Parallel Particle-In-Cell Solver for Curvilinear LES
Loading...
Searching...
No Matches
Boundaries.c
Go to the documentation of this file.
1#include "Boundaries.h" // The main header for our project
2#include <string.h> // For strcasecmp
3#include <ctype.h> // For isspace
4
5#undef __FUNCT__
6#define __FUNCT__ "CanRankServiceInletFace"
7/**
8 * @brief Internal helper implementation: `CanRankServiceInletFace()`.
9 * @details Local to this translation unit.
10 */
11PetscErrorCode CanRankServiceInletFace(UserCtx *user, const DMDALocalInfo *info,
12 PetscInt IM_nodes_global, PetscInt JM_nodes_global, PetscInt KM_nodes_global,
13 PetscBool *can_service_inlet_out)
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}
116
117#undef __FUNCT__
118#define __FUNCT__ "CanRankServiceFace"
119
120/**
121 * @brief Implementation of \ref CanRankServiceFace().
122 * @details Full API contract (arguments, ownership, side effects) is documented with
123 * the header declaration in `include/Boundaries.h`.
124 * @see CanRankServiceFace()
125 */
126PetscErrorCode CanRankServiceFace(const DMDALocalInfo *info, PetscInt IM_nodes_global, PetscInt JM_nodes_global, PetscInt KM_nodes_global,
127 BCFace face_id, PetscBool *can_service_out)
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}
204
205#undef __FUNCT__
206#define __FUNCT__ "GetDeterministicFaceGridLocation"
207
208/**
209 * @brief Internal helper implementation: `GetDeterministicFaceGridLocation()`.
210 * @details Local to this translation unit.
211 */
212PetscErrorCode GetDeterministicFaceGridLocation(
213 UserCtx *user, const DMDALocalInfo *info,
214 PetscInt xs_gnode_rank, PetscInt ys_gnode_rank, PetscInt zs_gnode_rank,
215 PetscInt IM_cells_global, PetscInt JM_cells_global, PetscInt KM_cells_global,
216 PetscInt64 particle_global_id,
217 PetscInt *ci_metric_lnode_out, PetscInt *cj_metric_lnode_out, PetscInt *ck_metric_lnode_out,
218 PetscReal *xi_metric_logic_out, PetscReal *eta_metric_logic_out, PetscReal *zta_metric_logic_out,
219 PetscBool *placement_successful_out)
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}
391
392#undef __FUNCT__
393#define __FUNCT__ "GetRandomFCellAndLogicOnInletFace"
394
395/**
396 * @brief Internal helper implementation: `GetRandomCellAndLogicalCoordsOnInletFace()`.
397 * @details Local to this translation unit.
398 */
399PetscErrorCode GetRandomCellAndLogicalCoordsOnInletFace(
400 UserCtx *user, const DMDALocalInfo *info,
401 PetscInt xs_gnode_rank, PetscInt ys_gnode_rank, PetscInt zs_gnode_rank, // Local starting node index (with ghosts) of the rank's DA patch
402 PetscInt IM_nodes_global, PetscInt JM_nodes_global, PetscInt KM_nodes_global,
403 PetscRandom *rand_logic_i_ptr, PetscRandom *rand_logic_j_ptr, PetscRandom *rand_logic_k_ptr,
404 PetscInt *ci_metric_lnode_out, PetscInt *cj_metric_lnode_out, PetscInt *ck_metric_lnode_out,
405 PetscReal *xi_metric_logic_out, PetscReal *eta_metric_logic_out, PetscReal *zta_metric_logic_out)
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}
582
583
584
585#undef __FUNCT__
586#define __FUNCT__ "EnforceRHSBoundaryConditions"
587/**
588 * @brief Internal helper implementation: `EnforceRHSBoundaryConditions()`.
589 * @details Local to this translation unit.
590 */
591PetscErrorCode EnforceRHSBoundaryConditions(UserCtx *user)
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}
736
737#undef __FUNCT__
738#define __FUNCT__ "BoundaryCondition_Create"
739/**
740 * @brief Internal helper implementation: `BoundaryCondition_Create()`.
741 * @details Local to this translation unit.
742 */
743
744PetscErrorCode BoundaryCondition_Create(BCHandlerType handler_type, BoundaryCondition **new_bc_ptr)
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
804 case BC_HANDLER_INLET_PROFILE_FROM_FILE:
805 LOG_ALLOW(LOCAL, LOG_DEBUG, "Dispatching to Create_InletProfileFromFile().\n");
806 ierr = Create_InletProfileFromFile(bc); CHKERRQ(ierr);
807 break;
808 //Add cases for other handlers here in future phases
809
810 default:
811 LOG_ALLOW(GLOBAL, LOG_ERROR, "Handler type (%s) is not recognized or implemented in the factory.\n", handler_name);
812 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Boundary handler type %d (%s) not recognized in factory.\n", handler_type, handler_name);
813 }
814
815 if(bc->priority < 0) {
816 LOG_ALLOW(GLOBAL, LOG_ERROR, "Handler type %d (%s) did not set a valid priority during creation.\n", handler_type, handler_name);
817 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);
818 }
819
820 LOG_ALLOW(LOCAL, LOG_DEBUG, "Successfully created and configured handler for %s.\n", handler_name);
821 PetscFunctionReturn(0);
822}
823
824#undef __FUNCT__
825#define __FUNCT__ "BoundarySystem_Validate"
826/**
827 * @brief Internal helper implementation: `BoundarySystem_Validate()`.
828 * @details Local to this translation unit.
829 */
830PetscErrorCode BoundarySystem_Validate(UserCtx *user)
831{
832 PetscErrorCode ierr;
833 PetscFunctionBeginUser;
834
835 LOG_ALLOW(GLOBAL, LOG_INFO, "Validating parsed boundary condition configuration...\n");
836
837 // --- Rule Set 1: Driven Flow Handler Consistency ---
838 // This specialized validator will check all rules related to driven flow handlers.
839 ierr = Validate_DrivenFlowConfiguration(user); CHKERRQ(ierr);
840
841 // --- Rule Set 2: (Future Extension) Overset Interface Consistency ---
842 // ierr = Validate_OversetConfiguration(user); CHKERRQ(ierr);
843
844 LOG_ALLOW(GLOBAL, LOG_INFO, "Boundary configuration is valid.\n");
845
846 PetscFunctionReturn(0);
847}
848
849//================================================================================
850//
851// PUBLIC MASTER SETUP FUNCTION
852//
853//================================================================================
854#undef __FUNCT__
855#define __FUNCT__ "BoundarySystem_Initialize"
856/**
857 * @brief Implementation of \ref BoundarySystem_Initialize().
858 * @details Full API contract (arguments, ownership, side effects) is documented with
859 * the header declaration in `include/Boundaries.h`.
860 * @see BoundarySystem_Initialize()
861 */
862PetscErrorCode BoundarySystem_Initialize(UserCtx *user, const char *bcs_filename)
863{
864 PetscErrorCode ierr;
865 PetscFunctionBeginUser;
866
867 LOG_ALLOW(GLOBAL, LOG_INFO, "Starting creation and initialization of all boundary handlers.\n");
868
869 // =========================================================================
870 // Step 0: Clear any existing boundary handlers (if re-initializing).
871 // This ensures no memory leaks if this function is called multiple times.
872 // =========================================================================
873 for (int i = 0; i < 6; i++) {
874 BoundaryFaceConfig *face_cfg = &user->boundary_faces[i];
875 if (face_cfg->handler) {
876 LOG_ALLOW(LOCAL, LOG_DEBUG, "Destroying existing handler on Face %s before re-initialization.\n", BCFaceToString((BCFace)i));
877 if (face_cfg->handler->Destroy) {
878 ierr = face_cfg->handler->Destroy(face_cfg->handler); CHKERRQ(ierr);
879 }
880 ierr = PetscFree(face_cfg->handler); CHKERRQ(ierr);
881 face_cfg->handler = NULL;
882 }
883 }
884 // =========================================================================
885
886 // Step 0.1: Initiate flux sums to zero
887 user->simCtx->FluxInSum = 0.0;
888 user->simCtx->FluxOutSum = 0.0;
889 user->simCtx->FarFluxInSum = 0.0;
890 user->simCtx->FarFluxOutSum = 0.0;
891 // =========================================================================
892
893 // Step 1: Parse the configuration file to determine user intent.
894 // This function, defined in io.c, populates the configuration enums and parameter
895 // lists within the user->boundary_faces array on all MPI ranks.
896 ierr = ParseAllBoundaryConditions(user, bcs_filename); CHKERRQ(ierr);
897 LOG_ALLOW(GLOBAL, LOG_INFO, "Configuration file '%s' parsed successfully.\n", bcs_filename);
898
899 // Step 1.1: Validate the parsed configuration to ensure there are no Boundary Condition conflicts
900 ierr = BoundarySystem_Validate(user); CHKERRQ(ierr);
901
902 // Step 2: Create and Initialize the handler object for each of the 6 faces.
903 for (int i = 0; i < 6; i++) {
904 BoundaryFaceConfig *face_cfg = &user->boundary_faces[i];
905
906 const char *face_name = BCFaceToString(face_cfg->face_id);
907 const char *type_name = BCTypeToString(face_cfg->mathematical_type);
908 const char *handler_name = BCHandlerTypeToString(face_cfg->handler_type);
909
910 LOG_ALLOW(LOCAL, LOG_DEBUG, "Creating handler for Face %s with Type %s and handler '%s'.\n", face_name, type_name,handler_name);
911
912 // Use the private factory to construct the correct handler object based on the parsed type.
913 // The factory returns a pointer to the new handler object, which we store in the config struct.
914 ierr = BoundaryCondition_Create(face_cfg->handler_type, &face_cfg->handler); CHKERRQ(ierr);
915
916 // Step 3: Call the specific Initialize() method for the newly created handler.
917 // This allows the handler to perform its own setup, like reading parameters from the
918 // face_cfg->params list and setting the initial field values on its face.
919 if (face_cfg->handler && face_cfg->handler->Initialize) {
920 LOG_ALLOW(LOCAL, LOG_DEBUG, "Calling Initialize() method for handler %s(%s) on Face %s.\n",type_name,handler_name,face_name);
921
922 // Prepare the context needed by the Initialize() function.
923 BCContext ctx = {
924 .user = user,
925 .face_id = face_cfg->face_id,
926 .global_inflow_sum = &user->simCtx->FluxInSum, // Global flux sums are not relevant during initialization.
927 .global_outflow_sum = &user->simCtx->FluxOutSum,
928 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
929 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
930 };
931
932 ierr = face_cfg->handler->Initialize(face_cfg->handler, &ctx); CHKERRQ(ierr);
933 } else {
934 LOG_ALLOW(LOCAL, LOG_DEBUG, "Handler %s(%s) for Face %s has no Initialize() method, skipping.\n", type_name,handler_name,face_name);
935 }
936 }
937 // =========================================================================
938 // NO SYNCHRONIZATION NEEDED HERE
939 // =========================================================================
940 // Initialize() only reads parameters and allocates memory.
941 // It does NOT modify field values (Ucat, Ucont, Ubcs).
942 // Field values are set by:
943 // 1. Initial conditions (before this function)
944 // 2. Apply() during timestepping (after this function)
945 // The first call to ApplyBoundaryConditions() will handle synchronization.
946 // =========================================================================
947
948 // ====================================================================================
949 // --- NEW: Step 4: Synchronize Vectors After Initialization ---
950 // This is the CRITICAL fix. The Initialize() calls have modified local vector
951 // arrays on some ranks but not others. We must now update the global vector state
952 // and then update all local ghost regions to be consistent.
953 // ====================================================================================
954
955 //LOG_ALLOW(GLOBAL, LOG_DEBUG, "Committing global boundary initializations to local vectors.\n");
956
957 // Commit changes from the global vectors (Ucat, Ucont) to the local vectors (lUcat, lUcont)
958 // NOTE: The Apply functions modified Ucat and Ucont via GetArray, which works on the global
959 // representation.
960 /*
961 ierr = DMGlobalToLocalBegin(user->fda, user->Ucat, INSERT_VALUES, user->lUcat); CHKERRQ(ierr);
962 ierr = DMGlobalToLocalEnd(user->fda, user->Ucat, INSERT_VALUES, user->lUcat); CHKERRQ(ierr);
963
964 ierr = DMGlobalToLocalBegin(user->fda, user->Ucont, INSERT_VALUES, user->lUcont); CHKERRQ(ierr);
965 ierr = DMGlobalToLocalEnd(user->fda, user->Ucont, INSERT_VALUES, user->lUcont); CHKERRQ(ierr);
966
967 // Now, update all local vectors (including ghost cells) from the newly consistent global vectors
968
969 ierr = DMLocalToGlobalBegin(user->fda, user->lUcat, INSERT_VALUES, user->Ucat); CHKERRQ(ierr);
970 ierr = DMLocalToGlobalEnd(user->fda, user->lUcat, INSERT_VALUES, user->Ucat); CHKERRQ(ierr);
971
972 ierr = DMLocalToGlobalBegin(user->fda, user->lUcont, INSERT_VALUES, user->Ucont); CHKERRQ(ierr);
973 ierr = DMLocalToGlobalEnd(user->fda, user->lUcont, INSERT_VALUES, user->Ucont); CHKERRQ(ierr);
974 */
975
976 LOG_ALLOW(GLOBAL, LOG_INFO, "All boundary handlers created and initialized successfully.\n");
977 PetscFunctionReturn(0);
978}
979
980
981#undef __FUNCT__
982#define __FUNCT__ "PropagateBoundaryConfigToCoarserLevels"
983/**
984 * @brief Internal helper implementation: `PropagateBoundaryConfigToCoarserLevels()`.
985 * @details Local to this translation unit.
986 */
987PetscErrorCode PropagateBoundaryConfigToCoarserLevels(SimCtx *simCtx)
988{
989 PetscErrorCode ierr;
990 UserMG *usermg = &simCtx->usermg;
991
992 PetscFunctionBeginUser;
993 PROFILE_FUNCTION_BEGIN;
994
995 LOG_ALLOW(GLOBAL, LOG_INFO, "Propagating BC configuration from finest to coarser multigrid levels...\n");
996
997 // Loop from second-finest down to coarsest
998 for (PetscInt level = usermg->mglevels - 2; level >= 0; level--) {
999 for (PetscInt bi = 0; bi < simCtx->block_number; bi++) {
1000 UserCtx *user_coarse = &usermg->mgctx[level].user[bi];
1001 UserCtx *user_fine = &usermg->mgctx[level + 1].user[bi];
1002
1003 LOG_ALLOW_SYNC(LOCAL, LOG_DEBUG, "Rank %d: Copying BC config from level %d to level %d, block %d\n",
1004 simCtx->rank, level + 1, level, bi);
1005
1006 // Copy the 6 boundary face configurations
1007 for (int face_i = 0; face_i < 6; face_i++) {
1008 user_coarse->boundary_faces[face_i].face_id = user_fine->boundary_faces[face_i].face_id;
1009 user_coarse->boundary_faces[face_i].mathematical_type = user_fine->boundary_faces[face_i].mathematical_type;
1010 user_coarse->boundary_faces[face_i].handler_type = user_fine->boundary_faces[face_i].handler_type;
1011
1012 // Copy parameter list (deep copy)
1013 FreeBC_ParamList(user_coarse->boundary_faces[face_i].params); // Clear any existing
1014 user_coarse->boundary_faces[face_i].params = NULL;
1015
1016 BC_Param **dst_next = &user_coarse->boundary_faces[face_i].params;
1017 for (BC_Param *src = user_fine->boundary_faces[face_i].params; src; src = src->next) {
1018 BC_Param *new_param;
1019 ierr = PetscMalloc1(1, &new_param); CHKERRQ(ierr);
1020 ierr = PetscStrallocpy(src->key, &new_param->key); CHKERRQ(ierr);
1021 ierr = PetscStrallocpy(src->value, &new_param->value); CHKERRQ(ierr);
1022 new_param->next = NULL;
1023 *dst_next = new_param;
1024 dst_next = &new_param->next;
1025 }
1026
1027 // IMPORTANT: Do NOT create handler objects for coarser levels
1028 // Handlers are only needed at finest level for timestepping Apply() calls
1029 user_coarse->boundary_faces[face_i].handler = NULL;
1030 }
1031
1032 // Propagate the particle inlet lookup fields to coarse levels as well.
1033 user_coarse->inletFaceDefined = user_fine->inletFaceDefined;
1034 user_coarse->identifiedInletBCFace = user_fine->identifiedInletBCFace;
1035 }
1036 }
1037
1038 LOG_ALLOW(GLOBAL, LOG_INFO, "BC configuration propagation complete.\n");
1039
1040 PROFILE_FUNCTION_END;
1041 PetscFunctionReturn(0);
1042}
1043
1044//================================================================================
1045//
1046// PUBLIC MASTER TIME-STEP FUNCTION
1047//
1048//================================================================================
1049
1050#undef __FUNCT__
1051#define __FUNCT__ "BoundarySystem_ExecuteStep"
1052/**
1053 * @brief Implementation of \ref BoundarySystem_ExecuteStep().
1054 * @details Full API contract (arguments, ownership, side effects) is documented with
1055 * the header declaration in `include/Boundaries.h`.
1056 * @see BoundarySystem_ExecuteStep()
1057 */
1058PetscErrorCode BoundarySystem_ExecuteStep(UserCtx *user)
1059{
1060 PetscErrorCode ierr;
1061 PetscFunctionBeginUser;
1062 PROFILE_FUNCTION_BEGIN;
1063
1064 LOG_ALLOW(LOCAL, LOG_DEBUG, "Starting.\n");
1065
1066 // =========================================================================
1067 // PRIORITY 0: INLETS
1068 // =========================================================================
1069
1070 PetscReal local_inflow_pre = 0.0;
1071 PetscReal local_inflow_post = 0.0;
1072 PetscReal global_inflow_pre = 0.0;
1073 PetscReal global_inflow_post = 0.0;
1074 PetscInt num_handlers[3] = {0,0,0};
1075
1076 LOG_ALLOW(LOCAL, LOG_TRACE, " (INLETS): Begin.\n");
1077
1078 // Phase 1: PreStep - Preparation (e.g., calculate profiles, read files)
1079 for (int i = 0; i < 6; i++) {
1080 BoundaryCondition *handler = user->boundary_faces[i].handler;
1081 if (!handler || handler->priority != BC_PRIORITY_INLET) continue;
1082 if (!handler->PreStep) continue;
1083
1084 num_handlers[0]++;
1085 BCContext ctx = {
1086 .user = user,
1087 .face_id = (BCFace)i,
1088 .global_inflow_sum = NULL,
1089 .global_outflow_sum = NULL,
1090 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1091 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1092 };
1093
1094 LOG_ALLOW(LOCAL, LOG_TRACE, " PreStep: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1095 ierr = handler->PreStep(handler, &ctx, &local_inflow_pre, NULL); CHKERRQ(ierr);
1096 }
1097
1098 // Optional: Global communication for PreStep (for debugging)
1099 if (local_inflow_pre != 0.0) {
1100 ierr = MPI_Allreduce(&local_inflow_pre, &global_inflow_pre, 1, MPIU_REAL,
1101 MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
1102 LOG_ALLOW(GLOBAL, LOG_TRACE, " PreStep predicted flux: %.6e\n", global_inflow_pre);
1103 }
1104
1105 // Phase 2: Apply - Set boundary conditions
1106 for (int i = 0; i < 6; i++) {
1107 BoundaryCondition *handler = user->boundary_faces[i].handler;
1108 if (!handler || handler->priority != BC_PRIORITY_INLET) continue;
1109 if(!handler->Apply) continue; // For example Periodic BCs
1110
1111 num_handlers[1]++;
1112
1113 BCContext ctx = {
1114 .user = user,
1115 .face_id = (BCFace)i,
1116 .global_inflow_sum = NULL,
1117 .global_outflow_sum = NULL,
1118 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1119 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1120 };
1121
1122 LOG_ALLOW(LOCAL, LOG_TRACE, " Apply: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1123 ierr = handler->Apply(handler, &ctx); CHKERRQ(ierr);
1124 }
1125
1126 // Phase 3: PostStep - Measure actual flux
1127 for (int i = 0; i < 6; i++) {
1128 BoundaryCondition *handler = user->boundary_faces[i].handler;
1129 if (!handler || handler->priority != BC_PRIORITY_INLET) continue;
1130 if (!handler->PostStep) continue;
1131
1132 num_handlers[2]++;
1133
1134 BCContext ctx = {
1135 .user = user,
1136 .face_id = (BCFace)i,
1137 .global_inflow_sum = NULL,
1138 .global_outflow_sum = NULL,
1139 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1140 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1141 };
1142
1143 LOG_ALLOW(LOCAL, LOG_TRACE, " PostStep: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1144 ierr = handler->PostStep(handler, &ctx, &local_inflow_post, NULL); CHKERRQ(ierr);
1145 }
1146
1147 // Phase 4: Global communication - Sum flux for other priorities to use
1148 ierr = MPI_Allreduce(&local_inflow_post, &global_inflow_post, 1, MPIU_REAL,
1149 MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
1150
1151 // Store for next priority levels
1152 user->simCtx->FluxInSum = global_inflow_post;
1153
1154 LOG_ALLOW(GLOBAL, LOG_INFO,
1155 " (INLETS): %d Prestep(s), %d Application(s), %d Poststep(s), FluxInSum = %.6e\n",
1156 num_handlers[0],num_handlers[1],num_handlers[2], global_inflow_post);
1157
1158 // =========================================================================
1159 // PRIORITY 1: FARFIELD
1160 // =========================================================================
1161
1162 PetscReal local_farfield_in_pre = 0.0;
1163 PetscReal local_farfield_out_pre = 0.0;
1164 PetscReal local_farfield_in_post = 0.0;
1165 PetscReal local_farfield_out_post = 0.0;
1166 PetscReal global_farfield_in_pre = 0.0;
1167 PetscReal global_farfield_out_pre = 0.0;
1168 PetscReal global_farfield_in_post = 0.0;
1169 PetscReal global_farfield_out_post = 0.0;
1170 memset(num_handlers,0,sizeof(num_handlers));
1171
1172 LOG_ALLOW(LOCAL, LOG_TRACE, " (FARFIELD): Begin.\n");
1173
1174 // Phase 1: PreStep - Analyze flow direction, measure initial flux
1175 for (int i = 0; i < 6; i++) {
1176 BoundaryCondition *handler = user->boundary_faces[i].handler;
1177 if (!handler || handler->priority != BC_PRIORITY_FARFIELD) continue;
1178 if (!handler->PreStep) continue;
1179
1180 num_handlers[0]++;
1181 BCContext ctx = {
1182 .user = user,
1183 .face_id = (BCFace)i,
1184 .global_inflow_sum = &user->simCtx->FluxInSum, // Available from Priority 0
1185 .global_outflow_sum = NULL,
1186 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1187 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1188 };
1189
1190 LOG_ALLOW(LOCAL, LOG_TRACE, " PreStep: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1191 ierr = handler->PreStep(handler, &ctx, &local_farfield_in_pre, &local_farfield_out_pre);
1192 CHKERRQ(ierr);
1193 }
1194
1195 // Phase 2: Global communication (optional, for debugging)
1196 if (local_farfield_in_pre != 0.0 || local_farfield_out_pre != 0.0) {
1197 ierr = MPI_Allreduce(&local_farfield_in_pre, &global_farfield_in_pre, 1, MPIU_REAL,
1198 MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
1199 ierr = MPI_Allreduce(&local_farfield_out_pre, &global_farfield_out_pre, 1, MPIU_REAL,
1200 MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
1201
1202 LOG_ALLOW(GLOBAL, LOG_DEBUG,
1203 " Farfield pre-analysis: In=%.6e, Out=%.6e\n",
1204 global_farfield_in_pre, global_farfield_out_pre);
1205 }
1206
1207 // Phase 3: Apply - Set farfield boundary conditions
1208 for (int i = 0; i < 6; i++) {
1209 BoundaryCondition *handler = user->boundary_faces[i].handler;
1210 if (!handler || handler->priority != BC_PRIORITY_FARFIELD) continue;
1211 if(!handler->Apply) continue; // For example Periodic BCs
1212
1213 num_handlers[1]++;
1214
1215 BCContext ctx = {
1216 .user = user,
1217 .face_id = (BCFace)i,
1218 .global_inflow_sum = &user->simCtx->FluxInSum,
1219 .global_outflow_sum = NULL,
1220 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1221 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1222 };
1223
1224 LOG_ALLOW(LOCAL, LOG_TRACE, " Apply: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1225 ierr = handler->Apply(handler, &ctx); CHKERRQ(ierr);
1226 }
1227
1228 // Phase 4: PostStep - Measure actual farfield fluxes
1229 for (int i = 0; i < 6; i++) {
1230 BoundaryCondition *handler = user->boundary_faces[i].handler;
1231 if (!handler || handler->priority != BC_PRIORITY_FARFIELD) continue;
1232 if (!handler->PostStep) continue;
1233
1234 num_handlers[2]++;
1235
1236 BCContext ctx = {
1237 .user = user,
1238 .face_id = (BCFace)i,
1239 .global_inflow_sum = &user->simCtx->FluxInSum,
1240 .global_outflow_sum = NULL,
1241 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1242 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1243 };
1244
1245 LOG_ALLOW(LOCAL, LOG_TRACE, " PostStep: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1246 ierr = handler->PostStep(handler, &ctx, &local_farfield_in_post, &local_farfield_out_post);
1247 CHKERRQ(ierr);
1248 }
1249
1250 // Phase 5: Global communication - Store for outlet priority
1251 if (num_handlers > 0) {
1252 ierr = MPI_Allreduce(&local_farfield_in_post, &global_farfield_in_post, 1, MPIU_REAL,
1253 MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
1254 ierr = MPI_Allreduce(&local_farfield_out_post, &global_farfield_out_post, 1, MPIU_REAL,
1255 MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
1256
1257 // Store for outlet handlers to use
1258 user->simCtx->FarFluxInSum = global_farfield_in_post;
1259 user->simCtx->FarFluxOutSum = global_farfield_out_post;
1260
1261 LOG_ALLOW(GLOBAL, LOG_INFO,
1262 " (FARFIELD): %d Prestep(s), %d Application(s), %d Poststep(s) , InFlux=%.6e, OutFlux=%.6e\n",
1263 num_handlers[0],num_handlers[1],num_handlers[2], global_farfield_in_post, global_farfield_out_post);
1264 } else {
1265 // No farfield handlers - zero out the fluxes
1266 user->simCtx->FarFluxInSum = 0.0;
1267 user->simCtx->FarFluxOutSum = 0.0;
1268 }
1269
1270
1271 // =========================================================================
1272 // PRIORITY 2: WALLS
1273 // =========================================================================
1274
1275 memset(num_handlers,0,sizeof(num_handlers));
1276
1277 LOG_ALLOW(LOCAL, LOG_TRACE, " (WALLS): Begin.\n");
1278
1279 // Phase 1: PreStep - Preparation (usually no-op for walls)
1280 for (int i = 0; i < 6; i++) {
1281 BoundaryCondition *handler = user->boundary_faces[i].handler;
1282 if (!handler || handler->priority != BC_PRIORITY_WALL) continue;
1283 if (!handler->PreStep) continue;
1284
1285 num_handlers[0]++;
1286 BCContext ctx = {
1287 .user = user,
1288 .face_id = (BCFace)i,
1289 .global_inflow_sum = &user->simCtx->FluxInSum,
1290 .global_outflow_sum = NULL,
1291 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1292 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1293 };
1294
1295 LOG_ALLOW(LOCAL, LOG_TRACE, " PreStep: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1296 ierr = handler->PreStep(handler, &ctx, NULL, NULL); CHKERRQ(ierr);
1297 }
1298
1299 // No global communication needed for walls
1300
1301 // Phase 2: Apply - Set boundary conditions
1302 for (int i = 0; i < 6; i++) {
1303 BoundaryCondition *handler = user->boundary_faces[i].handler;
1304 if (!handler || handler->priority != BC_PRIORITY_WALL) continue;
1305 if(!handler->Apply) continue; // For example Periodic BCs
1306
1307 num_handlers[1]++;
1308
1309 BCContext ctx = {
1310 .user = user,
1311 .face_id = (BCFace)i,
1312 .global_inflow_sum = &user->simCtx->FluxInSum,
1313 .global_outflow_sum = NULL,
1314 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1315 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1316 };
1317
1318 LOG_ALLOW(LOCAL, LOG_TRACE, " Apply: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1319 ierr = handler->Apply(handler, &ctx); CHKERRQ(ierr);
1320 }
1321
1322 // Phase 3: PostStep - Post-application processing (usually no-op for walls)
1323 for (int i = 0; i < 6; i++) {
1324 BoundaryCondition *handler = user->boundary_faces[i].handler;
1325 if (!handler || handler->priority != BC_PRIORITY_WALL) continue;
1326 if (!handler->PostStep) continue;
1327
1328 num_handlers[2]++;
1329
1330 BCContext ctx = {
1331 .user = user,
1332 .face_id = (BCFace)i,
1333 .global_inflow_sum = &user->simCtx->FluxInSum,
1334 .global_outflow_sum = NULL,
1335 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1336 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1337 };
1338
1339 LOG_ALLOW(LOCAL, LOG_TRACE, " PostStep: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1340 ierr = handler->PostStep(handler, &ctx, NULL, NULL); CHKERRQ(ierr);
1341 }
1342
1343 // No global communication needed for walls
1344
1345 LOG_ALLOW(GLOBAL, LOG_INFO, " (WALLS): %d Prestep(s), %d Application(s), %d Poststep(s) applied.\n",
1346 num_handlers[0],num_handlers[1],num_handlers[2]);
1347
1348
1349 // =========================================================================
1350 // PRIORITY 3: OUTLETS
1351 // =========================================================================
1352
1353 PetscReal local_outflow_pre = 0.0;
1354 PetscReal local_outflow_post = 0.0;
1355 PetscReal global_outflow_pre = 0.0;
1356 PetscReal global_outflow_post = 0.0;
1357 memset(num_handlers,0,sizeof(num_handlers));
1358
1359 LOG_ALLOW(LOCAL, LOG_TRACE, " (OUTLETS): Begin.\n");
1360
1361 // Phase 1: PreStep - Measure uncorrected outflow (from ucat)
1362 for (int i = 0; i < 6; i++) {
1363 BoundaryCondition *handler = user->boundary_faces[i].handler;
1364 if (!handler || handler->priority != BC_PRIORITY_OUTLET) continue;
1365 if (!handler->PreStep) continue;
1366
1367 num_handlers[0]++;
1368 BCContext ctx = {
1369 .user = user,
1370 .face_id = (BCFace)i,
1371 .global_inflow_sum = &user->simCtx->FluxInSum, // From Priority 0
1372 .global_outflow_sum = NULL,
1373 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1374 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1375 };
1376
1377 LOG_ALLOW(LOCAL, LOG_TRACE, " PreStep: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1378 ierr = handler->PreStep(handler, &ctx, NULL, &local_outflow_pre); CHKERRQ(ierr);
1379 }
1380
1381 // Phase 2: Global communication - Get uncorrected outflow sum
1382 ierr = MPI_Allreduce(&local_outflow_pre, &global_outflow_pre, 1, MPIU_REAL,
1383 MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
1384
1385 // Calculate total inflow (inlet + farfield inflow)
1386 PetscReal total_inflow = user->simCtx->FluxInSum + user->simCtx->FarFluxInSum;
1387
1388 LOG_ALLOW(GLOBAL, LOG_DEBUG,
1389 " Uncorrected outflow: %.6e, Total inflow: %.6e (Inlet: %.6e + Farfield: %.6e)\n",
1390 global_outflow_pre, total_inflow, user->simCtx->FluxInSum,
1391 user->simCtx->FarFluxInSum);
1392
1393 // Phase 3: Apply - Set corrected boundary conditions
1394 for (int i = 0; i < 6; i++) {
1395 BoundaryCondition *handler = user->boundary_faces[i].handler;
1396 if (!handler || handler->priority != BC_PRIORITY_OUTLET) continue;
1397 if(!handler->Apply) continue; // For example Periodic BCs
1398
1399 num_handlers[1]++;
1400
1401 BCContext ctx = {
1402 .user = user,
1403 .face_id = (BCFace)i,
1404 .global_inflow_sum = &user->simCtx->FluxInSum, // From Priority 0
1405 .global_outflow_sum = &global_outflow_pre, // From PreStep above
1406 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1407 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1408 };
1409
1410 LOG_ALLOW(LOCAL, LOG_TRACE, " Apply: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1411 ierr = handler->Apply(handler, &ctx); CHKERRQ(ierr);
1412 }
1413
1414 // Phase 4: PostStep - Measure corrected outflow (verification)
1415 for (int i = 0; i < 6; i++) {
1416 BoundaryCondition *handler = user->boundary_faces[i].handler;
1417 if (!handler || handler->priority != BC_PRIORITY_OUTLET) continue;
1418 if (!handler->PostStep) continue;
1419
1420 num_handlers[2]++;
1421
1422 BCContext ctx = {
1423 .user = user,
1424 .face_id = (BCFace)i,
1425 .global_inflow_sum = &user->simCtx->FluxInSum,
1426 .global_outflow_sum = &global_outflow_pre,
1427 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1428 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1429 };
1430
1431 LOG_ALLOW(LOCAL, LOG_TRACE, " PostStep: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1432 ierr = handler->PostStep(handler, &ctx, NULL, &local_outflow_post); CHKERRQ(ierr);
1433 }
1434
1435 // Phase 5: Global communication - Verify conservation
1436 ierr = MPI_Allreduce(&local_outflow_post, &global_outflow_post, 1, MPIU_REAL,
1437 MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
1438
1439 // Store for global reporting.
1440 user->simCtx->FluxOutSum = global_outflow_post;
1441
1442 // Conservation check (compare total outflow vs total inflow)
1443 PetscReal total_outflow = global_outflow_post + user->simCtx->FarFluxOutSum;
1444 PetscReal flux_error = PetscAbsReal(total_outflow - total_inflow);
1445 PetscReal relative_error = (total_inflow > 1e-16) ?
1446 flux_error / total_inflow : flux_error;
1447
1448 LOG_ALLOW(GLOBAL, LOG_INFO,
1449 " (OUTLETS): %d Prestep(s), %d Application(s), %d Poststep(s), FluxOutSum = %.6e\n",
1450 num_handlers[0],num_handlers[1],num_handlers[2], global_outflow_post);
1451 LOG_ALLOW(GLOBAL, LOG_INFO,
1452 " Conservation: Total In=%.6e, Total Out=%.6e, Error=%.3e (%.2e)%%)\n",
1453 total_inflow, total_outflow, flux_error, relative_error * 100.0);
1454
1455 if (relative_error > 1e-6) {
1456 LOG_ALLOW(GLOBAL, LOG_WARNING,
1457 " WARNING: Large mass conservation error (%.2e%%)!\n",
1458 relative_error * 100.0);
1459 }
1460
1461
1462 LOG_ALLOW(LOCAL, LOG_VERBOSE, "Complete.\n");
1463
1464 PROFILE_FUNCTION_END;
1465 PetscFunctionReturn(0);
1466}
1467
1468// =============================================================================
1469//
1470// PRIVATE "LIGHT" EXECUTION ENGINE
1471//
1472// =============================================================================
1473
1474#undef __FUNCT__
1475#define __FUNCT__ "BoundarySystem_RefreshUbcs"
1476/**
1477 * @brief Internal helper implementation: `BoundarySystem_RefreshUbcs()`.
1478 * @details Local to this translation unit.
1479 */
1480PetscErrorCode BoundarySystem_RefreshUbcs(UserCtx *user)
1481{
1482 PetscErrorCode ierr;
1483 PetscFunctionBeginUser;
1484
1485 LOG_ALLOW(GLOBAL, LOG_TRACE, "Refreshing `ubcs` targets for flow-dependent boundaries...\n");
1486
1487 // Loop through all 6 faces of the domain
1488 for (int i = 0; i < 6; i++) {
1489 BoundaryCondition *handler = user->boundary_faces[i].handler;
1490
1491 // THE FILTER:
1492 // This is the core logic. We only act if a handler exists for the face
1493 // AND that handler has explicitly implemented the `UpdateUbcs` method.
1494 if (handler && handler->UpdateUbcs) {
1495
1496 const char *face_name = BCFaceToString((BCFace)i);
1497 LOG_ALLOW(LOCAL, LOG_TRACE, " Calling UpdateUbcs() for handler on Face %s.\n", face_name);
1498
1499 // Prepare the context. For this refresh step, we don't need to pass flux sums.
1500 BCContext ctx = {
1501 .user = user,
1502 .face_id = (BCFace)i,
1503 .global_inflow_sum = NULL,
1504 .global_outflow_sum = NULL,
1505 .global_farfield_inflow_sum = NULL,
1506 .global_farfield_outflow_sum = NULL
1507 };
1508
1509 // Call the handler's specific UpdateUbcs function pointer.
1510 ierr = handler->UpdateUbcs(handler, &ctx); CHKERRQ(ierr);
1511 }
1512 }
1513
1514 PetscFunctionReturn(0);
1515}
1516
1517//================================================================================
1518//
1519// PUBLIC MASTER CLEANUP FUNCTION
1520//
1521//================================================================================
1522#undef __FUNCT__
1523#define __FUNCT__ "BoundarySystem_Destroy"
1524/**
1525 * @brief Implementation of \ref BoundarySystem_Destroy().
1526 * @details Full API contract (arguments, ownership, side effects) is documented with
1527 * the header declaration in `include/Boundaries.h`.
1528 * @see BoundarySystem_Destroy()
1529 */
1530PetscErrorCode BoundarySystem_Destroy(UserCtx *user)
1531{
1532 PetscErrorCode ierr;
1533 PetscFunctionBeginUser;
1534
1535
1536
1537 LOG_ALLOW(GLOBAL, LOG_INFO, "Starting destruction of all boundary handlers. \n");
1538
1539 for (int i = 0; i < 6; i++) {
1540 BoundaryFaceConfig *face_cfg = &user->boundary_faces[i];
1541 const char *face_name = BCFaceToString(face_cfg->face_id);
1542
1543 // --- Step 1: Free the parameter linked list associated with this face ---
1544 if (face_cfg->params) {
1545 LOG_ALLOW(LOCAL, LOG_DEBUG, " Freeing parameter list for Face %d (%s). \n", i, face_name);
1546 FreeBC_ParamList(face_cfg->params);
1547 face_cfg->params = NULL; // Good practice to nullify dangling pointers
1548 }
1549
1550 // --- Step 2: Destroy the handler object itself ---
1551 if (face_cfg->handler) {
1552 const char *handler_name = BCHandlerTypeToString(face_cfg->handler->type);
1553 LOG_ALLOW(LOCAL, LOG_DEBUG, " Destroying handler '%s' on Face %d (%s).\n", handler_name, i, face_name);
1554
1555 // Call the handler's specific cleanup function first, if it exists.
1556 // This will free any memory stored in the handler's private `data` pointer.
1557 if (face_cfg->handler->Destroy) {
1558 ierr = face_cfg->handler->Destroy(face_cfg->handler); CHKERRQ(ierr);
1559 }
1560
1561 // Finally, free the generic BoundaryCondition object itself.
1562 ierr = PetscFree(face_cfg->handler); CHKERRQ(ierr);
1563 face_cfg->handler = NULL;
1564 }
1565 }
1566
1567 LOG_ALLOW(GLOBAL, LOG_INFO, "Destruction complete.\n");
1568 PetscFunctionReturn(0);
1569}
1570
1571#undef __FUNCT__
1572#define __FUNCT__ "TransferPeriodicFieldByDirection"
1573/**
1574 * @brief Internal helper implementation: `TransferPeriodicFieldByDirection()`.
1575 * @details Local to this translation unit.
1576 */
1577PetscErrorCode TransferPeriodicFieldByDirection(UserCtx *user, const char *field_name, char direction)
1578{
1579 PetscErrorCode ierr;
1580 DMDALocalInfo info = user->info;
1581 PetscInt xs = info.xs, xe = info.xs + info.xm;
1582 PetscInt ys = info.ys, ye = info.ys + info.ym;
1583 PetscInt zs = info.zs, ze = info.zs + info.zm;
1584 PetscInt mx = info.mx, my = info.my, mz = info.mz;
1585
1586 // --- Dispatcher to get DM, Vecs, and DoF for the specified field ---
1587 DM dm;
1588 Vec global_vec;
1589 Vec local_vec;
1590 PetscInt dof;
1591 // (This dispatcher is identical to your TransferPeriodicField function)
1592 if (strcmp(field_name, "Ucat") == 0) {
1593 dm = user->fda; global_vec = user->Ucat; local_vec = user->lUcat; dof = 3;
1594 } else if (strcmp(field_name, "P") == 0) {
1595 dm = user->da; global_vec = user->P; local_vec = user->lP; dof = 1;
1596 } else if (strcmp(field_name, "Nvert") == 0) {
1597 dm = user->da; global_vec = user->Nvert; local_vec = user->lNvert; dof = 1;
1598 } else {
1599 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown field name '%s'", field_name);
1600 }
1601
1602 PetscFunctionBeginUser;
1603
1604 // --- Execute the copy logic based on DoF and Direction ---
1605 if (dof == 1) { // --- Handle SCALAR fields (PetscReal) ---
1606 PetscReal ***g_array, ***l_array;
1607 ierr = DMDAVecGetArray(dm, global_vec, &g_array); CHKERRQ(ierr);
1608 ierr = DMDAVecGetArrayRead(dm, local_vec, (void*)&l_array); CHKERRQ(ierr); // Use Read for safety
1609
1610 switch (direction) {
1611 case 'i':
1612 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];
1613 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];
1614 break;
1615 case 'j':
1616 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];
1617 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];
1618 break;
1619 case 'k':
1620 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];
1621 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];
1622 break;
1623 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid direction '%c'", direction);
1624 }
1625 ierr = DMDAVecRestoreArray(dm, global_vec, &g_array); CHKERRQ(ierr);
1626 ierr = DMDAVecRestoreArrayRead(dm, local_vec, (void*)&l_array); CHKERRQ(ierr);
1627
1628 } else if (dof == 3) { // --- Handle VECTOR fields (Cmpnts) ---
1629 Cmpnts ***g_array, ***l_array;
1630 ierr = DMDAVecGetArray(dm, global_vec, &g_array); CHKERRQ(ierr);
1631 ierr = DMDAVecGetArrayRead(dm, local_vec, (void*)&l_array); CHKERRQ(ierr);
1632
1633 switch (direction) {
1634 case 'i':
1635 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];
1636 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];
1637 break;
1638 case 'j':
1639 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];
1640 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];
1641 break;
1642 case 'k':
1643 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];
1644 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];
1645 break;
1646 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid direction '%c'", direction);
1647 }
1648 ierr = DMDAVecRestoreArray(dm, global_vec, &g_array); CHKERRQ(ierr);
1649 ierr = DMDAVecRestoreArrayRead(dm, local_vec, (void*)&l_array); CHKERRQ(ierr);
1650 }
1651
1652 PetscFunctionReturn(0);
1653}
1654
1655#undef __FUNCT__
1656#define __FUNCT__ "TransferPeriodicField"
1657/**
1658 * @brief Internal helper implementation: `TransferPeriodicField()`.
1659 * @details Local to this translation unit.
1660 */
1661PetscErrorCode TransferPeriodicField(UserCtx *user, const char *field_name)
1662{
1663 PetscErrorCode ierr;
1664 DMDALocalInfo info = user->info;
1665 PetscInt xs = info.xs, xe = info.xs + info.xm;
1666 PetscInt ys = info.ys, ye = info.ys + info.ym;
1667 PetscInt zs = info.zs, ze = info.zs + info.zm;
1668 PetscInt mx = info.mx, my = info.my, mz = info.mz;
1669
1670 // --- Local variables to hold the specific details of the chosen field ---
1671 DM dm;
1672 Vec global_vec;
1673 Vec local_vec;
1674 PetscInt dof;
1675
1676 PetscFunctionBeginUser;
1677 PROFILE_FUNCTION_BEGIN;
1678
1679 // --- STEP 1: Dispatcher - Set the specific DM, Vecs, and dof based on field_name ---
1680 // Field-extension note: to add a new periodic field, add one case here and then
1681 // add one call in ApplyPeriodicBCs()/UpdatePeriodicCornerNodes() where appropriate.
1682 if (strcmp(field_name, "Ucat") == 0) {
1683 dm = user->fda;
1684 global_vec = user->Ucat;
1685 local_vec = user->lUcat;
1686 dof = 3;
1687 } else if (strcmp(field_name, "P") == 0) {
1688 dm = user->da;
1689 global_vec = user->P;
1690 local_vec = user->lP;
1691 dof = 1;
1692 } else if (strcmp(field_name, "Nvert") == 0) {
1693 dm = user->da;
1694 global_vec = user->Nvert;
1695 local_vec = user->lNvert;
1696 dof = 1;
1697 } else if (strcmp(field_name, "Eddy Viscosity") == 0) {
1698 dm = user->da;
1699 global_vec = user->Nu_t;
1700 local_vec = user->lNu_t;
1701 dof = 1;
1702 }
1703 /*
1704 // Example for future extension:
1705 else if (strcmp(field_name, "Temperature") == 0) {
1706 dm = user->da; // Assuming Temperature is scalar
1707 global_vec = user->T;
1708 local_vec = user->lT;
1709 dof = 1;
1710 }
1711 */
1712 else {
1713 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown field name '%s' in TransferPeriodicFieldByName.", field_name);
1714 }
1715
1716 LOG_ALLOW(GLOBAL,LOG_TRACE,"Periodic Transform being performed for field: %s with %d DoF.\n",field_name,dof);
1717 // --- STEP 2: Execute the copy logic using the dispatched variables ---
1718 if (dof == 1) { // --- Handle SCALAR fields (PetscReal) ---
1719 PetscReal ***g_array, ***l_array;
1720 ierr = DMDAVecGetArray(dm, global_vec, &g_array); CHKERRQ(ierr);
1721 ierr = DMDAVecGetArray(dm, local_vec, &l_array); CHKERRQ(ierr);
1722
1723 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];
1724 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];
1725 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];
1726 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];
1727 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];
1728 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];
1729
1730 ierr = DMDAVecRestoreArray(dm, global_vec, &g_array); CHKERRQ(ierr);
1731 ierr = DMDAVecRestoreArray(dm, local_vec, &l_array); CHKERRQ(ierr);
1732
1733 } else if (dof == 3) { // --- Handle VECTOR fields (Cmpnts) ---
1734 Cmpnts ***g_array, ***l_array;
1735 ierr = DMDAVecGetArray(dm, global_vec, &g_array); CHKERRQ(ierr);
1736 ierr = DMDAVecGetArray(dm, local_vec, &l_array); CHKERRQ(ierr);
1737
1738 LOG_ALLOW(GLOBAL,LOG_VERBOSE,"Array %s read successfully (Global and Local).\n",field_name);
1739
1740 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];
1741 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];
1742 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];
1743 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];
1744 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];
1745 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];
1746
1747 ierr = DMDAVecRestoreArray(dm, global_vec, &g_array); CHKERRQ(ierr);
1748 ierr = DMDAVecRestoreArray(dm, local_vec, &l_array); CHKERRQ(ierr);
1749 }
1750 else{
1751 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "This function only accepts Fields with 1 or 3 DoF.");
1752 }
1753
1754 PROFILE_FUNCTION_END;
1755 PetscFunctionReturn(0);
1756}
1757
1758#undef __FUNCT__
1759#define __FUNCT__ "TransferPeriodicFaceField"
1760/**
1761 * @brief Internal helper implementation: `TransferPeriodicFaceField()`.
1762 * @details Local to this translation unit.
1763 */
1764PetscErrorCode TransferPeriodicFaceField(UserCtx *user, const char *field_name)
1765{
1766 PetscErrorCode ierr;
1767 DMDALocalInfo info = user->info;
1768 PetscInt gxs = info.gxs, gxe = info.gxs + info.gxm;
1769 PetscInt gys = info.gys, gye = info.gys + info.gym;
1770 PetscInt gzs = info.gzs, gze = info.gzs + info.gzm;
1771 PetscInt mx = info.mx, my = info.my, mz = info.mz;
1772
1773 // --- Dispatcher to get the correct DM, Vec, and DoF for the specified field ---
1774 // Field-extension note: add one case here when a new field needs periodic
1775 // ghost-face transfer (typically for stencil-heavy operators).
1776 DM dm;
1777 Vec local_vec;
1778 PetscInt dof;
1779 // (This dispatcher contains all 17 potential fields)
1780 if (strcmp(field_name, "Ucont") == 0) { dm = user->fda; local_vec = user->lUcont; dof = 3; }
1781 else if (strcmp(field_name, "Csi") == 0) { dm = user->fda; local_vec = user->lCsi; dof = 3; }
1782 else if (strcmp(field_name, "Eta") == 0) { dm = user->fda; local_vec = user->lEta; dof = 3; }
1783 else if (strcmp(field_name, "Zet") == 0) { dm = user->fda; local_vec = user->lZet; dof = 3; }
1784 else if (strcmp(field_name, "ICsi") == 0) { dm = user->fda; local_vec = user->lICsi; dof = 3; }
1785 else if (strcmp(field_name, "IEta") == 0) { dm = user->fda; local_vec = user->lIEta; dof = 3; }
1786 else if (strcmp(field_name, "IZet") == 0) { dm = user->fda; local_vec = user->lIZet; dof = 3; }
1787 else if (strcmp(field_name, "JCsi") == 0) { dm = user->fda; local_vec = user->lJCsi; dof = 3; }
1788 else if (strcmp(field_name, "JEta") == 0) { dm = user->fda; local_vec = user->lJEta; dof = 3; }
1789 else if (strcmp(field_name, "JZet") == 0) { dm = user->fda; local_vec = user->lJZet; dof = 3; }
1790 else if (strcmp(field_name, "KCsi") == 0) { dm = user->fda; local_vec = user->lKCsi; dof = 3; }
1791 else if (strcmp(field_name, "KEta") == 0) { dm = user->fda; local_vec = user->lKEta; dof = 3; }
1792 else if (strcmp(field_name, "KZet") == 0) { dm = user->fda; local_vec = user->lKZet; dof = 3; }
1793 else if (strcmp(field_name, "Aj") == 0) { dm = user->da; local_vec = user->lAj; dof = 1; }
1794 else if (strcmp(field_name, "IAj") == 0) { dm = user->da; local_vec = user->lIAj; dof = 1; }
1795 else if (strcmp(field_name, "JAj") == 0) { dm = user->da; local_vec = user->lJAj; dof = 1; }
1796 else if (strcmp(field_name, "KAj") == 0) { dm = user->da; local_vec = user->lKAj; dof = 1; }
1797 else {
1798 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown field name '%s' in TransferPeriodicFaceField.", field_name);
1799 }
1800
1801 PetscFunctionBeginUser;
1802
1803 void *l_array_ptr;
1804 ierr = DMDAVecGetArray(dm, local_vec, &l_array_ptr); CHKERRQ(ierr);
1805
1806 // --- I-DIRECTION ---
1807 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC) {
1808 for (PetscInt k=gzs; k<gze; k++) for (PetscInt j=gys; j<gye; j++) {
1809 if (dof == 1) {
1810 PetscReal ***arr = (PetscReal***)l_array_ptr;
1811 arr[k][j][0] = arr[k][j][mx-2];
1812 arr[k][j][-1] = arr[k][j][mx-3];
1813 } else {
1814 Cmpnts ***arr = (Cmpnts***)l_array_ptr;
1815 arr[k][j][0] = arr[k][j][mx-2];
1816 arr[k][j][-1] = arr[k][j][mx-3];
1817 }
1818 }
1819 }
1820 if (user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC) {
1821 for (PetscInt k=gzs; k<gze; k++) for (PetscInt j=gys; j<gye; j++) {
1822 if (dof == 1) {
1823 PetscReal ***arr = (PetscReal***)l_array_ptr;
1824 arr[k][j][mx-1] = arr[k][j][1];
1825 arr[k][j][mx] = arr[k][j][2];
1826 } else {
1827 Cmpnts ***arr = (Cmpnts***)l_array_ptr;
1828 arr[k][j][mx-1] = arr[k][j][1];
1829 arr[k][j][mx] = arr[k][j][2];
1830 }
1831 }
1832 }
1833
1834 // --- J-DIRECTION ---
1835 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC) {
1836 for (PetscInt k=gzs; k<gze; k++) for (PetscInt i=gxs; i<gxe; i++) {
1837 if (dof == 1) {
1838 PetscReal ***arr = (PetscReal***)l_array_ptr;
1839 arr[k][0][i] = arr[k][my-2][i];
1840 arr[k][-1][i] = arr[k][my-3][i];
1841 } else {
1842 Cmpnts ***arr = (Cmpnts***)l_array_ptr;
1843 arr[k][0][i] = arr[k][my-2][i];
1844 arr[k][-1][i] = arr[k][my-3][i];
1845 }
1846 }
1847 }
1848 if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC) {
1849 for (PetscInt k=gzs; k<gze; k++) for (PetscInt i=gxs; i<gxe; i++) {
1850 if (dof == 1) {
1851 PetscReal ***arr = (PetscReal***)l_array_ptr;
1852 arr[k][my-1][i] = arr[k][1][i];
1853 arr[k][my][i] = arr[k][2][i];
1854 } else {
1855 Cmpnts ***arr = (Cmpnts***)l_array_ptr;
1856 arr[k][my-1][i] = arr[k][1][i];
1857 arr[k][my][i] = arr[k][2][i];
1858 }
1859 }
1860 }
1861
1862 // --- K-DIRECTION ---
1863 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC) {
1864 for (PetscInt j=gys; j<gye; j++) for (PetscInt i=gxs; i<gxe; i++) {
1865 if (dof == 1) {
1866 PetscReal ***arr = (PetscReal***)l_array_ptr;
1867 arr[0][j][i] = arr[mz-2][j][i];
1868 arr[-1][j][i] = arr[mz-3][j][i];
1869 } else {
1870 Cmpnts ***arr = (Cmpnts***)l_array_ptr;
1871 arr[0][j][i] = arr[mz-2][j][i];
1872 arr[-1][j][i] = arr[mz-3][j][i];
1873 }
1874 }
1875 }
1876 if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC) {
1877 for (PetscInt j=gys; j<gye; j++) for (PetscInt i=gxs; i<gxe; i++) {
1878 if (dof == 1) {
1879 PetscReal ***arr = (PetscReal***)l_array_ptr;
1880 arr[mz-1][j][i] = arr[1][j][i];
1881 arr[mz][j][i] = arr[2][j][i];
1882 } else {
1883 Cmpnts ***arr = (Cmpnts***)l_array_ptr;
1884 arr[mz-1][j][i] = arr[1][j][i];
1885 arr[mz][j][i] = arr[2][j][i];
1886 }
1887 }
1888 }
1889
1890 ierr = DMDAVecRestoreArray(dm, local_vec, &l_array_ptr); CHKERRQ(ierr);
1891 PetscFunctionReturn(0);
1892}
1893
1894#undef __FUNCT__
1895#define __FUNCT__ "ApplyMetricsPeriodicBCs"
1896/**
1897 * @brief Internal helper implementation: `ApplyMetricsPeriodicBCs()`.
1898 * @details Local to this translation unit.
1899 */
1900PetscErrorCode ApplyMetricsPeriodicBCs(UserCtx *user)
1901{
1902 PetscErrorCode ierr;
1903 PetscFunctionBeginUser;
1904 PROFILE_FUNCTION_BEGIN;
1905
1906 const char* metric_fields[] = {
1907 "Csi", "Eta", "Zet", "ICsi", "JCsi", "KCsi", "IEta", "JEta", "KEta",
1908 "IZet", "JZet", "KZet", "Aj", "IAj", "JAj", "KAj"
1909 };
1910 PetscInt num_fields = sizeof(metric_fields) / sizeof(metric_fields[0]);
1911
1912 for (PetscInt i = 0; i < num_fields; i++) {
1913 ierr = TransferPeriodicFaceField(user, metric_fields[i]); CHKERRQ(ierr);
1914 LOG_ALLOW(GLOBAL,LOG_TRACE,"Periodic Transfer complete for %s.\n",metric_fields[i]);
1915 }
1916
1917 PROFILE_FUNCTION_END;
1918 PetscFunctionReturn(0);
1919}
1920
1921#undef __FUNCT__
1922#define __FUNCT__ "ApplyPeriodicBCs"
1923/**
1924 * @brief Internal helper implementation: `ApplyPeriodicBCs()`.
1925 * @details Local to this translation unit.
1926 */
1927PetscErrorCode ApplyPeriodicBCs(UserCtx *user)
1928{
1929 PetscErrorCode ierr;
1930 PetscBool is_any_periodic = PETSC_FALSE;
1931
1932 PetscFunctionBeginUser;
1933
1934 PROFILE_FUNCTION_BEGIN;
1935
1936 for (int i = 0; i < 6; i++) {
1937 if (user->boundary_faces[i].mathematical_type == PERIODIC) {
1938 is_any_periodic = PETSC_TRUE;
1939 break;
1940 }
1941 }
1942
1943 if (!is_any_periodic) {
1944 LOG_ALLOW(GLOBAL,LOG_TRACE, "No periodic boundaries defined; skipping ApplyPeriodicBCs.\n");
1945 PROFILE_FUNCTION_END;
1946 PetscFunctionReturn(0);
1947 }
1948
1949 LOG_ALLOW(GLOBAL, LOG_TRACE, "Applying periodic boundary conditions for all fields.\n");
1950
1951 // STEP 1: Perform the collective communication for periodic core fields.
1952 ierr = UpdateLocalGhosts(user, "Ucat"); CHKERRQ(ierr);
1953 ierr = UpdateLocalGhosts(user, "P"); CHKERRQ(ierr);
1954 ierr = UpdateLocalGhosts(user, "Nvert"); CHKERRQ(ierr);
1955 ierr = UpdateLocalGhosts(user, "Ucont"); CHKERRQ(ierr);
1956 /* if (user->solve_temperature) { ierr = UpdateLocalGhosts(user, "Temperature"); CHKERRQ(ierr); } */
1957
1958 // STEP 2: Call the generic copy routine for each face-centered/cell-centered field.
1959 ierr = TransferPeriodicField(user, "Ucat"); CHKERRQ(ierr);
1960 ierr = TransferPeriodicField(user, "P"); CHKERRQ(ierr);
1961 ierr = TransferPeriodicField(user, "Nvert"); CHKERRQ(ierr);
1962
1963 // STEP 3: Update contravariant flux periodic ghosts and enforce strict seam consistency.
1964 // This keeps the staggered Ucont field robust for QUICK-like stencils and periodic flux closure.
1965 ierr = ApplyUcontPeriodicBCs(user); CHKERRQ(ierr);
1966 ierr = EnforceUcontPeriodicity(user); CHKERRQ(ierr);
1967 ierr = UpdateLocalGhosts(user, "Ucont"); CHKERRQ(ierr);
1968
1969 // FUTURE EXTENSION: Adding a new scalar field like Temperature is now trivial.
1970 /*
1971 if (user->solve_temperature) {
1972 ierr = TransferPeriodicField(user, "Temperature"); CHKERRQ(ierr);
1973 }
1974 */
1975
1976 PROFILE_FUNCTION_END;
1977 PetscFunctionReturn(0);
1978}
1979
1980#undef __FUNCT__
1981#define __FUNCT__ "ApplyUcontPeriodicBCs"
1982/**
1983 * @brief Implementation of \ref ApplyUcontPeriodicBCs().
1984 * @details Full API contract (arguments, ownership, side effects) is documented with
1985 * the header declaration in `include/Boundaries.h`.
1986 * @see ApplyUcontPeriodicBCs()
1987 */
1988PetscErrorCode ApplyUcontPeriodicBCs(UserCtx *user)
1989{
1990 PetscErrorCode ierr;
1991 PetscFunctionBeginUser;
1992 PROFILE_FUNCTION_BEGIN;
1993
1994 // Update local periodic ghost faces (two-cells deep) for Ucont.
1995 ierr = TransferPeriodicFaceField(user, "Ucont"); CHKERRQ(ierr);
1996
1997 // Commit periodic-face updates back to global Ucont so follow-up routines
1998 // (including EnforceUcontPeriodicity) can safely refresh from global state.
1999 ierr = DMLocalToGlobalBegin(user->fda, user->lUcont, INSERT_VALUES, user->Ucont); CHKERRQ(ierr);
2000 ierr = DMLocalToGlobalEnd(user->fda, user->lUcont, INSERT_VALUES, user->Ucont); CHKERRQ(ierr);
2001
2002 PROFILE_FUNCTION_END;
2003 PetscFunctionReturn(0);
2004}
2005
2006#undef __FUNCT__
2007#define __FUNCT__ "EnforceUcontPeriodicity"
2008/**
2009 * @brief Internal helper implementation: `EnforceUcontPeriodicity()`.
2010 * @details Local to this translation unit.
2011 */
2012PetscErrorCode EnforceUcontPeriodicity(UserCtx *user)
2013{
2014 PetscErrorCode ierr;
2015 DMDALocalInfo info = user->info;
2016 PetscInt gxs = info.gxs, gxe = info.gxs + info.gxm;
2017 PetscInt gys = info.gys, gye = info.gys + info.gym;
2018 PetscInt gzs = info.gzs, gze = info.gzs + info.gzm;
2019 PetscInt mx = info.mx, my = info.my, mz = info.mz;
2020 Cmpnts ***ucont;
2021
2022 PetscFunctionBeginUser;
2023 PROFILE_FUNCTION_BEGIN;
2024
2025 // STEP 1: Ensure local ghost cells are up-to-date with current global state.
2026 ierr = DMGlobalToLocalBegin(user->fda, user->Ucont, INSERT_VALUES, user->lUcont); CHKERRQ(ierr);
2027 ierr = DMGlobalToLocalEnd(user->fda, user->Ucont, INSERT_VALUES, user->lUcont); CHKERRQ(ierr);
2028
2029 ierr = DMDAVecGetArray(user->fda, user->lUcont, &ucont); CHKERRQ(ierr);
2030
2031 // STEP 2: Perform the component-wise copy from ghost cells to the last interior faces.
2032 if (user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC) {
2033 for (PetscInt k=gzs; k<gze; k++) for (PetscInt j=gys; j<gye; j++) {
2034 ucont[k][j][mx-2].x = ucont[k][j][mx].x; // Note: ucont[mx] is ghost, gets value from ucont[0]
2035 }
2036 }
2037 if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC) {
2038 for (PetscInt k=gzs; k<gze; k++) for (PetscInt i=gxs; i<gxe; i++) {
2039 ucont[k][my-2][i].y = ucont[k][my][i].y; // Note: ucont[my] is ghost, gets value from ucont[0]
2040 }
2041 }
2042 if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC) {
2043 for (PetscInt j=gys; j<gye; j++) for (PetscInt i=gxs; i<gxe; i++) {
2044 ucont[mz-2][j][i].z = ucont[mz][j][i].z; // Note: ucont[mz] is ghost, gets value from ucont[0]
2045 }
2046 }
2047
2048 ierr = DMDAVecRestoreArray(user->fda, user->lUcont, &ucont); CHKERRQ(ierr);
2049
2050 // STEP 3: Communicate the changes made to the interior back to the global vector.
2051 ierr = DMLocalToGlobalBegin(user->fda, user->lUcont, INSERT_VALUES, user->Ucont); CHKERRQ(ierr);
2052 ierr = DMLocalToGlobalEnd(user->fda, user->lUcont, INSERT_VALUES, user->Ucont); CHKERRQ(ierr);
2053
2054 PROFILE_FUNCTION_END;
2055 PetscFunctionReturn(0);
2056}
2057
2058#undef __FUNCT__
2059#define __FUNCT__ "UpdateDummyCells"
2060/**
2061 * @brief Internal helper implementation: `UpdateDummyCells()`.
2062 * @details Local to this translation unit.
2063 */
2064PetscErrorCode UpdateDummyCells(UserCtx *user)
2065{
2066 PetscErrorCode ierr;
2067 DM fda = user->fda;
2068 DMDALocalInfo info = user->info;
2069 PetscInt xs = info.xs, xe = info.xs + info.xm;
2070 PetscInt ys = info.ys, ye = info.ys + info.ym;
2071 PetscInt zs = info.zs, ze = info.zs + info.zm;
2072 PetscInt mx = info.mx, my = info.my, mz = info.mz;
2073
2074 // --- Calculate shrunken loop ranges to avoid edges and corners ---
2075 PetscInt lxs = xs, lxe = xe;
2076 PetscInt lys = ys, lye = ye;
2077 PetscInt lzs = zs, lze = ze;
2078
2079 if (xs == 0) lxs = xs + 1;
2080 if (ys == 0) lys = ys + 1;
2081 if (zs == 0) lzs = zs + 1;
2082
2083 if (xe == mx) lxe = xe - 1;
2084 if (ye == my) lye = ye - 1;
2085 if (ze == mz) lze = ze - 1;
2086
2087 Cmpnts ***ucat, ***ubcs;
2088 PetscFunctionBeginUser;
2089
2090 ierr = DMDAVecGetArray(fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
2091 ierr = DMDAVecGetArray(fda, user->Ucat, &ucat); CHKERRQ(ierr);
2092
2093 // -X Face
2094 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC && xs == 0) {
2095 for (PetscInt k = lzs; k < lze; k++) for (PetscInt j = lys; j < lye; j++) {
2096 ucat[k][j][xs].x = 2.0 * ubcs[k][j][xs].x - ucat[k][j][xs + 1].x;
2097 ucat[k][j][xs].y = 2.0 * ubcs[k][j][xs].y - ucat[k][j][xs + 1].y;
2098 ucat[k][j][xs].z = 2.0 * ubcs[k][j][xs].z - ucat[k][j][xs + 1].z;
2099 }
2100 }
2101 // +X Face
2102 if (user->boundary_faces[BC_FACE_POS_X].mathematical_type != PERIODIC && xe == mx) {
2103 for (PetscInt k = lzs; k < lze; k++) for (PetscInt j = lys; j < lye; j++) {
2104 ucat[k][j][xe-1].x = 2.0 * ubcs[k][j][xe-1].x - ucat[k][j][xe - 2].x;
2105 ucat[k][j][xe-1].y = 2.0 * ubcs[k][j][xe-1].y - ucat[k][j][xe - 2].y;
2106 ucat[k][j][xe-1].z = 2.0 * ubcs[k][j][xe-1].z - ucat[k][j][xe - 2].z;
2107 }
2108 }
2109
2110 // -Y Face
2111 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC && ys == 0) {
2112 for (PetscInt k = lzs; k < lze; k++) for (PetscInt i = lxs; i < lxe; i++) {
2113 ucat[k][ys][i].x = 2.0 * ubcs[k][ys][i].x - ucat[k][ys + 1][i].x;
2114 ucat[k][ys][i].y = 2.0 * ubcs[k][ys][i].y - ucat[k][ys + 1][i].y;
2115 ucat[k][ys][i].z = 2.0 * ubcs[k][ys][i].z - ucat[k][ys + 1][i].z;
2116 }
2117 }
2118 // +Y Face
2119 if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type != PERIODIC && ye == my) {
2120 for (PetscInt k = lzs; k < lze; k++) for (PetscInt i = lxs; i < lxe; i++) {
2121 ucat[k][ye-1][i].x = 2.0 * ubcs[k][ye-1][i].x - ucat[k][ye-2][i].x;
2122 ucat[k][ye-1][i].y = 2.0 * ubcs[k][ye-1][i].y - ucat[k][ye-2][i].y;
2123 ucat[k][ye-1][i].z = 2.0 * ubcs[k][ye-1][i].z - ucat[k][ye-2][i].z;
2124 }
2125 }
2126
2127 // -Z Face
2128 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC && zs == 0) {
2129 for (PetscInt j = lys; j < lye; j++) for (PetscInt i = lxs; i < lxe; i++) {
2130 ucat[zs][j][i].x = 2.0 * ubcs[zs][j][i].x - ucat[zs + 1][j][i].x;
2131 ucat[zs][j][i].y = 2.0 * ubcs[zs][j][i].y - ucat[zs + 1][j][i].y;
2132 ucat[zs][j][i].z = 2.0 * ubcs[zs][j][i].z - ucat[zs + 1][j][i].z;
2133 }
2134 }
2135 // +Z Face
2136 if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type != PERIODIC && ze == mz) {
2137 for (PetscInt j = lys; j < lye; j++) for (PetscInt i = lxs; i < lxe; i++) {
2138 ucat[ze-1][j][i].x = 2.0 * ubcs[ze-1][j][i].x - ucat[ze-2][j][i].x;
2139 ucat[ze-1][j][i].y = 2.0 * ubcs[ze-1][j][i].y - ucat[ze-2][j][i].y;
2140 ucat[ze-1][j][i].z = 2.0 * ubcs[ze-1][j][i].z - ucat[ze-2][j][i].z;
2141 }
2142 }
2143
2144 ierr = DMDAVecRestoreArray(fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
2145 ierr = DMDAVecRestoreArray(fda, user->Ucat, &ucat); CHKERRQ(ierr);
2146
2147 PetscFunctionReturn(0);
2148}
2149
2150#undef __FUNCT__
2151#define __FUNCT__ "UpdateCornerNodes"
2152/**
2153 * @brief Internal helper implementation: `UpdateCornerNodes()`.
2154 * @details Local to this translation unit.
2155 */
2156PetscErrorCode UpdateCornerNodes(UserCtx *user)
2157{
2158 PetscErrorCode ierr;
2159 DM da = user->da, fda = user->fda;
2160 DMDALocalInfo info = user->info;
2161 PetscInt xs = info.xs, xe = info.xs + info.xm;
2162 PetscInt ys = info.ys, ye = info.ys + info.ym;
2163 PetscInt zs = info.zs, ze = info.zs + info.zm;
2164 PetscInt mx = info.mx, my = info.my, mz = info.mz;
2165
2166 Cmpnts ***ucat;
2167 PetscReal ***p;
2168
2169 PetscFunctionBeginUser;
2170
2171 ierr = DMDAVecGetArray(fda, user->Ucat, &ucat); CHKERRQ(ierr);
2172 ierr = DMDAVecGetArray(da, user->P, &p); CHKERRQ(ierr);
2173
2174 // --- Update Edges and Corners by Averaging ---
2175 // The order of these blocks ensures that corners (where 3 faces meet) are
2176 // computed using data from edges (where 2 faces meet), which are computed first.
2177// Edges connected to the -Z face (k=zs)
2178 if (zs == 0) {
2179 if (xs == 0) {
2180 for (PetscInt j = ys; j < ye; j++) {
2181 p[zs][j][xs] = 0.5 * (p[zs+1][j][xs] + p[zs][j][xs+1]);
2182 ucat[zs][j][xs].x = 0.5 * (ucat[zs+1][j][xs].x + ucat[zs][j][xs+1].x);
2183 ucat[zs][j][xs].y = 0.5 * (ucat[zs+1][j][xs].y + ucat[zs][j][xs+1].y);
2184 ucat[zs][j][xs].z = 0.5 * (ucat[zs+1][j][xs].z + ucat[zs][j][xs+1].z);
2185 }
2186 }
2187 if (xe == mx) {
2188 for (PetscInt j = ys; j < ye; j++) {
2189 p[zs][j][mx-1] = 0.5 * (p[zs+1][j][mx-1] + p[zs][j][mx-2]);
2190 ucat[zs][j][mx-1].x = 0.5 * (ucat[zs+1][j][mx-1].x + ucat[zs][j][mx-2].x);
2191 ucat[zs][j][mx-1].y = 0.5 * (ucat[zs+1][j][mx-1].y + ucat[zs][j][mx-2].y);
2192 ucat[zs][j][mx-1].z = 0.5 * (ucat[zs+1][j][mx-1].z + ucat[zs][j][mx-2].z);
2193 }
2194 }
2195 if (ys == 0) {
2196 for (PetscInt i = xs; i < xe; i++) {
2197 p[zs][ys][i] = 0.5 * (p[zs+1][ys][i] + p[zs][ys+1][i]);
2198 ucat[zs][ys][i].x = 0.5 * (ucat[zs+1][ys][i].x + ucat[zs][ys+1][i].x);
2199 ucat[zs][ys][i].y = 0.5 * (ucat[zs+1][ys][i].y + ucat[zs][ys+1][i].y);
2200 ucat[zs][ys][i].z = 0.5 * (ucat[zs+1][ys][i].z + ucat[zs][ys+1][i].z);
2201 }
2202 }
2203 if (ye == my) {
2204 for (PetscInt i = xs; i < xe; i++) {
2205 p[zs][my-1][i] = 0.5 * (p[zs+1][my-1][i] + p[zs][my-2][i]);
2206 ucat[zs][my-1][i].x = 0.5 * (ucat[zs+1][my-1][i].x + ucat[zs][my-2][i].x);
2207 ucat[zs][my-1][i].y = 0.5 * (ucat[zs+1][my-1][i].y + ucat[zs][my-2][i].y);
2208 ucat[zs][my-1][i].z = 0.5 * (ucat[zs+1][my-1][i].z + ucat[zs][my-2][i].z);
2209 }
2210 }
2211 }
2212
2213 // Edges connected to the +Z face (k=ze-1)
2214 if (ze == mz) {
2215 if (xs == 0) {
2216 for (PetscInt j = ys; j < ye; j++) {
2217 p[mz-1][j][xs] = 0.5 * (p[mz-2][j][xs] + p[mz-1][j][xs+1]);
2218 ucat[mz-1][j][xs].x = 0.5 * (ucat[mz-2][j][xs].x + ucat[mz-1][j][xs+1].x);
2219 ucat[mz-1][j][xs].y = 0.5 * (ucat[mz-2][j][xs].y + ucat[mz-1][j][xs+1].y);
2220 ucat[mz-1][j][xs].z = 0.5 * (ucat[mz-2][j][xs].z + ucat[mz-1][j][xs+1].z);
2221 }
2222 }
2223 if (xe == mx) {
2224 for (PetscInt j = ys; j < ye; j++) {
2225 p[mz-1][j][mx-1] = 0.5 * (p[mz-2][j][mx-1] + p[mz-1][j][mx-2]);
2226 ucat[mz-1][j][mx-1].x = 0.5 * (ucat[mz-2][j][mx-1].x + ucat[mz-1][j][mx-2].x);
2227 ucat[mz-1][j][mx-1].y = 0.5 * (ucat[mz-2][j][mx-1].y + ucat[mz-1][j][mx-2].y);
2228 ucat[mz-1][j][mx-1].z = 0.5 * (ucat[mz-2][j][mx-1].z + ucat[mz-1][j][mx-2].z);
2229 }
2230 }
2231 if (ys == 0) {
2232 for (PetscInt i = xs; i < xe; i++) {
2233 p[mz-1][ys][i] = 0.5 * (p[mz-2][ys][i] + p[mz-1][ys+1][i]);
2234 ucat[mz-1][ys][i].x = 0.5 * (ucat[mz-2][ys][i].x + ucat[mz-1][ys+1][i].x);
2235 ucat[mz-1][ys][i].y = 0.5 * (ucat[mz-2][ys][i].y + ucat[mz-1][ys+1][i].y);
2236 ucat[mz-1][ys][i].z = 0.5 * (ucat[mz-2][ys][i].z + ucat[mz-1][ys+1][i].z);
2237 }
2238 }
2239 if (ye == my) {
2240 for (PetscInt i = xs; i < xe; i++) {
2241 p[mz-1][my-1][i] = 0.5 * (p[mz-2][my-1][i] + p[mz-1][my-2][i]);
2242 ucat[mz-1][my-1][i].x = 0.5 * (ucat[mz-2][my-1][i].x + ucat[mz-1][my-2][i].x);
2243 ucat[mz-1][my-1][i].y = 0.5 * (ucat[mz-2][my-1][i].y + ucat[mz-1][my-2][i].y);
2244 ucat[mz-1][my-1][i].z = 0.5 * (ucat[mz-2][my-1][i].z + ucat[mz-1][my-2][i].z);
2245 }
2246 }
2247 }
2248
2249 // Remaining edges on the XY plane (that are not on Z faces)
2250 if (ys == 0) {
2251 if (xs == 0) {
2252 for (PetscInt k = zs; k < ze; k++) {
2253 p[k][ys][xs] = 0.5 * (p[k][ys+1][xs] + p[k][ys][xs+1]);
2254 ucat[k][ys][xs].x = 0.5 * (ucat[k][ys+1][xs].x + ucat[k][ys][xs+1].x);
2255 ucat[k][ys][xs].y = 0.5 * (ucat[k][ys+1][xs].y + ucat[k][ys][xs+1].y);
2256 ucat[k][ys][xs].z = 0.5 * (ucat[k][ys+1][xs].z + ucat[k][ys][xs+1].z);
2257 }
2258 }
2259 if (xe == mx) {
2260 for (PetscInt k = zs; k < ze; k++) {
2261 p[k][ys][mx-1] = 0.5 * (p[k][ys+1][mx-1] + p[k][ys][mx-2]);
2262 ucat[k][ys][mx-1].x = 0.5 * (ucat[k][ys+1][mx-1].x + ucat[k][ys][mx-2].x);
2263 ucat[k][ys][mx-1].y = 0.5 * (ucat[k][ys+1][mx-1].y + ucat[k][ys][mx-2].y);
2264 ucat[k][ys][mx-1].z = 0.5 * (ucat[k][ys+1][mx-1].z + ucat[k][ys][mx-2].z);
2265 }
2266 }
2267 }
2268
2269 if (ye == my) {
2270 if (xs == 0) {
2271 for (PetscInt k = zs; k < ze; k++) {
2272 p[k][my-1][xs] = 0.5 * (p[k][my-2][xs] + p[k][my-1][xs+1]);
2273 ucat[k][my-1][xs].x = 0.5 * (ucat[k][my-2][xs].x + ucat[k][my-1][xs+1].x);
2274 ucat[k][my-1][xs].y = 0.5 * (ucat[k][my-2][xs].y + ucat[k][my-1][xs+1].y);
2275 ucat[k][my-1][xs].z = 0.5 * (ucat[k][my-2][xs].z + ucat[k][my-1][xs+1].z);
2276 }
2277 }
2278 if (xe == mx) {
2279 for (PetscInt k = zs; k < ze; k++) {
2280 p[k][my-1][mx-1] = 0.5 * (p[k][my-2][mx-1] + p[k][my-1][mx-2]);
2281 ucat[k][my-1][mx-1].x = 0.5 * (ucat[k][my-2][mx-1].x + ucat[k][my-1][mx-2].x);
2282 ucat[k][my-1][mx-1].y = 0.5 * (ucat[k][my-2][mx-1].y + ucat[k][my-1][mx-2].y);
2283 ucat[k][my-1][mx-1].z = 0.5 * (ucat[k][my-2][mx-1].z + ucat[k][my-1][mx-2].z);
2284 }
2285 }
2286 }
2287
2288 ierr = DMDAVecRestoreArray(fda, user->Ucat, &ucat); CHKERRQ(ierr);
2289 ierr = DMDAVecRestoreArray(da, user->P, &p); CHKERRQ(ierr);
2290
2291 PetscFunctionReturn(0);
2292}
2293
2294#undef __FUNCT__
2295#define __FUNCT__ "UpdatePeriodicCornerNodes"
2296/**
2297 * @brief Internal helper implementation: `UpdatePeriodicCornerNodes()`.
2298 * @details Local to this translation unit.
2299 */
2300PetscErrorCode UpdatePeriodicCornerNodes(UserCtx *user, PetscInt num_fields, const char* field_names[])
2301{
2302 PetscErrorCode ierr;
2303 PetscFunctionBeginUser;
2304
2305 if (num_fields == 0) PetscFunctionReturn(0);
2306
2307 // --- I-DIRECTION ---
2308 for (PetscInt i = 0; i < num_fields; i++) {
2309 ierr = TransferPeriodicFieldByDirection(user, field_names[i], 'i'); CHKERRQ(ierr);
2310 }
2311 // --- SYNC ---
2312 for (PetscInt i = 0; i < num_fields; i++) {
2313 ierr = UpdateLocalGhosts(user, field_names[i]); CHKERRQ(ierr);
2314 }
2315
2316 // --- J-DIRECTION ---
2317 for (PetscInt i = 0; i < num_fields; i++) {
2318 ierr = TransferPeriodicFieldByDirection(user, field_names[i], 'j'); CHKERRQ(ierr);
2319 }
2320 // --- SYNC ---
2321 for (PetscInt i = 0; i < num_fields; i++) {
2322 ierr = UpdateLocalGhosts(user, field_names[i]); CHKERRQ(ierr);
2323 }
2324
2325 // --- K-DIRECTION ---
2326 for (PetscInt i = 0; i < num_fields; i++) {
2327 ierr = TransferPeriodicFieldByDirection(user, field_names[i], 'k'); CHKERRQ(ierr);
2328 }
2329 // --- FINAL SYNC ---
2330 for (PetscInt i = 0; i < num_fields; i++) {
2331 ierr = UpdateLocalGhosts(user, field_names[i]); CHKERRQ(ierr);
2332 }
2333
2334 PetscFunctionReturn(0);
2335}
2336
2337#undef __FUNCT__
2338#define __FUNCT__ "ApplyWallFunction"
2339/**
2340 * @brief Internal helper implementation: `ApplyWallFunction()`.
2341 * @details Local to this translation unit.
2342 */
2343PetscErrorCode ApplyWallFunction(UserCtx *user)
2344{
2345 PetscErrorCode ierr;
2346 SimCtx *simCtx = user->simCtx;
2347 DMDALocalInfo *info = &user->info;
2348
2349 PetscFunctionBeginUser;
2350
2351 // =========================================================================
2352 // STEP 0: Early exit if wall functions are disabled
2353 // =========================================================================
2354 if (!simCtx->wallfunction) {
2355 PetscFunctionReturn(0);
2356 }
2357
2358 LOG_ALLOW(LOCAL, LOG_DEBUG, "Processing wall function boundaries.\n");
2359
2360 // =========================================================================
2361 // STEP 1: Get read/write access to all necessary field arrays
2362 // =========================================================================
2363 Cmpnts ***velocity_cartesian; // Cartesian velocity (modified)
2364 Cmpnts ***velocity_contravariant; // Contravariant velocity (set to zero at walls)
2365 Cmpnts ***velocity_boundary; // Boundary condition velocity (kept at zero)
2366 Cmpnts ***csi, ***eta, ***zet; // Metric tensor components (face normals)
2367 PetscReal ***node_vertex_flag; // Fluid/solid indicator (0=fluid, 1=solid)
2368 PetscReal ***cell_jacobian; // Grid Jacobian (1/volume)
2369 PetscReal ***friction_velocity; // u_tau (friction velocity field)
2370
2371 ierr = DMDAVecGetArray(user->fda, user->Ucat, &velocity_cartesian); CHKERRQ(ierr);
2372 ierr = DMDAVecGetArray(user->fda, user->Ucont, &velocity_contravariant); CHKERRQ(ierr);
2373 ierr = DMDAVecGetArray(user->fda, user->Bcs.Ubcs, &velocity_boundary); CHKERRQ(ierr);
2374 ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
2375 ierr = DMDAVecGetArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
2376 ierr = DMDAVecGetArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
2377 ierr = DMDAVecGetArrayRead(user->da, user->lNvert, (const PetscReal***)&node_vertex_flag); CHKERRQ(ierr);
2378 ierr = DMDAVecGetArrayRead(user->da, user->lAj, (const PetscReal***)&cell_jacobian); CHKERRQ(ierr);
2379 ierr = DMDAVecGetArray(user->da, user->lFriction_Velocity, &friction_velocity); CHKERRQ(ierr);
2380
2381 // =========================================================================
2382 // STEP 2: Define loop bounds (owned portion of the grid for this MPI rank)
2383 // =========================================================================
2384 PetscInt grid_start_i = info->xs, grid_end_i = info->xs + info->xm;
2385 PetscInt grid_start_j = info->ys, grid_end_j = info->ys + info->ym;
2386 PetscInt grid_start_k = info->zs, grid_end_k = info->zs + info->zm;
2387 PetscInt grid_size_i = info->mx, grid_size_j = info->my, grid_size_k = info->mz;
2388
2389 // Shrunken loop bounds: exclude domain edges and corners to avoid double-counting
2390 PetscInt loop_start_i = grid_start_i, loop_end_i = grid_end_i;
2391 PetscInt loop_start_j = grid_start_j, loop_end_j = grid_end_j;
2392 PetscInt loop_start_k = grid_start_k, loop_end_k = grid_end_k;
2393
2394 if (grid_start_i == 0) loop_start_i = grid_start_i + 1;
2395 if (grid_end_i == grid_size_i) loop_end_i = grid_end_i - 1;
2396 if (grid_start_j == 0) loop_start_j = grid_start_j + 1;
2397 if (grid_end_j == grid_size_j) loop_end_j = grid_end_j - 1;
2398 if (grid_start_k == 0) loop_start_k = grid_start_k + 1;
2399 if (grid_end_k == grid_size_k) loop_end_k = grid_end_k - 1;
2400
2401 // Wall roughness parameter (currently smooth wall)
2402 const PetscReal wall_roughness_height = 1.e-16;
2403
2404 // =========================================================================
2405 // STEP 3: Process each of the 6 domain faces
2406 // =========================================================================
2407 for (int face_index = 0; face_index < 6; face_index++) {
2408 BCFace current_face_id = (BCFace)face_index;
2409 BoundaryFaceConfig *face_config = &user->boundary_faces[current_face_id];
2410
2411 // Only process faces that are mathematical walls (applies to no-slip, moving, slip, etc.)
2412 if (face_config->mathematical_type != WALL) {
2413 continue;
2414 }
2415
2416 // Check if this MPI rank owns part of this face
2417 PetscBool rank_owns_this_face;
2418 ierr = CanRankServiceFace(info, user->IM, user->JM, user->KM,
2419 current_face_id, &rank_owns_this_face); CHKERRQ(ierr);
2420
2421 if (!rank_owns_this_face) {
2422 continue;
2423 }
2424
2425 LOG_ALLOW(LOCAL, LOG_TRACE, "Processing Face %d (%s)\n",
2426 current_face_id, BCFaceToString(current_face_id));
2427
2428 // =====================================================================
2429 // Process each face with appropriate indexing
2430 // =====================================================================
2431 switch(current_face_id) {
2432
2433 // =================================================================
2434 // NEGATIVE X FACE (i = 0, normal points in +X direction)
2435 // =================================================================
2436 case BC_FACE_NEG_X: {
2437 if (grid_start_i == 0) {
2438 const PetscInt ghost_cell_index = grid_start_i;
2439 const PetscInt first_interior_cell = grid_start_i + 1;
2440 const PetscInt second_interior_cell = grid_start_i + 2;
2441
2442 for (PetscInt k = loop_start_k; k < loop_end_k; k++) {
2443 for (PetscInt j = loop_start_j; j < loop_end_j; j++) {
2444
2445 // Skip if this is a solid cell (embedded boundary)
2446 if (node_vertex_flag[k][j][first_interior_cell] < 0.1) {
2447
2448 // Calculate face area from contravariant metric tensor
2449 PetscReal face_area = sqrt(
2450 csi[k][j][ghost_cell_index].x * csi[k][j][ghost_cell_index].x +
2451 csi[k][j][ghost_cell_index].y * csi[k][j][ghost_cell_index].y +
2452 csi[k][j][ghost_cell_index].z * csi[k][j][ghost_cell_index].z
2453 );
2454
2455 // Compute wall-normal distances using cell Jacobians
2456 // sb = distance from wall to first interior cell center
2457 // sc = distance from wall to second interior cell center
2458 PetscReal distance_to_first_cell = 0.5 / cell_jacobian[k][j][first_interior_cell] / face_area;
2459 PetscReal distance_to_second_cell = 2.0 * distance_to_first_cell +
2460 0.5 / cell_jacobian[k][j][second_interior_cell] / face_area;
2461
2462 // Compute unit normal vector pointing INTO the domain
2463 PetscReal wall_normal[3];
2464 wall_normal[0] = csi[k][j][ghost_cell_index].x / face_area;
2465 wall_normal[1] = csi[k][j][ghost_cell_index].y / face_area;
2466 wall_normal[2] = csi[k][j][ghost_cell_index].z / face_area;
2467
2468 // Define velocities for wall function calculation
2469 Cmpnts wall_velocity; // Ua = velocity at wall (zero for stationary wall)
2470 Cmpnts reference_velocity; // Uc = velocity at second interior cell
2471
2472 wall_velocity.x = wall_velocity.y = wall_velocity.z = 0.0;
2473 reference_velocity = velocity_cartesian[k][j][second_interior_cell];
2474
2475 // Step 1: Linear interpolation (provides initial guess)
2476 noslip(user, distance_to_second_cell, distance_to_first_cell,
2477 wall_velocity, reference_velocity,
2478 &velocity_cartesian[k][j][first_interior_cell],
2479 wall_normal[0], wall_normal[1], wall_normal[2]);
2480
2481 // Step 2: Apply log-law correction (improves near-wall velocity)
2482 wall_function_loglaw(user, wall_roughness_height,
2483 distance_to_second_cell, distance_to_first_cell,
2484 wall_velocity, reference_velocity,
2485 &velocity_cartesian[k][j][first_interior_cell],
2486 &friction_velocity[k][j][first_interior_cell],
2487 wall_normal[0], wall_normal[1], wall_normal[2]);
2488
2489 // Ensure ghost cell BC remains zero (required for proper extrapolation)
2490 velocity_boundary[k][j][ghost_cell_index].x = 0.0;
2491 velocity_boundary[k][j][ghost_cell_index].y = 0.0;
2492 velocity_boundary[k][j][ghost_cell_index].z = 0.0;
2493 velocity_contravariant[k][j][ghost_cell_index].x = 0.0;
2494 }
2495 }
2496 }
2497 }
2498 } break;
2499
2500 // =================================================================
2501 // POSITIVE X FACE (i = mx-1, normal points in -X direction)
2502 // =================================================================
2503 case BC_FACE_POS_X: {
2504 if (grid_end_i == grid_size_i) {
2505 const PetscInt ghost_cell_index = grid_end_i - 1;
2506 const PetscInt first_interior_cell = grid_end_i - 2;
2507 const PetscInt second_interior_cell = grid_end_i - 3;
2508
2509 for (PetscInt k = loop_start_k; k < loop_end_k; k++) {
2510 for (PetscInt j = loop_start_j; j < loop_end_j; j++) {
2511
2512 if (node_vertex_flag[k][j][first_interior_cell] < 0.1) {
2513
2514 PetscReal face_area = sqrt(
2515 csi[k][j][first_interior_cell].x * csi[k][j][first_interior_cell].x +
2516 csi[k][j][first_interior_cell].y * csi[k][j][first_interior_cell].y +
2517 csi[k][j][first_interior_cell].z * csi[k][j][first_interior_cell].z
2518 );
2519
2520 PetscReal distance_to_first_cell = 0.5 / cell_jacobian[k][j][first_interior_cell] / face_area;
2521 PetscReal distance_to_second_cell = 2.0 * distance_to_first_cell +
2522 0.5 / cell_jacobian[k][j][second_interior_cell] / face_area;
2523
2524 // Note: Normal flipped for +X face to point INTO domain
2525 PetscReal wall_normal[3];
2526 wall_normal[0] = -csi[k][j][first_interior_cell].x / face_area;
2527 wall_normal[1] = -csi[k][j][first_interior_cell].y / face_area;
2528 wall_normal[2] = -csi[k][j][first_interior_cell].z / face_area;
2529
2530 Cmpnts wall_velocity, reference_velocity;
2531 wall_velocity.x = wall_velocity.y = wall_velocity.z = 0.0;
2532 reference_velocity = velocity_cartesian[k][j][second_interior_cell];
2533
2534 noslip(user, distance_to_second_cell, distance_to_first_cell,
2535 wall_velocity, reference_velocity,
2536 &velocity_cartesian[k][j][first_interior_cell],
2537 wall_normal[0], wall_normal[1], wall_normal[2]);
2538
2539 wall_function_loglaw(user, wall_roughness_height,
2540 distance_to_second_cell, distance_to_first_cell,
2541 wall_velocity, reference_velocity,
2542 &velocity_cartesian[k][j][first_interior_cell],
2543 &friction_velocity[k][j][first_interior_cell],
2544 wall_normal[0], wall_normal[1], wall_normal[2]);
2545
2546 velocity_boundary[k][j][ghost_cell_index].x = 0.0;
2547 velocity_boundary[k][j][ghost_cell_index].y = 0.0;
2548 velocity_boundary[k][j][ghost_cell_index].z = 0.0;
2549 velocity_contravariant[k][j][first_interior_cell].x = 0.0;
2550 }
2551 }
2552 }
2553 }
2554 } break;
2555
2556 // =================================================================
2557 // NEGATIVE Y FACE (j = 0, normal points in +Y direction)
2558 // =================================================================
2559 case BC_FACE_NEG_Y: {
2560 if (grid_start_j == 0) {
2561 const PetscInt ghost_cell_index = grid_start_j;
2562 const PetscInt first_interior_cell = grid_start_j + 1;
2563 const PetscInt second_interior_cell = grid_start_j + 2;
2564
2565 for (PetscInt k = loop_start_k; k < loop_end_k; k++) {
2566 for (PetscInt i = loop_start_i; i < loop_end_i; i++) {
2567
2568 if (node_vertex_flag[k][first_interior_cell][i] < 0.1) {
2569
2570 PetscReal face_area = sqrt(
2571 eta[k][ghost_cell_index][i].x * eta[k][ghost_cell_index][i].x +
2572 eta[k][ghost_cell_index][i].y * eta[k][ghost_cell_index][i].y +
2573 eta[k][ghost_cell_index][i].z * eta[k][ghost_cell_index][i].z
2574 );
2575
2576 PetscReal distance_to_first_cell = 0.5 / cell_jacobian[k][first_interior_cell][i] / face_area;
2577 PetscReal distance_to_second_cell = 2.0 * distance_to_first_cell +
2578 0.5 / cell_jacobian[k][second_interior_cell][i] / face_area;
2579
2580 PetscReal wall_normal[3];
2581 wall_normal[0] = eta[k][ghost_cell_index][i].x / face_area;
2582 wall_normal[1] = eta[k][ghost_cell_index][i].y / face_area;
2583 wall_normal[2] = eta[k][ghost_cell_index][i].z / face_area;
2584
2585 Cmpnts wall_velocity, reference_velocity;
2586 wall_velocity.x = wall_velocity.y = wall_velocity.z = 0.0;
2587 reference_velocity = velocity_cartesian[k][second_interior_cell][i];
2588
2589 noslip(user, distance_to_second_cell, distance_to_first_cell,
2590 wall_velocity, reference_velocity,
2591 &velocity_cartesian[k][first_interior_cell][i],
2592 wall_normal[0], wall_normal[1], wall_normal[2]);
2593
2594 wall_function_loglaw(user, wall_roughness_height,
2595 distance_to_second_cell, distance_to_first_cell,
2596 wall_velocity, reference_velocity,
2597 &velocity_cartesian[k][first_interior_cell][i],
2598 &friction_velocity[k][first_interior_cell][i],
2599 wall_normal[0], wall_normal[1], wall_normal[2]);
2600
2601 velocity_boundary[k][ghost_cell_index][i].x = 0.0;
2602 velocity_boundary[k][ghost_cell_index][i].y = 0.0;
2603 velocity_boundary[k][ghost_cell_index][i].z = 0.0;
2604 velocity_contravariant[k][ghost_cell_index][i].y = 0.0;
2605 }
2606 }
2607 }
2608 }
2609 } break;
2610
2611 // =================================================================
2612 // POSITIVE Y FACE (j = my-1, normal points in -Y direction)
2613 // =================================================================
2614 case BC_FACE_POS_Y: {
2615 if (grid_end_j == grid_size_j) {
2616 const PetscInt ghost_cell_index = grid_end_j - 1;
2617 const PetscInt first_interior_cell = grid_end_j - 2;
2618 const PetscInt second_interior_cell = grid_end_j - 3;
2619
2620 for (PetscInt k = loop_start_k; k < loop_end_k; k++) {
2621 for (PetscInt i = loop_start_i; i < loop_end_i; i++) {
2622
2623 if (node_vertex_flag[k][first_interior_cell][i] < 0.1) {
2624
2625 PetscReal face_area = sqrt(
2626 eta[k][first_interior_cell][i].x * eta[k][first_interior_cell][i].x +
2627 eta[k][first_interior_cell][i].y * eta[k][first_interior_cell][i].y +
2628 eta[k][first_interior_cell][i].z * eta[k][first_interior_cell][i].z
2629 );
2630
2631 PetscReal distance_to_first_cell = 0.5 / cell_jacobian[k][first_interior_cell][i] / face_area;
2632 PetscReal distance_to_second_cell = 2.0 * distance_to_first_cell +
2633 0.5 / cell_jacobian[k][second_interior_cell][i] / face_area;
2634
2635 PetscReal wall_normal[3];
2636 wall_normal[0] = -eta[k][first_interior_cell][i].x / face_area;
2637 wall_normal[1] = -eta[k][first_interior_cell][i].y / face_area;
2638 wall_normal[2] = -eta[k][first_interior_cell][i].z / face_area;
2639
2640 Cmpnts wall_velocity, reference_velocity;
2641 wall_velocity.x = wall_velocity.y = wall_velocity.z = 0.0;
2642 reference_velocity = velocity_cartesian[k][second_interior_cell][i];
2643
2644 noslip(user, distance_to_second_cell, distance_to_first_cell,
2645 wall_velocity, reference_velocity,
2646 &velocity_cartesian[k][first_interior_cell][i],
2647 wall_normal[0], wall_normal[1], wall_normal[2]);
2648
2649 wall_function_loglaw(user, wall_roughness_height,
2650 distance_to_second_cell, distance_to_first_cell,
2651 wall_velocity, reference_velocity,
2652 &velocity_cartesian[k][first_interior_cell][i],
2653 &friction_velocity[k][first_interior_cell][i],
2654 wall_normal[0], wall_normal[1], wall_normal[2]);
2655
2656 velocity_boundary[k][ghost_cell_index][i].x = 0.0;
2657 velocity_boundary[k][ghost_cell_index][i].y = 0.0;
2658 velocity_boundary[k][ghost_cell_index][i].z = 0.0;
2659 velocity_contravariant[k][first_interior_cell][i].y = 0.0;
2660 }
2661 }
2662 }
2663 }
2664 } break;
2665
2666 // =================================================================
2667 // NEGATIVE Z FACE (k = 0, normal points in +Z direction)
2668 // =================================================================
2669 case BC_FACE_NEG_Z: {
2670 if (grid_start_k == 0) {
2671 const PetscInt ghost_cell_index = grid_start_k;
2672 const PetscInt first_interior_cell = grid_start_k + 1;
2673 const PetscInt second_interior_cell = grid_start_k + 2;
2674
2675 for (PetscInt j = loop_start_j; j < loop_end_j; j++) {
2676 for (PetscInt i = loop_start_i; i < loop_end_i; i++) {
2677
2678 if (node_vertex_flag[first_interior_cell][j][i] < 0.1) {
2679
2680 PetscReal face_area = sqrt(
2681 zet[ghost_cell_index][j][i].x * zet[ghost_cell_index][j][i].x +
2682 zet[ghost_cell_index][j][i].y * zet[ghost_cell_index][j][i].y +
2683 zet[ghost_cell_index][j][i].z * zet[ghost_cell_index][j][i].z
2684 );
2685
2686 PetscReal distance_to_first_cell = 0.5 / cell_jacobian[first_interior_cell][j][i] / face_area;
2687 PetscReal distance_to_second_cell = 2.0 * distance_to_first_cell +
2688 0.5 / cell_jacobian[second_interior_cell][j][i] / face_area;
2689
2690 PetscReal wall_normal[3];
2691 wall_normal[0] = zet[ghost_cell_index][j][i].x / face_area;
2692 wall_normal[1] = zet[ghost_cell_index][j][i].y / face_area;
2693 wall_normal[2] = zet[ghost_cell_index][j][i].z / face_area;
2694
2695 Cmpnts wall_velocity, reference_velocity;
2696 wall_velocity.x = wall_velocity.y = wall_velocity.z = 0.0;
2697 reference_velocity = velocity_cartesian[second_interior_cell][j][i];
2698
2699 noslip(user, distance_to_second_cell, distance_to_first_cell,
2700 wall_velocity, reference_velocity,
2701 &velocity_cartesian[first_interior_cell][j][i],
2702 wall_normal[0], wall_normal[1], wall_normal[2]);
2703
2704 wall_function_loglaw(user, wall_roughness_height,
2705 distance_to_second_cell, distance_to_first_cell,
2706 wall_velocity, reference_velocity,
2707 &velocity_cartesian[first_interior_cell][j][i],
2708 &friction_velocity[first_interior_cell][j][i],
2709 wall_normal[0], wall_normal[1], wall_normal[2]);
2710
2711 velocity_boundary[ghost_cell_index][j][i].x = 0.0;
2712 velocity_boundary[ghost_cell_index][j][i].y = 0.0;
2713 velocity_boundary[ghost_cell_index][j][i].z = 0.0;
2714 velocity_contravariant[ghost_cell_index][j][i].z = 0.0;
2715 }
2716 }
2717 }
2718 }
2719 } break;
2720
2721 // =================================================================
2722 // POSITIVE Z FACE (k = mz-1, normal points in -Z direction)
2723 // =================================================================
2724 case BC_FACE_POS_Z: {
2725 if (grid_end_k == grid_size_k) {
2726 const PetscInt ghost_cell_index = grid_end_k - 1;
2727 const PetscInt first_interior_cell = grid_end_k - 2;
2728 const PetscInt second_interior_cell = grid_end_k - 3;
2729
2730 for (PetscInt j = loop_start_j; j < loop_end_j; j++) {
2731 for (PetscInt i = loop_start_i; i < loop_end_i; i++) {
2732
2733 if (node_vertex_flag[first_interior_cell][j][i] < 0.1) {
2734
2735 PetscReal face_area = sqrt(
2736 zet[first_interior_cell][j][i].x * zet[first_interior_cell][j][i].x +
2737 zet[first_interior_cell][j][i].y * zet[first_interior_cell][j][i].y +
2738 zet[first_interior_cell][j][i].z * zet[first_interior_cell][j][i].z
2739 );
2740
2741 PetscReal distance_to_first_cell = 0.5 / cell_jacobian[first_interior_cell][j][i] / face_area;
2742 PetscReal distance_to_second_cell = 2.0 * distance_to_first_cell +
2743 0.5 / cell_jacobian[second_interior_cell][j][i] / face_area;
2744
2745 PetscReal wall_normal[3];
2746 wall_normal[0] = -zet[first_interior_cell][j][i].x / face_area;
2747 wall_normal[1] = -zet[first_interior_cell][j][i].y / face_area;
2748 wall_normal[2] = -zet[first_interior_cell][j][i].z / face_area;
2749
2750 Cmpnts wall_velocity, reference_velocity;
2751 wall_velocity.x = wall_velocity.y = wall_velocity.z = 0.0;
2752 reference_velocity = velocity_cartesian[second_interior_cell][j][i];
2753
2754 noslip(user, distance_to_second_cell, distance_to_first_cell,
2755 wall_velocity, reference_velocity,
2756 &velocity_cartesian[first_interior_cell][j][i],
2757 wall_normal[0], wall_normal[1], wall_normal[2]);
2758
2759 wall_function_loglaw(user, wall_roughness_height,
2760 distance_to_second_cell, distance_to_first_cell,
2761 wall_velocity, reference_velocity,
2762 &velocity_cartesian[first_interior_cell][j][i],
2763 &friction_velocity[first_interior_cell][j][i],
2764 wall_normal[0], wall_normal[1], wall_normal[2]);
2765
2766 velocity_boundary[ghost_cell_index][j][i].x = 0.0;
2767 velocity_boundary[ghost_cell_index][j][i].y = 0.0;
2768 velocity_boundary[ghost_cell_index][j][i].z = 0.0;
2769 velocity_contravariant[first_interior_cell][j][i].z = 0.0;
2770 }
2771 }
2772 }
2773 }
2774 } break;
2775 }
2776 }
2777
2778 // =========================================================================
2779 // STEP 4: Restore all arrays and release memory
2780 // =========================================================================
2781 ierr = DMDAVecRestoreArray(user->fda, user->Ucat, &velocity_cartesian); CHKERRQ(ierr);
2782 ierr = DMDAVecRestoreArray(user->fda, user->Ucont, &velocity_contravariant); CHKERRQ(ierr);
2783 ierr = DMDAVecRestoreArray(user->fda, user->Bcs.Ubcs, &velocity_boundary); CHKERRQ(ierr);
2784 ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
2785 ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
2786 ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
2787 ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, (const PetscReal***)&node_vertex_flag); CHKERRQ(ierr);
2788 ierr = DMDAVecRestoreArrayRead(user->da, user->lAj, (const PetscReal***)&cell_jacobian); CHKERRQ(ierr);
2789 ierr = DMDAVecRestoreArray(user->da, user->lFriction_Velocity, &friction_velocity); CHKERRQ(ierr);
2790
2791 LOG_ALLOW(LOCAL, LOG_DEBUG, "Complete.\n");
2792
2793 PetscFunctionReturn(0);
2794}
2795
2796#undef __FUNCT__
2797#define __FUNCT__ "RefreshBoundaryGhostCells"
2798/**
2799 * @brief Internal helper implementation: `RefreshBoundaryGhostCells()`.
2800 * @details Local to this translation unit.
2801 */
2802PetscErrorCode RefreshBoundaryGhostCells(UserCtx *user)
2803{
2804 PetscErrorCode ierr;
2805 PetscFunctionBeginUser;
2806 PROFILE_FUNCTION_BEGIN;
2807
2808 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Starting post-projection refresh of boundary ghost cells.\n");
2809
2810 // -------------------------------------------------------------------------
2811 // STEP 1: Refresh Flow-Dependent Boundary Value Targets (`ubcs`)
2812 // -------------------------------------------------------------------------
2813 // This is the "physics" part of the refresh. It calls the lightweight engine
2814 // to update `ubcs` for any BCs (like Symmetry) that depend on the now-corrected
2815 // interior velocity field. This step does NOT touch `ucont`.
2816 ierr = BoundarySystem_RefreshUbcs(user); CHKERRQ(ierr);
2817 LOG_ALLOW(GLOBAL, LOG_VERBOSE, " `ubcs` targets refreshed.\n");
2818
2819 ierr = UpdateLocalGhosts(user,"Ucat");
2820
2821 // STEP 1.5: Apply Wall function if applicable.
2822 if(user->simCtx->wallfunction){
2823
2824 ierr = ApplyWallFunction(user); CHKERRQ(ierr);
2825
2826 }
2827 // -------------------------------------------------------------------------
2828 // STEP 2: Apply Geometric Ghost Cell Updates
2829 // -------------------------------------------------------------------------
2830 // With `ubcs` now fully up-to-date, we execute the purely geometric
2831 // operations to fill the entire ghost cell layer. The order is important.
2832
2833 // (a) Update the ghost cells on the faces of non-periodic boundaries.
2834 ierr = UpdateDummyCells(user); CHKERRQ(ierr);
2835 LOG_ALLOW(GLOBAL, LOG_VERBOSE, " Face ghost cells (dummy cells) updated.\n");
2836
2837 // (b) Handle periodic boundaries first. This is a direct data copy.
2838 // Ghost Update is done inside this function.s
2839 ierr = ApplyPeriodicBCs(user); CHKERRQ(ierr);
2840 LOG_ALLOW(GLOBAL, LOG_VERBOSE, " Periodic boundaries synchronized.\n");
2841
2842 // (c) Update the ghost cells at the edges and corners by averaging.
2843 ierr = UpdateCornerNodes(user); CHKERRQ(ierr);
2844 LOG_ALLOW(GLOBAL, LOG_VERBOSE, " Edge and corner ghost cells updated.\n");
2845
2846 // (d) Synchronize Ucat after setting corner nodes.
2847 ierr = UpdateLocalGhosts(user, "Ucat"); CHKERRQ(ierr);
2848
2849 // (e) Update the corners for periodic conditions sequentially
2850 // ensuring no race conditions are raised.
2851 const char* ucat_only[] = {"Ucat"};
2852 ierr = UpdatePeriodicCornerNodes(user,1,ucat_only);CHKERRQ(ierr);
2853
2854 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Boundary ghost cell refresh complete.\n");
2855 PROFILE_FUNCTION_END;
2856 PetscFunctionReturn(0);
2857}
2858
2859#undef __FUNCT__
2860#define __FUNCT__ "ApplyBoundaryConditions"
2861/**
2862 * @brief Implementation of \ref ApplyBoundaryConditions().
2863 * @details Full API contract (arguments, ownership, side effects) is documented with
2864 * the header declaration in `include/Boundaries.h`.
2865 * @see ApplyBoundaryConditions()
2866 */
2867PetscErrorCode ApplyBoundaryConditions(UserCtx *user)
2868{
2869 PetscErrorCode ierr;
2870 PetscFunctionBeginUser;
2871 PROFILE_FUNCTION_BEGIN;
2872
2873 LOG_ALLOW(GLOBAL,LOG_TRACE,"Boundary Condition Application begins.\n");
2874
2875 // STEP 1: Main iteration loop for applying and converging non-periodic BCs.
2876 // The number of iterations (e.g., 3) allows information to propagate
2877 // between coupled boundaries, like an inlet and a conserving outlet.
2878 for (PetscInt iter = 0; iter < 3; iter++) {
2879 // (a) Execute the boundary system. This phase calculates fluxes across
2880 // the domain and then applies the physical logic for each non-periodic
2881 // handler, setting the `ubcs` (boundary value) array.
2882 ierr = BoundarySystem_ExecuteStep(user); CHKERRQ(ierr);
2883
2884 LOG_ALLOW(GLOBAL,LOG_VERBOSE,"Boundary Condition Setup Executed.\n");
2885
2886 // (b) Synchronize the updated ghost cells across all processors to ensure
2887 // all ucont values are current before updating the dummy cells.
2888 ierr = UpdateLocalGhosts(user, "Ucont"); CHKERRQ(ierr);
2889
2890 // (c) Convert updated Contravariant velocities to Cartesian velocities.
2891 ierr = Contra2Cart(user); CHKERRQ(ierr);
2892
2893 // (d) Synchronize the updated Cartesian velocities across all processors
2894 // to ensure all ucat values are current before updating the dummy cells.
2895 ierr = UpdateLocalGhosts(user, "Ucat"); CHKERRQ(ierr);
2896
2897 // (e) If Wall functions are enabled, apply them now to adjust near-wall velocities.
2898 if(user->simCtx->wallfunction){
2899 // Apply wall function adjustments to the boundary velocities.
2900 ierr = ApplyWallFunction(user); CHKERRQ(ierr);
2901
2902 // Synchronize the updated Cartesian velocities after wall function adjustments.
2903 ierr = UpdateLocalGhosts(user, "Ucat"); CHKERRQ(ierr);
2904
2905 LOG_ALLOW(GLOBAL,LOG_VERBOSE,"Wall Function Applied at Walls.\n");
2906 }
2907
2908 // (f) Update the first layer of ghost cells for non-periodic faces using
2909 // the newly computed `ubcs` values.
2910 ierr = UpdateDummyCells(user); CHKERRQ(ierr);
2911
2912 LOG_ALLOW(GLOBAL,LOG_VERBOSE,"Dummy Cells/Ghost Cells Updated.\n");
2913
2914 // (g) Handle all periodic boundaries. This is a parallel direct copy
2915 // that sets the absolute constraints for the rest of the solve.
2916 // There is a Ghost update happening inside this function.
2917 ierr = ApplyPeriodicBCs(user); CHKERRQ(ierr);
2918
2919 // (h) Update the corner and edge ghost nodes. This routine calculates
2920 // values for corners/edges by averaging their neighbors, which have been
2921 // finalized in the steps above (both periodic and non-periodic).
2922 ierr = UpdateCornerNodes(user); CHKERRQ(ierr);
2923
2924 // (i) Synchronize the updated edge and corner cells across all processors to ensure
2925 // consistency before the next iteration or finalization.
2926 ierr = UpdateLocalGhosts(user, "P"); CHKERRQ(ierr);
2927 ierr = UpdateLocalGhosts(user, "Ucat"); CHKERRQ(ierr);
2928 ierr = UpdateLocalGhosts(user, "Ucont"); CHKERRQ(ierr);
2929
2930 // (j) Ensure All the corners are synchronized with a well defined protocol in case of Periodic boundary conditions
2931 // To avoid race conditions.
2932 const char* all_fields[] = {"Ucat", "P", "Nvert"};
2933 ierr = UpdatePeriodicCornerNodes(user,3,all_fields); CHKERRQ(ierr);
2934
2935 }
2936
2937 // STEP 3: Final ghost node synchronization. This ensures all changes made
2938 // to the global vectors are reflected in the local ghost regions of all
2939 // processors, making the state fully consistent before the next solver stage.
2940 ierr = UpdateLocalGhosts(user, "P"); CHKERRQ(ierr);
2941 ierr = UpdateLocalGhosts(user, "Ucat"); CHKERRQ(ierr);
2942 ierr = UpdateLocalGhosts(user, "Ucont"); CHKERRQ(ierr);
2943
2944 PROFILE_FUNCTION_END;
2945 PetscFunctionReturn(0);
2946}
Create_InletConstantVelocity
PetscErrorCode Create_InletConstantVelocity(BoundaryCondition *bc)
Configures a BoundaryCondition object to behave as a constant velocity inlet.
Definition BC_Handlers.c:323
Create_InletProfileFromFile
PetscErrorCode Create_InletProfileFromFile(BoundaryCondition *bc)
Configures a BoundaryCondition object for a file-prescribed inlet profile.
Definition BC_Handlers.c:1303
Create_PeriodicGeometric
PetscErrorCode Create_PeriodicGeometric(BoundaryCondition *bc)
Configures a BoundaryCondition object for geometric periodic coupling.
Definition BC_Handlers.c:2172
Validate_DrivenFlowConfiguration
PetscErrorCode Validate_DrivenFlowConfiguration(UserCtx *user)
(Private) Validates all consistency rules for a driven flow (channel/pipe) setup.
Definition BC_Handlers.c:15
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:2225
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:1655
ApplyPeriodicBCs
PetscErrorCode ApplyPeriodicBCs(UserCtx *user)
Internal helper implementation: ApplyPeriodicBCs().
Definition Boundaries.c:1927
BoundarySystem_Initialize
PetscErrorCode BoundarySystem_Initialize(UserCtx *user, const char *bcs_filename)
Implementation of BoundarySystem_Initialize().
Definition Boundaries.c:862
ApplyUcontPeriodicBCs
PetscErrorCode ApplyUcontPeriodicBCs(UserCtx *user)
Implementation of ApplyUcontPeriodicBCs().
Definition Boundaries.c:1988
TransferPeriodicFaceField
PetscErrorCode TransferPeriodicFaceField(UserCtx *user, const char *field_name)
Internal helper implementation: TransferPeriodicFaceField().
Definition Boundaries.c:1764
RefreshBoundaryGhostCells
PetscErrorCode RefreshBoundaryGhostCells(UserCtx *user)
Internal helper implementation: RefreshBoundaryGhostCells().
Definition Boundaries.c:2802
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().
Definition Boundaries.c:399
ApplyWallFunction
PetscErrorCode ApplyWallFunction(UserCtx *user)
Internal helper implementation: ApplyWallFunction().
Definition Boundaries.c:2343
PropagateBoundaryConfigToCoarserLevels
PetscErrorCode PropagateBoundaryConfigToCoarserLevels(SimCtx *simCtx)
Internal helper implementation: PropagateBoundaryConfigToCoarserLevels().
Definition Boundaries.c:987
ApplyMetricsPeriodicBCs
PetscErrorCode ApplyMetricsPeriodicBCs(UserCtx *user)
Internal helper implementation: ApplyMetricsPeriodicBCs().
Definition Boundaries.c:1900
EnforceRHSBoundaryConditions
PetscErrorCode EnforceRHSBoundaryConditions(UserCtx *user)
Internal helper implementation: EnforceRHSBoundaryConditions().
Definition Boundaries.c:591
EnforceUcontPeriodicity
PetscErrorCode EnforceUcontPeriodicity(UserCtx *user)
Internal helper implementation: EnforceUcontPeriodicity().
Definition Boundaries.c:2012
TransferPeriodicFieldByDirection
PetscErrorCode TransferPeriodicFieldByDirection(UserCtx *user, const char *field_name, char direction)
Internal helper implementation: TransferPeriodicFieldByDirection().
Definition Boundaries.c:1577
UpdateDummyCells
PetscErrorCode UpdateDummyCells(UserCtx *user)
Internal helper implementation: UpdateDummyCells().
Definition Boundaries.c:2064
BoundarySystem_RefreshUbcs
PetscErrorCode BoundarySystem_RefreshUbcs(UserCtx *user)
Internal helper implementation: BoundarySystem_RefreshUbcs().
Definition Boundaries.c:1480
BoundarySystem_Validate
PetscErrorCode BoundarySystem_Validate(UserCtx *user)
Internal helper implementation: BoundarySystem_Validate().
Definition Boundaries.c:830
UpdatePeriodicCornerNodes
PetscErrorCode UpdatePeriodicCornerNodes(UserCtx *user, PetscInt num_fields, const char *field_names[])
Internal helper implementation: UpdatePeriodicCornerNodes().
Definition Boundaries.c:2300
BoundaryCondition_Create
PetscErrorCode BoundaryCondition_Create(BCHandlerType handler_type, BoundaryCondition **new_bc_ptr)
Internal helper implementation: BoundaryCondition_Create().
Definition Boundaries.c:744
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().
Definition Boundaries.c:11
ApplyBoundaryConditions
PetscErrorCode ApplyBoundaryConditions(UserCtx *user)
Implementation of ApplyBoundaryConditions().
Definition Boundaries.c:2867
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
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().
Definition Boundaries.c:212
BoundarySystem_ExecuteStep
PetscErrorCode BoundarySystem_ExecuteStep(UserCtx *user)
Implementation of BoundarySystem_ExecuteStep().
Definition Boundaries.c:1058
BoundarySystem_Destroy
PetscErrorCode BoundarySystem_Destroy(UserCtx *user)
Implementation of BoundarySystem_Destroy().
Definition Boundaries.c:1530
UpdateCornerNodes
PetscErrorCode UpdateCornerNodes(UserCtx *user)
Internal helper implementation: UpdateCornerNodes().
Definition Boundaries.c:2156
TransferPeriodicField
PetscErrorCode TransferPeriodicField(UserCtx *user, const char *field_name)
Internal helper implementation: TransferPeriodicField().
Definition Boundaries.c:1661
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:450
FreeBC_ParamList
void FreeBC_ParamList(BC_Param *head)
Frees an entire linked list of boundary-condition parameters.
Definition io.c:302
BCHandlerTypeToString
const char * BCHandlerTypeToString(BCHandlerType handler_type)
Converts a BCHandlerType enum to its string representation.
Definition logging.c:771
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
LOCAL
#define LOCAL
Logging scope definitions for controlling message output.
Definition logging.h:44
GLOBAL
#define GLOBAL
Scope for global logging across all processes.
Definition logging.h:45
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:820
BCTypeToString
const char * BCTypeToString(BCType type)
Helper function to convert BCType enum to a string representation.
Definition logging.c:751
LOG_ERROR
@ LOG_ERROR
Critical errors that may halt the program.
Definition logging.h:28
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
LOG_VERBOSE
@ LOG_VERBOSE
Extremely detailed logs, typically for development use only.
Definition logging.h:33
PROFILE_FUNCTION_BEGIN
#define PROFILE_FUNCTION_BEGIN
Marks the beginning of a profiled code block (typically a function).
Definition logging.h:811
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:2032
Contra2Cart
PetscErrorCode Contra2Cart(UserCtx *user)
Reconstructs Cartesian velocity (Ucat) at cell centers from contravariant velocity (Ucont) defined on...
Definition setup.c:2396
UpdateLocalGhosts
PetscErrorCode UpdateLocalGhosts(UserCtx *user, const char *fieldName)
Updates the local vector (including ghost points) from its corresponding global vector.
Definition setup.c:1510
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
UserCtx::lFriction_Velocity
Vec lFriction_Velocity
Definition variables.h:852
SimCtx::FarFluxInSum
PetscReal FarFluxInSum
Definition variables.h:735
PERIODIC
@ PERIODIC
Definition variables.h:260
WALL
@ WALL
Definition variables.h:254
MGCtx::user
UserCtx * user
Definition variables.h:536
SimCtx::FarFluxOutSum
PetscReal FarFluxOutSum
Definition variables.h:735
UserCtx::inletFaceDefined
PetscBool inletFaceDefined
Definition variables.h:849
UserCtx::Rhs
Vec Rhs
Definition variables.h:864
SimCtx::rank
PetscMPIInt rank
Definition variables.h:654
UserCtx::boundary_faces
BoundaryFaceConfig boundary_faces[6]
Definition variables.h:848
SimCtx::block_number
PetscInt block_number
Definition variables.h:726
UserCtx::lIEta
Vec lIEta
Definition variables.h:881
UserCtx::identifiedInletBCFace
BCFace identifiedInletBCFace
Definition variables.h:850
UserCtx::lIZet
Vec lIZet
Definition variables.h:881
UserCtx::lNvert
Vec lNvert
Definition variables.h:856
UserCtx::simCtx
SimCtx * simCtx
Back-pointer to the master simulation context.
Definition variables.h:833
SimCtx::FluxOutSum
PetscReal FluxOutSum
Definition variables.h:735
BC_Param_s::next
struct BC_Param_s * next
Definition variables.h:307
BC_Param_s::key
char * key
Definition variables.h:305
UserCtx::KM
PetscInt KM
Definition variables.h:839
UserCtx::lZet
Vec lZet
Definition variables.h:879
SimCtx::usermg
UserMG usermg
Definition variables.h:778
BCHandlerType
BCHandlerType
Defines the specific computational "strategy" for a boundary handler.
Definition variables.h:271
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_INLET_PROFILE_FROM_FILE
@ BC_HANDLER_INLET_PROFILE_FROM_FILE
Definition variables.h:278
BC_HANDLER_WALL_NOSLIP
@ BC_HANDLER_WALL_NOSLIP
Definition variables.h:273
BC_HANDLER_OUTLET_CONSERVATION
@ BC_HANDLER_OUTLET_CONSERVATION
Definition variables.h:282
BoundaryFaceConfig::handler_type
BCHandlerType handler_type
Definition variables.h:337
UserCtx::lIAj
Vec lIAj
Definition variables.h:881
UserCtx::_this
PetscInt _this
Definition variables.h:843
UserCtx::lKEta
Vec lKEta
Definition variables.h:883
SimCtx::np
PetscInt np
Definition variables.h:753
UserCtx::lJCsi
Vec lJCsi
Definition variables.h:882
BC_Param_s::value
char * value
Definition variables.h:306
UserCtx::Ucont
Vec Ucont
Definition variables.h:856
BCS::Ubcs
Vec Ubcs
Physical Cartesian velocity at boundary faces. Full 3D array but only boundary-face entries are meani...
Definition variables.h:121
Cmpnts::x
PetscScalar x
Definition variables.h:101
UserCtx::Bcs
BCS Bcs
Definition variables.h:851
BCContext::user
UserCtx * user
Definition variables.h:312
UserCtx::lKZet
Vec lKZet
Definition variables.h:883
UserCtx::lNu_t
Vec lNu_t
Definition variables.h:886
SimCtx::FluxInSum
PetscReal FluxInSum
Definition variables.h:735
UserCtx::Nu_t
Vec Nu_t
Definition variables.h:886
UserCtx::lJEta
Vec lJEta
Definition variables.h:882
UserCtx::lCsi
Vec lCsi
Definition variables.h:879
BoundaryFaceConfig::params
BC_Param * params
Definition variables.h:338
Cmpnts::z
PetscScalar z
Definition variables.h:101
UserCtx::lKCsi
Vec lKCsi
Definition variables.h:883
UserCtx::Ucat
Vec Ucat
Definition variables.h:856
UserCtx::JM
PetscInt JM
Definition variables.h:839
SimCtx::wallfunction
PetscInt wallfunction
Definition variables.h:747
UserMG::mglevels
PetscInt mglevels
Definition variables.h:543
UserCtx::lJZet
Vec lJZet
Definition variables.h:882
UserCtx::lUcont
Vec lUcont
Definition variables.h:856
UserCtx::lAj
Vec lAj
Definition variables.h:879
UserCtx::lICsi
Vec lICsi
Definition variables.h:881
UserCtx::info
DMDALocalInfo info
Definition variables.h:837
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
UserCtx::lUcat
Vec lUcat
Definition variables.h:856
Cmpnts::y
PetscScalar y
Definition variables.h:101
UserCtx::IM
PetscInt IM
Definition variables.h:839
UserCtx::lEta
Vec lEta
Definition variables.h:879
UserCtx::Nvert
Vec Nvert
Definition variables.h:856
UserMG::mgctx
MGCtx * mgctx
Definition variables.h:546
BoundaryFaceConfig::mathematical_type
BCType mathematical_type
Definition variables.h:336
UserCtx::lJAj
Vec lJAj
Definition variables.h:882
UserCtx::lKAj
Vec lKAj
Definition variables.h:883
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
BoundaryFaceConfig::handler
BoundaryCondition * handler
Definition variables.h:339
BCContext
Provides execution context for a boundary condition handler.
Definition variables.h:311
BC_Param_s
A node in a linked list for storing key-value parameters from the bcs.dat file.
Definition variables.h:304
BoundaryFaceConfig
Holds the complete configuration for one of the six boundary faces.
Definition variables.h:334
Cmpnts
A 3D point or vector with PetscScalar components.
Definition variables.h:100
SimCtx
The master context for the entire simulation.
Definition variables.h:651
UserCtx
User-defined context containing data specific to a single computational grid level.
Definition variables.h:830
UserMG
User-level context for managing the entire multigrid hierarchy.
Definition variables.h:542
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
Generated by doxygen 1.9.8