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