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 PROFILE_FUNCTION_END;
27 PetscFunctionReturn(0);
28 }
29
30 // Get the range of cells owned by this rank in each dimension
31 PetscInt owned_start_cell_i, num_owned_cells_on_rank_i;
32 PetscInt owned_start_cell_j, num_owned_cells_on_rank_j;
33 PetscInt owned_start_cell_k, num_owned_cells_on_rank_k;
34
35 ierr = GetOwnedCellRange(info, 0, &owned_start_cell_i, &num_owned_cells_on_rank_i); CHKERRQ(ierr);
36 ierr = GetOwnedCellRange(info, 1, &owned_start_cell_j, &num_owned_cells_on_rank_j); CHKERRQ(ierr);
37 ierr = GetOwnedCellRange(info, 2, &owned_start_cell_k, &num_owned_cells_on_rank_k); CHKERRQ(ierr);
38
39 // Determine the global index of the last cell (0-indexed) in each direction.
40 // Example: If IM_nodes_global = 11 (nodes 0-10), there are 10 cells (0-9). Last cell index is 9.
41 // Formula: global_nodes - 1 (num cells) - 1 (0-indexed) = global_nodes - 2.
42 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)
43 PetscInt last_global_cell_idx_j = (JM_nodes_global > 1) ? (JM_nodes_global - 2) : -1;
44 PetscInt last_global_cell_idx_k = (KM_nodes_global > 1) ? (KM_nodes_global - 2) : -1;
45
46 switch (user->identifiedInletBCFace) {
47 case BC_FACE_NEG_X: // Inlet on the global I-minimum face (face of cell C_i=0)
48 // Rank services if its first owned node is global node 0 (info->xs == 0),
49 // and it owns cells in I, J, and K directions.
50 if (info->xs == 0 && num_owned_cells_on_rank_i > 0 &&
51 num_owned_cells_on_rank_j > 0 && num_owned_cells_on_rank_k > 0) {
52 *can_service_inlet_out = PETSC_TRUE;
53 }
54 break;
55 case BC_FACE_POS_X: // Inlet on the global I-maximum face (face of cell C_i=last_global_cell_idx_i)
56 // Rank services if it owns the last cell in I-direction,
57 // and has extent in J and K.
58 if (last_global_cell_idx_i >= 0 && /* Check for valid global domain */
59 (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 */
60 num_owned_cells_on_rank_j > 0 && num_owned_cells_on_rank_k > 0) {
61 *can_service_inlet_out = PETSC_TRUE;
62 }
63 break;
64 case BC_FACE_NEG_Y:
65 if (info->ys == 0 && num_owned_cells_on_rank_j > 0 &&
66 num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_k > 0) {
67 *can_service_inlet_out = PETSC_TRUE;
68 }
69 break;
70 case BC_FACE_POS_Y:
71 if (last_global_cell_idx_j >= 0 &&
72 (owned_start_cell_j + num_owned_cells_on_rank_j - 1) == last_global_cell_idx_j &&
73 num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_k > 0) {
74 *can_service_inlet_out = PETSC_TRUE;
75 }
76 break;
77 case BC_FACE_NEG_Z:
78 if (info->zs == 0 && num_owned_cells_on_rank_k > 0 &&
79 num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_j > 0) {
80 *can_service_inlet_out = PETSC_TRUE;
81 }
82 break;
83 case BC_FACE_POS_Z:
84 if (last_global_cell_idx_k >= 0 &&
85 (owned_start_cell_k + num_owned_cells_on_rank_k - 1) == last_global_cell_idx_k &&
86 num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_j > 0) {
87 *can_service_inlet_out = PETSC_TRUE;
88 }
89 break;
90 default:
91 LOG_ALLOW(LOCAL, LOG_WARNING, "[Rank %d]: Unknown inlet face %s.\n", rank_for_logging, BCFaceToString((BCFace)user->identifiedInletBCFace));
92 break;
93 }
94
95 LOG_ALLOW(LOCAL, LOG_TRACE,
96 "[Rank %d] Check Service for Inlet %s:\n"
97 " - Local Domain: starts at cell (%d,%d,%d), has (%d,%d,%d) cells.\n"
98 " - Global Domain: has (%d,%d,%d) nodes, so last cell is (%d,%d,%d).\n",
99 rank_for_logging,
100 BCFaceToString((BCFace)user->identifiedInletBCFace),
101 owned_start_cell_i, owned_start_cell_j, owned_start_cell_k,
102 num_owned_cells_on_rank_i, num_owned_cells_on_rank_j, num_owned_cells_on_rank_k,
103 IM_nodes_global, JM_nodes_global, KM_nodes_global,
104 last_global_cell_idx_i, last_global_cell_idx_j, last_global_cell_idx_k);
105
106 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",
107 rank_for_logging,
108 BCFaceToString((BCFace)user->identifiedInletBCFace),
109 (*can_service_inlet_out) ? "CAN SERVICE" : "CANNOT SERVICE",
110 num_owned_cells_on_rank_i, num_owned_cells_on_rank_j, num_owned_cells_on_rank_k,
111 owned_start_cell_i, owned_start_cell_j, owned_start_cell_k);
112
113 PROFILE_FUNCTION_END;
114
115 PetscFunctionReturn(0);
116}
117
118#undef __FUNCT__
119#define __FUNCT__ "CanRankServiceFace"
120
121/**
122 * @brief Implementation of \ref CanRankServiceFace().
123 * @details Full API contract (arguments, ownership, side effects) is documented with
124 * the header declaration in `include/Boundaries.h`.
125 * @see CanRankServiceFace()
126 */
127PetscErrorCode CanRankServiceFace(const DMDALocalInfo *info, PetscInt IM_nodes_global, PetscInt JM_nodes_global, PetscInt KM_nodes_global,
128 BCFace face_id, PetscBool *can_service_out)
129{
130 PetscErrorCode ierr;
131 PetscMPIInt rank_for_logging;
132 PetscFunctionBeginUser;
133
134 PROFILE_FUNCTION_BEGIN;
135
136 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank_for_logging); CHKERRQ(ierr);
137
138 *can_service_out = PETSC_FALSE; // Default to no service
139
140 // Get the range of cells owned by this rank
141 PetscInt owned_start_cell_i, num_owned_cells_on_rank_i;
142 PetscInt owned_start_cell_j, num_owned_cells_on_rank_j;
143 PetscInt owned_start_cell_k, num_owned_cells_on_rank_k;
144 ierr = GetOwnedCellRange(info, 0, &owned_start_cell_i, &num_owned_cells_on_rank_i); CHKERRQ(ierr);
145 ierr = GetOwnedCellRange(info, 1, &owned_start_cell_j, &num_owned_cells_on_rank_j); CHKERRQ(ierr);
146 ierr = GetOwnedCellRange(info, 2, &owned_start_cell_k, &num_owned_cells_on_rank_k); CHKERRQ(ierr);
147
148 // Determine the global index of the last cell (0-indexed) in each direction.
149 PetscInt last_global_cell_idx_i = (IM_nodes_global > 1) ? (IM_nodes_global - 2) : -1;
150 PetscInt last_global_cell_idx_j = (JM_nodes_global > 1) ? (JM_nodes_global - 2) : -1;
151 PetscInt last_global_cell_idx_k = (KM_nodes_global > 1) ? (KM_nodes_global - 2) : -1;
152
153 switch (face_id) {
154 case BC_FACE_NEG_X:
155 if (info->xs == 0 && num_owned_cells_on_rank_i > 0 &&
156 num_owned_cells_on_rank_j > 0 && num_owned_cells_on_rank_k > 0) {
157 *can_service_out = PETSC_TRUE;
158 }
159 break;
160 case BC_FACE_POS_X:
161 if (last_global_cell_idx_i >= 0 &&
162 (owned_start_cell_i + num_owned_cells_on_rank_i - 1) == last_global_cell_idx_i &&
163 num_owned_cells_on_rank_j > 0 && num_owned_cells_on_rank_k > 0) {
164 *can_service_out = PETSC_TRUE;
165 }
166 break;
167 case BC_FACE_NEG_Y:
168 if (info->ys == 0 && num_owned_cells_on_rank_j > 0 &&
169 num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_k > 0) {
170 *can_service_out = PETSC_TRUE;
171 }
172 break;
173 case BC_FACE_POS_Y:
174 if (last_global_cell_idx_j >= 0 &&
175 (owned_start_cell_j + num_owned_cells_on_rank_j - 1) == last_global_cell_idx_j &&
176 num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_k > 0) {
177 *can_service_out = PETSC_TRUE;
178 }
179 break;
180 case BC_FACE_NEG_Z:
181 if (info->zs == 0 && num_owned_cells_on_rank_k > 0 &&
182 num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_j > 0) {
183 *can_service_out = PETSC_TRUE;
184 }
185 break;
186 case BC_FACE_POS_Z:
187 if (last_global_cell_idx_k >= 0 &&
188 (owned_start_cell_k + num_owned_cells_on_rank_k - 1) == last_global_cell_idx_k &&
189 num_owned_cells_on_rank_i > 0 && num_owned_cells_on_rank_j > 0) {
190 *can_service_out = PETSC_TRUE;
191 }
192 break;
193 default:
194 LOG_ALLOW(LOCAL, LOG_WARNING, "Rank %d: Unknown face enum %d. \n", rank_for_logging, face_id);
195 break;
196 }
197
198 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d check for face %s: Result=%s. \n",
199 rank_for_logging, BCFaceToString((BCFace)face_id), (*can_service_out ? "TRUE" : "FALSE"));
200
201 PROFILE_FUNCTION_END;
202
203 PetscFunctionReturn(0);
204}
205
206#undef __FUNCT__
207#define __FUNCT__ "GetDeterministicFaceGridLocation"
208
209/**
210 * @brief Internal helper implementation: `GetDeterministicFaceGridLocation()`.
211 * @details Local to this translation unit.
212 */
213PetscErrorCode GetDeterministicFaceGridLocation(
214 UserCtx *user, const DMDALocalInfo *info,
215 PetscInt xs_gnode_rank, PetscInt ys_gnode_rank, PetscInt zs_gnode_rank,
216 PetscInt IM_cells_global, PetscInt JM_cells_global, PetscInt KM_cells_global,
217 PetscInt64 particle_global_id,
218 PetscInt *ci_metric_lnode_out, PetscInt *cj_metric_lnode_out, PetscInt *ck_metric_lnode_out,
219 PetscReal *xi_metric_logic_out, PetscReal *eta_metric_logic_out, PetscReal *zta_metric_logic_out,
220 PetscBool *placement_successful_out)
221{
222 SimCtx *simCtx = user->simCtx;
223 PetscReal global_logic_i = 0.0, global_logic_j = 0.0, global_logic_k = 0.0;
224 PetscErrorCode ierr;
225 PetscMPIInt rank_for_logging;
226
227 PetscFunctionBeginUser;
228 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank_for_logging); CHKERRQ(ierr);
229
230 *placement_successful_out = PETSC_FALSE; // Default to failure
231
232 // --- Step 1: Configuration and Input Validation ---
233
234 // *** Hardcoded number of grid layers. Change this value to alter the pattern. ***
235 const PetscInt grid_layers = 2;
236
237 LOG_ALLOW(LOCAL, LOG_DEBUG,
238 "[Rank %d] Placing particle %lld on face %s with grid_layers=%d in global domain (%d,%d,%d) cells.\n",
239 rank_for_logging, (long long)particle_global_id, BCFaceToString(user->identifiedInletBCFace), grid_layers,
240 IM_cells_global, JM_cells_global, KM_cells_global);
241
242 const char *face_name = BCFaceToString(user->identifiedInletBCFace);
243
244 // Fatal Error Checks: Ensure the requested grid is geometrically possible.
245 // The total layers from opposite faces (2 * grid_layers) must be less than the domain size.
246 switch (user->identifiedInletBCFace) {
247 case BC_FACE_NEG_X: case BC_FACE_POS_X:
248 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);
249 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);
250 break;
251 case BC_FACE_NEG_Y: case BC_FACE_POS_Y:
252 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);
253 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);
254 break;
255 case BC_FACE_NEG_Z: case BC_FACE_POS_Z:
256 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);
257 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);
258 break;
259 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid identifiedInletBCFace specified: %d", user->identifiedInletBCFace);
260 }
261
262 const PetscInt num_lines_total = 4 * grid_layers;
263 if (simCtx->np < num_lines_total) {
264 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);
265 }
266 if (simCtx->np > 0 && simCtx->np % num_lines_total != 0) {
267 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);
268 }
269
270 // --- Step 2: Map global particle ID to a line and a point on that line ---
271 if (simCtx->np == 0) PetscFunctionReturn(0); // Nothing to do
272
273 LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Distributing %lld particles over %d lines on face %s.\n",
274 rank_for_logging, (long long)simCtx->np, num_lines_total, face_name);
275
276 const PetscInt points_per_line = PetscMax(1, simCtx->np / num_lines_total);
277 PetscInt line_index = particle_global_id / points_per_line;
278 PetscInt point_index_on_line = particle_global_id % points_per_line;
279 line_index = PetscMin(line_index, num_lines_total - 1); // Clamp to handle uneven division
280
281 // Decode the line_index into an edge group (0-3) and a layer within that group (0 to grid_layers-1)
282 const PetscInt edge_group = line_index / grid_layers;
283 const PetscInt layer_index = line_index % grid_layers;
284
285 // --- Step 3: Calculate placement coordinates based on the decoded indices ---
286 const PetscReal layer_spacing_norm_i = (IM_cells_global > 0) ? 1.0 / (PetscReal)IM_cells_global : 0.0;
287 const PetscReal layer_spacing_norm_j = (JM_cells_global > 0) ? 1.0 / (PetscReal)JM_cells_global : 0.0;
288 const PetscReal layer_spacing_norm_k = (KM_cells_global > 0) ? 1.0 / (PetscReal)KM_cells_global : 0.0;
289
290 // Grid-aware epsilon: scale with minimum cell size to keep particles away from rank boundaries
291 const PetscReal min_layer_spacing = PetscMin(layer_spacing_norm_i, PetscMin(layer_spacing_norm_j, layer_spacing_norm_k));
292 const PetscReal epsilon = 0.5 * min_layer_spacing; // Keep particles 10% of cell width from boundaries
293
294 PetscReal variable_coord; // The coordinate that varies along a line
295 if (points_per_line <= 1) {
296 variable_coord = 0.5; // Place single point in the middle
297 } else {
298 variable_coord = ((PetscReal)point_index_on_line + 0.5)/ (PetscReal)(points_per_line);
299 }
300 variable_coord = PetscMin(1.0 - epsilon, PetscMax(epsilon, variable_coord)); // Clamp within [eps, 1-eps]
301
302 // Main logic switch to determine the three global logical coordinates
303 switch (user->identifiedInletBCFace) {
304 case BC_FACE_NEG_X:
305 global_logic_i = 0.5 * layer_spacing_norm_i; // Place near the face, in the middle of the first cell
306 if (edge_group == 0) { global_logic_j = (PetscReal)layer_index * layer_spacing_norm_j + epsilon; global_logic_k = variable_coord; }
307 else if (edge_group == 1) { global_logic_j = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_j) - epsilon; global_logic_k = variable_coord; }
308 else if (edge_group == 2) { global_logic_k = (PetscReal)layer_index * layer_spacing_norm_k + epsilon; global_logic_j = variable_coord; }
309 else /* edge_group == 3 */ { global_logic_k = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_k) - epsilon; global_logic_j = variable_coord; }
310 break;
311 case BC_FACE_POS_X:
312 global_logic_i = 1.0 - (0.5 * layer_spacing_norm_i); // Place near the face, in the middle of the last cell
313 if (edge_group == 0) { global_logic_j = (PetscReal)layer_index * layer_spacing_norm_j + epsilon; global_logic_k = variable_coord; }
314 else if (edge_group == 1) { global_logic_j = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_j) - epsilon; global_logic_k = variable_coord; }
315 else if (edge_group == 2) { global_logic_k = (PetscReal)layer_index * layer_spacing_norm_k + epsilon; global_logic_j = variable_coord; }
316 else /* edge_group == 3 */ { global_logic_k = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_k) - epsilon; global_logic_j = variable_coord; }
317 break;
318 case BC_FACE_NEG_Y:
319 global_logic_j = 0.5 * layer_spacing_norm_j;
320 if (edge_group == 0) { global_logic_i = (PetscReal)layer_index * layer_spacing_norm_i + epsilon; global_logic_k = variable_coord; }
321 else if (edge_group == 1) { global_logic_i = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_i) - epsilon; global_logic_k = variable_coord; }
322 else if (edge_group == 2) { global_logic_k = (PetscReal)layer_index * layer_spacing_norm_k + epsilon; global_logic_i = variable_coord; }
323 else /* edge_group == 3 */ { global_logic_k = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_k) - epsilon; global_logic_i = variable_coord; }
324 break;
325 case BC_FACE_POS_Y:
326 global_logic_j = 1.0 - (0.5 * layer_spacing_norm_j);
327 if (edge_group == 0) { global_logic_i = (PetscReal)layer_index * layer_spacing_norm_i + epsilon; global_logic_k = variable_coord; }
328 else if (edge_group == 1) { global_logic_i = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_i) - epsilon; global_logic_k = variable_coord; }
329 else if (edge_group == 2) { global_logic_k = (PetscReal)layer_index * layer_spacing_norm_k + epsilon; global_logic_i = variable_coord; }
330 else /* edge_group == 3 */ { global_logic_k = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_k) - epsilon; global_logic_i = variable_coord; }
331 break;
332 case BC_FACE_NEG_Z:
333 global_logic_k = 0.5 * layer_spacing_norm_k;
334 if (edge_group == 0) { global_logic_i = (PetscReal)layer_index * layer_spacing_norm_i + epsilon; global_logic_j = variable_coord; }
335 else if (edge_group == 1) { global_logic_i = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_i) - epsilon; global_logic_j = variable_coord; }
336 else if (edge_group == 2) { global_logic_j = (PetscReal)layer_index * layer_spacing_norm_j + epsilon; global_logic_i = variable_coord; }
337 else /* edge_group == 3 */ { global_logic_j = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_j) - epsilon; global_logic_i = variable_coord; }
338 break;
339 case BC_FACE_POS_Z:
340 global_logic_k = 1.0 - (0.5 * layer_spacing_norm_k);
341 if (edge_group == 0) { global_logic_i = (PetscReal)layer_index * layer_spacing_norm_i + epsilon; global_logic_j = variable_coord; }
342 else if (edge_group == 1) { global_logic_i = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_i) - epsilon; global_logic_j = variable_coord; }
343 else if (edge_group == 2) { global_logic_j = (PetscReal)layer_index * layer_spacing_norm_j + epsilon; global_logic_i = variable_coord; }
344 else /* edge_group == 3 */ { global_logic_j = 1.0 - ((PetscReal)layer_index * layer_spacing_norm_j) - epsilon; global_logic_i = variable_coord; }
345 break;
346 }
347
348 LOG_ALLOW(LOCAL, LOG_TRACE,
349 "[Rank %d] Particle %lld assigned to line %d (edge group %d, layer %d) with variable_coord=%.4f.\n"
350 " -> Global logical coords: (i,j,k) = (%.6f, %.6f, %.6f)\n",
351 rank_for_logging, (long long)particle_global_id, line_index, edge_group, layer_index, variable_coord,
352 global_logic_i, global_logic_j, global_logic_k);
353
354 // --- Step 4: Convert global logical coordinate to global cell index and intra-cell logicals ---
355 PetscReal global_cell_coord_i = global_logic_i * IM_cells_global;
356 PetscInt I_g = (PetscInt)global_cell_coord_i;
357 *xi_metric_logic_out = global_cell_coord_i - I_g;
358
359 PetscReal global_cell_coord_j = global_logic_j * JM_cells_global;
360 PetscInt J_g = (PetscInt)global_cell_coord_j;
361 *eta_metric_logic_out = global_cell_coord_j - J_g;
362
363 PetscReal global_cell_coord_k = global_logic_k * KM_cells_global;
364 PetscInt K_g = (PetscInt)global_cell_coord_k;
365 *zta_metric_logic_out = global_cell_coord_k - K_g;
366
367 // --- Step 5: Check if this rank owns the target cell and finalize outputs ---
368 if ((I_g >= info->xs && I_g < info->xs + info->xm) &&
369 (J_g >= info->ys && J_g < info->ys + info->ym) &&
370 (K_g >= info->zs && K_g < info->zs + info->zm))
371 {
372 // Convert global cell index to the local node index for this rank's DA patch
373 *ci_metric_lnode_out = (I_g - info->xs) + xs_gnode_rank;
374 *cj_metric_lnode_out = (J_g - info->ys) + ys_gnode_rank;
375 *ck_metric_lnode_out = (K_g - info->zs) + zs_gnode_rank;
376 *placement_successful_out = PETSC_TRUE;
377 }
378
379 LOG_ALLOW(LOCAL, LOG_VERBOSE,
380 "[Rank %d] Particle %lld placement %s.\n",
381 rank_for_logging, (long long)particle_global_id,
382 (*placement_successful_out ? "SUCCESSFUL" : "NOT ON THIS RANK"));
383
384 if(*placement_successful_out){
385 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",
386 *ci_metric_lnode_out, *cj_metric_lnode_out, *ck_metric_lnode_out,
387 *xi_metric_logic_out, *eta_metric_logic_out, *zta_metric_logic_out);
388 }
389
390 PetscFunctionReturn(0);
391}
392
393#undef __FUNCT__
394#define __FUNCT__ "GetRandomFCellAndLogicOnInletFace"
395
396/**
397 * @brief Internal helper implementation: `GetRandomCellAndLogicalCoordsOnInletFace()`.
398 * @details Local to this translation unit.
399 */
400PetscErrorCode GetRandomCellAndLogicalCoordsOnInletFace(
401 UserCtx *user, const DMDALocalInfo *info,
402 PetscInt xs_gnode_rank, PetscInt ys_gnode_rank, PetscInt zs_gnode_rank, // Local starting node index (with ghosts) of the rank's DA patch
403 PetscInt IM_nodes_global, PetscInt JM_nodes_global, PetscInt KM_nodes_global,
404 PetscRandom *rand_logic_i_ptr, PetscRandom *rand_logic_j_ptr, PetscRandom *rand_logic_k_ptr,
405 PetscInt *ci_metric_lnode_out, PetscInt *cj_metric_lnode_out, PetscInt *ck_metric_lnode_out,
406 PetscReal *xi_metric_logic_out, PetscReal *eta_metric_logic_out, PetscReal *zta_metric_logic_out)
407{
408 PetscErrorCode ierr = 0;
409 PetscReal r_val_i_sel, r_val_j_sel, r_val_k_sel;
410 PetscInt local_cell_idx_on_face_dim1 = 0; // 0-indexed relative to owned cells on face
411 PetscInt local_cell_idx_on_face_dim2 = 0;
412 PetscMPIInt rank_for_logging;
413
414 PetscFunctionBeginUser;
415
416 PROFILE_FUNCTION_BEGIN;
417
418 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank_for_logging); CHKERRQ(ierr);
419
420 // Get number of cells this rank owns in each dimension (tangential to the face mainly)
421 PetscInt owned_start_cell_i, num_owned_cells_on_rank_i;
422 PetscInt owned_start_cell_j, num_owned_cells_on_rank_j;
423 PetscInt owned_start_cell_k, num_owned_cells_on_rank_k;
424
425 ierr = GetOwnedCellRange(info, 0, &owned_start_cell_i, &num_owned_cells_on_rank_i); CHKERRQ(ierr);
426 ierr = GetOwnedCellRange(info, 1, &owned_start_cell_j, &num_owned_cells_on_rank_j); CHKERRQ(ierr);
427 ierr = GetOwnedCellRange(info, 2, &owned_start_cell_k, &num_owned_cells_on_rank_k); CHKERRQ(ierr);
428
429 // Defaults for cell origin node (local index for the rank's DA patch, including ghosts)
430 *ci_metric_lnode_out = xs_gnode_rank; *cj_metric_lnode_out = ys_gnode_rank; *ck_metric_lnode_out = zs_gnode_rank;
431 // Defaults for logical coordinates
432 *xi_metric_logic_out = 0.5; *eta_metric_logic_out = 0.5; *zta_metric_logic_out = 0.5;
433
434 // Index of the last cell (0-indexed) in each global direction
435 PetscInt last_global_cell_idx_i = (IM_nodes_global > 1) ? (IM_nodes_global - 2) : -1;
436 PetscInt last_global_cell_idx_j = (JM_nodes_global > 1) ? (JM_nodes_global - 2) : -1;
437 PetscInt last_global_cell_idx_k = (KM_nodes_global > 1) ? (KM_nodes_global - 2) : -1;
438
439 LOG_ALLOW(LOCAL, LOG_INFO, "PARTICLE_INIT_DEBUG Rank %d: Inlet face %s.\n"
440 " Owned cells (i,j,k): (%d,%d,%d)\n"
441 " Global nodes (I,J,K): (%d,%d,%d)\n"
442 " info->xs,ys,zs (first owned node GLOBAL): (%d,%d,%d)\n"
443 " info->xm,ym,zm (num owned nodes GLOBAL): (%d,%d,%d)\n"
444 " xs_gnode_rank,ys_gnode_rank,zs_gnode_rank (DMDAGetCorners): (%d,%d,%d)\n"
445 " owned_start_cell (i,j,k) GLOBAL: (%d,%d,%d)\n"
446 " last_global_cell_idx (i,j,k): (%d,%d,%d)\n",
447 rank_for_logging, BCFaceToString((BCFace)user->identifiedInletBCFace),
448 num_owned_cells_on_rank_i,num_owned_cells_on_rank_j,num_owned_cells_on_rank_k,
449 IM_nodes_global,JM_nodes_global,KM_nodes_global,
450 info->xs, info->ys, info->zs,
451 info->xm, info->ym, info->zm,
452 xs_gnode_rank,ys_gnode_rank,zs_gnode_rank,
453 owned_start_cell_i, owned_start_cell_j, owned_start_cell_k,
454 last_global_cell_idx_i, last_global_cell_idx_j, last_global_cell_idx_k);
455
456
457 switch (user->identifiedInletBCFace) {
458 case BC_FACE_NEG_X: // Particle on -X face of cell C_0 (origin node N_0)
459 // Cell origin node is the first owned node in I by this rank (global index info->xs).
460 // Its local index within the rank's DA (incl ghosts) is xs_gnode_rank.
461 *ci_metric_lnode_out = xs_gnode_rank;
462 *xi_metric_logic_out = 1.0e-6;
463
464 // Tangential dimensions are J and K. Select an owned cell randomly on this face.
465 // num_owned_cells_on_rank_j/k must be > 0 (checked by CanRankServiceInletFace)
466 ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, &r_val_j_sel); CHKERRQ(ierr);
467 local_cell_idx_on_face_dim1 = (PetscInt)(r_val_j_sel * num_owned_cells_on_rank_j); // Index among owned J-cells
468 local_cell_idx_on_face_dim1 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim1), num_owned_cells_on_rank_j - 1);
469 *cj_metric_lnode_out = ys_gnode_rank + local_cell_idx_on_face_dim1; // Offset from start of rank's J-nodes
470
471 ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, &r_val_k_sel); CHKERRQ(ierr);
472 local_cell_idx_on_face_dim2 = (PetscInt)(r_val_k_sel * num_owned_cells_on_rank_k);
473 local_cell_idx_on_face_dim2 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim2), num_owned_cells_on_rank_k - 1);
474 *ck_metric_lnode_out = zs_gnode_rank + local_cell_idx_on_face_dim2;
475
476 ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, eta_metric_logic_out); CHKERRQ(ierr);
477 ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, zta_metric_logic_out); CHKERRQ(ierr);
478 break;
479
480 case BC_FACE_POS_X: // Particle on +X face of cell C_last_I (origin node N_last_I_origin)
481 // Origin node of the last I-cell is global_node_idx = last_global_cell_idx_i.
482 // Its local index in rank's DA: (last_global_cell_idx_i - info->xs) + xs_gnode_rank
483 *ci_metric_lnode_out = xs_gnode_rank + (last_global_cell_idx_i - info->xs);
484 *xi_metric_logic_out = 1.0 - 1.0e-6;
485
486 ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, &r_val_j_sel); CHKERRQ(ierr);
487 local_cell_idx_on_face_dim1 = (PetscInt)(r_val_j_sel * num_owned_cells_on_rank_j);
488 local_cell_idx_on_face_dim1 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim1), num_owned_cells_on_rank_j - 1);
489 *cj_metric_lnode_out = ys_gnode_rank + local_cell_idx_on_face_dim1;
490
491 ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, &r_val_k_sel); CHKERRQ(ierr);
492 local_cell_idx_on_face_dim2 = (PetscInt)(r_val_k_sel * num_owned_cells_on_rank_k);
493 local_cell_idx_on_face_dim2 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim2), num_owned_cells_on_rank_k - 1);
494 *ck_metric_lnode_out = zs_gnode_rank + local_cell_idx_on_face_dim2;
495
496 ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, eta_metric_logic_out); CHKERRQ(ierr);
497 ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, zta_metric_logic_out); CHKERRQ(ierr);
498 break;
499 // ... (Cases for Y and Z faces, following the same pattern) ...
500 case BC_FACE_NEG_Y:
501 *cj_metric_lnode_out = ys_gnode_rank;
502 *eta_metric_logic_out = 1.0e-6;
503 ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, &r_val_i_sel); CHKERRQ(ierr);
504 local_cell_idx_on_face_dim1 = (PetscInt)(r_val_i_sel * num_owned_cells_on_rank_i);
505 local_cell_idx_on_face_dim1 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim1), num_owned_cells_on_rank_i - 1);
506 *ci_metric_lnode_out = xs_gnode_rank + local_cell_idx_on_face_dim1;
507 ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, &r_val_k_sel); CHKERRQ(ierr);
508 local_cell_idx_on_face_dim2 = (PetscInt)(r_val_k_sel * num_owned_cells_on_rank_k);
509 local_cell_idx_on_face_dim2 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim2), num_owned_cells_on_rank_k - 1);
510 *ck_metric_lnode_out = zs_gnode_rank + local_cell_idx_on_face_dim2;
511 ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, xi_metric_logic_out); CHKERRQ(ierr);
512 ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, zta_metric_logic_out); CHKERRQ(ierr);
513 break;
514 case BC_FACE_POS_Y:
515 *cj_metric_lnode_out = ys_gnode_rank + (last_global_cell_idx_j - info->ys);
516 *eta_metric_logic_out = 1.0 - 1.0e-6;
517 ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, &r_val_i_sel); CHKERRQ(ierr);
518 local_cell_idx_on_face_dim1 = (PetscInt)(r_val_i_sel * num_owned_cells_on_rank_i);
519 local_cell_idx_on_face_dim1 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim1), num_owned_cells_on_rank_i - 1);
520 *ci_metric_lnode_out = xs_gnode_rank + local_cell_idx_on_face_dim1;
521 ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, &r_val_k_sel); CHKERRQ(ierr);
522 local_cell_idx_on_face_dim2 = (PetscInt)(r_val_k_sel * num_owned_cells_on_rank_k);
523 local_cell_idx_on_face_dim2 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim2), num_owned_cells_on_rank_k - 1);
524 *ck_metric_lnode_out = zs_gnode_rank + local_cell_idx_on_face_dim2;
525 ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, xi_metric_logic_out); CHKERRQ(ierr);
526 ierr = PetscRandomGetValueReal(*rand_logic_k_ptr, zta_metric_logic_out); CHKERRQ(ierr);
527 break;
528 case BC_FACE_NEG_Z: // Your example case
529 *ck_metric_lnode_out = zs_gnode_rank; // Cell origin is the first owned node in K by this rank
530 *zta_metric_logic_out = 1.0e-6; // Place particle slightly inside this cell from its -Z face
531 // Tangential dimensions are I and J
532 ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, &r_val_i_sel); CHKERRQ(ierr);
533 local_cell_idx_on_face_dim1 = (PetscInt)(r_val_i_sel * num_owned_cells_on_rank_i);
534 local_cell_idx_on_face_dim1 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim1), num_owned_cells_on_rank_i - 1);
535 *ci_metric_lnode_out = xs_gnode_rank + local_cell_idx_on_face_dim1;
536
537 ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, &r_val_j_sel); CHKERRQ(ierr);
538 local_cell_idx_on_face_dim2 = (PetscInt)(r_val_j_sel * num_owned_cells_on_rank_j);
539 local_cell_idx_on_face_dim2 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim2), num_owned_cells_on_rank_j - 1);
540 *cj_metric_lnode_out = ys_gnode_rank + local_cell_idx_on_face_dim2;
541
542 ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, xi_metric_logic_out); CHKERRQ(ierr); // Intra-cell logical for I
543 ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, eta_metric_logic_out); CHKERRQ(ierr); // Intra-cell logical for J
544 break;
545 case BC_FACE_POS_Z:
546 *ck_metric_lnode_out = zs_gnode_rank + (last_global_cell_idx_k - info->zs);
547 *zta_metric_logic_out = 1.0 - 1.0e-6;
548 ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, &r_val_i_sel); CHKERRQ(ierr);
549 local_cell_idx_on_face_dim1 = (PetscInt)(r_val_i_sel * num_owned_cells_on_rank_i);
550 local_cell_idx_on_face_dim1 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim1), num_owned_cells_on_rank_i - 1);
551 *ci_metric_lnode_out = xs_gnode_rank + local_cell_idx_on_face_dim1;
552 ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, &r_val_j_sel); CHKERRQ(ierr);
553 local_cell_idx_on_face_dim2 = (PetscInt)(r_val_j_sel * num_owned_cells_on_rank_j);
554 local_cell_idx_on_face_dim2 = PetscMin(PetscMax(0, local_cell_idx_on_face_dim2), num_owned_cells_on_rank_j - 1);
555 *cj_metric_lnode_out = ys_gnode_rank + local_cell_idx_on_face_dim2;
556 ierr = PetscRandomGetValueReal(*rand_logic_i_ptr, xi_metric_logic_out); CHKERRQ(ierr);
557 ierr = PetscRandomGetValueReal(*rand_logic_j_ptr, eta_metric_logic_out); CHKERRQ(ierr);
558 break;
559 default:
560 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "GetRandomCellAndLogicOnInletFace: Invalid user->identifiedInletBCFace %d. \n", user->identifiedInletBCFace);
561 }
562
563 PetscReal eps = 1.0e-7;
564 if (user->identifiedInletBCFace == BC_FACE_NEG_X || user->identifiedInletBCFace == BC_FACE_POS_X) {
565 *eta_metric_logic_out = PetscMin(PetscMax(0.0, *eta_metric_logic_out), 1.0 - eps);
566 *zta_metric_logic_out = PetscMin(PetscMax(0.0, *zta_metric_logic_out), 1.0 - eps);
567 } else if (user->identifiedInletBCFace == BC_FACE_NEG_Y || user->identifiedInletBCFace == BC_FACE_POS_Y) {
568 *xi_metric_logic_out = PetscMin(PetscMax(0.0, *xi_metric_logic_out), 1.0 - eps);
569 *zta_metric_logic_out = PetscMin(PetscMax(0.0, *zta_metric_logic_out), 1.0 - eps);
570 } else {
571 *xi_metric_logic_out = PetscMin(PetscMax(0.0, *xi_metric_logic_out), 1.0 - eps);
572 *eta_metric_logic_out = PetscMin(PetscMax(0.0, *eta_metric_logic_out), 1.0 - eps);
573 }
574
575 LOG_ALLOW(LOCAL, LOG_VERBOSE, "Rank %d: Target Cell Node =(%d,%d,%d). (xi,et,zt)=(%.2e,%.2f,%.2f). \n",
576 rank_for_logging, *ci_metric_lnode_out, *cj_metric_lnode_out, *ck_metric_lnode_out,
577 *xi_metric_logic_out, *eta_metric_logic_out, *zta_metric_logic_out);
578
579 PROFILE_FUNCTION_END;
580
581 PetscFunctionReturn(0);
582}
583
584
585
586#undef __FUNCT__
587#define __FUNCT__ "ClassifyMomentumRow"
588/**
589 * @brief Implementation of \ref ClassifyMomentumRow().
590 *
591 * Pure index/boundary-type arithmetic: no field reads, no communication.
592 * Precedence matters. A conditioned row is reported first because its explicit
593 * Dirichlet value is more specific than the homogeneous fallback, and a periodic
594 * duplicate is reported before the homogeneous case because the Newton path needs
595 * its representative index to build `F = X_dup - X_rep`.
596 */
597MomentumRowType ClassifyMomentumRow(UserCtx *user, PetscInt i, PetscInt j, PetscInt k,
598 PetscInt component, PetscInt *ri, PetscInt *rj, PetscInt *rk)
599{
600 const PetscInt mx = user->info.mx, my = user->info.my, mz = user->info.mz;
601 const PetscInt coord[3] = {i, j, k};
602 const PetscInt size[3] = {mx, my, mz};
603 const BCFace neg_face[3] = {BC_FACE_NEG_X, BC_FACE_NEG_Y, BC_FACE_NEG_Z};
604 PetscBool periodic[3], periodic_duplicate = PETSC_FALSE;
605 PetscBool residual_zeroed = PETSC_FALSE, conditioned = PETSC_FALSE;
606
607 *ri = i; *rj = j; *rk = k;
608 for (PetscInt axis = 0; axis < 3; ++axis) {
609 periodic[axis] = (PetscBool)(
610 user->boundary_faces[neg_face[axis]].mathematical_type == PERIODIC);
611 if (periodic[axis] && coord[axis] == 0) {
612 periodic_duplicate = PETSC_TRUE;
613 if (axis == 0) *ri = -2;
614 else if (axis == 1) *rj = -2;
615 else *rk = -2;
616 }
617 if (periodic[axis] && coord[axis] == size[axis] - 1) {
618 periodic_duplicate = PETSC_TRUE;
619 if (axis == 0) *ri = mx + 1;
620 else if (axis == 1) *rj = my + 1;
621 else *rk = mz + 1;
622 }
623
624 if (!periodic[axis] && coord[axis] == 0) residual_zeroed = PETSC_TRUE;
625 if (coord[axis] == size[axis] - 1) residual_zeroed = PETSC_TRUE;
626 if (!periodic[axis] && coord[axis] == size[axis] - 2 && component == axis)
627 residual_zeroed = PETSC_TRUE;
628 }
629
630 if (!periodic[component] &&
631 (coord[component] == 0 || coord[component] == size[component] - 2)) {
632 PetscBool tangential_interior = PETSC_TRUE;
633 for (PetscInt axis = 0; axis < 3; ++axis) {
634 if (axis == component) continue;
635 if (coord[axis] < 1 || coord[axis] > size[axis] - 2)
636 tangential_interior = PETSC_FALSE;
637 }
638 conditioned = tangential_interior;
639 }
640
641 if (conditioned) return MOM_ROW_FIXED_CONDITIONED;
642 if (periodic_duplicate) return MOM_ROW_PERIODIC_DUPLICATE;
643 if (residual_zeroed) return MOM_ROW_FIXED_HOMOGENEOUS;
644 return MOM_ROW_PHYSICAL;
645}
646
647#undef __FUNCT__
648#define __FUNCT__ "EnforceRHSBoundaryConditions"
649/**
650 * @brief Implementation of \ref EnforceRHSBoundaryConditions().
651 *
652 * The sweep is deliberately expressed over every owned location rather than over
653 * the six boundary slabs: restating "which indices can be non-physical" here is
654 * exactly the duplication that let the periodic duplicate column go unzeroed.
655 * ClassifyMomentumRow() is a handful of integer comparisons and the walk is a
656 * single pass with no stencil access, which is negligible next to the several
657 * ghosted stencil passes ComputeRHS() has already made over the same range.
658 * MomentumNewtonKrylov_ApplyConstraints() walks the same range the same way.
659 */
660PetscErrorCode EnforceRHSBoundaryConditions(UserCtx *user)
661{
662 PetscErrorCode ierr;
663 DMDALocalInfo info = user->info;
664 Cmpnts ***rhs;
665
666 PetscFunctionBeginUser;
667 PROFILE_FUNCTION_BEGIN;
668
669 // Get a writable pointer to the local data of the global RHS vector.
670 ierr = DMDAVecGetArray(user->fda, user->Rhs, &rhs); CHKERRQ(ierr);
671
672 for (PetscInt k = info.zs; k < info.zs + info.zm; k++) {
673 for (PetscInt j = info.ys; j < info.ys + info.ym; j++) {
674 for (PetscInt i = info.xs; i < info.xs + info.xm; i++) {
675 PetscScalar *row = &rhs[k][j][i].x; /* .x/.y/.z are contiguous */
676 for (PetscInt component = 0; component < 3; component++) {
677 PetscInt ri, rj, rk;
678 if (ClassifyMomentumRow(user, i, j, k, component, &ri, &rj, &rk) != MOM_ROW_PHYSICAL)
679 row[component] = 0.0;
680 }
681 }
682 }
683 }
684
685 // --- Release the pointer to the local data ---
686 ierr = DMDAVecRestoreArray(user->fda, user->Rhs, &rhs); CHKERRQ(ierr);
687
688 LOG_ALLOW(LOCAL, LOG_TRACE, "Rank %d, Block %d: Finished enforcing RHS boundary conditions.\n",
689 user->simCtx->rank, user->_this);
690
691 PROFILE_FUNCTION_END;
692
693 PetscFunctionReturn(0);
694}
695
696#undef __FUNCT__
697#define __FUNCT__ "BoundaryCondition_Create"
698/**
699 * @brief Internal helper implementation: `BoundaryCondition_Create()`.
700 * @details Local to this translation unit.
701 */
702
703PetscErrorCode BoundaryCondition_Create(BCHandlerType handler_type, BoundaryCondition **new_bc_ptr)
704{
705 PetscErrorCode ierr;
706 PetscFunctionBeginUser;
707
708 const char* handler_name = BCHandlerTypeToString(handler_type);
709 LOG_ALLOW(LOCAL, LOG_DEBUG, "Factory called for handler type %s. \n", handler_name);
710
711 ierr = PetscMalloc1(1, new_bc_ptr); CHKERRQ(ierr);
712 BoundaryCondition *bc = *new_bc_ptr;
713
714 bc->type = handler_type;
715 bc->priority = -1; // Default priority; can be overridden in specific handlers
716 bc->data = NULL;
717 bc->Initialize = NULL;
718 bc->PreStep = NULL;
719 bc->Apply = NULL;
720 bc->PostStep = NULL;
721 bc->UpdateUbcs = NULL;
722 bc->Destroy = NULL;
723
724 LOG_ALLOW(LOCAL, LOG_DEBUG, "Allocated generic handler object at address %p.\n", (void*)bc);
725
726 switch (handler_type) {
727
728 case BC_HANDLER_OUTLET_CONSERVATION:
729 LOG_ALLOW(LOCAL, LOG_DEBUG, "Dispatching to Create_OutletConservation().\n");
730 ierr = Create_OutletConservation(bc); CHKERRQ(ierr);
731 break;
732
733 case BC_HANDLER_WALL_NOSLIP:
734 LOG_ALLOW(LOCAL, LOG_DEBUG, "Dispatching to Create_WallNoSlip().\n");
735 ierr = Create_WallNoSlip(bc); CHKERRQ(ierr);
736 break;
737
738 case BC_HANDLER_INLET_CONSTANT_VELOCITY:
739 LOG_ALLOW(LOCAL, LOG_DEBUG, "Dispatching to Create_InletConstantVelocity().\n");
740 ierr = Create_InletConstantVelocity(bc); CHKERRQ(ierr);
741 break;
742
743 case BC_HANDLER_PERIODIC_GEOMETRIC:
744 LOG_ALLOW(LOCAL,LOG_DEBUG,"Dispatching to Create_PeriodicGeometric().\n");
745 ierr = Create_PeriodicGeometric(bc);
746 break;
747
748 case BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX:
749 LOG_ALLOW(LOCAL,LOG_DEBUG,"Dispatching to Create_PeriodicDrivenConstant().\n");
750 ierr = Create_PeriodicDrivenConstant(bc);
751 break;
752
753 case BC_HANDLER_PERIODIC_DRIVEN_INITIAL_FLUX:
754 LOG_ALLOW(LOCAL,LOG_DEBUG,"Dispatching to Create_PeriodicDrivenInitial().\n");
755 ierr = Create_PeriodicDrivenInitial(bc);
756 break;
757
758 case BC_HANDLER_INLET_PARABOLIC:
759 LOG_ALLOW(LOCAL, LOG_DEBUG, "Dispatching to Create_InletParabolicProfile().\n");
760 ierr = Create_InletParabolicProfile(bc); CHKERRQ(ierr);
761 break;
762
763 case BC_HANDLER_INLET_PROFILE_FROM_FILE:
764 LOG_ALLOW(LOCAL, LOG_DEBUG, "Dispatching to Create_InletProfileFromFile().\n");
765 ierr = Create_InletProfileFromFile(bc); CHKERRQ(ierr);
766 break;
767 //Add cases for other handlers here in future phases
768
769 default:
770 LOG_ALLOW(GLOBAL, LOG_ERROR, "Handler type (%s) is not recognized or implemented in the factory.\n", handler_name);
771 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Boundary handler type %d (%s) not recognized in factory.\n", handler_type, handler_name);
772 }
773
774 if(bc->priority < 0) {
775 LOG_ALLOW(GLOBAL, LOG_ERROR, "Handler type %d (%s) did not set a valid priority during creation.\n", handler_type, handler_name);
776 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);
777 }
778
779 LOG_ALLOW(LOCAL, LOG_DEBUG, "Successfully created and configured handler for %s.\n", handler_name);
780 PetscFunctionReturn(0);
781}
782
783#undef __FUNCT__
784#define __FUNCT__ "BoundarySystem_Validate"
785/**
786 * @brief Internal helper implementation: `BoundarySystem_Validate()`.
787 * @details Local to this translation unit.
788 */
789PetscErrorCode BoundarySystem_Validate(UserCtx *user)
790{
791 PetscErrorCode ierr;
792 const BCFace neg_faces[3] = {BC_FACE_NEG_X, BC_FACE_NEG_Y, BC_FACE_NEG_Z};
793 const BCFace pos_faces[3] = {BC_FACE_POS_X, BC_FACE_POS_Y, BC_FACE_POS_Z};
794 const char axis_names[3] = {'X', 'Y', 'Z'};
795 DMBoundaryType bx, by, bz;
796 PetscBool dm_periodic[3];
797 PetscFunctionBeginUser;
798
799 LOG_ALLOW(GLOBAL, LOG_INFO, "Validating parsed boundary condition configuration...\n");
800 ierr = DMDAGetInfo(user->da, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
801 &bx, &by, &bz, NULL); CHKERRQ(ierr);
802 dm_periodic[0] = (PetscBool)(bx == DM_BOUNDARY_PERIODIC);
803 dm_periodic[1] = (PetscBool)(by == DM_BOUNDARY_PERIODIC);
804 dm_periodic[2] = (PetscBool)(bz == DM_BOUNDARY_PERIODIC);
805
806 // --- Rule Set 1: Geometric periodic faces must be paired and match the DM topology. ---
807 for (PetscInt axis = 0; axis < 3; axis++) {
808 const PetscBool neg_periodic =
809 user->boundary_faces[neg_faces[axis]].mathematical_type == PERIODIC;
810 const PetscBool pos_periodic =
811 user->boundary_faces[pos_faces[axis]].mathematical_type == PERIODIC;
812
813 PetscCheck(neg_periodic == pos_periodic, PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT,
814 "Configuration Error: Periodic boundaries in the %c direction must be paired; "
815 "%s is %s while %s is %s.",
816 axis_names[axis],
817 BCFaceToString(neg_faces[axis]), neg_periodic ? "PERIODIC" : "not periodic",
818 BCFaceToString(pos_faces[axis]), pos_periodic ? "PERIODIC" : "not periodic");
819 PetscCheck(dm_periodic[axis] == neg_periodic, PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT,
820 "Configuration Error: The %c-direction DM periodic flag (%d) does not match "
821 "the paired boundary configuration (%s).",
822 axis_names[axis], (int)dm_periodic[axis], neg_periodic ? "PERIODIC" : "not periodic");
823 }
824
825 // --- Rule Set 2: Driven Flow Handler Consistency ---
826 // This specialized validator will check all rules related to driven flow handlers.
827 ierr = Validate_DrivenFlowConfiguration(user); CHKERRQ(ierr);
828
829 // --- Rule Set 3: (Future Extension) Overset Interface Consistency ---
830 // ierr = Validate_OversetConfiguration(user); CHKERRQ(ierr);
831
832 LOG_ALLOW(GLOBAL, LOG_INFO, "Boundary configuration is valid.\n");
833
834 PetscFunctionReturn(0);
835}
836
837//================================================================================
838//
839// PUBLIC MASTER SETUP FUNCTION
840//
841//================================================================================
842#undef __FUNCT__
843#define __FUNCT__ "BoundarySystem_Initialize"
844/**
845 * @brief Implementation of \ref BoundarySystem_Initialize().
846 * @details Full API contract (arguments, ownership, side effects) is documented with
847 * the header declaration in `include/Boundaries.h`.
848 * @see BoundarySystem_Initialize()
849 */
850PetscErrorCode BoundarySystem_Initialize(UserCtx *user, const char *bcs_filename)
851{
852 PetscErrorCode ierr;
853 PetscFunctionBeginUser;
854
855 LOG_ALLOW(GLOBAL, LOG_INFO, "Starting creation and initialization of all boundary handlers.\n");
856
857 // =========================================================================
858 // Step 0: Clear any existing boundary handlers (if re-initializing).
859 // This ensures no memory leaks if this function is called multiple times.
860 // =========================================================================
861 for (int i = 0; i < 6; i++) {
862 BoundaryFaceConfig *face_cfg = &user->boundary_faces[i];
863 if (face_cfg->handler) {
864 LOG_ALLOW(LOCAL, LOG_DEBUG, "Destroying existing handler on Face %s before re-initialization.\n", BCFaceToString((BCFace)i));
865 if (face_cfg->handler->Destroy) {
866 ierr = face_cfg->handler->Destroy(face_cfg->handler); CHKERRQ(ierr);
867 }
868 ierr = PetscFree(face_cfg->handler); CHKERRQ(ierr);
869 face_cfg->handler = NULL;
870 }
871 }
872 // =========================================================================
873
874 // Step 0.1: Initiate flux sums to zero
875 user->simCtx->FluxInSum = 0.0;
876 user->simCtx->FluxOutSum = 0.0;
877 user->simCtx->FarFluxInSum = 0.0;
878 user->simCtx->FarFluxOutSum = 0.0;
879 // =========================================================================
880
881 // Step 1: Parse the configuration file to determine user intent.
882 // This function, defined in io.c, populates the configuration enums and parameter
883 // lists within the user->boundary_faces array on all MPI ranks.
884 ierr = ParseAllBoundaryConditions(user, bcs_filename); CHKERRQ(ierr);
885 LOG_ALLOW(GLOBAL, LOG_INFO, "Configuration file '%s' parsed successfully.\n", bcs_filename);
886
887 // Step 1.1: Validate the parsed configuration to ensure there are no Boundary Condition conflicts
888 ierr = BoundarySystem_Validate(user); CHKERRQ(ierr);
889
890 // Step 2: Create and Initialize the handler object for each of the 6 faces.
891 for (int i = 0; i < 6; i++) {
892 BoundaryFaceConfig *face_cfg = &user->boundary_faces[i];
893
894 const char *face_name = BCFaceToString(face_cfg->face_id);
895 const char *type_name = BCTypeToString(face_cfg->mathematical_type);
896 const char *handler_name = BCHandlerTypeToString(face_cfg->handler_type);
897
898 LOG_ALLOW(LOCAL, LOG_DEBUG, "Creating handler for Face %s with Type %s and handler '%s'.\n", face_name, type_name,handler_name);
899
900 // Use the private factory to construct the correct handler object based on the parsed type.
901 // The factory returns a pointer to the new handler object, which we store in the config struct.
902 ierr = BoundaryCondition_Create(face_cfg->handler_type, &face_cfg->handler); CHKERRQ(ierr);
903
904 // Step 3: Call the specific Initialize() method for the newly created handler.
905 // This allows the handler to perform its own setup, like reading parameters from the
906 // face_cfg->params list and setting the initial field values on its face.
907 if (face_cfg->handler && face_cfg->handler->Initialize) {
908 LOG_ALLOW(LOCAL, LOG_DEBUG, "Calling Initialize() method for handler %s(%s) on Face %s.\n",type_name,handler_name,face_name);
909
910 // Prepare the context needed by the Initialize() function.
911 BCContext ctx = {
912 .user = user,
913 .face_id = face_cfg->face_id,
914 .global_inflow_sum = &user->simCtx->FluxInSum, // Global flux sums are not relevant during initialization.
915 .global_outflow_sum = &user->simCtx->FluxOutSum,
916 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
917 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
918 };
919
920 ierr = face_cfg->handler->Initialize(face_cfg->handler, &ctx); CHKERRQ(ierr);
921 } else {
922 LOG_ALLOW(LOCAL, LOG_DEBUG, "Handler %s(%s) for Face %s has no Initialize() method, skipping.\n", type_name,handler_name,face_name);
923 }
924 }
925 // =========================================================================
926 // NO SYNCHRONIZATION NEEDED HERE
927 // =========================================================================
928 // Initialize() only reads parameters and allocates memory.
929 // It does NOT modify field values (Ucat, Ucont, Ubcs).
930 // Field values are set by:
931 // 1. Initial conditions (before this function)
932 // 2. Apply() during timestepping (after this function)
933 // The first call to ApplyBoundaryConditions() will handle synchronization.
934 // =========================================================================
935
936 LOG_ALLOW(GLOBAL, LOG_INFO, "All boundary handlers created and initialized successfully.\n");
937 PetscFunctionReturn(0);
938}
939
940
941#undef __FUNCT__
942#define __FUNCT__ "PropagateBoundaryConfigToCoarserLevels"
943/**
944 * @brief Internal helper implementation: `PropagateBoundaryConfigToCoarserLevels()`.
945 * @details Local to this translation unit.
946 */
947PetscErrorCode PropagateBoundaryConfigToCoarserLevels(SimCtx *simCtx)
948{
949 PetscErrorCode ierr;
950 UserMG *usermg = &simCtx->usermg;
951
952 PetscFunctionBeginUser;
953 PROFILE_FUNCTION_BEGIN;
954
955 LOG_ALLOW(GLOBAL, LOG_INFO, "Propagating BC configuration from finest to coarser multigrid levels...\n");
956
957 // Loop from second-finest down to coarsest
958 for (PetscInt level = usermg->mglevels - 2; level >= 0; level--) {
959 for (PetscInt bi = 0; bi < simCtx->block_number; bi++) {
960 UserCtx *user_coarse = &usermg->mgctx[level].user[bi];
961 UserCtx *user_fine = &usermg->mgctx[level + 1].user[bi];
962
963 LOG_ALLOW_SYNC(LOCAL, LOG_DEBUG, "Rank %d: Copying BC config from level %d to level %d, block %d\n",
964 simCtx->rank, level + 1, level, bi);
965
966 // Copy the 6 boundary face configurations
967 for (int face_i = 0; face_i < 6; face_i++) {
968 user_coarse->boundary_faces[face_i].face_id = user_fine->boundary_faces[face_i].face_id;
969 user_coarse->boundary_faces[face_i].mathematical_type = user_fine->boundary_faces[face_i].mathematical_type;
970 user_coarse->boundary_faces[face_i].handler_type = user_fine->boundary_faces[face_i].handler_type;
971
972 // Copy parameter list (deep copy)
973 FreeBC_ParamList(user_coarse->boundary_faces[face_i].params); // Clear any existing
974 user_coarse->boundary_faces[face_i].params = NULL;
975
976 BC_Param **dst_next = &user_coarse->boundary_faces[face_i].params;
977 for (BC_Param *src = user_fine->boundary_faces[face_i].params; src; src = src->next) {
978 BC_Param *new_param;
979 ierr = PetscMalloc1(1, &new_param); CHKERRQ(ierr);
980 ierr = PetscStrallocpy(src->key, &new_param->key); CHKERRQ(ierr);
981 ierr = PetscStrallocpy(src->value, &new_param->value); CHKERRQ(ierr);
982 new_param->next = NULL;
983 *dst_next = new_param;
984 dst_next = &new_param->next;
985 }
986
987 // IMPORTANT: Do NOT create handler objects for coarser levels
988 // Handlers are only needed at finest level for timestepping Apply() calls
989 user_coarse->boundary_faces[face_i].handler = NULL;
990 }
991
992 // Propagate the particle inlet lookup fields to coarse levels as well.
993 user_coarse->inletFaceDefined = user_fine->inletFaceDefined;
994 user_coarse->identifiedInletBCFace = user_fine->identifiedInletBCFace;
995 }
996 }
997
998 LOG_ALLOW(GLOBAL, LOG_INFO, "BC configuration propagation complete.\n");
999
1000 PROFILE_FUNCTION_END;
1001 PetscFunctionReturn(0);
1002}
1003
1004//================================================================================
1005//
1006// PUBLIC MASTER TIME-STEP FUNCTION
1007//
1008//================================================================================
1009
1010#undef __FUNCT__
1011#define __FUNCT__ "BoundarySystem_ExecuteStep"
1012/**
1013 * @brief Implementation of \ref BoundarySystem_ExecuteStep().
1014 * @details Full API contract (arguments, ownership, side effects) is documented with
1015 * the header declaration in `include/Boundaries.h`.
1016 * @see BoundarySystem_ExecuteStep()
1017 */
1018PetscErrorCode BoundarySystem_ExecuteStep(UserCtx *user)
1019{
1020 PetscErrorCode ierr;
1021 PetscFunctionBeginUser;
1022 PROFILE_FUNCTION_BEGIN;
1023
1024 LOG_ALLOW(LOCAL, LOG_DEBUG, "Starting.\n");
1025
1026 // =========================================================================
1027 // PRIORITY 0: INLETS
1028 // =========================================================================
1029
1030 PetscReal local_inflow_pre = 0.0;
1031 PetscReal local_inflow_post = 0.0;
1032 PetscReal global_inflow_pre = 0.0;
1033 PetscReal global_inflow_post = 0.0;
1034 PetscInt num_handlers[3] = {0,0,0};
1035
1036 LOG_ALLOW(LOCAL, LOG_TRACE, " (INLETS): Begin.\n");
1037
1038 // Phase 1: PreStep - Preparation (e.g., calculate profiles, read files)
1039 for (int i = 0; i < 6; i++) {
1040 BoundaryCondition *handler = user->boundary_faces[i].handler;
1041 if (!handler || handler->priority != BC_PRIORITY_INLET) continue;
1042 if (!handler->PreStep) continue;
1043
1044 num_handlers[0]++;
1045 BCContext ctx = {
1046 .user = user,
1047 .face_id = (BCFace)i,
1048 .global_inflow_sum = NULL,
1049 .global_outflow_sum = NULL,
1050 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1051 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1052 };
1053
1054 LOG_ALLOW(LOCAL, LOG_TRACE, " PreStep: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1055 ierr = handler->PreStep(handler, &ctx, &local_inflow_pre, NULL); CHKERRQ(ierr);
1056 }
1057
1058 // Optional: Global communication for PreStep (for debugging)
1059 if (local_inflow_pre != 0.0) {
1060 ierr = MPI_Allreduce(&local_inflow_pre, &global_inflow_pre, 1, MPIU_REAL,
1061 MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
1062 LOG_ALLOW(GLOBAL, LOG_TRACE, " PreStep predicted flux: %.6e\n", global_inflow_pre);
1063 }
1064
1065 // Phase 2: Apply - Set boundary conditions
1066 for (int i = 0; i < 6; i++) {
1067 BoundaryCondition *handler = user->boundary_faces[i].handler;
1068 if (!handler || handler->priority != BC_PRIORITY_INLET) continue;
1069 if(!handler->Apply) continue; // For example Periodic BCs
1070
1071 num_handlers[1]++;
1072
1073 BCContext ctx = {
1074 .user = user,
1075 .face_id = (BCFace)i,
1076 .global_inflow_sum = NULL,
1077 .global_outflow_sum = NULL,
1078 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1079 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1080 };
1081
1082 LOG_ALLOW(LOCAL, LOG_TRACE, " Apply: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1083 ierr = handler->Apply(handler, &ctx); CHKERRQ(ierr);
1084 }
1085
1086 // Phase 3: PostStep - Measure actual flux
1087 for (int i = 0; i < 6; i++) {
1088 BoundaryCondition *handler = user->boundary_faces[i].handler;
1089 if (!handler || handler->priority != BC_PRIORITY_INLET) continue;
1090 if (!handler->PostStep) continue;
1091
1092 num_handlers[2]++;
1093
1094 BCContext ctx = {
1095 .user = user,
1096 .face_id = (BCFace)i,
1097 .global_inflow_sum = NULL,
1098 .global_outflow_sum = NULL,
1099 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1100 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1101 };
1102
1103 LOG_ALLOW(LOCAL, LOG_TRACE, " PostStep: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1104 ierr = handler->PostStep(handler, &ctx, &local_inflow_post, NULL); CHKERRQ(ierr);
1105 }
1106
1107 // Phase 4: Global communication - Sum flux for other priorities to use
1108 ierr = MPI_Allreduce(&local_inflow_post, &global_inflow_post, 1, MPIU_REAL,
1109 MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
1110
1111 // Store for next priority levels
1112 user->simCtx->FluxInSum = global_inflow_post;
1113
1114 LOG_ALLOW(GLOBAL, LOG_INFO,
1115 " (INLETS): %d Prestep(s), %d Application(s), %d Poststep(s), FluxInSum = %.6e\n",
1116 num_handlers[0],num_handlers[1],num_handlers[2], global_inflow_post);
1117
1118 // =========================================================================
1119 // PRIORITY 1: FARFIELD
1120 // =========================================================================
1121
1122 PetscReal local_farfield_in_pre = 0.0;
1123 PetscReal local_farfield_out_pre = 0.0;
1124 PetscReal local_farfield_in_post = 0.0;
1125 PetscReal local_farfield_out_post = 0.0;
1126 PetscReal global_farfield_in_pre = 0.0;
1127 PetscReal global_farfield_out_pre = 0.0;
1128 PetscReal global_farfield_in_post = 0.0;
1129 PetscReal global_farfield_out_post = 0.0;
1130 memset(num_handlers,0,sizeof(num_handlers));
1131
1132 LOG_ALLOW(LOCAL, LOG_TRACE, " (FARFIELD): Begin.\n");
1133
1134 // Phase 1: PreStep - Analyze flow direction, measure initial flux
1135 for (int i = 0; i < 6; i++) {
1136 BoundaryCondition *handler = user->boundary_faces[i].handler;
1137 if (!handler || handler->priority != BC_PRIORITY_FARFIELD) continue;
1138 if (!handler->PreStep) continue;
1139
1140 num_handlers[0]++;
1141 BCContext ctx = {
1142 .user = user,
1143 .face_id = (BCFace)i,
1144 .global_inflow_sum = &user->simCtx->FluxInSum, // Available from Priority 0
1145 .global_outflow_sum = NULL,
1146 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1147 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1148 };
1149
1150 LOG_ALLOW(LOCAL, LOG_TRACE, " PreStep: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1151 ierr = handler->PreStep(handler, &ctx, &local_farfield_in_pre, &local_farfield_out_pre);
1152 CHKERRQ(ierr);
1153 }
1154
1155 // Phase 2: Global communication (optional, for debugging)
1156 if (local_farfield_in_pre != 0.0 || local_farfield_out_pre != 0.0) {
1157 ierr = MPI_Allreduce(&local_farfield_in_pre, &global_farfield_in_pre, 1, MPIU_REAL,
1158 MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
1159 ierr = MPI_Allreduce(&local_farfield_out_pre, &global_farfield_out_pre, 1, MPIU_REAL,
1160 MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
1161
1162 LOG_ALLOW(GLOBAL, LOG_DEBUG,
1163 " Farfield pre-analysis: In=%.6e, Out=%.6e\n",
1164 global_farfield_in_pre, global_farfield_out_pre);
1165 }
1166
1167 // Phase 3: Apply - Set farfield boundary conditions
1168 for (int i = 0; i < 6; i++) {
1169 BoundaryCondition *handler = user->boundary_faces[i].handler;
1170 if (!handler || handler->priority != BC_PRIORITY_FARFIELD) continue;
1171 if(!handler->Apply) continue; // For example Periodic BCs
1172
1173 num_handlers[1]++;
1174
1175 BCContext ctx = {
1176 .user = user,
1177 .face_id = (BCFace)i,
1178 .global_inflow_sum = &user->simCtx->FluxInSum,
1179 .global_outflow_sum = NULL,
1180 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1181 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1182 };
1183
1184 LOG_ALLOW(LOCAL, LOG_TRACE, " Apply: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1185 ierr = handler->Apply(handler, &ctx); CHKERRQ(ierr);
1186 }
1187
1188 // Phase 4: PostStep - Measure actual farfield fluxes
1189 for (int i = 0; i < 6; i++) {
1190 BoundaryCondition *handler = user->boundary_faces[i].handler;
1191 if (!handler || handler->priority != BC_PRIORITY_FARFIELD) continue;
1192 if (!handler->PostStep) continue;
1193
1194 num_handlers[2]++;
1195
1196 BCContext ctx = {
1197 .user = user,
1198 .face_id = (BCFace)i,
1199 .global_inflow_sum = &user->simCtx->FluxInSum,
1200 .global_outflow_sum = NULL,
1201 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1202 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1203 };
1204
1205 LOG_ALLOW(LOCAL, LOG_TRACE, " PostStep: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1206 ierr = handler->PostStep(handler, &ctx, &local_farfield_in_post, &local_farfield_out_post);
1207 CHKERRQ(ierr);
1208 }
1209
1210 // Phase 5: Global communication - Store for outlet priority
1211 if (num_handlers > 0) {
1212 ierr = MPI_Allreduce(&local_farfield_in_post, &global_farfield_in_post, 1, MPIU_REAL,
1213 MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
1214 ierr = MPI_Allreduce(&local_farfield_out_post, &global_farfield_out_post, 1, MPIU_REAL,
1215 MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
1216
1217 // Store for outlet handlers to use
1218 user->simCtx->FarFluxInSum = global_farfield_in_post;
1219 user->simCtx->FarFluxOutSum = global_farfield_out_post;
1220
1221 LOG_ALLOW(GLOBAL, LOG_INFO,
1222 " (FARFIELD): %d Prestep(s), %d Application(s), %d Poststep(s) , InFlux=%.6e, OutFlux=%.6e\n",
1223 num_handlers[0],num_handlers[1],num_handlers[2], global_farfield_in_post, global_farfield_out_post);
1224 } else {
1225 // No farfield handlers - zero out the fluxes
1226 user->simCtx->FarFluxInSum = 0.0;
1227 user->simCtx->FarFluxOutSum = 0.0;
1228 }
1229
1230
1231 // =========================================================================
1232 // PRIORITY 2: WALLS
1233 // =========================================================================
1234
1235 memset(num_handlers,0,sizeof(num_handlers));
1236
1237 LOG_ALLOW(LOCAL, LOG_TRACE, " (WALLS): Begin.\n");
1238
1239 // Phase 1: PreStep - Preparation (usually no-op for walls)
1240 for (int i = 0; i < 6; i++) {
1241 BoundaryCondition *handler = user->boundary_faces[i].handler;
1242 if (!handler || handler->priority != BC_PRIORITY_WALL) continue;
1243 if (!handler->PreStep) continue;
1244
1245 num_handlers[0]++;
1246 BCContext ctx = {
1247 .user = user,
1248 .face_id = (BCFace)i,
1249 .global_inflow_sum = &user->simCtx->FluxInSum,
1250 .global_outflow_sum = NULL,
1251 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1252 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1253 };
1254
1255 LOG_ALLOW(LOCAL, LOG_TRACE, " PreStep: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1256 ierr = handler->PreStep(handler, &ctx, NULL, NULL); CHKERRQ(ierr);
1257 }
1258
1259 // No global communication needed for walls
1260
1261 // Phase 2: Apply - Set boundary conditions
1262 for (int i = 0; i < 6; i++) {
1263 BoundaryCondition *handler = user->boundary_faces[i].handler;
1264 if (!handler || handler->priority != BC_PRIORITY_WALL) continue;
1265 if(!handler->Apply) continue; // For example Periodic BCs
1266
1267 num_handlers[1]++;
1268
1269 BCContext ctx = {
1270 .user = user,
1271 .face_id = (BCFace)i,
1272 .global_inflow_sum = &user->simCtx->FluxInSum,
1273 .global_outflow_sum = NULL,
1274 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1275 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1276 };
1277
1278 LOG_ALLOW(LOCAL, LOG_TRACE, " Apply: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1279 ierr = handler->Apply(handler, &ctx); CHKERRQ(ierr);
1280 }
1281
1282 // Phase 3: PostStep - Post-application processing (usually no-op for walls)
1283 for (int i = 0; i < 6; i++) {
1284 BoundaryCondition *handler = user->boundary_faces[i].handler;
1285 if (!handler || handler->priority != BC_PRIORITY_WALL) continue;
1286 if (!handler->PostStep) continue;
1287
1288 num_handlers[2]++;
1289
1290 BCContext ctx = {
1291 .user = user,
1292 .face_id = (BCFace)i,
1293 .global_inflow_sum = &user->simCtx->FluxInSum,
1294 .global_outflow_sum = NULL,
1295 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1296 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1297 };
1298
1299 LOG_ALLOW(LOCAL, LOG_TRACE, " PostStep: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1300 ierr = handler->PostStep(handler, &ctx, NULL, NULL); CHKERRQ(ierr);
1301 }
1302
1303 // No global communication needed for walls
1304
1305 LOG_ALLOW(GLOBAL, LOG_INFO, " (WALLS): %d Prestep(s), %d Application(s), %d Poststep(s) applied.\n",
1306 num_handlers[0],num_handlers[1],num_handlers[2]);
1307
1308
1309 // =========================================================================
1310 // PRIORITY 3: OUTLETS
1311 // =========================================================================
1312
1313 PetscReal local_outflow_pre = 0.0;
1314 PetscReal local_outflow_post = 0.0;
1315 PetscReal global_outflow_pre = 0.0;
1316 PetscReal global_outflow_post = 0.0;
1317 memset(num_handlers,0,sizeof(num_handlers));
1318
1319 LOG_ALLOW(LOCAL, LOG_TRACE, " (OUTLETS): Begin.\n");
1320
1321 // Phase 1: PreStep - Measure uncorrected outflow (from ucat)
1322 for (int i = 0; i < 6; i++) {
1323 BoundaryCondition *handler = user->boundary_faces[i].handler;
1324 if (!handler || handler->priority != BC_PRIORITY_OUTLET) continue;
1325 if (!handler->PreStep) continue;
1326
1327 num_handlers[0]++;
1328 BCContext ctx = {
1329 .user = user,
1330 .face_id = (BCFace)i,
1331 .global_inflow_sum = &user->simCtx->FluxInSum, // From Priority 0
1332 .global_outflow_sum = NULL,
1333 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1334 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1335 };
1336
1337 LOG_ALLOW(LOCAL, LOG_TRACE, " PreStep: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1338 ierr = handler->PreStep(handler, &ctx, NULL, &local_outflow_pre); CHKERRQ(ierr);
1339 }
1340
1341 // Phase 2: Global communication - Get uncorrected outflow sum
1342 ierr = MPI_Allreduce(&local_outflow_pre, &global_outflow_pre, 1, MPIU_REAL,
1343 MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
1344
1345 // Calculate total inflow (inlet + farfield inflow)
1346 PetscReal total_inflow = user->simCtx->FluxInSum + user->simCtx->FarFluxInSum;
1347
1348 LOG_ALLOW(GLOBAL, LOG_DEBUG,
1349 " Uncorrected outflow: %.6e, Total inflow: %.6e (Inlet: %.6e + Farfield: %.6e)\n",
1350 global_outflow_pre, total_inflow, user->simCtx->FluxInSum,
1351 user->simCtx->FarFluxInSum);
1352
1353 // Phase 3: Apply - Set corrected boundary conditions
1354 for (int i = 0; i < 6; i++) {
1355 BoundaryCondition *handler = user->boundary_faces[i].handler;
1356 if (!handler || handler->priority != BC_PRIORITY_OUTLET) continue;
1357 if(!handler->Apply) continue; // For example Periodic BCs
1358
1359 num_handlers[1]++;
1360
1361 BCContext ctx = {
1362 .user = user,
1363 .face_id = (BCFace)i,
1364 .global_inflow_sum = &user->simCtx->FluxInSum, // From Priority 0
1365 .global_outflow_sum = &global_outflow_pre, // From PreStep above
1366 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1367 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1368 };
1369
1370 LOG_ALLOW(LOCAL, LOG_TRACE, " Apply: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1371 ierr = handler->Apply(handler, &ctx); CHKERRQ(ierr);
1372 }
1373
1374 // Phase 4: PostStep - Measure corrected outflow (verification)
1375 for (int i = 0; i < 6; i++) {
1376 BoundaryCondition *handler = user->boundary_faces[i].handler;
1377 if (!handler || handler->priority != BC_PRIORITY_OUTLET) continue;
1378 if (!handler->PostStep) continue;
1379
1380 num_handlers[2]++;
1381
1382 BCContext ctx = {
1383 .user = user,
1384 .face_id = (BCFace)i,
1385 .global_inflow_sum = &user->simCtx->FluxInSum,
1386 .global_outflow_sum = &global_outflow_pre,
1387 .global_farfield_inflow_sum = &user->simCtx->FarFluxInSum,
1388 .global_farfield_outflow_sum = &user->simCtx->FarFluxOutSum
1389 };
1390
1391 LOG_ALLOW(LOCAL, LOG_TRACE, " PostStep: Face %d (%s)\n", i, BCFaceToString((BCFace)i));
1392 ierr = handler->PostStep(handler, &ctx, NULL, &local_outflow_post); CHKERRQ(ierr);
1393 }
1394
1395 // Phase 5: Global communication - Verify conservation
1396 ierr = MPI_Allreduce(&local_outflow_post, &global_outflow_post, 1, MPIU_REAL,
1397 MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
1398
1399 // Store for global reporting.
1400 user->simCtx->FluxOutSum = global_outflow_post;
1401
1402 // Conservation check (compare total outflow vs total inflow)
1403 PetscReal total_outflow = global_outflow_post + user->simCtx->FarFluxOutSum;
1404 PetscReal flux_error = PetscAbsReal(total_outflow - total_inflow);
1405 PetscReal relative_error = (total_inflow > 1e-16) ?
1406 flux_error / total_inflow : flux_error;
1407
1408 LOG_ALLOW(GLOBAL, LOG_INFO,
1409 " (OUTLETS): %d Prestep(s), %d Application(s), %d Poststep(s), FluxOutSum = %.6e\n",
1410 num_handlers[0],num_handlers[1],num_handlers[2], global_outflow_post);
1411 LOG_ALLOW(GLOBAL, LOG_INFO,
1412 " Conservation: Total In=%.6e, Total Out=%.6e, Error=%.3e (%.2e)%%)\n",
1413 total_inflow, total_outflow, flux_error, relative_error * 100.0);
1414
1415 if (relative_error > 1e-6) {
1416 LOG_ALLOW(GLOBAL, LOG_WARNING,
1417 " WARNING: Large mass conservation error (%.2e%%)!\n",
1418 relative_error * 100.0);
1419 }
1420
1421
1422 LOG_ALLOW(LOCAL, LOG_VERBOSE, "Complete.\n");
1423
1424 PROFILE_FUNCTION_END;
1425 PetscFunctionReturn(0);
1426}
1427
1428// =============================================================================
1429//
1430// PRIVATE "LIGHT" EXECUTION ENGINE
1431//
1432// =============================================================================
1433
1434#undef __FUNCT__
1435#define __FUNCT__ "BoundarySystem_RefreshUbcs"
1436/**
1437 * @brief Internal helper implementation: `BoundarySystem_RefreshUbcs()`.
1438 * @details Local to this translation unit.
1439 */
1440PetscErrorCode BoundarySystem_RefreshUbcs(UserCtx *user)
1441{
1442 PetscErrorCode ierr;
1443 PetscFunctionBeginUser;
1444
1445 LOG_ALLOW(GLOBAL, LOG_TRACE, "Refreshing `ubcs` targets for flow-dependent boundaries...\n");
1446
1447 // Loop through all 6 faces of the domain
1448 for (int i = 0; i < 6; i++) {
1449 BoundaryCondition *handler = user->boundary_faces[i].handler;
1450
1451 // THE FILTER:
1452 // This is the core logic. We only act if a handler exists for the face
1453 // AND that handler has explicitly implemented the `UpdateUbcs` method.
1454 if (handler && handler->UpdateUbcs) {
1455
1456 const char *face_name = BCFaceToString((BCFace)i);
1457 LOG_ALLOW(LOCAL, LOG_TRACE, " Calling UpdateUbcs() for handler on Face %s.\n", face_name);
1458
1459 // Prepare the context. For this refresh step, we don't need to pass flux sums.
1460 BCContext ctx = {
1461 .user = user,
1462 .face_id = (BCFace)i,
1463 .global_inflow_sum = NULL,
1464 .global_outflow_sum = NULL,
1465 .global_farfield_inflow_sum = NULL,
1466 .global_farfield_outflow_sum = NULL
1467 };
1468
1469 // Call the handler's specific UpdateUbcs function pointer.
1470 ierr = handler->UpdateUbcs(handler, &ctx); CHKERRQ(ierr);
1471 }
1472 }
1473
1474 PetscFunctionReturn(0);
1475}
1476
1477//================================================================================
1478//
1479// PUBLIC MASTER CLEANUP FUNCTION
1480//
1481//================================================================================
1482#undef __FUNCT__
1483#define __FUNCT__ "BoundarySystem_Destroy"
1484/**
1485 * @brief Implementation of \ref BoundarySystem_Destroy().
1486 * @details Full API contract (arguments, ownership, side effects) is documented with
1487 * the header declaration in `include/Boundaries.h`.
1488 * @see BoundarySystem_Destroy()
1489 */
1490PetscErrorCode BoundarySystem_Destroy(UserCtx *user)
1491{
1492 PetscErrorCode ierr;
1493 PetscFunctionBeginUser;
1494
1495
1496
1497 LOG_ALLOW(GLOBAL, LOG_INFO, "Starting destruction of all boundary handlers. \n");
1498
1499 for (int i = 0; i < 6; i++) {
1500 BoundaryFaceConfig *face_cfg = &user->boundary_faces[i];
1501 const char *face_name = BCFaceToString(face_cfg->face_id);
1502
1503 // --- Step 1: Free the parameter linked list associated with this face ---
1504 if (face_cfg->params) {
1505 LOG_ALLOW(LOCAL, LOG_DEBUG, " Freeing parameter list for Face %d (%s). \n", i, face_name);
1506 FreeBC_ParamList(face_cfg->params);
1507 face_cfg->params = NULL; // Good practice to nullify dangling pointers
1508 }
1509
1510 // --- Step 2: Destroy the handler object itself ---
1511 if (face_cfg->handler) {
1512 const char *handler_name = BCHandlerTypeToString(face_cfg->handler->type);
1513 LOG_ALLOW(LOCAL, LOG_DEBUG, " Destroying handler '%s' on Face %d (%s).\n", handler_name, i, face_name);
1514
1515 // Call the handler's specific cleanup function first, if it exists.
1516 // This will free any memory stored in the handler's private `data` pointer.
1517 if (face_cfg->handler->Destroy) {
1518 ierr = face_cfg->handler->Destroy(face_cfg->handler); CHKERRQ(ierr);
1519 }
1520
1521 // Finally, free the generic BoundaryCondition object itself.
1522 ierr = PetscFree(face_cfg->handler); CHKERRQ(ierr);
1523 face_cfg->handler = NULL;
1524 }
1525 }
1526
1527 LOG_ALLOW(GLOBAL, LOG_INFO, "Destruction complete.\n");
1528 PetscFunctionReturn(0);
1529}
1530
1531#undef __FUNCT__
1532#define __FUNCT__ "TransferPeriodicFieldByDirection"
1533/**
1534 * @brief Copies one cell field's wrapped local values onto the owned periodic duplicate plane.
1535 * @details Handles scalar and three-component storage for one selected logical axis.
1536 */
1537static PetscErrorCode TransferPeriodicFieldByDirection(UserCtx *user, FieldId field_id, char direction)
1538{
1539 PetscErrorCode ierr;
1540 DMDALocalInfo info = user->info;
1541 PetscInt xs = info.xs, xe = info.xs + info.xm;
1542 PetscInt ys = info.ys, ye = info.ys + info.ym;
1543 PetscInt zs = info.zs, ze = info.zs + info.zm;
1544 PetscInt mx = info.mx, my = info.my, mz = info.mz;
1545
1546 FieldView field_view;
1547 DM dm;
1548 Vec global_vec;
1549 Vec local_vec;
1550 PetscInt dof;
1551
1552 PetscFunctionBeginUser;
1553 PetscCall(FieldGetView(user, field_id, &field_view));
1554 PetscCheck((field_view.descriptor->capabilities & FIELD_CAPABILITY_PERIODIC_CELL_SYNC) != 0u,
1555 PETSC_COMM_SELF, PETSC_ERR_SUP,
1556 "Field '%s' is not registered for periodic cell synchronization.",
1557 field_view.descriptor->canonical_name);
1558 dm = field_view.dm;
1559 global_vec = field_view.global_vec;
1560 local_vec = field_view.local_vec;
1561 dof = field_view.descriptor->dof;
1562
1563 // --- Execute the copy logic based on DoF and Direction ---
1564 if (dof == 1) { // --- Handle SCALAR fields (PetscReal) ---
1565 PetscReal ***g_array, ***l_array;
1566 ierr = DMDAVecGetArray(dm, global_vec, &g_array); CHKERRQ(ierr);
1567 ierr = DMDAVecGetArrayRead(dm, local_vec, (void*)&l_array); CHKERRQ(ierr); // Use Read for safety
1568
1569 switch (direction) {
1570 case 'i':
1571 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];
1572 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];
1573 break;
1574 case 'j':
1575 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];
1576 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];
1577 break;
1578 case 'k':
1579 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];
1580 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];
1581 break;
1582 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid direction '%c'", direction);
1583 }
1584 ierr = DMDAVecRestoreArray(dm, global_vec, &g_array); CHKERRQ(ierr);
1585 ierr = DMDAVecRestoreArrayRead(dm, local_vec, (void*)&l_array); CHKERRQ(ierr);
1586
1587 } else if (dof == 3) { // --- Handle VECTOR fields (Cmpnts) ---
1588 Cmpnts ***g_array, ***l_array;
1589 ierr = DMDAVecGetArray(dm, global_vec, &g_array); CHKERRQ(ierr);
1590 ierr = DMDAVecGetArrayRead(dm, local_vec, (void*)&l_array); CHKERRQ(ierr);
1591
1592 switch (direction) {
1593 case 'i':
1594 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];
1595 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];
1596 break;
1597 case 'j':
1598 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];
1599 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];
1600 break;
1601 case 'k':
1602 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];
1603 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];
1604 break;
1605 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid direction '%c'", direction);
1606 }
1607 ierr = DMDAVecRestoreArray(dm, global_vec, &g_array); CHKERRQ(ierr);
1608 ierr = DMDAVecRestoreArrayRead(dm, local_vec, (void*)&l_array); CHKERRQ(ierr);
1609 }
1610
1611 PetscFunctionReturn(0);
1612}
1613
1614#undef __FUNCT__
1615#define __FUNCT__ "SynchronizePeriodicCellFields"
1616/**
1617 * @brief Implementation of \ref SynchronizePeriodicCellFields().
1618 * @details Full API contract is documented with the header declaration in
1619 * `include/Boundaries.h`.
1620 */
1621PetscErrorCode SynchronizePeriodicCellFields(UserCtx *user, PetscInt num_fields, const FieldId field_ids[])
1622{
1623 PetscErrorCode ierr;
1624 PetscBool periodic_i;
1625 PetscBool periodic_j;
1626 PetscBool periodic_k;
1627
1628 PetscFunctionBeginUser;
1629
1630 PetscCheck(num_fields >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE,
1631 "Number of cell fields cannot be negative.");
1632 if (num_fields == 0) PetscFunctionReturn(0);
1633 PetscCheck(field_ids != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
1634 "Cell field-ID array cannot be NULL.");
1635
1636 periodic_i =
1637 user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC ||
1638 user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC;
1639 periodic_j =
1640 user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC ||
1641 user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC;
1642 periodic_k =
1643 user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC ||
1644 user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC;
1645
1646 if (!periodic_i && !periodic_j && !periodic_k) PetscFunctionReturn(0);
1647
1648 for (PetscInt field = 0; field < num_fields; field++) {
1649 ierr = UpdateLocalGhosts(user, field_ids[field]); CHKERRQ(ierr);
1650 }
1651
1652 if (periodic_i) {
1653 for (PetscInt field = 0; field < num_fields; field++) {
1654 ierr = TransferPeriodicFieldByDirection(user, field_ids[field], 'i'); CHKERRQ(ierr);
1655 }
1656 for (PetscInt field = 0; field < num_fields; field++) {
1657 ierr = UpdateLocalGhosts(user, field_ids[field]); CHKERRQ(ierr);
1658 }
1659 }
1660
1661 if (periodic_j) {
1662 for (PetscInt field = 0; field < num_fields; field++) {
1663 ierr = TransferPeriodicFieldByDirection(user, field_ids[field], 'j'); CHKERRQ(ierr);
1664 }
1665 for (PetscInt field = 0; field < num_fields; field++) {
1666 ierr = UpdateLocalGhosts(user, field_ids[field]); CHKERRQ(ierr);
1667 }
1668 }
1669
1670 if (periodic_k) {
1671 for (PetscInt field = 0; field < num_fields; field++) {
1672 ierr = TransferPeriodicFieldByDirection(user, field_ids[field], 'k'); CHKERRQ(ierr);
1673 }
1674 for (PetscInt field = 0; field < num_fields; field++) {
1675 ierr = UpdateLocalGhosts(user, field_ids[field]); CHKERRQ(ierr);
1676 }
1677 }
1678
1679 PetscFunctionReturn(0);
1680}
1681
1682#undef __FUNCT__
1683#define __FUNCT__ "GetPersistentFaceField"
1684/**
1685 * @brief Resolves one registered persistent single-face-family field.
1686 */
1687static PetscErrorCode GetPersistentFaceField(UserCtx *user, FieldId field_id,
1688 char face_direction, DM *dm,
1689 Vec *global_vec, Vec *local_vec,
1690 PetscInt *dof)
1691{
1692 FieldView field_view;
1693 FieldLayout expected_layout;
1694
1695 PetscFunctionBeginUser;
1696 PetscCheck(face_direction == 'i' || face_direction == 'j' || face_direction == 'k',
1697 PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
1698 "Invalid face direction '%c'; expected 'i', 'j', or 'k'.", face_direction);
1699
1700 expected_layout = face_direction == 'i' ? FIELD_LAYOUT_I_FACE :
1701 (face_direction == 'j' ? FIELD_LAYOUT_J_FACE : FIELD_LAYOUT_K_FACE);
1702 PetscCall(FieldGetView(user, field_id, &field_view));
1703 PetscCheck((field_view.descriptor->capabilities & FIELD_CAPABILITY_PERIODIC_FACE_SYNC) != 0u &&
1704 field_view.descriptor->layout == expected_layout,
1705 PETSC_COMM_SELF, PETSC_ERR_SUP,
1706 "Field '%s' is not registered for %c-face periodic synchronization.",
1707 field_view.descriptor->canonical_name, face_direction);
1708
1709 *dm = field_view.dm;
1710 *global_vec = field_view.global_vec;
1711 *local_vec = field_view.local_vec;
1712 *dof = field_view.descriptor->dof;
1713 PetscFunctionReturn(0);
1714}
1715
1716/**
1717 * @brief Returns whether a registered face field stores physical coordinates.
1718 */
1719static PetscErrorCode IsFaceCenterCoordinateField(FieldId field_id, PetscBool *is_coordinate)
1720{
1721 const FieldDescriptor *descriptor = NULL;
1722
1723 PetscFunctionBeginUser;
1724 PetscCheck(is_coordinate != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
1725 "Face-coordinate result cannot be NULL.");
1726 PetscCall(FieldGetDescriptor(field_id, &descriptor));
1727 *is_coordinate = (PetscBool)((descriptor->capabilities & FIELD_CAPABILITY_PERIODIC_GEOMETRY_SHIFT) != 0u);
1728 PetscFunctionReturn(0);
1729}
1730
1731/**
1732 * @brief Applies geometric translations to wrapped face-center ghost coordinates.
1733 */
1734static PetscErrorCode TranslatePeriodicFaceCenterGhosts(UserCtx *user, Vec local_vec)
1735{
1736 DMDALocalInfo info;
1737 Cmpnts ***array;
1738 const BCFace negative_faces[3] = {BC_FACE_NEG_X, BC_FACE_NEG_Y, BC_FACE_NEG_Z};
1739 const BCFace positive_faces[3] = {BC_FACE_POS_X, BC_FACE_POS_Y, BC_FACE_POS_Z};
1740
1741 PetscFunctionBeginUser;
1742 PetscCall(DMDAGetLocalInfo(user->fda, &info));
1743 PetscCall(DMDAVecGetArray(user->fda, local_vec, &array));
1744
1745 for (PetscInt axis = 0; axis < 3; axis++) {
1746 const PetscBool active =
1747 user->boundary_faces[negative_faces[axis]].mathematical_type == PERIODIC ||
1748 user->boundary_faces[positive_faces[axis]].mathematical_type == PERIODIC;
1749 Cmpnts translation;
1750
1751 if (!active) continue;
1752 PetscCheck(user->periodic_translation_valid[axis], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE,
1753 "Periodic face-center synchronization requires validated %c-direction geometry.",
1754 "XYZ"[axis]);
1755 translation = user->periodic_translation[axis];
1756
1757 for (PetscInt k = info.gzs; k < info.gzs + info.gzm; k++) {
1758 for (PetscInt j = info.gys; j < info.gys + info.gym; j++) {
1759 for (PetscInt i = info.gxs; i < info.gxs + info.gxm; i++) {
1760 PetscReal scale = 0.0;
1761 const PetscInt index = axis == 0 ? i : (axis == 1 ? j : k);
1762 const PetscInt size = axis == 0 ? info.mx : (axis == 1 ? info.my : info.mz);
1763 if (index < 0) scale = -1.0;
1764 else if (index >= size) scale = 1.0;
1765 if (scale == 0.0) continue;
1766 array[k][j][i].x += scale * translation.x;
1767 array[k][j][i].y += scale * translation.y;
1768 array[k][j][i].z += scale * translation.z;
1769 }
1770 }
1771 }
1772 }
1773
1774 PetscCall(DMDAVecRestoreArray(user->fda, local_vec, &array));
1775 PetscFunctionReturn(0);
1776}
1777
1778#undef __FUNCT__
1779#define __FUNCT__ "TransferPeriodicFaceFieldByDirection"
1780/** @brief Transfers one registered face-family field along one periodic axis. */
1781static PetscErrorCode TransferPeriodicFaceFieldByDirection(UserCtx *user, FieldId field_id,
1782 char face_direction, char periodic_direction)
1783{
1784 PetscErrorCode ierr;
1785 DMDALocalInfo info = user->info;
1786 PetscInt xs = info.xs, xe = info.xs + info.xm;
1787 PetscInt ys = info.ys, ye = info.ys + info.ym;
1788 PetscInt zs = info.zs, ze = info.zs + info.zm;
1789 PetscInt mx = info.mx, my = info.my, mz = info.mz;
1790 DM dm;
1791 Vec global_vec, local_vec;
1792 PetscInt dof;
1793
1794 PetscFunctionBeginUser;
1795 PetscCheck(periodic_direction == 'i' || periodic_direction == 'j' || periodic_direction == 'k',
1796 PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
1797 "Invalid periodic direction '%c'; expected 'i', 'j', or 'k'.", periodic_direction);
1798 PetscCall(GetPersistentFaceField(user, field_id, face_direction, &dm, &global_vec, &local_vec, &dof));
1799
1800 if (dof == 1) {
1801 PetscReal ***global_array, ***local_array;
1802 ierr = DMDAVecGetArray(dm, global_vec, &global_array); CHKERRQ(ierr);
1803 ierr = DMDAVecGetArrayRead(dm, local_vec, &local_array); CHKERRQ(ierr);
1804
1805 if (periodic_direction == 'i') {
1806 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && xs == 0)
1807 for (PetscInt k=zs; k<ze; k++) for (PetscInt j=ys; j<ye; j++) global_array[k][j][0] = local_array[k][j][-2];
1808 if (user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC && xe == mx)
1809 for (PetscInt k=zs; k<ze; k++) for (PetscInt j=ys; j<ye; j++) global_array[k][j][mx-1] = local_array[k][j][mx+1];
1810 } else if (periodic_direction == 'j') {
1811 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && ys == 0)
1812 for (PetscInt k=zs; k<ze; k++) for (PetscInt i=xs; i<xe; i++) global_array[k][0][i] = local_array[k][-2][i];
1813 if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC && ye == my)
1814 for (PetscInt k=zs; k<ze; k++) for (PetscInt i=xs; i<xe; i++) global_array[k][my-1][i] = local_array[k][my+1][i];
1815 } else {
1816 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && zs == 0)
1817 for (PetscInt j=ys; j<ye; j++) for (PetscInt i=xs; i<xe; i++) global_array[0][j][i] = local_array[-2][j][i];
1818 if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC && ze == mz)
1819 for (PetscInt j=ys; j<ye; j++) for (PetscInt i=xs; i<xe; i++) global_array[mz-1][j][i] = local_array[mz+1][j][i];
1820 }
1821 ierr = DMDAVecRestoreArrayRead(dm, local_vec, &local_array); CHKERRQ(ierr);
1822 ierr = DMDAVecRestoreArray(dm, global_vec, &global_array); CHKERRQ(ierr);
1823 } else {
1824 Cmpnts ***global_array, ***local_array;
1825 ierr = DMDAVecGetArray(dm, global_vec, &global_array); CHKERRQ(ierr);
1826 ierr = DMDAVecGetArrayRead(dm, local_vec, &local_array); CHKERRQ(ierr);
1827
1828 if (periodic_direction == 'i') {
1829 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && xs == 0)
1830 for (PetscInt k=zs; k<ze; k++) for (PetscInt j=ys; j<ye; j++) global_array[k][j][0] = local_array[k][j][-2];
1831 if (user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC && xe == mx)
1832 for (PetscInt k=zs; k<ze; k++) for (PetscInt j=ys; j<ye; j++) global_array[k][j][mx-1] = local_array[k][j][mx+1];
1833 } else if (periodic_direction == 'j') {
1834 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && ys == 0)
1835 for (PetscInt k=zs; k<ze; k++) for (PetscInt i=xs; i<xe; i++) global_array[k][0][i] = local_array[k][-2][i];
1836 if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC && ye == my)
1837 for (PetscInt k=zs; k<ze; k++) for (PetscInt i=xs; i<xe; i++) global_array[k][my-1][i] = local_array[k][my+1][i];
1838 } else {
1839 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && zs == 0)
1840 for (PetscInt j=ys; j<ye; j++) for (PetscInt i=xs; i<xe; i++) global_array[0][j][i] = local_array[-2][j][i];
1841 if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC && ze == mz)
1842 for (PetscInt j=ys; j<ye; j++) for (PetscInt i=xs; i<xe; i++) global_array[mz-1][j][i] = local_array[mz+1][j][i];
1843 }
1844 ierr = DMDAVecRestoreArrayRead(dm, local_vec, &local_array); CHKERRQ(ierr);
1845 ierr = DMDAVecRestoreArray(dm, global_vec, &global_array); CHKERRQ(ierr);
1846 }
1847
1848 PetscFunctionReturn(0);
1849}
1850
1851#undef __FUNCT__
1852#define __FUNCT__ "SynchronizePeriodicFaceFields"
1853// Implements SynchronizePeriodicFaceFields(); the public header owns the
1854// rendered API contract.
1855PetscErrorCode SynchronizePeriodicFaceFields(UserCtx *user, char face_direction,
1856 PetscInt num_fields, const FieldId field_ids[])
1857{
1858 PetscErrorCode ierr;
1859 const char periodic_directions[3] = {'i', 'j', 'k'};
1860 const BCFace negative_faces[3] = {BC_FACE_NEG_X, BC_FACE_NEG_Y, BC_FACE_NEG_Z};
1861 const BCFace positive_faces[3] = {BC_FACE_POS_X, BC_FACE_POS_Y, BC_FACE_POS_Z};
1862
1863 PetscFunctionBeginUser;
1864 PetscCheck(num_fields >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE,
1865 "Number of face fields cannot be negative.");
1866 if (num_fields == 0) PetscFunctionReturn(0);
1867 PetscCheck(field_ids != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
1868 "Face field-ID array cannot be NULL.");
1869
1870 for (PetscInt field = 0; field < num_fields; field++) {
1871 DM dm;
1872 Vec global_vec, local_vec;
1873 PetscInt dof;
1874 PetscBool is_coordinate;
1875 PetscCall(GetPersistentFaceField(user, field_ids[field], face_direction,
1876 &dm, &global_vec, &local_vec, &dof));
1877 ierr = UpdateLocalGhosts(user, field_ids[field]); CHKERRQ(ierr);
1878 ierr = IsFaceCenterCoordinateField(field_ids[field], &is_coordinate); CHKERRQ(ierr);
1879 if (is_coordinate) {
1880 ierr = TranslatePeriodicFaceCenterGhosts(user, local_vec); CHKERRQ(ierr);
1881 }
1882 }
1883
1884 for (PetscInt direction = 0; direction < 3; direction++) {
1885 const PetscBool active =
1886 user->boundary_faces[negative_faces[direction]].mathematical_type == PERIODIC ||
1887 user->boundary_faces[positive_faces[direction]].mathematical_type == PERIODIC;
1888 if (!active) continue;
1889
1890 for (PetscInt field = 0; field < num_fields; field++) {
1891 ierr = TransferPeriodicFaceFieldByDirection(user, field_ids[field], face_direction,
1892 periodic_directions[direction]); CHKERRQ(ierr);
1893 }
1894 for (PetscInt field = 0; field < num_fields; field++) {
1895 DM dm;
1896 Vec global_vec, local_vec;
1897 PetscInt dof;
1898 PetscBool is_coordinate;
1899 ierr = UpdateLocalGhosts(user, field_ids[field]); CHKERRQ(ierr);
1900 ierr = IsFaceCenterCoordinateField(field_ids[field], &is_coordinate); CHKERRQ(ierr);
1901 if (is_coordinate) {
1902 PetscCall(GetPersistentFaceField(user, field_ids[field], face_direction,
1903 &dm, &global_vec, &local_vec, &dof));
1904 ierr = TranslatePeriodicFaceCenterGhosts(user, local_vec); CHKERRQ(ierr);
1905 }
1906 }
1907 }
1908
1909 PetscFunctionReturn(0);
1910}
1911
1912#undef __FUNCT__
1913#define __FUNCT__ "GetPersistentStaggeredField"
1914/**
1915 * @brief Resolves one registered persistent component-staggered field.
1916 */
1917static PetscErrorCode GetPersistentStaggeredField(UserCtx *user, FieldId field_id,
1918 DM *dm, Vec *global_vec, Vec *local_vec)
1919{
1920 FieldView field_view;
1921
1922 PetscFunctionBeginUser;
1923 PetscCall(FieldGetView(user, field_id, &field_view));
1924 PetscCheck((field_view.descriptor->capabilities & FIELD_CAPABILITY_PERIODIC_STAGGERED_SYNC) != 0u &&
1925 field_view.descriptor->layout == FIELD_LAYOUT_COMPONENT_STAGGERED,
1926 PETSC_COMM_SELF, PETSC_ERR_SUP,
1927 "Field '%s' is not registered for component-staggered periodic synchronization.",
1928 field_view.descriptor->canonical_name);
1929 *dm = field_view.dm;
1930 *global_vec = field_view.global_vec;
1931 *local_vec = field_view.local_vec;
1932 PetscFunctionReturn(0);
1933}
1934
1935#undef __FUNCT__
1936#define __FUNCT__ "TransferPeriodicStaggeredFieldByDirection"
1937/**
1938 * @brief Transfers one component-staggered field along one periodic axis.
1939 */
1940static PetscErrorCode TransferPeriodicStaggeredFieldByDirection(UserCtx *user, FieldId field_id,
1941 char periodic_direction)
1942{
1943 PetscErrorCode ierr;
1944 DMDALocalInfo info = user->info;
1945 PetscInt xs = info.xs, xe = info.xs + info.xm;
1946 PetscInt ys = info.ys, ye = info.ys + info.ym;
1947 PetscInt zs = info.zs, ze = info.zs + info.zm;
1948 PetscInt mx = info.mx, my = info.my, mz = info.mz;
1949 DM dm;
1950 Vec global_vec, local_vec;
1951 Cmpnts ***global_array, ***local_array;
1952
1953 PetscFunctionBeginUser;
1954 PetscCheck(periodic_direction == 'i' || periodic_direction == 'j' || periodic_direction == 'k',
1955 PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
1956 "Invalid periodic direction '%c'; expected 'i', 'j', or 'k'.", periodic_direction);
1957 PetscCall(GetPersistentStaggeredField(user, field_id, &dm, &global_vec, &local_vec));
1958
1959 ierr = DMDAVecGetArray(dm, global_vec, &global_array); CHKERRQ(ierr);
1960 ierr = DMDAVecGetArrayRead(dm, local_vec, &local_array); CHKERRQ(ierr);
1961
1962 if (periodic_direction == 'i') {
1963 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && xs == 0)
1964 for (PetscInt k=zs; k<ze; k++) for (PetscInt j=ys; j<ye; j++) global_array[k][j][0] = local_array[k][j][-2];
1965 if (user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC && xe == mx)
1966 for (PetscInt k=zs; k<ze; k++) for (PetscInt j=ys; j<ye; j++) global_array[k][j][mx-1] = local_array[k][j][mx+1];
1967 } else if (periodic_direction == 'j') {
1968 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && ys == 0)
1969 for (PetscInt k=zs; k<ze; k++) for (PetscInt i=xs; i<xe; i++) global_array[k][0][i] = local_array[k][-2][i];
1970 if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC && ye == my)
1971 for (PetscInt k=zs; k<ze; k++) for (PetscInt i=xs; i<xe; i++) global_array[k][my-1][i] = local_array[k][my+1][i];
1972 } else {
1973 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && zs == 0)
1974 for (PetscInt j=ys; j<ye; j++) for (PetscInt i=xs; i<xe; i++) global_array[0][j][i] = local_array[-2][j][i];
1975 if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC && ze == mz)
1976 for (PetscInt j=ys; j<ye; j++) for (PetscInt i=xs; i<xe; i++) global_array[mz-1][j][i] = local_array[mz+1][j][i];
1977 }
1978
1979 ierr = DMDAVecRestoreArrayRead(dm, local_vec, &local_array); CHKERRQ(ierr);
1980 ierr = DMDAVecRestoreArray(dm, global_vec, &global_array); CHKERRQ(ierr);
1981 PetscFunctionReturn(0);
1982}
1983
1984#undef __FUNCT__
1985#define __FUNCT__ "SynchronizePeriodicStaggeredFields"
1986/**
1987 * @brief Implementation of \ref SynchronizePeriodicStaggeredFields().
1988 */
1989PetscErrorCode SynchronizePeriodicStaggeredFields(UserCtx *user, PetscInt num_fields,
1990 const FieldId field_ids[])
1991{
1992 PetscErrorCode ierr;
1993 const char periodic_directions[3] = {'i', 'j', 'k'};
1994 const BCFace negative_faces[3] = {BC_FACE_NEG_X, BC_FACE_NEG_Y, BC_FACE_NEG_Z};
1995 const BCFace positive_faces[3] = {BC_FACE_POS_X, BC_FACE_POS_Y, BC_FACE_POS_Z};
1996
1997 PetscFunctionBeginUser;
1998 PetscCheck(num_fields >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE,
1999 "Number of staggered fields cannot be negative.");
2000 if (num_fields == 0) PetscFunctionReturn(0);
2001 PetscCheck(field_ids != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
2002 "Staggered field-ID array cannot be NULL.");
2003
2004 for (PetscInt field = 0; field < num_fields; field++) {
2005 DM dm;
2006 Vec global_vec, local_vec;
2007 PetscCall(GetPersistentStaggeredField(user, field_ids[field], &dm, &global_vec, &local_vec));
2008 ierr = UpdateLocalGhosts(user, field_ids[field]); CHKERRQ(ierr);
2009 }
2010
2011 for (PetscInt direction = 0; direction < 3; direction++) {
2012 const PetscBool active =
2013 user->boundary_faces[negative_faces[direction]].mathematical_type == PERIODIC ||
2014 user->boundary_faces[positive_faces[direction]].mathematical_type == PERIODIC;
2015 if (!active) continue;
2016
2017 for (PetscInt field = 0; field < num_fields; field++) {
2018 ierr = TransferPeriodicStaggeredFieldByDirection(user, field_ids[field],
2019 periodic_directions[direction]); CHKERRQ(ierr);
2020 }
2021 for (PetscInt field = 0; field < num_fields; field++) {
2022 ierr = UpdateLocalGhosts(user, field_ids[field]); CHKERRQ(ierr);
2023 }
2024 }
2025
2026 PetscFunctionReturn(0);
2027}
2028
2029#undef __FUNCT__
2030#define __FUNCT__ "PreparePeriodicQuickStencilFields"
2031/**
2032 * @brief Implementation of \ref PreparePeriodicQuickStencilFields().
2033 */
2034PetscErrorCode PreparePeriodicQuickStencilFields(UserCtx *user, Vec local_vector_field,
2035 Vec local_scalar_field)
2036{
2037 DMDALocalInfo info = user->info;
2038 Cmpnts ***vector_array;
2039 PetscReal ***scalar_array;
2040 const PetscInt xs = info.xs, xe = info.xs + info.xm;
2041 const PetscInt ys = info.ys, ye = info.ys + info.ym;
2042 const PetscInt zs = info.zs, ze = info.zs + info.zm;
2043 const PetscInt gxs = info.gxs, gxe = info.gxs + info.gxm;
2044 const PetscInt gys = info.gys, gye = info.gys + info.gym;
2045 const PetscInt gzs = info.gzs, gze = info.gzs + info.gzm;
2046
2047 PetscFunctionBeginUser;
2048 PetscCheck(local_vector_field && local_scalar_field, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
2049 "QUICK stencil repair requires both local vector and scalar fields.");
2050 PetscCall(DMDAVecGetArray(user->fda, local_vector_field, &vector_array));
2051 PetscCall(DMDAVecGetArray(user->da, local_scalar_field, &scalar_array));
2052
2053 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && xs == 0) {
2054 for (PetscInt k = gzs; k < gze; k++) for (PetscInt j = gys; j < gye; j++) {
2055 vector_array[k][j][-1] = vector_array[k][j][-3];
2056 scalar_array[k][j][-1] = scalar_array[k][j][-3];
2057 }
2058 }
2059 if (user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC && xe == info.mx) {
2060 for (PetscInt k = gzs; k < gze; k++) for (PetscInt j = gys; j < gye; j++) {
2061 vector_array[k][j][info.mx] = vector_array[k][j][info.mx + 2];
2062 scalar_array[k][j][info.mx] = scalar_array[k][j][info.mx + 2];
2063 }
2064 }
2065 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && ys == 0) {
2066 for (PetscInt k = gzs; k < gze; k++) for (PetscInt i = gxs; i < gxe; i++) {
2067 vector_array[k][-1][i] = vector_array[k][-3][i];
2068 scalar_array[k][-1][i] = scalar_array[k][-3][i];
2069 }
2070 }
2071 if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC && ye == info.my) {
2072 for (PetscInt k = gzs; k < gze; k++) for (PetscInt i = gxs; i < gxe; i++) {
2073 vector_array[k][info.my][i] = vector_array[k][info.my + 2][i];
2074 scalar_array[k][info.my][i] = scalar_array[k][info.my + 2][i];
2075 }
2076 }
2077 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && zs == 0) {
2078 for (PetscInt j = gys; j < gye; j++) for (PetscInt i = gxs; i < gxe; i++) {
2079 vector_array[-1][j][i] = vector_array[-3][j][i];
2080 scalar_array[-1][j][i] = scalar_array[-3][j][i];
2081 }
2082 }
2083 if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC && ze == info.mz) {
2084 for (PetscInt j = gys; j < gye; j++) for (PetscInt i = gxs; i < gxe; i++) {
2085 vector_array[info.mz][j][i] = vector_array[info.mz + 2][j][i];
2086 scalar_array[info.mz][j][i] = scalar_array[info.mz + 2][j][i];
2087 }
2088 }
2089
2090 PetscCall(DMDAVecRestoreArray(user->da, local_scalar_field, &scalar_array));
2091 PetscCall(DMDAVecRestoreArray(user->fda, local_vector_field, &vector_array));
2092 PetscFunctionReturn(0);
2093}
2094
2095#undef __FUNCT__
2096#define __FUNCT__ "SynchronizePeriodicLocalStaggeredField"
2097/**
2098 * @brief Implementation of \ref SynchronizePeriodicLocalStaggeredField().
2099 */
2100PetscErrorCode SynchronizePeriodicLocalStaggeredField(UserCtx *user, Vec local_field)
2101{
2102 DMDALocalInfo info = user->info;
2103 Cmpnts ***array;
2104 const PetscInt xs = info.xs, xe = info.xs + info.xm;
2105 const PetscInt ys = info.ys, ye = info.ys + info.ym;
2106 const PetscInt zs = info.zs, ze = info.zs + info.zm;
2107
2108 PetscFunctionBeginUser;
2109 PetscCheck(local_field, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
2110 "Local staggered field cannot be NULL.");
2111 PetscCall(DMLocalToLocalBegin(user->fda, local_field, INSERT_VALUES, local_field));
2112 PetscCall(DMLocalToLocalEnd(user->fda, local_field, INSERT_VALUES, local_field));
2113 PetscCall(DMDAVecGetArray(user->fda, local_field, &array));
2114
2115 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && xs == 0)
2116 for (PetscInt k = zs; k < ze; k++) for (PetscInt j = ys; j < ye; j++) array[k][j][0].x = array[k][j][-2].x;
2117 if (user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC && xe == info.mx)
2118 for (PetscInt k = zs; k < ze; k++) for (PetscInt j = ys; j < ye; j++) array[k][j][info.mx - 1].x = array[k][j][info.mx + 1].x;
2119 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && ys == 0)
2120 for (PetscInt k = zs; k < ze; k++) for (PetscInt i = xs; i < xe; i++) array[k][0][i].y = array[k][-2][i].y;
2121 if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC && ye == info.my)
2122 for (PetscInt k = zs; k < ze; k++) for (PetscInt i = xs; i < xe; i++) array[k][info.my - 1][i].y = array[k][info.my + 1][i].y;
2123 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && zs == 0)
2124 for (PetscInt j = ys; j < ye; j++) for (PetscInt i = xs; i < xe; i++) array[0][j][i].z = array[-2][j][i].z;
2125 if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC && ze == info.mz)
2126 for (PetscInt j = ys; j < ye; j++) for (PetscInt i = xs; i < xe; i++) array[info.mz - 1][j][i].z = array[info.mz + 1][j][i].z;
2127
2128 PetscCall(DMDAVecRestoreArray(user->fda, local_field, &array));
2129 PetscCall(DMLocalToLocalBegin(user->fda, local_field, INSERT_VALUES, local_field));
2130 PetscCall(DMLocalToLocalEnd(user->fda, local_field, INSERT_VALUES, local_field));
2131 PetscFunctionReturn(0);
2132}
2133
2134#undef __FUNCT__
2135#define __FUNCT__ "ApplyMetricsPeriodicBCs"
2136/**
2137 * @brief Internal helper implementation: `ApplyMetricsPeriodicBCs()`.
2138 * @details Local to this translation unit.
2139 */
2140PetscErrorCode ApplyMetricsPeriodicBCs(UserCtx *user)
2141{
2142 PetscErrorCode ierr;
2143 PetscFunctionBeginUser;
2144 PROFILE_FUNCTION_BEGIN;
2145
2146 const FieldId cell_fields[] = {FIELD_ID_AJ};
2147 const FieldId i_face_fields[] = {FIELD_ID_CENTX, FIELD_ID_CSI, FIELD_ID_ICSI,
2148 FIELD_ID_IETA, FIELD_ID_IZET, FIELD_ID_IAJ};
2149 const FieldId j_face_fields[] = {FIELD_ID_CENTY, FIELD_ID_ETA, FIELD_ID_JCSI,
2150 FIELD_ID_JETA, FIELD_ID_JZET, FIELD_ID_JAJ};
2151 const FieldId k_face_fields[] = {FIELD_ID_CENTZ, FIELD_ID_ZET, FIELD_ID_KCSI,
2152 FIELD_ID_KETA, FIELD_ID_KZET, FIELD_ID_KAJ};
2153
2154 ierr = SynchronizePeriodicCellFields(user, 1, cell_fields); CHKERRQ(ierr);
2155 ierr = SynchronizePeriodicFaceFields(user, 'i', 6, i_face_fields); CHKERRQ(ierr);
2156 ierr = SynchronizePeriodicFaceFields(user, 'j', 6, j_face_fields); CHKERRQ(ierr);
2157 ierr = SynchronizePeriodicFaceFields(user, 'k', 6, k_face_fields); CHKERRQ(ierr);
2158
2159 PROFILE_FUNCTION_END;
2160 PetscFunctionReturn(0);
2161}
2162
2163#undef __FUNCT__
2164#define __FUNCT__ "ApplyPeriodicBCs"
2165/**
2166 * @brief Internal helper implementation: `ApplyPeriodicBCs()`.
2167 * @details Local to this translation unit.
2168 */
2169PetscErrorCode ApplyPeriodicBCs(UserCtx *user)
2170{
2171 PetscErrorCode ierr;
2172 PetscBool is_any_periodic = PETSC_FALSE;
2173
2174 PetscFunctionBeginUser;
2175
2176 PROFILE_FUNCTION_BEGIN;
2177
2178 for (int i = 0; i < 6; i++) {
2179 if (user->boundary_faces[i].mathematical_type == PERIODIC) {
2180 is_any_periodic = PETSC_TRUE;
2181 break;
2182 }
2183 }
2184
2185 if (!is_any_periodic) {
2186 LOG_ALLOW(GLOBAL,LOG_TRACE, "No periodic boundaries defined; skipping ApplyPeriodicBCs.\n");
2187 PROFILE_FUNCTION_END;
2188 PetscFunctionReturn(0);
2189 }
2190
2191 LOG_ALLOW(GLOBAL, LOG_TRACE, "Applying periodic boundary conditions for all fields.\n");
2192
2193 // STEP 1: Synchronize periodic cell-centered fields in deterministic direction order.
2194 const FieldId cell_fields[] = {FIELD_ID_UCAT, FIELD_ID_P, FIELD_ID_NVERT};
2195 ierr = SynchronizePeriodicCellFields(user, 3, cell_fields); CHKERRQ(ierr);
2196
2197 /* A future temperature field must be catalogued before requesting its typed ghost update. */
2198
2199 // STEP 2: Synchronize persistent staggered endpoints and repair local
2200 // component-normal ghosts through UpdateLocalGhosts().
2201 const FieldId staggered_fields[] = {FIELD_ID_UCONT};
2202 ierr = SynchronizePeriodicStaggeredFields(user, 1, staggered_fields); CHKERRQ(ierr);
2203
2204 // FUTURE EXTENSION: Add new cell fields through SynchronizePeriodicCellFields().
2205 /*
2206 if (user->solve_temperature) {
2207 const char *temperature_field[] = {"Temperature"};
2208 ierr = SynchronizePeriodicCellFields(user, 1, temperature_field); CHKERRQ(ierr);
2209 }
2210 */
2211
2212 PROFILE_FUNCTION_END;
2213 PetscFunctionReturn(0);
2214}
2215
2216#undef __FUNCT__
2217#define __FUNCT__ "UpdateDummyCells"
2218/**
2219 * @brief Internal helper implementation: `UpdateDummyCells()`.
2220 * @details Local to this translation unit.
2221 */
2222PetscErrorCode UpdateDummyCells(UserCtx *user)
2223{
2224 PetscErrorCode ierr;
2225 DM fda = user->fda;
2226 DMDALocalInfo info = user->info;
2227 PetscInt xs = info.xs, xe = info.xs + info.xm;
2228 PetscInt ys = info.ys, ye = info.ys + info.ym;
2229 PetscInt zs = info.zs, ze = info.zs + info.zm;
2230 PetscInt mx = info.mx, my = info.my, mz = info.mz;
2231
2232 // --- Calculate shrunken loop ranges to avoid edges and corners ---
2233 PetscInt lxs = xs, lxe = xe;
2234 PetscInt lys = ys, lye = ye;
2235 PetscInt lzs = zs, lze = ze;
2236
2237 if (xs == 0) lxs = xs + 1;
2238 if (ys == 0) lys = ys + 1;
2239 if (zs == 0) lzs = zs + 1;
2240
2241 if (xe == mx) lxe = xe - 1;
2242 if (ye == my) lye = ye - 1;
2243 if (ze == mz) lze = ze - 1;
2244
2245 Cmpnts ***ucat, ***ubcs;
2246 PetscFunctionBeginUser;
2247
2248 ierr = DMDAVecGetArray(fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
2249 ierr = DMDAVecGetArray(fda, user->Ucat, &ucat); CHKERRQ(ierr);
2250
2251 // -X Face
2252 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC && xs == 0) {
2253 for (PetscInt k = lzs; k < lze; k++) for (PetscInt j = lys; j < lye; j++) {
2254 ucat[k][j][xs].x = 2.0 * ubcs[k][j][xs].x - ucat[k][j][xs + 1].x;
2255 ucat[k][j][xs].y = 2.0 * ubcs[k][j][xs].y - ucat[k][j][xs + 1].y;
2256 ucat[k][j][xs].z = 2.0 * ubcs[k][j][xs].z - ucat[k][j][xs + 1].z;
2257 }
2258 }
2259 // +X Face
2260 if (user->boundary_faces[BC_FACE_POS_X].mathematical_type != PERIODIC && xe == mx) {
2261 for (PetscInt k = lzs; k < lze; k++) for (PetscInt j = lys; j < lye; j++) {
2262 ucat[k][j][xe-1].x = 2.0 * ubcs[k][j][xe-1].x - ucat[k][j][xe - 2].x;
2263 ucat[k][j][xe-1].y = 2.0 * ubcs[k][j][xe-1].y - ucat[k][j][xe - 2].y;
2264 ucat[k][j][xe-1].z = 2.0 * ubcs[k][j][xe-1].z - ucat[k][j][xe - 2].z;
2265 }
2266 }
2267
2268 // -Y Face
2269 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC && ys == 0) {
2270 for (PetscInt k = lzs; k < lze; k++) for (PetscInt i = lxs; i < lxe; i++) {
2271 ucat[k][ys][i].x = 2.0 * ubcs[k][ys][i].x - ucat[k][ys + 1][i].x;
2272 ucat[k][ys][i].y = 2.0 * ubcs[k][ys][i].y - ucat[k][ys + 1][i].y;
2273 ucat[k][ys][i].z = 2.0 * ubcs[k][ys][i].z - ucat[k][ys + 1][i].z;
2274 }
2275 }
2276 // +Y Face
2277 if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type != PERIODIC && ye == my) {
2278 for (PetscInt k = lzs; k < lze; k++) for (PetscInt i = lxs; i < lxe; i++) {
2279 ucat[k][ye-1][i].x = 2.0 * ubcs[k][ye-1][i].x - ucat[k][ye-2][i].x;
2280 ucat[k][ye-1][i].y = 2.0 * ubcs[k][ye-1][i].y - ucat[k][ye-2][i].y;
2281 ucat[k][ye-1][i].z = 2.0 * ubcs[k][ye-1][i].z - ucat[k][ye-2][i].z;
2282 }
2283 }
2284
2285 // -Z Face
2286 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC && zs == 0) {
2287 for (PetscInt j = lys; j < lye; j++) for (PetscInt i = lxs; i < lxe; i++) {
2288 ucat[zs][j][i].x = 2.0 * ubcs[zs][j][i].x - ucat[zs + 1][j][i].x;
2289 ucat[zs][j][i].y = 2.0 * ubcs[zs][j][i].y - ucat[zs + 1][j][i].y;
2290 ucat[zs][j][i].z = 2.0 * ubcs[zs][j][i].z - ucat[zs + 1][j][i].z;
2291 }
2292 }
2293 // +Z Face
2294 if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type != PERIODIC && ze == mz) {
2295 for (PetscInt j = lys; j < lye; j++) for (PetscInt i = lxs; i < lxe; i++) {
2296 ucat[ze-1][j][i].x = 2.0 * ubcs[ze-1][j][i].x - ucat[ze-2][j][i].x;
2297 ucat[ze-1][j][i].y = 2.0 * ubcs[ze-1][j][i].y - ucat[ze-2][j][i].y;
2298 ucat[ze-1][j][i].z = 2.0 * ubcs[ze-1][j][i].z - ucat[ze-2][j][i].z;
2299 }
2300 }
2301
2302 ierr = DMDAVecRestoreArray(fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
2303 ierr = DMDAVecRestoreArray(fda, user->Ucat, &ucat); CHKERRQ(ierr);
2304
2305 PetscFunctionReturn(0);
2306}
2307
2308#undef __FUNCT__
2309#define __FUNCT__ "UpdateCornerNodes"
2310/**
2311 * @brief Internal helper implementation: `UpdateCornerNodes()`.
2312 * @details Local to this translation unit.
2313 */
2314PetscErrorCode UpdateCornerNodes(UserCtx *user)
2315{
2316 PetscErrorCode ierr;
2317 DM da = user->da, fda = user->fda;
2318 DMDALocalInfo info = user->info;
2319 PetscInt xs = info.xs, xe = info.xs + info.xm;
2320 PetscInt ys = info.ys, ye = info.ys + info.ym;
2321 PetscInt zs = info.zs, ze = info.zs + info.zm;
2322 PetscInt mx = info.mx, my = info.my, mz = info.mz;
2323
2324 Cmpnts ***ucat;
2325 PetscReal ***p;
2326
2327 PetscFunctionBeginUser;
2328
2329 ierr = DMDAVecGetArray(fda, user->Ucat, &ucat); CHKERRQ(ierr);
2330 ierr = DMDAVecGetArray(da, user->P, &p); CHKERRQ(ierr);
2331
2332 // --- Update Edges and Corners by Averaging ---
2333 // The order of these blocks ensures that corners (where 3 faces meet) are
2334 // computed using data from edges (where 2 faces meet), which are computed first.
2335// Edges connected to the -Z face (k=zs)
2336 if (zs == 0) {
2337 if (xs == 0) {
2338 for (PetscInt j = ys; j < ye; j++) {
2339 p[zs][j][xs] = 0.5 * (p[zs+1][j][xs] + p[zs][j][xs+1]);
2340 ucat[zs][j][xs].x = 0.5 * (ucat[zs+1][j][xs].x + ucat[zs][j][xs+1].x);
2341 ucat[zs][j][xs].y = 0.5 * (ucat[zs+1][j][xs].y + ucat[zs][j][xs+1].y);
2342 ucat[zs][j][xs].z = 0.5 * (ucat[zs+1][j][xs].z + ucat[zs][j][xs+1].z);
2343 }
2344 }
2345 if (xe == mx) {
2346 for (PetscInt j = ys; j < ye; j++) {
2347 p[zs][j][mx-1] = 0.5 * (p[zs+1][j][mx-1] + p[zs][j][mx-2]);
2348 ucat[zs][j][mx-1].x = 0.5 * (ucat[zs+1][j][mx-1].x + ucat[zs][j][mx-2].x);
2349 ucat[zs][j][mx-1].y = 0.5 * (ucat[zs+1][j][mx-1].y + ucat[zs][j][mx-2].y);
2350 ucat[zs][j][mx-1].z = 0.5 * (ucat[zs+1][j][mx-1].z + ucat[zs][j][mx-2].z);
2351 }
2352 }
2353 if (ys == 0) {
2354 for (PetscInt i = xs; i < xe; i++) {
2355 p[zs][ys][i] = 0.5 * (p[zs+1][ys][i] + p[zs][ys+1][i]);
2356 ucat[zs][ys][i].x = 0.5 * (ucat[zs+1][ys][i].x + ucat[zs][ys+1][i].x);
2357 ucat[zs][ys][i].y = 0.5 * (ucat[zs+1][ys][i].y + ucat[zs][ys+1][i].y);
2358 ucat[zs][ys][i].z = 0.5 * (ucat[zs+1][ys][i].z + ucat[zs][ys+1][i].z);
2359 }
2360 }
2361 if (ye == my) {
2362 for (PetscInt i = xs; i < xe; i++) {
2363 p[zs][my-1][i] = 0.5 * (p[zs+1][my-1][i] + p[zs][my-2][i]);
2364 ucat[zs][my-1][i].x = 0.5 * (ucat[zs+1][my-1][i].x + ucat[zs][my-2][i].x);
2365 ucat[zs][my-1][i].y = 0.5 * (ucat[zs+1][my-1][i].y + ucat[zs][my-2][i].y);
2366 ucat[zs][my-1][i].z = 0.5 * (ucat[zs+1][my-1][i].z + ucat[zs][my-2][i].z);
2367 }
2368 }
2369 }
2370
2371 // Edges connected to the +Z face (k=ze-1)
2372 if (ze == mz) {
2373 if (xs == 0) {
2374 for (PetscInt j = ys; j < ye; j++) {
2375 p[mz-1][j][xs] = 0.5 * (p[mz-2][j][xs] + p[mz-1][j][xs+1]);
2376 ucat[mz-1][j][xs].x = 0.5 * (ucat[mz-2][j][xs].x + ucat[mz-1][j][xs+1].x);
2377 ucat[mz-1][j][xs].y = 0.5 * (ucat[mz-2][j][xs].y + ucat[mz-1][j][xs+1].y);
2378 ucat[mz-1][j][xs].z = 0.5 * (ucat[mz-2][j][xs].z + ucat[mz-1][j][xs+1].z);
2379 }
2380 }
2381 if (xe == mx) {
2382 for (PetscInt j = ys; j < ye; j++) {
2383 p[mz-1][j][mx-1] = 0.5 * (p[mz-2][j][mx-1] + p[mz-1][j][mx-2]);
2384 ucat[mz-1][j][mx-1].x = 0.5 * (ucat[mz-2][j][mx-1].x + ucat[mz-1][j][mx-2].x);
2385 ucat[mz-1][j][mx-1].y = 0.5 * (ucat[mz-2][j][mx-1].y + ucat[mz-1][j][mx-2].y);
2386 ucat[mz-1][j][mx-1].z = 0.5 * (ucat[mz-2][j][mx-1].z + ucat[mz-1][j][mx-2].z);
2387 }
2388 }
2389 if (ys == 0) {
2390 for (PetscInt i = xs; i < xe; i++) {
2391 p[mz-1][ys][i] = 0.5 * (p[mz-2][ys][i] + p[mz-1][ys+1][i]);
2392 ucat[mz-1][ys][i].x = 0.5 * (ucat[mz-2][ys][i].x + ucat[mz-1][ys+1][i].x);
2393 ucat[mz-1][ys][i].y = 0.5 * (ucat[mz-2][ys][i].y + ucat[mz-1][ys+1][i].y);
2394 ucat[mz-1][ys][i].z = 0.5 * (ucat[mz-2][ys][i].z + ucat[mz-1][ys+1][i].z);
2395 }
2396 }
2397 if (ye == my) {
2398 for (PetscInt i = xs; i < xe; i++) {
2399 p[mz-1][my-1][i] = 0.5 * (p[mz-2][my-1][i] + p[mz-1][my-2][i]);
2400 ucat[mz-1][my-1][i].x = 0.5 * (ucat[mz-2][my-1][i].x + ucat[mz-1][my-2][i].x);
2401 ucat[mz-1][my-1][i].y = 0.5 * (ucat[mz-2][my-1][i].y + ucat[mz-1][my-2][i].y);
2402 ucat[mz-1][my-1][i].z = 0.5 * (ucat[mz-2][my-1][i].z + ucat[mz-1][my-2][i].z);
2403 }
2404 }
2405 }
2406
2407 // Remaining edges on the XY plane (that are not on Z faces)
2408 if (ys == 0) {
2409 if (xs == 0) {
2410 for (PetscInt k = zs; k < ze; k++) {
2411 p[k][ys][xs] = 0.5 * (p[k][ys+1][xs] + p[k][ys][xs+1]);
2412 ucat[k][ys][xs].x = 0.5 * (ucat[k][ys+1][xs].x + ucat[k][ys][xs+1].x);
2413 ucat[k][ys][xs].y = 0.5 * (ucat[k][ys+1][xs].y + ucat[k][ys][xs+1].y);
2414 ucat[k][ys][xs].z = 0.5 * (ucat[k][ys+1][xs].z + ucat[k][ys][xs+1].z);
2415 }
2416 }
2417 if (xe == mx) {
2418 for (PetscInt k = zs; k < ze; k++) {
2419 p[k][ys][mx-1] = 0.5 * (p[k][ys+1][mx-1] + p[k][ys][mx-2]);
2420 ucat[k][ys][mx-1].x = 0.5 * (ucat[k][ys+1][mx-1].x + ucat[k][ys][mx-2].x);
2421 ucat[k][ys][mx-1].y = 0.5 * (ucat[k][ys+1][mx-1].y + ucat[k][ys][mx-2].y);
2422 ucat[k][ys][mx-1].z = 0.5 * (ucat[k][ys+1][mx-1].z + ucat[k][ys][mx-2].z);
2423 }
2424 }
2425 }
2426
2427 if (ye == my) {
2428 if (xs == 0) {
2429 for (PetscInt k = zs; k < ze; k++) {
2430 p[k][my-1][xs] = 0.5 * (p[k][my-2][xs] + p[k][my-1][xs+1]);
2431 ucat[k][my-1][xs].x = 0.5 * (ucat[k][my-2][xs].x + ucat[k][my-1][xs+1].x);
2432 ucat[k][my-1][xs].y = 0.5 * (ucat[k][my-2][xs].y + ucat[k][my-1][xs+1].y);
2433 ucat[k][my-1][xs].z = 0.5 * (ucat[k][my-2][xs].z + ucat[k][my-1][xs+1].z);
2434 }
2435 }
2436 if (xe == mx) {
2437 for (PetscInt k = zs; k < ze; k++) {
2438 p[k][my-1][mx-1] = 0.5 * (p[k][my-2][mx-1] + p[k][my-1][mx-2]);
2439 ucat[k][my-1][mx-1].x = 0.5 * (ucat[k][my-2][mx-1].x + ucat[k][my-1][mx-2].x);
2440 ucat[k][my-1][mx-1].y = 0.5 * (ucat[k][my-2][mx-1].y + ucat[k][my-1][mx-2].y);
2441 ucat[k][my-1][mx-1].z = 0.5 * (ucat[k][my-2][mx-1].z + ucat[k][my-1][mx-2].z);
2442 }
2443 }
2444 }
2445
2446 ierr = DMDAVecRestoreArray(fda, user->Ucat, &ucat); CHKERRQ(ierr);
2447 ierr = DMDAVecRestoreArray(da, user->P, &p); CHKERRQ(ierr);
2448
2449 PetscFunctionReturn(0);
2450}
2451
2452#undef __FUNCT__
2453#define __FUNCT__ "ApplyWallFunction"
2454/**
2455 * @brief Internal helper implementation: `ApplyWallFunction()`.
2456 * @details Local to this translation unit.
2457 */
2458PetscErrorCode ApplyWallFunction(UserCtx *user)
2459{
2460 PetscErrorCode ierr;
2461 SimCtx *simCtx = user->simCtx;
2462 DMDALocalInfo *info = &user->info;
2463
2464 PetscFunctionBeginUser;
2465
2466 // =========================================================================
2467 // STEP 0: Early exit if wall functions are disabled
2468 // =========================================================================
2469 if (!simCtx->wallfunction) {
2470 PetscFunctionReturn(0);
2471 }
2472
2473 LOG_ALLOW(LOCAL, LOG_DEBUG, "Processing wall function boundaries.\n");
2474
2475 // =========================================================================
2476 // STEP 1: Get read/write access to all necessary field arrays
2477 // =========================================================================
2478 Cmpnts ***velocity_cartesian; // Cartesian velocity (modified)
2479 Cmpnts ***velocity_contravariant; // Contravariant velocity (set to zero at walls)
2480 Cmpnts ***velocity_boundary; // Boundary condition velocity (kept at zero)
2481 Cmpnts ***csi, ***eta, ***zet; // Metric tensor components (face normals)
2482 PetscReal ***node_vertex_flag; // Fluid/solid indicator (0=fluid, 1=solid)
2483 PetscReal ***cell_jacobian; // Grid Jacobian (1/volume)
2484 PetscReal ***friction_velocity; // u_tau (friction velocity field)
2485
2486 ierr = DMDAVecGetArray(user->fda, user->Ucat, &velocity_cartesian); CHKERRQ(ierr);
2487 ierr = DMDAVecGetArray(user->fda, user->Ucont, &velocity_contravariant); CHKERRQ(ierr);
2488 ierr = DMDAVecGetArray(user->fda, user->Bcs.Ubcs, &velocity_boundary); CHKERRQ(ierr);
2489 ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
2490 ierr = DMDAVecGetArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
2491 ierr = DMDAVecGetArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
2492 ierr = DMDAVecGetArrayRead(user->da, user->lNvert, (const PetscReal***)&node_vertex_flag); CHKERRQ(ierr);
2493 ierr = DMDAVecGetArrayRead(user->da, user->lAj, (const PetscReal***)&cell_jacobian); CHKERRQ(ierr);
2494 ierr = DMDAVecGetArray(user->da, user->lFriction_Velocity, &friction_velocity); CHKERRQ(ierr);
2495
2496 // =========================================================================
2497 // STEP 2: Define loop bounds (owned portion of the grid for this MPI rank)
2498 // =========================================================================
2499 PetscInt grid_start_i = info->xs, grid_end_i = info->xs + info->xm;
2500 PetscInt grid_start_j = info->ys, grid_end_j = info->ys + info->ym;
2501 PetscInt grid_start_k = info->zs, grid_end_k = info->zs + info->zm;
2502 PetscInt grid_size_i = info->mx, grid_size_j = info->my, grid_size_k = info->mz;
2503
2504 // Shrunken loop bounds: exclude domain edges and corners to avoid double-counting
2505 PetscInt loop_start_i = grid_start_i, loop_end_i = grid_end_i;
2506 PetscInt loop_start_j = grid_start_j, loop_end_j = grid_end_j;
2507 PetscInt loop_start_k = grid_start_k, loop_end_k = grid_end_k;
2508
2509 if (grid_start_i == 0) loop_start_i = grid_start_i + 1;
2510 if (grid_end_i == grid_size_i) loop_end_i = grid_end_i - 1;
2511 if (grid_start_j == 0) loop_start_j = grid_start_j + 1;
2512 if (grid_end_j == grid_size_j) loop_end_j = grid_end_j - 1;
2513 if (grid_start_k == 0) loop_start_k = grid_start_k + 1;
2514 if (grid_end_k == grid_size_k) loop_end_k = grid_end_k - 1;
2515
2516 // Wall roughness parameter (smooth wall by default, configurable via -wall_roughness).
2517 const PetscReal wall_roughness_height = user->simCtx->wall_roughness_height;
2518
2519 // =========================================================================
2520 // STEP 3: Process each of the 6 domain faces
2521 // =========================================================================
2522 for (int face_index = 0; face_index < 6; face_index++) {
2523 BCFace current_face_id = (BCFace)face_index;
2524 BoundaryFaceConfig *face_config = &user->boundary_faces[current_face_id];
2525
2526 // Only process faces that are mathematical walls (applies to no-slip, moving, slip, etc.)
2527 if (face_config->mathematical_type != WALL) {
2528 continue;
2529 }
2530
2531 // Check if this MPI rank owns part of this face
2532 PetscBool rank_owns_this_face;
2533 ierr = CanRankServiceFace(info, user->IM, user->JM, user->KM,
2534 current_face_id, &rank_owns_this_face); CHKERRQ(ierr);
2535
2536 if (!rank_owns_this_face) {
2537 continue;
2538 }
2539
2540 LOG_ALLOW(LOCAL, LOG_TRACE, "Processing Face %d (%s)\n",
2541 current_face_id, BCFaceToString(current_face_id));
2542
2543 // =====================================================================
2544 // Process each face with appropriate indexing
2545 // =====================================================================
2546 switch(current_face_id) {
2547
2548 // =================================================================
2549 // NEGATIVE X FACE (i = 0, normal points in +X direction)
2550 // =================================================================
2551 case BC_FACE_NEG_X: {
2552 if (grid_start_i == 0) {
2553 const PetscInt ghost_cell_index = grid_start_i;
2554 const PetscInt first_interior_cell = grid_start_i + 1;
2555 const PetscInt second_interior_cell = grid_start_i + 2;
2556
2557 for (PetscInt k = loop_start_k; k < loop_end_k; k++) {
2558 for (PetscInt j = loop_start_j; j < loop_end_j; j++) {
2559
2560 // Skip if this is a solid cell (embedded boundary)
2561 if (node_vertex_flag[k][j][first_interior_cell] < 0.1) {
2562
2563 // Calculate face area from contravariant metric tensor
2564 PetscReal face_area = sqrt(
2565 csi[k][j][ghost_cell_index].x * csi[k][j][ghost_cell_index].x +
2566 csi[k][j][ghost_cell_index].y * csi[k][j][ghost_cell_index].y +
2567 csi[k][j][ghost_cell_index].z * csi[k][j][ghost_cell_index].z
2568 );
2569
2570 // Compute wall-normal distances using cell Jacobians
2571 // sb = distance from wall to first interior cell center
2572 // sc = distance from wall to second interior cell center
2573 PetscReal distance_to_first_cell = 0.5 / cell_jacobian[k][j][first_interior_cell] / face_area;
2574 PetscReal distance_to_second_cell = 2.0 * distance_to_first_cell +
2575 0.5 / cell_jacobian[k][j][second_interior_cell] / face_area;
2576
2577 // Compute unit normal vector pointing INTO the domain
2578 PetscReal wall_normal[3];
2579 wall_normal[0] = csi[k][j][ghost_cell_index].x / face_area;
2580 wall_normal[1] = csi[k][j][ghost_cell_index].y / face_area;
2581 wall_normal[2] = csi[k][j][ghost_cell_index].z / face_area;
2582
2583 // Define velocities for wall function calculation
2584 Cmpnts wall_velocity; // Ua = velocity at wall (zero for stationary wall)
2585 Cmpnts reference_velocity; // Uc = velocity at second interior cell
2586
2587 wall_velocity.x = wall_velocity.y = wall_velocity.z = 0.0;
2588 reference_velocity = velocity_cartesian[k][j][second_interior_cell];
2589
2590 // Step 1: Linear interpolation (provides initial guess)
2591 noslip(user, distance_to_second_cell, distance_to_first_cell,
2592 wall_velocity, reference_velocity,
2593 &velocity_cartesian[k][j][first_interior_cell],
2594 wall_normal[0], wall_normal[1], wall_normal[2]);
2595
2596 // Step 2: Apply log-law correction (improves near-wall velocity)
2597 wall_function_loglaw(user, wall_roughness_height,
2598 distance_to_second_cell, distance_to_first_cell,
2599 wall_velocity, reference_velocity,
2600 &velocity_cartesian[k][j][first_interior_cell],
2601 &friction_velocity[k][j][first_interior_cell],
2602 wall_normal[0], wall_normal[1], wall_normal[2]);
2603
2604 // Ensure ghost cell BC remains zero (required for proper extrapolation)
2605 velocity_boundary[k][j][ghost_cell_index].x = 0.0;
2606 velocity_boundary[k][j][ghost_cell_index].y = 0.0;
2607 velocity_boundary[k][j][ghost_cell_index].z = 0.0;
2608 velocity_contravariant[k][j][ghost_cell_index].x = 0.0;
2609 }
2610 }
2611 }
2612 }
2613 } break;
2614
2615 // =================================================================
2616 // POSITIVE X FACE (i = mx-1, normal points in -X direction)
2617 // =================================================================
2618 case BC_FACE_POS_X: {
2619 if (grid_end_i == grid_size_i) {
2620 const PetscInt ghost_cell_index = grid_end_i - 1;
2621 const PetscInt first_interior_cell = grid_end_i - 2;
2622 const PetscInt second_interior_cell = grid_end_i - 3;
2623
2624 for (PetscInt k = loop_start_k; k < loop_end_k; k++) {
2625 for (PetscInt j = loop_start_j; j < loop_end_j; j++) {
2626
2627 if (node_vertex_flag[k][j][first_interior_cell] < 0.1) {
2628
2629 PetscReal face_area = sqrt(
2630 csi[k][j][first_interior_cell].x * csi[k][j][first_interior_cell].x +
2631 csi[k][j][first_interior_cell].y * csi[k][j][first_interior_cell].y +
2632 csi[k][j][first_interior_cell].z * csi[k][j][first_interior_cell].z
2633 );
2634
2635 PetscReal distance_to_first_cell = 0.5 / cell_jacobian[k][j][first_interior_cell] / face_area;
2636 PetscReal distance_to_second_cell = 2.0 * distance_to_first_cell +
2637 0.5 / cell_jacobian[k][j][second_interior_cell] / face_area;
2638
2639 // Note: Normal flipped for +X face to point INTO domain
2640 PetscReal wall_normal[3];
2641 wall_normal[0] = -csi[k][j][first_interior_cell].x / face_area;
2642 wall_normal[1] = -csi[k][j][first_interior_cell].y / face_area;
2643 wall_normal[2] = -csi[k][j][first_interior_cell].z / face_area;
2644
2645 Cmpnts wall_velocity, reference_velocity;
2646 wall_velocity.x = wall_velocity.y = wall_velocity.z = 0.0;
2647 reference_velocity = velocity_cartesian[k][j][second_interior_cell];
2648
2649 noslip(user, distance_to_second_cell, distance_to_first_cell,
2650 wall_velocity, reference_velocity,
2651 &velocity_cartesian[k][j][first_interior_cell],
2652 wall_normal[0], wall_normal[1], wall_normal[2]);
2653
2654 wall_function_loglaw(user, wall_roughness_height,
2655 distance_to_second_cell, distance_to_first_cell,
2656 wall_velocity, reference_velocity,
2657 &velocity_cartesian[k][j][first_interior_cell],
2658 &friction_velocity[k][j][first_interior_cell],
2659 wall_normal[0], wall_normal[1], wall_normal[2]);
2660
2661 velocity_boundary[k][j][ghost_cell_index].x = 0.0;
2662 velocity_boundary[k][j][ghost_cell_index].y = 0.0;
2663 velocity_boundary[k][j][ghost_cell_index].z = 0.0;
2664 velocity_contravariant[k][j][first_interior_cell].x = 0.0;
2665 }
2666 }
2667 }
2668 }
2669 } break;
2670
2671 // =================================================================
2672 // NEGATIVE Y FACE (j = 0, normal points in +Y direction)
2673 // =================================================================
2674 case BC_FACE_NEG_Y: {
2675 if (grid_start_j == 0) {
2676 const PetscInt ghost_cell_index = grid_start_j;
2677 const PetscInt first_interior_cell = grid_start_j + 1;
2678 const PetscInt second_interior_cell = grid_start_j + 2;
2679
2680 for (PetscInt k = loop_start_k; k < loop_end_k; k++) {
2681 for (PetscInt i = loop_start_i; i < loop_end_i; i++) {
2682
2683 if (node_vertex_flag[k][first_interior_cell][i] < 0.1) {
2684
2685 PetscReal face_area = sqrt(
2686 eta[k][ghost_cell_index][i].x * eta[k][ghost_cell_index][i].x +
2687 eta[k][ghost_cell_index][i].y * eta[k][ghost_cell_index][i].y +
2688 eta[k][ghost_cell_index][i].z * eta[k][ghost_cell_index][i].z
2689 );
2690
2691 PetscReal distance_to_first_cell = 0.5 / cell_jacobian[k][first_interior_cell][i] / face_area;
2692 PetscReal distance_to_second_cell = 2.0 * distance_to_first_cell +
2693 0.5 / cell_jacobian[k][second_interior_cell][i] / face_area;
2694
2695 PetscReal wall_normal[3];
2696 wall_normal[0] = eta[k][ghost_cell_index][i].x / face_area;
2697 wall_normal[1] = eta[k][ghost_cell_index][i].y / face_area;
2698 wall_normal[2] = eta[k][ghost_cell_index][i].z / face_area;
2699
2700 Cmpnts wall_velocity, reference_velocity;
2701 wall_velocity.x = wall_velocity.y = wall_velocity.z = 0.0;
2702 reference_velocity = velocity_cartesian[k][second_interior_cell][i];
2703
2704 noslip(user, distance_to_second_cell, distance_to_first_cell,
2705 wall_velocity, reference_velocity,
2706 &velocity_cartesian[k][first_interior_cell][i],
2707 wall_normal[0], wall_normal[1], wall_normal[2]);
2708
2709 wall_function_loglaw(user, wall_roughness_height,
2710 distance_to_second_cell, distance_to_first_cell,
2711 wall_velocity, reference_velocity,
2712 &velocity_cartesian[k][first_interior_cell][i],
2713 &friction_velocity[k][first_interior_cell][i],
2714 wall_normal[0], wall_normal[1], wall_normal[2]);
2715
2716 velocity_boundary[k][ghost_cell_index][i].x = 0.0;
2717 velocity_boundary[k][ghost_cell_index][i].y = 0.0;
2718 velocity_boundary[k][ghost_cell_index][i].z = 0.0;
2719 velocity_contravariant[k][ghost_cell_index][i].y = 0.0;
2720 }
2721 }
2722 }
2723 }
2724 } break;
2725
2726 // =================================================================
2727 // POSITIVE Y FACE (j = my-1, normal points in -Y direction)
2728 // =================================================================
2729 case BC_FACE_POS_Y: {
2730 if (grid_end_j == grid_size_j) {
2731 const PetscInt ghost_cell_index = grid_end_j - 1;
2732 const PetscInt first_interior_cell = grid_end_j - 2;
2733 const PetscInt second_interior_cell = grid_end_j - 3;
2734
2735 for (PetscInt k = loop_start_k; k < loop_end_k; k++) {
2736 for (PetscInt i = loop_start_i; i < loop_end_i; i++) {
2737
2738 if (node_vertex_flag[k][first_interior_cell][i] < 0.1) {
2739
2740 PetscReal face_area = sqrt(
2741 eta[k][first_interior_cell][i].x * eta[k][first_interior_cell][i].x +
2742 eta[k][first_interior_cell][i].y * eta[k][first_interior_cell][i].y +
2743 eta[k][first_interior_cell][i].z * eta[k][first_interior_cell][i].z
2744 );
2745
2746 PetscReal distance_to_first_cell = 0.5 / cell_jacobian[k][first_interior_cell][i] / face_area;
2747 PetscReal distance_to_second_cell = 2.0 * distance_to_first_cell +
2748 0.5 / cell_jacobian[k][second_interior_cell][i] / face_area;
2749
2750 PetscReal wall_normal[3];
2751 wall_normal[0] = -eta[k][first_interior_cell][i].x / face_area;
2752 wall_normal[1] = -eta[k][first_interior_cell][i].y / face_area;
2753 wall_normal[2] = -eta[k][first_interior_cell][i].z / face_area;
2754
2755 Cmpnts wall_velocity, reference_velocity;
2756 wall_velocity.x = wall_velocity.y = wall_velocity.z = 0.0;
2757 reference_velocity = velocity_cartesian[k][second_interior_cell][i];
2758
2759 noslip(user, distance_to_second_cell, distance_to_first_cell,
2760 wall_velocity, reference_velocity,
2761 &velocity_cartesian[k][first_interior_cell][i],
2762 wall_normal[0], wall_normal[1], wall_normal[2]);
2763
2764 wall_function_loglaw(user, wall_roughness_height,
2765 distance_to_second_cell, distance_to_first_cell,
2766 wall_velocity, reference_velocity,
2767 &velocity_cartesian[k][first_interior_cell][i],
2768 &friction_velocity[k][first_interior_cell][i],
2769 wall_normal[0], wall_normal[1], wall_normal[2]);
2770
2771 velocity_boundary[k][ghost_cell_index][i].x = 0.0;
2772 velocity_boundary[k][ghost_cell_index][i].y = 0.0;
2773 velocity_boundary[k][ghost_cell_index][i].z = 0.0;
2774 velocity_contravariant[k][first_interior_cell][i].y = 0.0;
2775 }
2776 }
2777 }
2778 }
2779 } break;
2780
2781 // =================================================================
2782 // NEGATIVE Z FACE (k = 0, normal points in +Z direction)
2783 // =================================================================
2784 case BC_FACE_NEG_Z: {
2785 if (grid_start_k == 0) {
2786 const PetscInt ghost_cell_index = grid_start_k;
2787 const PetscInt first_interior_cell = grid_start_k + 1;
2788 const PetscInt second_interior_cell = grid_start_k + 2;
2789
2790 for (PetscInt j = loop_start_j; j < loop_end_j; j++) {
2791 for (PetscInt i = loop_start_i; i < loop_end_i; i++) {
2792
2793 if (node_vertex_flag[first_interior_cell][j][i] < 0.1) {
2794
2795 PetscReal face_area = sqrt(
2796 zet[ghost_cell_index][j][i].x * zet[ghost_cell_index][j][i].x +
2797 zet[ghost_cell_index][j][i].y * zet[ghost_cell_index][j][i].y +
2798 zet[ghost_cell_index][j][i].z * zet[ghost_cell_index][j][i].z
2799 );
2800
2801 PetscReal distance_to_first_cell = 0.5 / cell_jacobian[first_interior_cell][j][i] / face_area;
2802 PetscReal distance_to_second_cell = 2.0 * distance_to_first_cell +
2803 0.5 / cell_jacobian[second_interior_cell][j][i] / face_area;
2804
2805 PetscReal wall_normal[3];
2806 wall_normal[0] = zet[ghost_cell_index][j][i].x / face_area;
2807 wall_normal[1] = zet[ghost_cell_index][j][i].y / face_area;
2808 wall_normal[2] = zet[ghost_cell_index][j][i].z / face_area;
2809
2810 Cmpnts wall_velocity, reference_velocity;
2811 wall_velocity.x = wall_velocity.y = wall_velocity.z = 0.0;
2812 reference_velocity = velocity_cartesian[second_interior_cell][j][i];
2813
2814 noslip(user, distance_to_second_cell, distance_to_first_cell,
2815 wall_velocity, reference_velocity,
2816 &velocity_cartesian[first_interior_cell][j][i],
2817 wall_normal[0], wall_normal[1], wall_normal[2]);
2818
2819 wall_function_loglaw(user, wall_roughness_height,
2820 distance_to_second_cell, distance_to_first_cell,
2821 wall_velocity, reference_velocity,
2822 &velocity_cartesian[first_interior_cell][j][i],
2823 &friction_velocity[first_interior_cell][j][i],
2824 wall_normal[0], wall_normal[1], wall_normal[2]);
2825
2826 velocity_boundary[ghost_cell_index][j][i].x = 0.0;
2827 velocity_boundary[ghost_cell_index][j][i].y = 0.0;
2828 velocity_boundary[ghost_cell_index][j][i].z = 0.0;
2829 velocity_contravariant[ghost_cell_index][j][i].z = 0.0;
2830 }
2831 }
2832 }
2833 }
2834 } break;
2835
2836 // =================================================================
2837 // POSITIVE Z FACE (k = mz-1, normal points in -Z direction)
2838 // =================================================================
2839 case BC_FACE_POS_Z: {
2840 if (grid_end_k == grid_size_k) {
2841 const PetscInt ghost_cell_index = grid_end_k - 1;
2842 const PetscInt first_interior_cell = grid_end_k - 2;
2843 const PetscInt second_interior_cell = grid_end_k - 3;
2844
2845 for (PetscInt j = loop_start_j; j < loop_end_j; j++) {
2846 for (PetscInt i = loop_start_i; i < loop_end_i; i++) {
2847
2848 if (node_vertex_flag[first_interior_cell][j][i] < 0.1) {
2849
2850 PetscReal face_area = sqrt(
2851 zet[first_interior_cell][j][i].x * zet[first_interior_cell][j][i].x +
2852 zet[first_interior_cell][j][i].y * zet[first_interior_cell][j][i].y +
2853 zet[first_interior_cell][j][i].z * zet[first_interior_cell][j][i].z
2854 );
2855
2856 PetscReal distance_to_first_cell = 0.5 / cell_jacobian[first_interior_cell][j][i] / face_area;
2857 PetscReal distance_to_second_cell = 2.0 * distance_to_first_cell +
2858 0.5 / cell_jacobian[second_interior_cell][j][i] / face_area;
2859
2860 PetscReal wall_normal[3];
2861 wall_normal[0] = -zet[first_interior_cell][j][i].x / face_area;
2862 wall_normal[1] = -zet[first_interior_cell][j][i].y / face_area;
2863 wall_normal[2] = -zet[first_interior_cell][j][i].z / face_area;
2864
2865 Cmpnts wall_velocity, reference_velocity;
2866 wall_velocity.x = wall_velocity.y = wall_velocity.z = 0.0;
2867 reference_velocity = velocity_cartesian[second_interior_cell][j][i];
2868
2869 noslip(user, distance_to_second_cell, distance_to_first_cell,
2870 wall_velocity, reference_velocity,
2871 &velocity_cartesian[first_interior_cell][j][i],
2872 wall_normal[0], wall_normal[1], wall_normal[2]);
2873
2874 wall_function_loglaw(user, wall_roughness_height,
2875 distance_to_second_cell, distance_to_first_cell,
2876 wall_velocity, reference_velocity,
2877 &velocity_cartesian[first_interior_cell][j][i],
2878 &friction_velocity[first_interior_cell][j][i],
2879 wall_normal[0], wall_normal[1], wall_normal[2]);
2880
2881 velocity_boundary[ghost_cell_index][j][i].x = 0.0;
2882 velocity_boundary[ghost_cell_index][j][i].y = 0.0;
2883 velocity_boundary[ghost_cell_index][j][i].z = 0.0;
2884 velocity_contravariant[first_interior_cell][j][i].z = 0.0;
2885 }
2886 }
2887 }
2888 }
2889 } break;
2890 }
2891 }
2892
2893 // =========================================================================
2894 // STEP 4: Restore all arrays and release memory
2895 // =========================================================================
2896 ierr = DMDAVecRestoreArray(user->fda, user->Ucat, &velocity_cartesian); CHKERRQ(ierr);
2897 ierr = DMDAVecRestoreArray(user->fda, user->Ucont, &velocity_contravariant); CHKERRQ(ierr);
2898 ierr = DMDAVecRestoreArray(user->fda, user->Bcs.Ubcs, &velocity_boundary); CHKERRQ(ierr);
2899 ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
2900 ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
2901 ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
2902 ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, (const PetscReal***)&node_vertex_flag); CHKERRQ(ierr);
2903 ierr = DMDAVecRestoreArrayRead(user->da, user->lAj, (const PetscReal***)&cell_jacobian); CHKERRQ(ierr);
2904 ierr = DMDAVecRestoreArray(user->da, user->lFriction_Velocity, &friction_velocity); CHKERRQ(ierr);
2905
2906 LOG_ALLOW(LOCAL, LOG_DEBUG, "Complete.\n");
2907
2908 PetscFunctionReturn(0);
2909}
2910
2911#undef __FUNCT__
2912#define __FUNCT__ "FinalizePostProjectionCellFields"
2913/**
2914 * @brief Implementation of \ref FinalizePostProjectionCellFields().
2915 * @details Full API contract is documented with the header declaration in
2916 * `include/Boundaries.h`.
2917 */
2918PetscErrorCode FinalizePostProjectionCellFields(UserCtx *user)
2919{
2920 PetscErrorCode ierr;
2921 const FieldId cell_fields[] = {FIELD_ID_UCAT, FIELD_ID_P};
2922
2923 PetscFunctionBeginUser;
2924 PROFILE_FUNCTION_BEGIN;
2925
2926 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Finalizing post-projection cell-centered fields.\n");
2927
2928 // Ensure flow-dependent Ubcs handlers see the newly reconstructed Ucat.
2929 ierr = UpdateLocalGhosts(user, FIELD_ID_UCAT); CHKERRQ(ierr);
2930 ierr = BoundarySystem_RefreshUbcs(user); CHKERRQ(ierr);
2931
2932 // Establish flat non-periodic faces and periodic endpoints before corners.
2933 ierr = UpdateDummyCells(user); CHKERRQ(ierr);
2934 ierr = SynchronizePeriodicCellFields(user, 2, cell_fields); CHKERRQ(ierr);
2935
2936 // Corner averaging can overwrite periodic endpoints, so restore them after.
2937 ierr = UpdateCornerNodes(user); CHKERRQ(ierr);
2938 ierr = SynchronizePeriodicCellFields(user, 2, cell_fields); CHKERRQ(ierr);
2939
2940 // Synchronize explicitly because the periodic helper is a no-op when every
2941 // direction is non-periodic.
2942 ierr = UpdateLocalGhosts(user, FIELD_ID_UCAT); CHKERRQ(ierr);
2943 ierr = UpdateLocalGhosts(user, FIELD_ID_P); CHKERRQ(ierr);
2944
2945 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Post-projection cell-centered fields finalized.\n");
2946 PROFILE_FUNCTION_END;
2947 PetscFunctionReturn(0);
2948}
2949
2950#undef __FUNCT__
2951#define __FUNCT__ "ApplyBoundaryConditions"
2952/**
2953 * @brief Implementation of \ref ApplyBoundaryConditions().
2954 * @details Full API contract (arguments, ownership, side effects) is documented with
2955 * the header declaration in `include/Boundaries.h`.
2956 * @see ApplyBoundaryConditions()
2957 */
2958PetscErrorCode ApplyBoundaryConditions(UserCtx *user)
2959{
2960 PetscErrorCode ierr;
2961 const FieldId staggered_fields[] = {FIELD_ID_UCONT};
2962 PetscFunctionBeginUser;
2963 PROFILE_FUNCTION_BEGIN;
2964
2965 LOG_ALLOW(GLOBAL,LOG_TRACE,"Boundary Condition Application begins.\n");
2966
2967 // STEP 1: Main iteration loop for applying and converging non-periodic BCs.
2968 // The number of iterations (e.g., 3) allows information to propagate
2969 // between coupled boundaries, like an inlet and a conserving outlet.
2970 for (PetscInt iter = 0; iter < 3; iter++) {
2971 // (a) Execute the boundary system. This phase calculates fluxes across
2972 // the domain and then applies the physical logic for each non-periodic
2973 // handler, setting the `ubcs` (boundary value) array.
2974 ierr = BoundarySystem_ExecuteStep(user); CHKERRQ(ierr);
2975
2976 LOG_ALLOW(GLOBAL,LOG_VERBOSE,"Boundary Condition Setup Executed.\n");
2977
2978 // (b) Synchronize the updated ghost cells across all processors to ensure
2979 // all ucont values are current before updating the dummy cells.
2980 ierr = SynchronizePeriodicStaggeredFields(user, 1, staggered_fields); CHKERRQ(ierr);
2981
2982 // (c) Convert updated Contravariant velocities to Cartesian velocities.
2983 ierr = Contra2Cart(user); CHKERRQ(ierr);
2984
2985 // (d) Synchronize the updated Cartesian velocities across all processors
2986 // to ensure all ucat values are current before updating the dummy cells.
2987 ierr = UpdateLocalGhosts(user, FIELD_ID_UCAT); CHKERRQ(ierr);
2988
2989 // (e) If Wall functions are enabled, apply them now to adjust near-wall velocities.
2990 if(user->simCtx->wallfunction){
2991 // Apply wall function adjustments to the boundary velocities.
2992 ierr = ApplyWallFunction(user); CHKERRQ(ierr);
2993
2994 // Synchronize the updated Cartesian velocities after wall function adjustments.
2995 ierr = UpdateLocalGhosts(user, FIELD_ID_UCAT); CHKERRQ(ierr);
2996
2997 LOG_ALLOW(GLOBAL,LOG_VERBOSE,"Wall Function Applied at Walls.\n");
2998 }
2999
3000 // (f) Update the first layer of ghost cells for non-periodic faces using
3001 // the newly computed `ubcs` values.
3002 ierr = UpdateDummyCells(user); CHKERRQ(ierr);
3003
3004 LOG_ALLOW(GLOBAL,LOG_VERBOSE,"Dummy Cells/Ghost Cells Updated.\n");
3005
3006 // (g) Handle all periodic boundaries. This is a parallel direct copy
3007 // that sets the absolute constraints for the rest of the solve.
3008 // There is a Ghost update happening inside this function.
3009 ierr = ApplyPeriodicBCs(user); CHKERRQ(ierr);
3010
3011 // (h) Update the corner and edge ghost nodes. This routine calculates
3012 // values for corners/edges by averaging their neighbors, which have been
3013 // finalized in the steps above (both periodic and non-periodic).
3014 ierr = UpdateCornerNodes(user); CHKERRQ(ierr);
3015
3016 // (i) Synchronize the updated edge and corner cells across all processors to ensure
3017 // consistency before the next iteration or finalization.
3018 ierr = UpdateLocalGhosts(user, FIELD_ID_P); CHKERRQ(ierr);
3019 ierr = UpdateLocalGhosts(user, FIELD_ID_UCAT); CHKERRQ(ierr);
3020 ierr = SynchronizePeriodicStaggeredFields(user, 1, staggered_fields); CHKERRQ(ierr);
3021
3022 // (j) Ensure All the corners are synchronized with a well defined protocol in case of Periodic boundary conditions
3023 // To avoid race conditions.
3024 const FieldId all_fields[] = {FIELD_ID_UCAT, FIELD_ID_P, FIELD_ID_NVERT};
3025 ierr = SynchronizePeriodicCellFields(user, 3, all_fields); CHKERRQ(ierr);
3026
3027 }
3028
3029 // STEP 3: Final ghost node synchronization. This ensures all changes made
3030 // to the global vectors are reflected in the local ghost regions of all
3031 // processors, making the state fully consistent before the next solver stage.
3032 ierr = UpdateLocalGhosts(user, FIELD_ID_P); CHKERRQ(ierr);
3033 ierr = UpdateLocalGhosts(user, FIELD_ID_UCAT); CHKERRQ(ierr);
3034 ierr = SynchronizePeriodicStaggeredFields(user, 1, staggered_fields); CHKERRQ(ierr);
3035
3036 PROFILE_FUNCTION_END;
3037 PetscFunctionReturn(0);
3038}
Create_InletConstantVelocity
PetscErrorCode Create_InletConstantVelocity(BoundaryCondition *bc)
Configures a BoundaryCondition object to behave as a constant velocity inlet.
Definition BC_Handlers.c:322
Create_InletProfileFromFile
PetscErrorCode Create_InletProfileFromFile(BoundaryCondition *bc)
Configures a BoundaryCondition object for a file-prescribed inlet profile.
Definition BC_Handlers.c:1292
Create_PeriodicGeometric
PetscErrorCode Create_PeriodicGeometric(BoundaryCondition *bc)
Configures a BoundaryCondition object for geometric periodic coupling.
Definition BC_Handlers.c:2159
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:715
Create_PeriodicDrivenInitial
PetscErrorCode Create_PeriodicDrivenInitial(BoundaryCondition *bc)
Configures a BoundaryCondition object for initial-flux periodic driving.
Definition BC_Handlers.c:2669
Create_PeriodicDrivenConstant
PetscErrorCode Create_PeriodicDrivenConstant(BoundaryCondition *bc)
Configures a BoundaryCondition object for periodic driven-flow forcing.
Definition BC_Handlers.c:2325
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:1644
ApplyPeriodicBCs
PetscErrorCode ApplyPeriodicBCs(UserCtx *user)
Internal helper implementation: ApplyPeriodicBCs().
Definition Boundaries.c:2169
PreparePeriodicQuickStencilFields
PetscErrorCode PreparePeriodicQuickStencilFields(UserCtx *user, Vec local_vector_field, Vec local_scalar_field)
Implementation of PreparePeriodicQuickStencilFields().
Definition Boundaries.c:2034
BoundarySystem_Initialize
PetscErrorCode BoundarySystem_Initialize(UserCtx *user, const char *bcs_filename)
Implementation of BoundarySystem_Initialize().
Definition Boundaries.c:850
TransferPeriodicFieldByDirection
static PetscErrorCode TransferPeriodicFieldByDirection(UserCtx *user, FieldId field_id, char direction)
Copies one cell field's wrapped local values onto the owned periodic duplicate plane.
Definition Boundaries.c:1537
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:400
ApplyWallFunction
PetscErrorCode ApplyWallFunction(UserCtx *user)
Internal helper implementation: ApplyWallFunction().
Definition Boundaries.c:2458
GetPersistentFaceField
static PetscErrorCode GetPersistentFaceField(UserCtx *user, FieldId field_id, char face_direction, DM *dm, Vec *global_vec, Vec *local_vec, PetscInt *dof)
Resolves one registered persistent single-face-family field.
Definition Boundaries.c:1687
ClassifyMomentumRow
MomentumRowType ClassifyMomentumRow(UserCtx *user, PetscInt i, PetscInt j, PetscInt k, PetscInt component, PetscInt *ri, PetscInt *rj, PetscInt *rk)
Implementation of ClassifyMomentumRow().
Definition Boundaries.c:597
PropagateBoundaryConfigToCoarserLevels
PetscErrorCode PropagateBoundaryConfigToCoarserLevels(SimCtx *simCtx)
Internal helper implementation: PropagateBoundaryConfigToCoarserLevels().
Definition Boundaries.c:947
GetPersistentStaggeredField
static PetscErrorCode GetPersistentStaggeredField(UserCtx *user, FieldId field_id, DM *dm, Vec *global_vec, Vec *local_vec)
Resolves one registered persistent component-staggered field.
Definition Boundaries.c:1917
ApplyMetricsPeriodicBCs
PetscErrorCode ApplyMetricsPeriodicBCs(UserCtx *user)
Internal helper implementation: ApplyMetricsPeriodicBCs().
Definition Boundaries.c:2140
EnforceRHSBoundaryConditions
PetscErrorCode EnforceRHSBoundaryConditions(UserCtx *user)
Implementation of EnforceRHSBoundaryConditions().
Definition Boundaries.c:660
TransferPeriodicStaggeredFieldByDirection
static PetscErrorCode TransferPeriodicStaggeredFieldByDirection(UserCtx *user, FieldId field_id, char periodic_direction)
Transfers one component-staggered field along one periodic axis.
Definition Boundaries.c:1940
UpdateDummyCells
PetscErrorCode UpdateDummyCells(UserCtx *user)
Internal helper implementation: UpdateDummyCells().
Definition Boundaries.c:2222
BoundarySystem_RefreshUbcs
PetscErrorCode BoundarySystem_RefreshUbcs(UserCtx *user)
Internal helper implementation: BoundarySystem_RefreshUbcs().
Definition Boundaries.c:1440
SynchronizePeriodicLocalStaggeredField
PetscErrorCode SynchronizePeriodicLocalStaggeredField(UserCtx *user, Vec local_field)
Implementation of SynchronizePeriodicLocalStaggeredField().
Definition Boundaries.c:2100
TranslatePeriodicFaceCenterGhosts
static PetscErrorCode TranslatePeriodicFaceCenterGhosts(UserCtx *user, Vec local_vec)
Applies geometric translations to wrapped face-center ghost coordinates.
Definition Boundaries.c:1734
BoundarySystem_Validate
PetscErrorCode BoundarySystem_Validate(UserCtx *user)
Internal helper implementation: BoundarySystem_Validate().
Definition Boundaries.c:789
BoundaryCondition_Create
PetscErrorCode BoundaryCondition_Create(BCHandlerType handler_type, BoundaryCondition **new_bc_ptr)
Internal helper implementation: BoundaryCondition_Create().
Definition Boundaries.c:703
IsFaceCenterCoordinateField
static PetscErrorCode IsFaceCenterCoordinateField(FieldId field_id, PetscBool *is_coordinate)
Returns whether a registered face field stores physical coordinates.
Definition Boundaries.c:1719
SynchronizePeriodicStaggeredFields
PetscErrorCode SynchronizePeriodicStaggeredFields(UserCtx *user, PetscInt num_fields, const FieldId field_ids[])
Implementation of SynchronizePeriodicStaggeredFields().
Definition Boundaries.c:1989
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:2958
FinalizePostProjectionCellFields
PetscErrorCode FinalizePostProjectionCellFields(UserCtx *user)
Implementation of FinalizePostProjectionCellFields().
Definition Boundaries.c:2918
TransferPeriodicFaceFieldByDirection
static PetscErrorCode TransferPeriodicFaceFieldByDirection(UserCtx *user, FieldId field_id, char face_direction, char periodic_direction)
Transfers one registered face-family field along one periodic axis.
Definition Boundaries.c:1781
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:127
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:213
BoundarySystem_ExecuteStep
PetscErrorCode BoundarySystem_ExecuteStep(UserCtx *user)
Implementation of BoundarySystem_ExecuteStep().
Definition Boundaries.c:1018
SynchronizePeriodicFaceFields
PetscErrorCode SynchronizePeriodicFaceFields(UserCtx *user, char face_direction, PetscInt num_fields, const FieldId field_ids[])
Synchronizes persistent fields belonging to one face family.
Definition Boundaries.c:1855
BoundarySystem_Destroy
PetscErrorCode BoundarySystem_Destroy(UserCtx *user)
Implementation of BoundarySystem_Destroy().
Definition Boundaries.c:1490
SynchronizePeriodicCellFields
PetscErrorCode SynchronizePeriodicCellFields(UserCtx *user, PetscInt num_fields, const FieldId field_ids[])
Implementation of SynchronizePeriodicCellFields().
Definition Boundaries.c:1621
UpdateCornerNodes
PetscErrorCode UpdateCornerNodes(UserCtx *user)
Internal helper implementation: UpdateCornerNodes().
Definition Boundaries.c:2314
MomentumRowType
MomentumRowType
Classification of one staggered momentum row (location + component).
Definition Boundaries.h:252
MOM_ROW_FIXED_HOMOGENEOUS
@ MOM_ROW_FIXED_HOMOGENEOUS
Dummy/tangential row carrying no unknown at all.
Definition Boundaries.h:255
MOM_ROW_PHYSICAL
@ MOM_ROW_PHYSICAL
Independent unknown governed by the momentum equation.
Definition Boundaries.h:253
MOM_ROW_PERIODIC_DUPLICATE
@ MOM_ROW_PERIODIC_DUPLICATE
Duplicate of a wrapped representative row (see ri, rj, rk).
Definition Boundaries.h:256
MOM_ROW_FIXED_CONDITIONED
@ MOM_ROW_FIXED_CONDITIONED
Strong Dirichlet row; the value comes from ApplyBoundaryConditions().
Definition Boundaries.h:254
FieldDescriptor::layout
FieldLayout layout
Definition field_catalog.h:126
FIELD_CAPABILITY_PERIODIC_GEOMETRY_SHIFT
@ FIELD_CAPABILITY_PERIODIC_GEOMETRY_SHIFT
Definition field_catalog.h:114
FIELD_CAPABILITY_PERIODIC_CELL_SYNC
@ FIELD_CAPABILITY_PERIODIC_CELL_SYNC
Definition field_catalog.h:111
FIELD_CAPABILITY_PERIODIC_FACE_SYNC
@ FIELD_CAPABILITY_PERIODIC_FACE_SYNC
Definition field_catalog.h:112
FIELD_CAPABILITY_PERIODIC_STAGGERED_SYNC
@ FIELD_CAPABILITY_PERIODIC_STAGGERED_SYNC
Definition field_catalog.h:113
FieldDescriptor::capabilities
unsigned int capabilities
Definition field_catalog.h:129
FieldView::descriptor
const FieldDescriptor * descriptor
Definition field_catalog.h:136
FieldGetView
PetscErrorCode FieldGetView(UserCtx *user, FieldId field_id, FieldView *view)
Resolve the existing DM and global/local vectors for one field.
Definition field_catalog.c:259
FieldView::local_vec
Vec local_vec
Definition field_catalog.h:139
FieldView::dm
DM dm
Definition field_catalog.h:137
FieldDescriptor::dof
PetscInt dof
Definition field_catalog.h:124
FieldView::global_vec
Vec global_vec
Definition field_catalog.h:138
FieldLayout
FieldLayout
Logical storage topology of a field.
Definition field_catalog.h:71
FIELD_LAYOUT_K_FACE
@ FIELD_LAYOUT_K_FACE
Definition field_catalog.h:76
FIELD_LAYOUT_I_FACE
@ FIELD_LAYOUT_I_FACE
Definition field_catalog.h:74
FIELD_LAYOUT_COMPONENT_STAGGERED
@ FIELD_LAYOUT_COMPONENT_STAGGERED
Definition field_catalog.h:77
FIELD_LAYOUT_J_FACE
@ FIELD_LAYOUT_J_FACE
Definition field_catalog.h:75
FieldDescriptor::canonical_name
const char * canonical_name
Definition field_catalog.h:121
FieldGetDescriptor
PetscErrorCode FieldGetDescriptor(FieldId field_id, const FieldDescriptor **descriptor)
Return immutable metadata for a valid field identifier.
Definition field_catalog.c:178
FieldId
FieldId
Compile-time identity for a catalogued Eulerian field.
Definition field_catalog.h:23
FIELD_ID_JETA
@ FIELD_ID_JETA
Definition field_catalog.h:49
FIELD_ID_CENTZ
@ FIELD_ID_CENTZ
Definition field_catalog.h:44
FIELD_ID_CSI
@ FIELD_ID_CSI
Definition field_catalog.h:35
FIELD_ID_IAJ
@ FIELD_ID_IAJ
Definition field_catalog.h:54
FIELD_ID_NVERT
@ FIELD_ID_NVERT
Definition field_catalog.h:38
FIELD_ID_UCAT
@ FIELD_ID_UCAT
Definition field_catalog.h:26
FIELD_ID_KETA
@ FIELD_ID_KETA
Definition field_catalog.h:52
FIELD_ID_JAJ
@ FIELD_ID_JAJ
Definition field_catalog.h:55
FIELD_ID_KAJ
@ FIELD_ID_KAJ
Definition field_catalog.h:56
FIELD_ID_AJ
@ FIELD_ID_AJ
Definition field_catalog.h:39
FIELD_ID_CENTY
@ FIELD_ID_CENTY
Definition field_catalog.h:43
FIELD_ID_KZET
@ FIELD_ID_KZET
Definition field_catalog.h:53
FIELD_ID_IETA
@ FIELD_ID_IETA
Definition field_catalog.h:46
FIELD_ID_UCONT
@ FIELD_ID_UCONT
Definition field_catalog.h:27
FIELD_ID_ICSI
@ FIELD_ID_ICSI
Definition field_catalog.h:45
FIELD_ID_ETA
@ FIELD_ID_ETA
Definition field_catalog.h:36
FIELD_ID_JCSI
@ FIELD_ID_JCSI
Definition field_catalog.h:48
FIELD_ID_JZET
@ FIELD_ID_JZET
Definition field_catalog.h:50
FIELD_ID_P
@ FIELD_ID_P
Definition field_catalog.h:30
FIELD_ID_IZET
@ FIELD_ID_IZET
Definition field_catalog.h:47
FIELD_ID_ZET
@ FIELD_ID_ZET
Definition field_catalog.h:37
FIELD_ID_KCSI
@ FIELD_ID_KCSI
Definition field_catalog.h:51
FIELD_ID_CENTX
@ FIELD_ID_CENTX
Definition field_catalog.h:42
FieldDescriptor
Immutable metadata for one field identity.
Definition field_catalog.h:119
FieldView
Non-owning runtime objects resolved for one field and UserCtx.
Definition field_catalog.h:135
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:837
FreeBC_ParamList
void FreeBC_ParamList(BC_Param *head)
Frees an entire linked list of boundary-condition parameters.
Definition io.c:664
BCHandlerTypeToString
const char * BCHandlerTypeToString(BCHandlerType handler_type)
Converts a BCHandlerType enum to its string representation.
Definition logging.c:793
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:253
LOCAL
#define LOCAL
Logging scope definitions for controlling message output.
Definition logging.h:45
GLOBAL
#define GLOBAL
Scope for global logging across all processes.
Definition logging.h:46
BCFaceToString
const char * BCFaceToString(BCFace face)
Returns the canonical log token for a boundary-face enum value.
Definition logging.c:671
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:200
PROFILE_FUNCTION_END
#define PROFILE_FUNCTION_END
Marks the end of a profiled code block.
Definition logging.h:859
BCTypeToString
const char * BCTypeToString(BCType type)
Returns the canonical log token for a boundary mathematical type.
Definition logging.c:773
LOG_ERROR
@ LOG_ERROR
Critical errors that may halt the program.
Definition logging.h:29
LOG_TRACE
@ LOG_TRACE
Very fine-grained tracing information for in-depth debugging.
Definition logging.h:33
LOG_INFO
@ LOG_INFO
Informational messages about program execution.
Definition logging.h:31
LOG_WARNING
@ LOG_WARNING
Non-critical issues that warrant attention.
Definition logging.h:30
LOG_DEBUG
@ LOG_DEBUG
Detailed debugging information.
Definition logging.h:32
LOG_VERBOSE
@ LOG_VERBOSE
Extremely detailed logs, typically for development use only.
Definition logging.h:34
PROFILE_FUNCTION_BEGIN
#define PROFILE_FUNCTION_BEGIN
Marks the beginning of a profiled code block (typically a function).
Definition logging.h:850
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:2285
Contra2Cart
PetscErrorCode Contra2Cart(UserCtx *user)
Reconstructs Cartesian velocity (Ucat) at cell centers from contravariant velocity (Ucont) defined on...
Definition setup.c:2649
UpdateLocalGhosts
PetscErrorCode UpdateLocalGhosts(UserCtx *user, FieldId field_id)
Updates the local vector (including ghost points) from its corresponding global vector.
Definition setup.c:1838
BoundaryCondition
The "virtual table" struct for a boundary condition handler object.
Definition variables.h:353
BoundaryCondition::PostStep
PetscErrorCode(* PostStep)(BoundaryCondition *self, BCContext *ctx, PetscReal *local_inflow, PetscReal *local_outflow)
Definition variables.h:360
BoundaryCondition::PreStep
PetscErrorCode(* PreStep)(BoundaryCondition *self, BCContext *ctx, PetscReal *local_inflow, PetscReal *local_outflow)
Definition variables.h:358
BoundaryCondition::type
BCHandlerType type
Definition variables.h:354
BoundaryCondition::Destroy
PetscErrorCode(* Destroy)(BoundaryCondition *self)
Definition variables.h:362
BoundaryCondition::Initialize
PetscErrorCode(* Initialize)(BoundaryCondition *self, BCContext *ctx)
Definition variables.h:357
BoundaryCondition::UpdateUbcs
PetscErrorCode(* UpdateUbcs)(BoundaryCondition *self, BCContext *ctx)
Definition variables.h:361
BoundaryCondition::Apply
PetscErrorCode(* Apply)(BoundaryCondition *self, BCContext *ctx)
Definition variables.h:359
BoundaryCondition::priority
BCPriorityType priority
Definition variables.h:355
UserCtx::lFriction_Velocity
Vec lFriction_Velocity
Definition variables.h:935
SimCtx::FarFluxInSum
PetscReal FarFluxInSum
Definition variables.h:799
PERIODIC
@ PERIODIC
Definition variables.h:292
WALL
@ WALL
Definition variables.h:286
MGCtx::user
UserCtx * user
Definition variables.h:571
SimCtx::FarFluxOutSum
PetscReal FarFluxOutSum
Definition variables.h:799
UserCtx::inletFaceDefined
PetscBool inletFaceDefined
Definition variables.h:932
UserCtx::Rhs
Vec Rhs
Definition variables.h:947
SimCtx::rank
PetscMPIInt rank
Definition variables.h:698
UserCtx::boundary_faces
BoundaryFaceConfig boundary_faces[6]
Definition variables.h:931
SimCtx::block_number
PetscInt block_number
Definition variables.h:790
UserCtx::identifiedInletBCFace
BCFace identifiedInletBCFace
Definition variables.h:933
UserCtx::lNvert
Vec lNvert
Definition variables.h:939
UserCtx::simCtx
SimCtx * simCtx
Back-pointer to the master simulation context.
Definition variables.h:909
SimCtx::FluxOutSum
PetscReal FluxOutSum
Definition variables.h:799
BC_Param_s::next
struct BC_Param_s * next
Definition variables.h:339
BC_Param_s::key
char * key
Definition variables.h:337
UserCtx::KM
PetscInt KM
Definition variables.h:920
UserCtx::lZet
Vec lZet
Definition variables.h:974
SimCtx::usermg
UserMG usermg
Definition variables.h:852
BCHandlerType
BCHandlerType
Defines the specific computational "strategy" for a boundary handler.
Definition variables.h:303
BC_HANDLER_PERIODIC_GEOMETRIC
@ BC_HANDLER_PERIODIC_GEOMETRIC
Definition variables.h:316
BC_HANDLER_INLET_PARABOLIC
@ BC_HANDLER_INLET_PARABOLIC
Definition variables.h:309
BC_HANDLER_INLET_CONSTANT_VELOCITY
@ BC_HANDLER_INLET_CONSTANT_VELOCITY
Definition variables.h:308
BC_HANDLER_PERIODIC_DRIVEN_INITIAL_FLUX
@ BC_HANDLER_PERIODIC_DRIVEN_INITIAL_FLUX
Definition variables.h:319
BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX
@ BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX
Definition variables.h:318
BC_HANDLER_INLET_PROFILE_FROM_FILE
@ BC_HANDLER_INLET_PROFILE_FROM_FILE
Definition variables.h:310
BC_HANDLER_WALL_NOSLIP
@ BC_HANDLER_WALL_NOSLIP
Definition variables.h:305
BC_HANDLER_OUTLET_CONSERVATION
@ BC_HANDLER_OUTLET_CONSERVATION
Definition variables.h:314
BoundaryFaceConfig::handler_type
BCHandlerType handler_type
Definition variables.h:369
UserCtx::_this
PetscInt _this
Definition variables.h:924
SimCtx::np
PetscInt np
Definition variables.h:827
BC_Param_s::value
char * value
Definition variables.h:338
UserCtx::Ucont
Vec Ucont
Definition variables.h:939
BCS::Ubcs
Vec Ubcs
Physical Cartesian velocity at boundary faces. Full 3D array but only boundary-face entries are meani...
Definition variables.h:123
Cmpnts::x
PetscScalar x
Definition variables.h:103
UserCtx::Bcs
BCS Bcs
Definition variables.h:934
BCContext::user
UserCtx * user
Definition variables.h:344
SimCtx::FluxInSum
PetscReal FluxInSum
Definition variables.h:799
UserCtx::lCsi
Vec lCsi
Definition variables.h:974
BoundaryFaceConfig::params
BC_Param * params
Definition variables.h:370
SimCtx::wall_roughness_height
PetscReal wall_roughness_height
Definition variables.h:786
Cmpnts::z
PetscScalar z
Definition variables.h:103
UserCtx::Ucat
Vec Ucat
Definition variables.h:939
UserCtx::JM
PetscInt JM
Definition variables.h:920
SimCtx::wallfunction
PetscInt wallfunction
Definition variables.h:822
UserMG::mglevels
PetscInt mglevels
Definition variables.h:578
UserCtx::lAj
Vec lAj
Definition variables.h:974
UserCtx::info
DMDALocalInfo info
Definition variables.h:918
BC_PRIORITY_OUTLET
@ BC_PRIORITY_OUTLET
Definition variables.h:328
BC_PRIORITY_FARFIELD
@ BC_PRIORITY_FARFIELD
Definition variables.h:326
BC_PRIORITY_WALL
@ BC_PRIORITY_WALL
Definition variables.h:327
BC_PRIORITY_INLET
@ BC_PRIORITY_INLET
Definition variables.h:325
Cmpnts::y
PetscScalar y
Definition variables.h:103
UserCtx::IM
PetscInt IM
Definition variables.h:920
UserCtx::periodic_translation
Cmpnts periodic_translation[3]
Definition variables.h:927
UserCtx::lEta
Vec lEta
Definition variables.h:974
UserMG::mgctx
MGCtx * mgctx
Definition variables.h:581
UserCtx::periodic_translation_valid
PetscBool periodic_translation_valid[3]
Definition variables.h:928
BoundaryFaceConfig::mathematical_type
BCType mathematical_type
Definition variables.h:368
BCFace
BCFace
Identifies the six logical faces of a structured computational block.
Definition variables.h:261
BC_FACE_NEG_X
@ BC_FACE_NEG_X
Definition variables.h:262
BC_FACE_POS_Z
@ BC_FACE_POS_Z
Definition variables.h:264
BC_FACE_POS_Y
@ BC_FACE_POS_Y
Definition variables.h:263
BC_FACE_NEG_Z
@ BC_FACE_NEG_Z
Definition variables.h:264
BC_FACE_POS_X
@ BC_FACE_POS_X
Definition variables.h:262
BC_FACE_NEG_Y
@ BC_FACE_NEG_Y
Definition variables.h:263
BoundaryFaceConfig::handler
BoundaryCondition * handler
Definition variables.h:371
BCContext
Provides execution context for a boundary condition handler.
Definition variables.h:343
BC_Param_s
A node in a linked list for storing key-value parameters from the bcs.dat file.
Definition variables.h:336
BoundaryFaceConfig
Holds the complete configuration for one of the six boundary faces.
Definition variables.h:366
Cmpnts
A 3D point or vector with PetscScalar components.
Definition variables.h:102
SimCtx
The master context for the entire simulation.
Definition variables.h:695
UserCtx
User-defined context containing data specific to a single computational grid level.
Definition variables.h:906
UserMG
User-level context for managing the entire multigrid hierarchy.
Definition variables.h:577
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:706
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