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