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