PICurv 0.1.0
A Parallel Particle-In-Cell Solver for Curvilinear LES
Loading...
Searching...
No Matches
logging.c
Go to the documentation of this file.
1// logging.c
2#include "logging.h"
3#include "AnalyticalSolutions.h"
4#include "verification_sources.h"
5
6/* Maximum temporary buffer size for converting numbers to strings */
7#define TMP_BUF_SIZE 128
8
9// --------------------- Static Variable for Log Level ---------------------
10
11/**
12 * @brief Static variable to cache the current logging level.
13 *
14 * Initialized to -1 to indicate that the log level has not been set yet.
15 */
16static LogLevel current_log_level = -1;
17
18// --------------------- Static Variables for Allow-List -------------------
19
20/**
21 * @brief Global/static array of function names allowed to log.
22 */
23static char** gAllowedFunctions = NULL;
24
25/**
26 * @brief Number of entries in the gAllowedFunctions array.
27 */
28static int gNumAllowed = 0;
29
30enum {
31 SEARCH_METRIC_SUM_SEARCH_ATTEMPTS = 0,
32 SEARCH_METRIC_SUM_SEARCH_POPULATION,
33 SEARCH_METRIC_SUM_SEARCH_LOCATED,
34 SEARCH_METRIC_SUM_SEARCH_LOST,
35 SEARCH_METRIC_SUM_TRAVERSAL_STEPS,
36 SEARCH_METRIC_SUM_RESEARCH,
37 SEARCH_METRIC_SUM_TIE_BREAKS,
38 SEARCH_METRIC_SUM_BOUNDARY_CLAMPS,
39 SEARCH_METRIC_SUM_BBOX_GUESS_SUCCESS,
40 SEARCH_METRIC_SUM_BBOX_GUESS_FALLBACK,
41 SEARCH_METRIC_SUM_MAX_TRAVERSAL_FAILS,
42 SEARCH_METRIC_MAX_TRAVERSAL_STEPS,
43 SEARCH_METRIC_MAX_PASS_DEPTH,
44 SEARCH_METRIC_REDUCTION_LEN
45};
46
47/**
48 * @brief Internal reduction callback for packed search metrics.
49 * @details Local to this translation unit.
50 */
51static void SearchMetricsReduceOp(void *invec, void *inoutvec, int *len, MPI_Datatype *datatype)
52{
53 PetscReal *in = (PetscReal *)invec;
54 PetscReal *inout = (PetscReal *)inoutvec;
55 (void)datatype;
56
57 for (int idx = 0; idx < *len; idx += SEARCH_METRIC_REDUCTION_LEN) {
58 inout[idx + SEARCH_METRIC_SUM_SEARCH_ATTEMPTS] += in[idx + SEARCH_METRIC_SUM_SEARCH_ATTEMPTS];
59 inout[idx + SEARCH_METRIC_SUM_SEARCH_POPULATION] += in[idx + SEARCH_METRIC_SUM_SEARCH_POPULATION];
60 inout[idx + SEARCH_METRIC_SUM_SEARCH_LOCATED] += in[idx + SEARCH_METRIC_SUM_SEARCH_LOCATED];
61 inout[idx + SEARCH_METRIC_SUM_SEARCH_LOST] += in[idx + SEARCH_METRIC_SUM_SEARCH_LOST];
62 inout[idx + SEARCH_METRIC_SUM_TRAVERSAL_STEPS] += in[idx + SEARCH_METRIC_SUM_TRAVERSAL_STEPS];
63 inout[idx + SEARCH_METRIC_SUM_RESEARCH] += in[idx + SEARCH_METRIC_SUM_RESEARCH];
64 inout[idx + SEARCH_METRIC_SUM_TIE_BREAKS] += in[idx + SEARCH_METRIC_SUM_TIE_BREAKS];
65 inout[idx + SEARCH_METRIC_SUM_BOUNDARY_CLAMPS] += in[idx + SEARCH_METRIC_SUM_BOUNDARY_CLAMPS];
66 inout[idx + SEARCH_METRIC_SUM_BBOX_GUESS_SUCCESS] += in[idx + SEARCH_METRIC_SUM_BBOX_GUESS_SUCCESS];
67 inout[idx + SEARCH_METRIC_SUM_BBOX_GUESS_FALLBACK] += in[idx + SEARCH_METRIC_SUM_BBOX_GUESS_FALLBACK];
68 inout[idx + SEARCH_METRIC_SUM_MAX_TRAVERSAL_FAILS] += in[idx + SEARCH_METRIC_SUM_MAX_TRAVERSAL_FAILS];
69 inout[idx + SEARCH_METRIC_MAX_TRAVERSAL_STEPS] = PetscMax(inout[idx + SEARCH_METRIC_MAX_TRAVERSAL_STEPS],
70 in[idx + SEARCH_METRIC_MAX_TRAVERSAL_STEPS]);
71 inout[idx + SEARCH_METRIC_MAX_PASS_DEPTH] = PetscMax(inout[idx + SEARCH_METRIC_MAX_PASS_DEPTH],
72 in[idx + SEARCH_METRIC_MAX_PASS_DEPTH]);
73 }
74}
75
76// --------------------- Function Implementations ---------------------
77
78/**
79 * @brief Implementation of \ref get_log_level().
80 * @details Full API contract (arguments, ownership, side effects) is documented with
81 * the header declaration in `include/logging.h`.
82 * @see get_log_level()
83 */
84LogLevel get_log_level() {
85 if (current_log_level == -1) { // Log level not set yet
86 const char *env = getenv("LOG_LEVEL");
87 if (!env) {
88 current_log_level = LOG_ERROR; // Default level
89 }
90 else if (strcmp(env, "DEBUG") == 0) {
91 current_log_level = LOG_DEBUG;
92 }
93 else if (strcmp(env, "INFO") == 0) {
94 current_log_level = LOG_INFO;
95 }
96 else if (strcmp(env, "WARNING") == 0) {
97 current_log_level = LOG_WARNING;
98 }
99 else if (strcmp(env, "VERBOSE") == 0) {
100 current_log_level = LOG_VERBOSE;
101 }
102 else if (strcmp(env, "TRACE") == 0) {
103 current_log_level = LOG_TRACE;
104 }
105 else {
106 current_log_level = LOG_ERROR; // Default if unrecognized
107 }
108 }
109 return current_log_level;
110}
111
112/**
113 * @brief Internal helper implementation: `print_log_level()`.
114 * @details Local to this translation unit.
115 */
116PetscErrorCode print_log_level(void)
117{
118 PetscMPIInt rank;
119 PetscErrorCode ierr;
120 int level;
121 const char *level_name;
122
123 PetscFunctionBeginUser;
124 /* get MPI rank */
125 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRMPI(ierr);
126
127 /* decide level name */
128 level = get_log_level();
129 level_name = (level == LOG_ERROR) ? "ERROR" :
130 (level == LOG_WARNING) ? "WARNING" :
131 (level == LOG_INFO) ? "INFO" :
132 (level == LOG_DEBUG) ? "DEBUG" :
133 (level == LOG_VERBOSE) ? "VERBOSE" :
134 (level == LOG_TRACE) ? "TRACE" :
135 "UNKNOWN";
136
137 /* print it out */
138 ierr = PetscPrintf(PETSC_COMM_SELF,
139 "Current log level: %s (%d) | rank: %d\n",
140 level_name, level, (int)rank);
141 CHKERRMPI(ierr);
142
143 PetscFunctionReturn(PETSC_SUCCESS);
144}
145
146/**
147 * @brief Implementation of \ref set_allowed_functions().
148 * @details Full API contract (arguments, ownership, side effects) is documented with
149 * the header declaration in `include/logging.h`.
150 * @see set_allowed_functions()
151 */
152void set_allowed_functions(const char** functionList, int count)
153{
154 // 1. Free any existing entries
155 if (gAllowedFunctions) {
156 for (int i = 0; i < gNumAllowed; ++i) {
157 free(gAllowedFunctions[i]); // each was strdup'ed
158 }
159 free(gAllowedFunctions);
160 gAllowedFunctions = NULL;
161 gNumAllowed = 0;
162 }
163
164 // 2. Allocate new array
165 if (count > 0) {
166 gAllowedFunctions = (char**)malloc(sizeof(char*) * count);
167 }
168
169 // 3. Copy the new entries
170 for (int i = 0; i < count; ++i) {
171 // strdup is a POSIX function. If not available, implement your own string copy.
172 gAllowedFunctions[i] = strdup(functionList[i]);
173 }
174 gNumAllowed = count;
175}
176
177/**
178 * @brief Implementation of \ref is_function_allowed().
179 * @details Full API contract (arguments, ownership, side effects) is documented with
180 * the header declaration in `include/logging.h`.
181 * @see is_function_allowed()
182 */
183PetscBool is_function_allowed(const char* functionName)
184{
185 /* no list ⇒ allow all */
186 if (gNumAllowed == 0) {
187 return PETSC_TRUE;
188 }
189
190 /* otherwise only the listed functions are allowed */
191 for (int i = 0; i < gNumAllowed; ++i) {
192 if (strcmp(gAllowedFunctions[i], functionName) == 0) {
193 return PETSC_TRUE;
194 }
195 }
196 return PETSC_FALSE;
197}
198
199/**
200 * @brief Implementation of \ref LOG_CELL_VERTICES().
201 * @details Full API contract (arguments, ownership, side effects) is documented with
202 * the header declaration in `include/logging.h`.
203 * @see LOG_CELL_VERTICES()
204 */
205PetscErrorCode LOG_CELL_VERTICES(const Cell *cell, PetscMPIInt rank)
206{
207
208 // Validate input pointers
209 if (cell == NULL) {
210 LOG_ALLOW(LOCAL,LOG_ERROR, "'cell' is NULL.\n");
211 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "LOG_CELL_VERTICES - Input parameter 'cell' is NULL.");
212 }
213
214 LOG_ALLOW(LOCAL,LOG_VERBOSE, "Rank %d, Cell Vertices:\n", rank);
215 for(int i = 0; i < 8; i++){
216 LOG_ALLOW(LOCAL,LOG_VERBOSE, " Vertex[%d]: (%.2f, %.2f, %.2f)\n",
217 i, cell->vertices[i].x, cell->vertices[i].y, cell->vertices[i].z);
218 }
219
220 return 0; // Indicate successful execution
221}
222
223
224/**
225 * @brief Implementation of \ref LOG_FACE_DISTANCES().
226 * @details Full API contract (arguments, ownership, side effects) is documented with
227 * the header declaration in `include/logging.h`.
228 * @see LOG_FACE_DISTANCES()
229 */
230PetscErrorCode LOG_FACE_DISTANCES(PetscReal* d)
231{
232
233 // Validate input array
234 if (d == NULL) {
235 LOG_ALLOW(LOCAL,LOG_ERROR, " 'd' is NULL.\n");
236 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, " Input array 'd' is NULL.");
237 }
238
239 PetscPrintf(PETSC_COMM_SELF, " Face Distances:\n");
240 PetscPrintf(PETSC_COMM_SELF, " LEFT(%d): %.15f\n", LEFT, d[LEFT]);
241 PetscPrintf(PETSC_COMM_SELF, " RIGHT(%d): %.15f\n", RIGHT, d[RIGHT]);
242 PetscPrintf(PETSC_COMM_SELF, " BOTTOM(%d): %.15f\n", BOTTOM, d[BOTTOM]);
243 PetscPrintf(PETSC_COMM_SELF, " TOP(%d): %.15f\n", TOP, d[TOP]);
244 PetscPrintf(PETSC_COMM_SELF, " FRONT(%d): %.15f\n", FRONT, d[FRONT]);
245 PetscPrintf(PETSC_COMM_SELF, " BACK(%d): %.15f\n", BACK, d[BACK]);
246
247 return 0; // Indicate successful execution
248}
249
250/*
251 * Helper function: Converts an integer (of type int) to a string.
252 */
253static void IntToStr(int value, char *buf, size_t bufsize)
254{
255 snprintf(buf, bufsize, "%d", value);
256}
257
258/*
259 * Helper function: Converts a 64‐bit integer to a string.
260 */
261static void Int64ToStr(PetscInt64 value, char *buf, size_t bufsize)
262{
263 snprintf(buf, bufsize, "%ld", value);
264}
265
266/*
267 * Helper function: Converts three integers into a formatted string "(i, j, k)".
268 */
269static void CellToStr(const PetscInt *cell, char *buf, size_t bufsize)
270{
271 snprintf(buf, bufsize, "(%d, %d, %d)", cell[0], cell[1], cell[2]);
272}
273
274/*
275 * Helper function: Converts three PetscReal values into a formatted string "(x, y, z)".
276 */
277static void TripleRealToStr(const PetscReal *arr, char *buf, size_t bufsize)
278{
279 snprintf(buf, bufsize, "(%.4f, %.4f, %.4f)", arr[0], arr[1], arr[2]);
280}
281
282/*
283 * Helper function: Computes the maximum string length for each column (across all particles).
284 *
285 * The function examines every particle (from 0 to nParticles-1) and converts the value to a
286 * string using the helper functions above. The maximum length is stored in the pointers provided.
287 *
288 * @param nParticles Number of particles.
289 * @param ranks Array of particle MPI ranks.
290 * @param pids Array of particle IDs.
291 * @param cellIDs Array of cell IDs (stored consecutively, 3 per particle).
292 * @param positions Array of positions (3 per particle).
293 * @param velocities Array of velocities (3 per particle).
294 * @param weights Array of weights (3 per particle).
295 * @param wRank [out] Maximum width for Rank column.
296 * @param wPID [out] Maximum width for PID column.
297 * @param wCell [out] Maximum width for Cell column.
298 * @param wPos [out] Maximum width for Position column.
299 * @param wVel [out] Maximum width for Velocity column.
300 * @param wWt [out] Maximum width for Weights column.
301 */
302static PetscErrorCode ComputeMaxColumnWidths(PetscInt nParticles,
303 const PetscMPIInt *ranks,
304 const PetscInt64 *pids,
305 const PetscInt *cellIDs,
306 const PetscReal *positions,
307 const PetscReal *velocities,
308 const PetscReal *weights,
309 int *wRank, int *wPID, int *wCell,
310 int *wPos, int *wVel, int *wWt)
311{
312 char tmp[TMP_BUF_SIZE];
313
314 *wRank = strlen("Rank"); /* Start with the header label lengths */
315 *wPID = strlen("PID");
316 *wCell = strlen("Cell (i,j,k)");
317 *wPos = strlen("Position (x,y,z)");
318 *wVel = strlen("Velocity (x,y,z)");
319 *wWt = strlen("Weights (a1,a2,a3)");
320
321 for (PetscInt i = 0; i < nParticles; i++) {
322 /* Rank */
323 IntToStr(ranks[i], tmp, TMP_BUF_SIZE);
324 if ((int)strlen(tmp) > *wRank) *wRank = (int)strlen(tmp);
325
326 /* PID */
327 Int64ToStr(pids[i], tmp, TMP_BUF_SIZE);
328 if ((int)strlen(tmp) > *wPID) *wPID = (int)strlen(tmp);
329
330 /* Cell: use the three consecutive values */
331 CellToStr(&cellIDs[3 * i], tmp, TMP_BUF_SIZE);
332 if ((int)strlen(tmp) > *wCell) *wCell = (int)strlen(tmp);
333
334 /* Position */
335 TripleRealToStr(&positions[3 * i], tmp, TMP_BUF_SIZE);
336 if ((int)strlen(tmp) > *wPos) *wPos = (int)strlen(tmp);
337
338 /* Velocity */
339 TripleRealToStr(&velocities[3 * i], tmp, TMP_BUF_SIZE);
340 if ((int)strlen(tmp) > *wVel) *wVel = (int)strlen(tmp);
341
342 /* Weights */
343 TripleRealToStr(&weights[3 * i], tmp, TMP_BUF_SIZE);
344 if ((int)strlen(tmp) > *wWt) *wWt = (int)strlen(tmp);
345 }
346 return 0;
347}
348
349/*
350 * Helper function: Builds a format string for a table row.
351 *
352 * The format string will include proper width specifiers for each column.
353 * For example, it might create something like:
354 *
355 * "| %-6s | %-8s | %-20s | %-25s | %-25s | %-25s |\n"
356 *
357 * @param wRank Maximum width for the Rank column.
358 * @param wPID Maximum width for the PID column.
359 * @param wCell Maximum width for the Cell column.
360 * @param wPos Maximum width for the Position column.
361 * @param wVel Maximum width for the Velocity column.
362 * @param wWt Maximum width for the Weights column.
363 * @param fmtStr Buffer in which to build the format string.
364 * @param bufSize Size of fmtStr.
365 */
366static void BuildRowFormatString(PetscMPIInt wRank, PetscInt wPID, PetscInt wCell, PetscInt wPos, PetscInt wVel, PetscInt wWt, char *fmtStr, size_t bufSize)
367{
368 // Build a format string using snprintf.
369 // We assume that the Rank is an int (%d), PID is a 64-bit int (%ld)
370 // and the remaining columns are strings (which have been formatted already).
371 snprintf(fmtStr, bufSize,
372 "| %%-%dd | %%-%dd | %%-%ds | %%-%ds | %%-%ds | %%-%ds |\n",
373 wRank, wPID, wCell, wPos, wVel, wWt);
374}
375
376/*
377 * Helper function: Builds a header string for the table using column titles.
378 */
379static void BuildHeaderString(char *headerStr, size_t bufSize, PetscMPIInt wRank, PetscInt wPID, PetscInt wCell, PetscInt wPos, PetscInt wVel, PetscInt wWt)
380{
381 snprintf(headerStr, bufSize,
382 "| %-*s | %-*s | %-*s | %-*s | %-*s | %-*s |\n",
383 (int)wRank, "Rank",
384 (int)wPID, "PID",
385 (int)wCell, "Cell (i,j,k)",
386 (int)wPos, "Position (x,y,z)",
387 (int)wVel, "Velocity (x,y,z)",
388 (int)wWt, "Weights (a1,a2,a3)");
389}
390
391/**
392 * @brief Implementation of \ref LOG_PARTICLE_FIELDS().
393 * @details Full API contract (arguments, ownership, side effects) is documented with
394 * the header declaration in `include/logging.h`.
395 * @see LOG_PARTICLE_FIELDS()
396 */
397PetscErrorCode LOG_PARTICLE_FIELDS(UserCtx* user, PetscInt printInterval)
398{
399 DM swarm = user->swarm;
400 PetscErrorCode ierr;
401 PetscInt localNumParticles;
402 PetscReal *positions = NULL;
403 PetscInt64 *particleIDs = NULL;
404 PetscMPIInt *particleRanks = NULL;
405 PetscInt *cellIDs = NULL;
406 PetscReal *weights = NULL;
407 PetscReal *velocities = NULL;
408 PetscMPIInt rank;
409
410 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
411 LOG_ALLOW(LOCAL,LOG_INFO, "Rank %d is retrieving particle data.\n", rank);
412
413 ierr = DMSwarmGetLocalSize(swarm, &localNumParticles); CHKERRQ(ierr);
414 LOG_ALLOW(LOCAL,LOG_DEBUG,"Rank %d has %d particles.\n", rank, localNumParticles);
415
416 ierr = DMSwarmGetField(swarm, "position", NULL, NULL, (void**)&positions); CHKERRQ(ierr);
417 ierr = DMSwarmGetField(swarm, "DMSwarm_pid", NULL, NULL, (void**)&particleIDs); CHKERRQ(ierr);
418 ierr = DMSwarmGetField(swarm, "DMSwarm_rank", NULL, NULL, (void**)&particleRanks); CHKERRQ(ierr);
419 ierr = DMSwarmGetField(swarm, "DMSwarm_CellID", NULL, NULL, (void**)&cellIDs); CHKERRQ(ierr);
420 ierr = DMSwarmGetField(swarm, "weight", NULL, NULL, (void**)&weights); CHKERRQ(ierr);
421 ierr = DMSwarmGetField(swarm, "velocity", NULL, NULL, (void**)&velocities); CHKERRQ(ierr);
422
423 /* Compute maximum column widths. */
424 int wRank, wPID, wCell, wPos, wVel, wWt;
425 wRank = wPID = wCell = wPos = wVel = wWt = 0;
426 ierr = ComputeMaxColumnWidths(localNumParticles, particleRanks, particleIDs, cellIDs,
427 positions, velocities, weights,
428 &wRank, &wPID, &wCell, &wPos, &wVel, &wWt); CHKERRQ(ierr);
429
430 /* Build a header string and a row format string. */
431 char headerFmt[256];
432 char rowFmt[256];
433 BuildHeaderString(headerFmt, sizeof(headerFmt), wRank, wPID, wCell, wPos, wVel, wWt);
434 BuildRowFormatString(wRank, wPID, wCell, wPos, wVel, wWt, rowFmt, sizeof(rowFmt));
435
436 /* Print header (using synchronized printing for parallel output). */
437 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "--------------------------------------------------------------------------------------------------------------\n"); CHKERRQ(ierr);
438 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%s", headerFmt); CHKERRQ(ierr);
439 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "--------------------------------------------------------------------------------------------------------------\n"); CHKERRQ(ierr);
440
441 /* Loop over particles and print every printInterval-th row. */
442 char rowStr[256];
443 for (PetscInt i = 0; i < localNumParticles; i++) {
444 if (i % printInterval == 0) {
445 // ------- DEBUG
446 //char cellStr[TMP_BUF_SIZE], posStr[TMP_BUF_SIZE], velStr[TMP_BUF_SIZE], wtStr[TMP_BUF_SIZE];
447 //CellToStr(&cellIDs[3*i], cellStr, TMP_BUF_SIZE);
448 //TripleRealToStr(&positions[3*i], posStr, TMP_BUF_SIZE);
449 //TripleRealToStr(&velocities[3*i], velStr, TMP_BUF_SIZE);
450 // TripleRealToStr(&weights[3*i], wtStr, TMP_BUF_SIZE);
451
452 // if (rank == 0) { // Or whatever rank is Rank 0
453 //PetscPrintf(PETSC_COMM_SELF, "[Rank 0 DEBUG LPF] Particle %lld: PID=%lld, Rank=%d\n", (long long)i, (long long)particleIDs[i], particleRanks[i]);
454 //PetscPrintf(PETSC_COMM_SELF, "[Rank 0 DEBUG LPF] Raw Pos: (%.10e, %.10e, %.10e)\n", positions[3*i+0], positions[3*i+1], positions[3*i+2]);
455 //PetscPrintf(PETSC_COMM_SELF, "[Rank 0 DEBUG LPF] Str Pos: %s\n", posStr);
456 //PetscPrintf(PETSC_COMM_SELF, "[Rank 0 DEBUG LPF] Raw Vel: (%.10e, %.10e, %.10e)\n", velocities[3*i+0], velocities[3*i+1], velocities[3*i+2]);
457 // PetscPrintf(PETSC_COMM_SELF, "[Rank 0 DEBUG LPF] Str Vel: %s\n", velStr);
458 // Add similar for cell, weights
459 // PetscPrintf(PETSC_COMM_SELF, "[Rank 0 DEBUG LPF] About to build rowStr for particle %lld\n", (long long)i);
460 // fflush(stdout);
461 // }
462
463 // snprintf(rowStr, sizeof(rowStr), rowFmt,
464 // particleRanks[i],
465 // particleIDs[i],
466 // cellStr,
467 // posStr,
468 // velStr,
469 // wtStr);
470
471
472 // ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%s", rowStr); CHKERRQ(ierr);
473
474 // ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%s", rowStr); CHKERRQ(ierr);
475
476 // -------- DEBUG
477 /* Format the row by converting each field to a string first.
478 * We use temporary buffers and then build the row string.
479 */
480
481 char cellStr[TMP_BUF_SIZE], posStr[TMP_BUF_SIZE], velStr[TMP_BUF_SIZE], wtStr[TMP_BUF_SIZE];
482 CellToStr(&cellIDs[3*i], cellStr, TMP_BUF_SIZE);
483 TripleRealToStr(&positions[3*i], posStr, TMP_BUF_SIZE);
484 TripleRealToStr(&velocities[3*i], velStr, TMP_BUF_SIZE);
485 TripleRealToStr(&weights[3*i], wtStr, TMP_BUF_SIZE);
486
487 /* Build the row string. Note that for the integer fields we can use the row format string. */
488 snprintf(rowStr, sizeof(rowStr), rowFmt,
489 particleRanks[i],
490 particleIDs[i],
491 cellStr,
492 posStr,
493 velStr,
494 wtStr);
495 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%s", rowStr); CHKERRQ(ierr);
496 }
497 }
498
499
500 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "--------------------------------------------------------------------------------------------------------------\n"); CHKERRQ(ierr);
501 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\n"); CHKERRQ(ierr);
502 ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT); CHKERRQ(ierr);
503
504 LOG_ALLOW_SYNC(GLOBAL,LOG_DEBUG,"Completed printing on Rank %d.\n", rank);
505
506 /* Restore fields */
507 ierr = DMSwarmRestoreField(swarm, "position", NULL, NULL, (void**)&positions); CHKERRQ(ierr);
508 ierr = DMSwarmRestoreField(swarm, "DMSwarm_pid", NULL, NULL, (void**)&particleIDs); CHKERRQ(ierr);
509 ierr = DMSwarmRestoreField(swarm, "DMSwarm_rank", NULL, NULL, (void**)&particleRanks); CHKERRQ(ierr);
510 ierr = DMSwarmRestoreField(swarm, "DMSwarm_CellID", NULL, NULL, (void**)&cellIDs); CHKERRQ(ierr);
511 ierr = DMSwarmRestoreField(swarm, "weight", NULL, NULL, (void**)&weights); CHKERRQ(ierr);
512 ierr = DMSwarmRestoreField(swarm, "velocity", NULL, NULL, (void**)&velocities); CHKERRQ(ierr);
513
514 LOG_ALLOW(LOCAL,LOG_DEBUG, "Restored all particle fields.\n");
515 return 0;
516}
517
518/**
519 * @brief Implementation of \ref IsParticleConsoleSnapshotEnabled().
520 * @details Full API contract (arguments, ownership, side effects) is documented with
521 * the header declaration in `include/logging.h`.
522 * @see IsParticleConsoleSnapshotEnabled()
523 */
524
525PetscBool IsParticleConsoleSnapshotEnabled(const SimCtx *simCtx)
526{
527 if (!simCtx) {
528 return PETSC_FALSE;
529 }
530 return (PetscBool)(simCtx->np > 0 &&
531 simCtx->particleConsoleOutputFreq > 0 &&
532 get_log_level() >= LOG_INFO);
533}
534
535/**
536 * @brief Implementation of \ref ShouldEmitPeriodicParticleConsoleSnapshot().
537 * @details Full API contract (arguments, ownership, side effects) is documented with
538 * the header declaration in `include/logging.h`.
539 * @see ShouldEmitPeriodicParticleConsoleSnapshot()
540 */
541
542PetscBool ShouldEmitPeriodicParticleConsoleSnapshot(const SimCtx *simCtx, PetscInt completed_step)
543{
544 return (PetscBool)(IsParticleConsoleSnapshotEnabled(simCtx) &&
545 completed_step > 0 &&
546 completed_step % simCtx->particleConsoleOutputFreq == 0);
547}
548
549/**
550 * @brief Implementation of \ref EmitParticleConsoleSnapshot().
551 * @details Full API contract (arguments, ownership, side effects) is documented with
552 * the header declaration in `include/logging.h`.
553 * @see EmitParticleConsoleSnapshot()
554 */
555
556PetscErrorCode EmitParticleConsoleSnapshot(UserCtx *user, SimCtx *simCtx, PetscInt step)
557{
558 PetscErrorCode ierr;
559
560 PetscFunctionBeginUser;
561 LOG(GLOBAL, LOG_INFO, "Particle states at step %d:\n", step);
562 ierr = LOG_PARTICLE_FIELDS(user, simCtx->LoggingFrequency); CHKERRQ(ierr);
563 PetscFunctionReturn(0);
564}
565
566
567/**
568 * @brief Internal helper implementation: `trim()`.
569 * @details Local to this translation unit.
570 */
571static void trim(char *s)
572{
573 if (!s) return;
574
575 /* ---- 1. strip leading blanks ----------------------------------- */
576 char *p = s;
577 while (*p && isspace((unsigned char)*p))
578 ++p;
579
580 if (p != s) /* move the trimmed text forward */
581 memmove(s, p, strlen(p) + 1); /* +1 to copy the final NUL */
582
583 /* ---- 2. strip trailing blanks ---------------------------------- */
584 size_t len = strlen(s);
585 while (len > 0 && isspace((unsigned char)s[len - 1]))
586 s[--len] = '\0';
587}
588
589/* ------------------------------------------------------------------------- */
590/**
591 * @brief Implementation of \ref LoadAllowedFunctionsFromFile().
592 * @details Full API contract (arguments, ownership, side effects) is documented with
593 * the header declaration in `include/logging.h`.
594 * @see LoadAllowedFunctionsFromFile()
595 */
596PetscErrorCode LoadAllowedFunctionsFromFile(const char filename[],
597 char ***funcsOut,
598 PetscInt *nOut)
599{
600 FILE *fp = NULL;
601 char **funcs = NULL;
602 size_t cap = 16; /* initial capacity */
603 size_t n = 0; /* number of names */
604 char line[PETSC_MAX_PATH_LEN];
605 PetscErrorCode ierr;
606
607 PetscFunctionBegin;
608
609 /* ---------------------------------------------------------------------- */
610 /* 1. Open file */
611 fp = fopen(filename, "r");
612 if (!fp) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN,
613 "Cannot open %s", filename);
614
615 /* 2. Allocate initial pointer array */
616 ierr = PetscMalloc1(cap, &funcs); CHKERRQ(ierr);
617
618 /* 3. Read file line by line */
619 while (fgets(line, sizeof line, fp)) {
620 /* Strip everything after a comment character '#'. */
621 char *hash = strchr(line, '#');
622 if (hash) *hash = '\0';
623
624 trim(line); /* remove leading/trailing blanks */
625 if (!*line) continue; /* skip if empty */
626
627 /* Grow the array if necessary */
628 if (n == cap) {
629 cap *= 2;
630 ierr = PetscRealloc(cap * sizeof(*funcs), (void **)&funcs); CHKERRQ(ierr);
631 }
632
633 /* Deep‑copy the cleaned identifier */
634 ierr = PetscStrallocpy(line, &funcs[n++]); CHKERRQ(ierr);
635 }
636 fclose(fp);
637
638 /* 4. Return results to caller */
639 *funcsOut = funcs;
640 *nOut = (PetscInt)n;
641
642 PetscFunctionReturn(0);
643}
644
645/* ------------------------------------------------------------------------- */
646/**
647 * @brief Internal helper implementation: `FreeAllowedFunctions()`.
648 * @details Local to this translation unit.
649 */
650PetscErrorCode FreeAllowedFunctions(char **funcs, PetscInt n)
651{
652 PetscErrorCode ierr;
653 PetscFunctionBegin;
654 if (funcs) {
655 for (PetscInt i = 0; i < n; ++i) {
656 ierr = PetscFree(funcs[i]); CHKERRQ(ierr);
657 }
658 ierr = PetscFree(funcs); CHKERRQ(ierr);
659 }
660 PetscFunctionReturn(0);
661}
662
663/**
664 * @brief Implementation of \ref BCFaceToString().
665 * @details Full API contract (arguments, ownership, side effects) is documented with
666 * the header declaration in `include/logging.h`.
667 * @see BCFaceToString()
668 */
669const char* BCFaceToString(BCFace face) {
670 switch (face) {
671 case BC_FACE_NEG_X: return "-Xi (I-Min)";
672 case BC_FACE_POS_X: return "+Xi (I-Max)";
673 case BC_FACE_NEG_Y: return "-Eta (J-Min)";
674 case BC_FACE_POS_Y: return "+Eta (J-Max)";
675 case BC_FACE_NEG_Z: return "-Zeta (K-Min)";
676 case BC_FACE_POS_Z: return "+Zeta (K-Max)";
677 default: return "Unknown Face";
678 }
679}
680
681/**
682 * @brief Implementation of \ref InitialConditionModeToString().
683 * @details Full API contract (arguments, ownership, side effects) is documented with
684 * the header declaration in `include/logging.h`.
685 * @see InitialConditionModeToString()
686 */
687const char* InitialConditionModeToString(InitialConditionMode mode)
688{
689 switch(mode){
690 case IC_MODE_ZERO: return "Zero";
691 case IC_MODE_CONSTANT_CARTESIAN: return "Cartesian Constant";
692 case IC_MODE_POISEUILLE: return "Poiseuille";
693 case IC_MODE_CONSTANT_STREAMWISE: return "Streamwise Constant";
694 case IC_MODE_FILE: return "File";
695 default: return "Unknown Initial Condition";
696 }
697}
698
699/**
700 * @brief Convert a FlowDirection enum value to its token string.
701 * @param[in] fd FlowDirection value.
702 * @return Constant string token (e.g. "+Zeta") or "from INLET" when FLOW_DIR_UNSET.
703 */
704const char* FlowDirectionToString(FlowDirection fd)
705{
706 switch ((int)fd) {
707 case FLOW_DIR_POS_XI: return "+Xi";
708 case FLOW_DIR_NEG_XI: return "-Xi";
709 case FLOW_DIR_POS_ETA: return "+Eta";
710 case FLOW_DIR_NEG_ETA: return "-Eta";
711 case FLOW_DIR_POS_ZETA: return "+Zeta";
712 case FLOW_DIR_NEG_ZETA: return "-Zeta";
713 default: return "from INLET";
714 }
715}
716
717/**
718 * @brief Implementation of \ref ParticleInitializationToString().
719 * @details Full API contract (arguments, ownership, side effects) is documented with
720 * the header declaration in `include/logging.h`.
721 * @see ParticleInitializationToString()
722 */
723const char* ParticleInitializationToString(ParticleInitializationType ParticleInitialization)
724{
725 switch(ParticleInitialization){
726 case PARTICLE_INIT_SURFACE_RANDOM: return "Surface: Random";
727 case PARTICLE_INIT_VOLUME: return "Volume";
728 case PARTICLE_INIT_POINT_SOURCE: return "Point Source";
729 case PARTICLE_INIT_SURFACE_EDGES: return "Surface: At edges";
730 default: return "Unknown Particle Initialization";
731 }
732}
733
734/**
735 * @brief Implementation of \ref LESModelToString().
736 * @details Full API contract (arguments, ownership, side effects) is documented with
737 * the header declaration in `include/logging.h`.
738 * @see LESModelToString()
739 */
740const char* LESModelToString(LESModelType LESFlag)
741{
742 switch(LESFlag){
743 case NO_LES_MODEL: return "No LES";
744 case CONSTANT_SMAGORINSKY: return "Constant Smagorinsky";
745 case DYNAMIC_SMAGORINSKY: return "Dynamic Smagorinsky";
746 default: return "Unknown LES Flag";
747 }
748}
749
750/**
751 * @brief Implementation of \ref MomentumSolverTypeToString().
752 * @details Full API contract (arguments, ownership, side effects) is documented with
753 * the header declaration in `include/logging.h`.
754 * @see MomentumSolverTypeToString()
755 */
756const char* MomentumSolverTypeToString(MomentumSolverType SolverFlag)
757{
758 switch(SolverFlag){
759 case MOMENTUM_SOLVER_EXPLICIT_RK: return "Explicit 4 stage Runge-Kutta ";
760 case MOMENTUM_SOLVER_DUALTIME_PICARD_JAMESON_RK: return "Dual Time Picard with 4-stage Jameson RK Smoothing";
761 case MOMENTUM_SOLVER_NEWTON_KRYLOV: return "Newton Krylov";
762 default: return "Unknown Momentum Solver Type";
763 }
764}
765
766/**
767 * @brief Implementation of \ref BCTypeToString().
768 * @details Full API contract (arguments, ownership, side effects) is documented with
769 * the header declaration in `include/logging.h`.
770 * @see BCTypeToString()
771 */
772const char* BCTypeToString(BCType type) {
773 switch (type) {
774 // case DIRICHLET: return "DIRICHLET";
775 // case NEUMANN: return "NEUMANN";
776 case WALL: return "WALL";
777 case INLET: return "INLET";
778 case OUTLET: return "OUTLET";
779 case FARFIELD: return "FARFIELD";
780 case PERIODIC: return "PERIODIC";
781 case INTERFACE: return "INTERFACE";
782
783 // case CUSTOM: return "CUSTOM";
784 default: return "Unknown BC Type";
785 }
786}
787
788/**
789 * @brief Internal helper implementation: `BCHandlerTypeToString()`.
790 * @details Local to this translation unit.
791 */
792const char* BCHandlerTypeToString(BCHandlerType handler_type) {
793 switch (handler_type) {
794 // Wall & Symmetry Handlers
795 case BC_HANDLER_WALL_NOSLIP: return "noslip";
796 case BC_HANDLER_WALL_MOVING: return "moving";
797 case BC_HANDLER_SYMMETRY_PLANE: return "symmetry_plane";
798
799 // Inlet Handlers
800 case BC_HANDLER_INLET_CONSTANT_VELOCITY: return "constant_velocity";
801 case BC_HANDLER_INLET_PULSATILE_FLUX: return "pulsatile_flux";
802 case BC_HANDLER_INLET_PARABOLIC: return "parabolic";
803 case BC_HANDLER_INLET_PROFILE_FROM_FILE: return "prescribed_flow";
804
805 // Outlet Handlers
806 case BC_HANDLER_OUTLET_CONSERVATION: return "conservation";
807 case BC_HANDLER_OUTLET_PRESSURE: return "pressure";
808
809 // Other Physical Handlers
810 case BC_HANDLER_FARFIELD_NONREFLECTING: return "nonreflecting";
811
812 // Multi-Block / Interface Handlers
813 case BC_HANDLER_PERIODIC_GEOMETRIC: return "geometric";
814 case BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX: return "constant flux";
815 case BC_HANDLER_PERIODIC_DRIVEN_INITIAL_FLUX: return "initial flux";
816 case BC_HANDLER_INTERFACE_OVERSET: return "overset";
817
818 // Default case
819 case BC_HANDLER_UNDEFINED:
820 default: return "UNKNOWN_HANDLER";
821 }
822}
823
824/**
825 * @brief Implementation of \ref DualMonitorDestroy().
826 * @details Full API contract (arguments, ownership, side effects) is documented with
827 * the header declaration in `include/logging.h`.
828 * @see DualMonitorDestroy()
829 */
830PetscErrorCode DualMonitorDestroy(void **ctx)
831{
832 DualMonitorCtx *monctx = (DualMonitorCtx*)*ctx;
833 PetscErrorCode ierr;
834 PetscMPIInt rank;
835
836 PetscFunctionBeginUser;
837 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank); CHKERRQ(ierr);
838 if(!rank && monctx->file_handle){
839 fclose(monctx->file_handle);
840 }
841
842 ierr = PetscFree(monctx); CHKERRQ(ierr);
843 *ctx = NULL;
844 PetscFunctionReturn(0);
845}
846
847/**
848 * @brief A custom KSP monitor that logs the true residual to a file and optionally to the console.
849 *
850 * This function replicates the behavior of KSPMonitorTrueResidualNorm by calculating
851 * the true residual norm ||b - Ax|| itself. It unconditionally logs to a file
852 * viewer and conditionally logs to the console based on a flag in the context.
853 *
854 * @param ksp The Krylov subspace context.
855 * @param it The current iteration number.
856 * @param rnorm The preconditioned residual norm (ignored, we compute our own).
857 * @param ctx A pointer to the DualMonitorCtx structure.
858 * @return PetscErrorCode 0 on success.
859 */
860#undef __FUNCT__
861#define __FUNCT__ "DualKSPMonitor"
862/**
863 * @brief Implementation of \ref DualKSPMonitor().
864 * @details Full API contract (arguments, ownership, side effects) is documented with
865 * the header declaration in `include/logging.h`.
866 * @see DualKSPMonitor()
867 */
868
869PetscErrorCode DualKSPMonitor(KSP ksp, PetscInt it, PetscReal rnorm, void *ctx)
870{
871 DualMonitorCtx *monctx = (DualMonitorCtx*)ctx;
872 PetscErrorCode ierr;
873 PetscReal trnorm, relnorm;
874 Vec r;
875 char norm_buf[256];
876 PetscMPIInt rank;
877
878 PetscFunctionBeginUser;
879 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank); CHKERRQ(ierr);
880
881 // 1. Calculate the true residual norm.
882 ierr = KSPBuildResidual(ksp, NULL, NULL, &r); CHKERRQ(ierr);
883 ierr = VecNorm(r, NORM_2, &trnorm); CHKERRQ(ierr);
884 ierr = VecDestroy(&r); CHKERRQ(ierr);
885
886 // 2. On the first iteration, compute and store the norm of the RHS vector `b`.
887 if (it == 0) {
888 Vec b;
889 ierr = KSPGetRhs(ksp, &b); CHKERRQ(ierr);
890 ierr = VecNorm(b, NORM_2, &monctx->bnorm); CHKERRQ(ierr);
891 }
892
893 if(!rank){
894 // 3. Compute the relative norm and format the output string.
895 if (monctx->bnorm > 1.e-15) {
896 relnorm = trnorm / monctx->bnorm;
897 sprintf(norm_buf, "ts: %-5d | block: %-2d | iter: %-3d | Unprecond Norm: %12.5e | True Norm: %12.5e | Rel Norm: %12.5e",(int)monctx->step, (int)monctx->block_id, (int)it, (double)rnorm, (double)trnorm, (double)relnorm);
898 } else {
899 sprintf(norm_buf,"ts: %-5d | block: %-2d | iter: %-3d | Unprecond Norm: %12.5e | True Norm: %12.5e",(int)monctx->step, (int)monctx->block_id, (int)it, (double)rnorm, (double)trnorm);
900 }
901
902 // 4. Log to the file viewer (unconditionally).
903 if(monctx->file_handle){
904 ierr = PetscFPrintf(PETSC_COMM_SELF,monctx->file_handle,"%s\n", norm_buf); CHKERRQ(ierr);
905 }
906 // 5. Log to the console (conditionally).
907 if (monctx->log_to_console) {
908 PetscFPrintf(PETSC_COMM_SELF,stdout, "%s\n", norm_buf); CHKERRQ(ierr);
909 }
910
911 } //rank
912
913 PetscFunctionReturn(0);
914}
915
916#define SOLUTION_CONVERGENCE_FLUID_THRESHOLD 0.1
917#define SOLUTION_CONVERGENCE_REL_EPS 1.0e-30
918
930
935
936/**
937 * @brief Forms a guarded relative metric for solution-convergence logging.
938 *
939 * This helper centralizes the divide-by-nearly-zero protection used by the
940 * solution-convergence logger when turning an absolute drift into a relative
941 * one. The denominator is clamped away from zero so warmup rows, quiescent
942 * fields, and statistically small observables do not generate infinities.
943 *
944 * @param[in] numerator Absolute quantity or drift magnitude.
945 * @param[in] denominator Reference magnitude used for normalization.
946 * @return Guarded relative value `numerator / max(|denominator|, eps)`.
947 */
948static PetscReal SolutionConvergenceSafeRelative(PetscReal numerator, PetscReal denominator)
949{
950 return numerator / PetscMax(PetscAbsReal(denominator), SOLUTION_CONVERGENCE_REL_EPS);
951}
952
953/**
954 * @brief Computes instantaneous global flow observables for statistical mode.
955 *
956 * The statistical solution-convergence path does not compare full Eulerian
957 * fields. Instead, it tracks a compact history of global observables derived
958 * from the completed Eulerian state. This helper computes the current
959 * volume-weighted fluid-domain mean speed and mean kinetic energy from `Ucat`.
960 *
961 * Only physical fluid cells contribute:
962 * - solid/immersed cells are excluded using `Nvert`
963 * - cells with near-zero metric Jacobian are ignored to avoid invalid volume
964 * weights
965 *
966 * Each MPI rank accumulates local partial sums and the routine reduces them to
967 * one global pair of observables.
968 *
969 * @param[in] simCtx Simulation context owning the finest-level flow
970 * fields.
971 * @param[out] mean_speed_out Volume-weighted domain mean of `|u|`.
972 * @param[out] mean_ke_out Volume-weighted domain mean of `0.5 |u|^2`.
973 * @return PetscErrorCode 0 on success.
974 */
975static PetscErrorCode ComputeCurrentFlowObservables(SimCtx *simCtx, PetscReal *mean_speed_out, PetscReal *mean_ke_out)
976{
977 PetscReal local[3] = {0.0, 0.0, 0.0};
978 PetscReal global[3] = {0.0, 0.0, 0.0};
979 UserCtx *user = NULL;
980
981 PetscFunctionBeginUser;
982 if (!simCtx || !mean_speed_out || !mean_ke_out) {
983 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "ComputeCurrentFlowObservables received a NULL argument.");
984 }
985
986 user = simCtx->usermg.mgctx[simCtx->usermg.mglevels - 1].user;
987
988 for (PetscInt bi = 0; bi < simCtx->block_number; ++bi) {
989 const DMDALocalInfo info = user[bi].info;
990 const PetscBool x_per = (PetscBool)(simCtx->i_periodic != 0);
991 const PetscBool y_per = (PetscBool)(simCtx->j_periodic != 0);
992 const PetscBool z_per = (PetscBool)(simCtx->k_periodic != 0);
993 const PetscInt i_end = (x_per && (info.xs + info.xm == info.mx)) ? info.mx - 1 : info.xs + info.xm;
994 const PetscInt j_end = (y_per && (info.ys + info.ym == info.my)) ? info.my - 1 : info.ys + info.ym;
995 const PetscInt k_end = (z_per && (info.zs + info.zm == info.mz)) ? info.mz - 1 : info.zs + info.zm;
996 Cmpnts ***ucat = NULL;
997 PetscReal ***aj = NULL;
998 PetscReal ***nvert = NULL;
999
1000 PetscCall(DMDAVecGetArrayRead(user[bi].fda, user[bi].Ucat, &ucat));
1001 PetscCall(DMDAVecGetArrayRead(user[bi].da, user[bi].Aj, &aj));
1002 PetscCall(DMDAVecGetArrayRead(user[bi].da, user[bi].Nvert, &nvert));
1003
1004 for (PetscInt k = info.zs; k < k_end; ++k) {
1005 for (PetscInt j = info.ys; j < j_end; ++j) {
1006 for (PetscInt i = info.xs; i < i_end; ++i) {
1007 PetscReal jac = aj[k][j][i];
1008 PetscReal cell_volume = 0.0;
1009 PetscReal speed = 0.0;
1010 PetscReal ke = 0.0;
1011
1012 if (nvert[k][j][i] > SOLUTION_CONVERGENCE_FLUID_THRESHOLD) continue;
1013 if (PetscAbsReal(jac) <= 1.0e-14) continue;
1014
1015 cell_volume = 1.0 / jac;
1016 speed = PetscSqrtReal(ucat[k][j][i].x * ucat[k][j][i].x +
1017 ucat[k][j][i].y * ucat[k][j][i].y +
1018 ucat[k][j][i].z * ucat[k][j][i].z);
1019 ke = 0.5 * speed * speed;
1020
1021 local[0] += cell_volume;
1022 local[1] += speed * cell_volume;
1023 local[2] += ke * cell_volume;
1024 }
1025 }
1026 }
1027
1028 PetscCall(DMDAVecRestoreArrayRead(user[bi].da, user[bi].Nvert, &nvert));
1029 PetscCall(DMDAVecRestoreArrayRead(user[bi].da, user[bi].Aj, &aj));
1030 PetscCall(DMDAVecRestoreArrayRead(user[bi].fda, user[bi].Ucat, &ucat));
1031 }
1032
1033 PetscCallMPI(MPI_Allreduce(local, global, 3, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD));
1034
1035 if (global[0] <= 0.0) {
1036 *mean_speed_out = 0.0;
1037 *mean_ke_out = 0.0;
1038 } else {
1039 *mean_speed_out = global[1] / global[0];
1040 *mean_ke_out = global[2] / global[0];
1041 }
1042
1043 PetscFunctionReturn(0);
1044}
1045
1046/**
1047 * @brief Computes deterministic solution-drift metrics for the current step.
1048 *
1049 * This helper powers both `steady_deterministic` and
1050 * `periodic_deterministic` solution-convergence modes. It compares the
1051 * completed Eulerian state against either:
1052 * - the previous physical timestep (`Ucat_o`, `P_o`) for steady/transient use
1053 * - the stored phase-aligned snapshot ring for periodic use
1054 *
1055 * The routine performs two passes over the fluid cells:
1056 * 1. velocity/observable pass
1057 * - computes current and reference mean speed / mean KE
1058 * - computes absolute/relative velocity L2 drift
1059 * - accumulates pressure means needed for gauge removal
1060 * 2. pressure-only pass
1061 * - subtracts the volume-weighted mean pressure from both states
1062 * - computes gauge-invariant pressure L2 drift
1063 *
1064 * Warmup behavior is handled here. If no valid reference exists yet, all drift
1065 * outputs are left at zero and `has_reference_out` is set to `PETSC_FALSE`,
1066 * while current observables are still reported.
1067 *
1068 * @param[in] simCtx Simulation context owning the current state.
1069 * @param[in] periodic_mode `PETSC_TRUE` when comparing against
1070 * phase-aligned periodic storage.
1071 * @param[in] phase_step Active phase slot for periodic mode, or `-1`
1072 * when unused.
1073 * @param[in] samples_before Number of solution-convergence samples
1074 * already recorded before this timestep.
1075 * @param[out] has_reference_out Whether a valid comparison state existed.
1076 * @param[out] u_abs_l2_out Absolute L2 drift of Cartesian velocity.
1077 * @param[out] u_rel_l2_out Relative L2 drift of Cartesian velocity.
1078 * @param[out] p_abs_l2_out Gauge-invariant absolute L2 pressure drift.
1079 * @param[out] p_rel_l2_out Gauge-invariant relative L2 pressure drift.
1080 * @param[out] mean_speed_out Current volume-weighted mean speed.
1081 * @param[out] mean_speed_ref_out Reference volume-weighted mean speed.
1082 * @param[out] mean_speed_abs_out Absolute drift of mean speed.
1083 * @param[out] mean_speed_rel_out Relative drift of mean speed.
1084 * @param[out] mean_ke_out Current volume-weighted mean kinetic energy.
1085 * @param[out] mean_ke_ref_out Reference volume-weighted mean kinetic
1086 * energy.
1087 * @param[out] mean_ke_abs_out Absolute drift of mean kinetic energy.
1088 * @param[out] mean_ke_rel_out Relative drift of mean kinetic energy.
1089 * @return PetscErrorCode 0 on success.
1090 */
1091static PetscErrorCode ComputeDeterministicSolutionMetrics(SimCtx *simCtx,
1092 PetscBool periodic_mode,
1093 PetscInt phase_step,
1094 PetscInt samples_before,
1095 PetscBool *has_reference_out,
1096 PetscReal *u_abs_l2_out,
1097 PetscReal *u_rel_l2_out,
1098 PetscReal *p_abs_l2_out,
1099 PetscReal *p_rel_l2_out,
1100 PetscReal *mean_speed_out,
1101 PetscReal *mean_speed_ref_out,
1102 PetscReal *mean_speed_abs_out,
1103 PetscReal *mean_speed_rel_out,
1104 PetscReal *mean_ke_out,
1105 PetscReal *mean_ke_ref_out,
1106 PetscReal *mean_ke_abs_out,
1107 PetscReal *mean_ke_rel_out)
1108{
1109 SolutionConvergenceDeterministicPass1 local_pass1 = {0};
1110 SolutionConvergenceDeterministicPass1 global_pass1 = {0};
1111 SolutionConvergenceDeterministicPass2 local_pass2 = {0};
1112 SolutionConvergenceDeterministicPass2 global_pass2 = {0};
1113 PetscReal current_pressure_mean = 0.0;
1114 PetscReal reference_pressure_mean = 0.0;
1115 UserCtx *user = NULL;
1116
1117 PetscFunctionBeginUser;
1118 if (!simCtx || !has_reference_out || !u_abs_l2_out || !u_rel_l2_out || !p_abs_l2_out || !p_rel_l2_out ||
1119 !mean_speed_out || !mean_speed_ref_out || !mean_speed_abs_out || !mean_speed_rel_out ||
1120 !mean_ke_out || !mean_ke_ref_out || !mean_ke_abs_out || !mean_ke_rel_out) {
1121 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "ComputeDeterministicSolutionMetrics received a NULL output pointer.");
1122 }
1123
1124 *has_reference_out = periodic_mode
1125 ? (PetscBool)(simCtx->solutionConvergencePeriodSteps > 0 &&
1126 phase_step >= 0 &&
1127 phase_step < simCtx->solutionConvergencePeriodSteps &&
1128 samples_before >= simCtx->solutionConvergencePeriodSteps)
1129 : (PetscBool)(samples_before > 0);
1130
1131 *u_abs_l2_out = 0.0;
1132 *u_rel_l2_out = 0.0;
1133 *p_abs_l2_out = 0.0;
1134 *p_rel_l2_out = 0.0;
1135 *mean_speed_out = 0.0;
1136 *mean_speed_ref_out = 0.0;
1137 *mean_speed_abs_out = 0.0;
1138 *mean_speed_rel_out = 0.0;
1139 *mean_ke_out = 0.0;
1140 *mean_ke_ref_out = 0.0;
1141 *mean_ke_abs_out = 0.0;
1142 *mean_ke_rel_out = 0.0;
1143
1144 user = simCtx->usermg.mgctx[simCtx->usermg.mglevels - 1].user;
1145
1146 for (PetscInt bi = 0; bi < simCtx->block_number; ++bi) {
1147 const DMDALocalInfo info = user[bi].info;
1148 const PetscBool x_per = (PetscBool)(simCtx->i_periodic != 0);
1149 const PetscBool y_per = (PetscBool)(simCtx->j_periodic != 0);
1150 const PetscBool z_per = (PetscBool)(simCtx->k_periodic != 0);
1151 const PetscInt i_end = (x_per && (info.xs + info.xm == info.mx)) ? info.mx - 1 : info.xs + info.xm;
1152 const PetscInt j_end = (y_per && (info.ys + info.ym == info.my)) ? info.my - 1 : info.ys + info.ym;
1153 const PetscInt k_end = (z_per && (info.zs + info.zm == info.mz)) ? info.mz - 1 : info.zs + info.zm;
1154 Cmpnts ***ucat = NULL;
1155 Cmpnts ***ucat_ref = NULL;
1156 PetscReal ***pressure = NULL;
1157 PetscReal ***pressure_ref = NULL;
1158 PetscReal ***aj = NULL;
1159 PetscReal ***nvert = NULL;
1160 Vec ucat_reference_vec = NULL;
1161 Vec pressure_reference_vec = NULL;
1162
1163 if (*has_reference_out) {
1164 if (periodic_mode) {
1165 ucat_reference_vec = user[bi].solutionConvergencePeriodicUcatRef[phase_step];
1166 pressure_reference_vec = user[bi].solutionConvergencePeriodicPRef[phase_step];
1167 } else {
1168 ucat_reference_vec = user[bi].Ucat_o;
1169 pressure_reference_vec = user[bi].P_o;
1170 }
1171 }
1172
1173 PetscCall(DMDAVecGetArrayRead(user[bi].fda, user[bi].Ucat, &ucat));
1174 PetscCall(DMDAVecGetArrayRead(user[bi].da, user[bi].P, &pressure));
1175 PetscCall(DMDAVecGetArrayRead(user[bi].da, user[bi].Aj, &aj));
1176 PetscCall(DMDAVecGetArrayRead(user[bi].da, user[bi].Nvert, &nvert));
1177 if (*has_reference_out) {
1178 PetscCall(DMDAVecGetArrayRead(user[bi].fda, ucat_reference_vec, &ucat_ref));
1179 PetscCall(DMDAVecGetArrayRead(user[bi].da, pressure_reference_vec, &pressure_ref));
1180 }
1181
1182 for (PetscInt k = info.zs; k < k_end; ++k) {
1183 for (PetscInt j = info.ys; j < j_end; ++j) {
1184 for (PetscInt i = info.xs; i < i_end; ++i) {
1185 PetscReal jac = aj[k][j][i];
1186 PetscReal cell_volume = 0.0;
1187 PetscReal speed = 0.0;
1188 PetscReal ke = 0.0;
1189
1190 if (nvert[k][j][i] > SOLUTION_CONVERGENCE_FLUID_THRESHOLD) continue;
1191 if (PetscAbsReal(jac) <= 1.0e-14) continue;
1192
1193 cell_volume = 1.0 / jac;
1194 speed = PetscSqrtReal(ucat[k][j][i].x * ucat[k][j][i].x +
1195 ucat[k][j][i].y * ucat[k][j][i].y +
1196 ucat[k][j][i].z * ucat[k][j][i].z);
1197 ke = 0.5 * speed * speed;
1198
1199 local_pass1.fluid_volume += cell_volume;
1200 local_pass1.current_speed_sum += speed * cell_volume;
1201 local_pass1.current_ke_sum += ke * cell_volume;
1202 local_pass1.current_u_norm_sq += (ucat[k][j][i].x * ucat[k][j][i].x +
1203 ucat[k][j][i].y * ucat[k][j][i].y +
1204 ucat[k][j][i].z * ucat[k][j][i].z) * cell_volume;
1205
1206 if (*has_reference_out) {
1207 PetscReal ref_speed = PetscSqrtReal(ucat_ref[k][j][i].x * ucat_ref[k][j][i].x +
1208 ucat_ref[k][j][i].y * ucat_ref[k][j][i].y +
1209 ucat_ref[k][j][i].z * ucat_ref[k][j][i].z);
1210 PetscReal ref_ke = 0.5 * ref_speed * ref_speed;
1211 PetscReal dux = ucat[k][j][i].x - ucat_ref[k][j][i].x;
1212 PetscReal duy = ucat[k][j][i].y - ucat_ref[k][j][i].y;
1213 PetscReal duz = ucat[k][j][i].z - ucat_ref[k][j][i].z;
1214
1215 local_pass1.reference_speed_sum += ref_speed * cell_volume;
1216 local_pass1.reference_ke_sum += ref_ke * cell_volume;
1217 local_pass1.delta_u_norm_sq += (dux * dux + duy * duy + duz * duz) * cell_volume;
1218 local_pass1.current_pressure_sum += pressure[k][j][i] * cell_volume;
1219 local_pass1.reference_pressure_sum += pressure_ref[k][j][i] * cell_volume;
1220 }
1221 }
1222 }
1223 }
1224
1225 if (*has_reference_out) {
1226 PetscCall(DMDAVecRestoreArrayRead(user[bi].da, pressure_reference_vec, &pressure_ref));
1227 PetscCall(DMDAVecRestoreArrayRead(user[bi].fda, ucat_reference_vec, &ucat_ref));
1228 }
1229 PetscCall(DMDAVecRestoreArrayRead(user[bi].da, user[bi].Nvert, &nvert));
1230 PetscCall(DMDAVecRestoreArrayRead(user[bi].da, user[bi].Aj, &aj));
1231 PetscCall(DMDAVecRestoreArrayRead(user[bi].da, user[bi].P, &pressure));
1232 PetscCall(DMDAVecRestoreArrayRead(user[bi].fda, user[bi].Ucat, &ucat));
1233 }
1234
1235 PetscCallMPI(MPI_Allreduce(&local_pass1, &global_pass1,
1236 sizeof(SolutionConvergenceDeterministicPass1) / sizeof(PetscReal),
1237 MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD));
1238
1239 if (global_pass1.fluid_volume <= 0.0) PetscFunctionReturn(0);
1240
1241 *mean_speed_out = global_pass1.current_speed_sum / global_pass1.fluid_volume;
1242 *mean_ke_out = global_pass1.current_ke_sum / global_pass1.fluid_volume;
1243
1244 if (!*has_reference_out) PetscFunctionReturn(0);
1245
1246 *mean_speed_ref_out = global_pass1.reference_speed_sum / global_pass1.fluid_volume;
1247 *mean_speed_abs_out = PetscAbsReal(*mean_speed_out - *mean_speed_ref_out);
1248 *mean_speed_rel_out = SolutionConvergenceSafeRelative(*mean_speed_abs_out, *mean_speed_out);
1249 *mean_ke_ref_out = global_pass1.reference_ke_sum / global_pass1.fluid_volume;
1250 *mean_ke_abs_out = PetscAbsReal(*mean_ke_out - *mean_ke_ref_out);
1251 *mean_ke_rel_out = SolutionConvergenceSafeRelative(*mean_ke_abs_out, *mean_ke_out);
1252
1253 current_pressure_mean = global_pass1.current_pressure_sum / global_pass1.fluid_volume;
1254 reference_pressure_mean = global_pass1.reference_pressure_sum / global_pass1.fluid_volume;
1255
1256 for (PetscInt bi = 0; bi < simCtx->block_number; ++bi) {
1257 const DMDALocalInfo info = user[bi].info;
1258 const PetscBool x_per = (PetscBool)(simCtx->i_periodic != 0);
1259 const PetscBool y_per = (PetscBool)(simCtx->j_periodic != 0);
1260 const PetscBool z_per = (PetscBool)(simCtx->k_periodic != 0);
1261 const PetscInt i_end = (x_per && (info.xs + info.xm == info.mx)) ? info.mx - 1 : info.xs + info.xm;
1262 const PetscInt j_end = (y_per && (info.ys + info.ym == info.my)) ? info.my - 1 : info.ys + info.ym;
1263 const PetscInt k_end = (z_per && (info.zs + info.zm == info.mz)) ? info.mz - 1 : info.zs + info.zm;
1264 PetscReal ***pressure = NULL;
1265 PetscReal ***pressure_ref = NULL;
1266 PetscReal ***aj = NULL;
1267 PetscReal ***nvert = NULL;
1268 Vec pressure_reference_vec = periodic_mode ? user[bi].solutionConvergencePeriodicPRef[phase_step] : user[bi].P_o;
1269
1270 PetscCall(DMDAVecGetArrayRead(user[bi].da, user[bi].P, &pressure));
1271 PetscCall(DMDAVecGetArrayRead(user[bi].da, pressure_reference_vec, &pressure_ref));
1272 PetscCall(DMDAVecGetArrayRead(user[bi].da, user[bi].Aj, &aj));
1273 PetscCall(DMDAVecGetArrayRead(user[bi].da, user[bi].Nvert, &nvert));
1274
1275 for (PetscInt k = info.zs; k < k_end; ++k) {
1276 for (PetscInt j = info.ys; j < j_end; ++j) {
1277 for (PetscInt i = info.xs; i < i_end; ++i) {
1278 PetscReal jac = aj[k][j][i];
1279 PetscReal cell_volume = 0.0;
1280 PetscReal current_pressure = 0.0;
1281 PetscReal reference_pressure = 0.0;
1282 PetscReal delta_pressure = 0.0;
1283
1284 if (nvert[k][j][i] > SOLUTION_CONVERGENCE_FLUID_THRESHOLD) continue;
1285 if (PetscAbsReal(jac) <= 1.0e-14) continue;
1286
1287 cell_volume = 1.0 / jac;
1288 current_pressure = pressure[k][j][i] - current_pressure_mean;
1289 reference_pressure = pressure_ref[k][j][i] - reference_pressure_mean;
1290 delta_pressure = current_pressure - reference_pressure;
1291
1292 local_pass2.current_pressure_norm_sq += current_pressure * current_pressure * cell_volume;
1293 local_pass2.delta_pressure_norm_sq += delta_pressure * delta_pressure * cell_volume;
1294 }
1295 }
1296 }
1297
1298 PetscCall(DMDAVecRestoreArrayRead(user[bi].da, user[bi].Nvert, &nvert));
1299 PetscCall(DMDAVecRestoreArrayRead(user[bi].da, user[bi].Aj, &aj));
1300 PetscCall(DMDAVecRestoreArrayRead(user[bi].da, pressure_reference_vec, &pressure_ref));
1301 PetscCall(DMDAVecRestoreArrayRead(user[bi].da, user[bi].P, &pressure));
1302 }
1303
1304 PetscCallMPI(MPI_Allreduce(&local_pass2, &global_pass2,
1305 sizeof(SolutionConvergenceDeterministicPass2) / sizeof(PetscReal),
1306 MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD));
1307
1308 *u_abs_l2_out = PetscSqrtReal(global_pass1.delta_u_norm_sq);
1309 *u_rel_l2_out = SolutionConvergenceSafeRelative(*u_abs_l2_out, PetscSqrtReal(global_pass1.current_u_norm_sq));
1310 *p_abs_l2_out = PetscSqrtReal(global_pass2.delta_pressure_norm_sq);
1311 *p_rel_l2_out = SolutionConvergenceSafeRelative(*p_abs_l2_out, PetscSqrtReal(global_pass2.current_pressure_norm_sq));
1312
1313 PetscFunctionReturn(0);
1314}
1315
1316/**
1317 * @brief Reads one sample from the statistical ring buffer by age.
1318 *
1319 * The statistical logger stores scalar observables in a compact circular
1320 * buffer. This helper interprets the buffer using `samples_available` as the
1321 * logical end of the history and returns the entry `offset_from_latest` steps
1322 * back from the newest stored sample.
1323 *
1324 * Out-of-range requests return zero so warmup handling can remain simple and
1325 * deterministic.
1326 *
1327 * @param[in] history Ring-buffer storage array.
1328 * @param[in] capacity Total ring-buffer capacity.
1329 * @param[in] samples_available Number of logical samples available to read.
1330 * @param[in] offset_from_latest `0` means newest sample, `1` previous sample,
1331 * and so on.
1332 * @return Requested historical sample, or `0.0` if unavailable.
1333 */
1334static PetscReal SolutionConvergenceHistoryGet(const PetscReal *history,
1335 PetscInt capacity,
1336 PetscInt samples_available,
1337 PetscInt offset_from_latest)
1338{
1339 PetscInt count = 0;
1340 PetscInt index = 0;
1341
1342 if (!history || capacity <= 0 || samples_available <= 0 || offset_from_latest < 0) {
1343 return 0.0;
1344 }
1345
1346 count = PetscMin(samples_available, capacity);
1347 if (offset_from_latest >= count) return 0.0;
1348
1349 index = (samples_available - 1 - offset_from_latest) % capacity;
1350 if (index < 0) index += capacity;
1351 return history[index];
1352}
1353
1354/**
1355 * @brief Appends one timestep's scalar observables to the statistical history.
1356 *
1357 * Statistical solution-convergence compares adjacent windows of scalar
1358 * observables rather than full fields. This helper writes the current
1359 * `mean_speed` and `mean_ke` into the rolling history arrays using the current
1360 * sample count to choose the circular-buffer slot.
1361 *
1362 * @param[in,out] simCtx Simulation context owning the history arrays.
1363 * @param[in] samples_before Number of samples present before appending the
1364 * current timestep.
1365 * @param[in] mean_speed Current timestep mean-speed observable.
1366 * @param[in] mean_ke Current timestep mean-KE observable.
1367 * @return PetscErrorCode 0 on success.
1368 */
1369static PetscErrorCode AppendStatisticalObservableSample(SimCtx *simCtx,
1370 PetscInt samples_before,
1371 PetscReal mean_speed,
1372 PetscReal mean_ke)
1373{
1374 PetscInt history_capacity = 0;
1375 PetscInt slot = 0;
1376
1377 PetscFunctionBeginUser;
1378 if (!simCtx || !simCtx->solutionConvergenceMeanSpeedHistory || !simCtx->solutionConvergenceMeanKEHistory) {
1379 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Statistical solution-convergence history is not allocated.");
1380 }
1381
1382 history_capacity = 2 * simCtx->solutionConvergenceWindowSteps;
1383 if (history_capacity <= 0) {
1384 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Statistical solution-convergence history capacity must be positive.");
1385 }
1386
1387 slot = samples_before % history_capacity;
1388 simCtx->solutionConvergenceMeanSpeedHistory[slot] = mean_speed;
1389 simCtx->solutionConvergenceMeanKEHistory[slot] = mean_ke;
1390
1391 PetscFunctionReturn(0);
1392}
1393
1394/**
1395 * @brief Computes adjacent-window drift metrics for statistical steady mode.
1396 *
1397 * Once enough samples have been accumulated, this helper forms:
1398 * - the current window over the most recent `window_steps` samples
1399 * - the previous adjacent window over the preceding `window_steps` samples
1400 *
1401 * From those windows it computes means, RMS values, and absolute/relative
1402 * drift for both tracked observables (`mean_speed` and `mean_ke`). When the
1403 * history is still warming up:
1404 * - fewer than `window_steps` samples: no window metrics are available
1405 * - between `window_steps` and `2*window_steps - 1` samples: current-window
1406 * metrics are available, but no reference window exists yet
1407 *
1408 * @param[in] simCtx Simulation context owning the
1409 * statistical history.
1410 * @param[in] samples_available Number of samples available after
1411 * appending the current timestep.
1412 * @param[out] has_reference_out Whether both adjacent windows
1413 * exist and drift metrics are
1414 * meaningful.
1415 * @param[out] mean_speed_window_out Mean speed over the current
1416 * window.
1417 * @param[out] mean_speed_window_prev_out Mean speed over the previous
1418 * window.
1419 * @param[out] mean_speed_window_abs_out Absolute drift between current
1420 * and previous window means.
1421 * @param[out] mean_speed_window_rel_out Relative drift between current
1422 * and previous window means.
1423 * @param[out] mean_speed_rms_window_out RMS of mean-speed samples in the
1424 * current window.
1425 * @param[out] mean_speed_rms_window_prev_out RMS of mean-speed samples in the
1426 * previous window.
1427 * @param[out] mean_speed_rms_window_abs_out Absolute drift between window RMS
1428 * values.
1429 * @param[out] mean_speed_rms_window_rel_out Relative drift between window RMS
1430 * values.
1431 * @param[out] mean_ke_window_out Mean kinetic energy over the
1432 * current window.
1433 * @param[out] mean_ke_window_prev_out Mean kinetic energy over the
1434 * previous window.
1435 * @param[out] mean_ke_window_abs_out Absolute drift between current
1436 * and previous KE-window means.
1437 * @param[out] mean_ke_window_rel_out Relative drift between current
1438 * and previous KE-window means.
1439 * @param[out] mean_ke_rms_window_out RMS of mean-KE samples in the
1440 * current window.
1441 * @param[out] mean_ke_rms_window_prev_out RMS of mean-KE samples in the
1442 * previous window.
1443 * @param[out] mean_ke_rms_window_abs_out Absolute drift between KE-window
1444 * RMS values.
1445 * @param[out] mean_ke_rms_window_rel_out Relative drift between KE-window
1446 * RMS values.
1447 * @return PetscErrorCode 0 on success.
1448 */
1449static PetscErrorCode ComputeStatisticalWindowMetrics(const SimCtx *simCtx,
1450 PetscInt samples_available,
1451 PetscBool *has_reference_out,
1452 PetscReal *mean_speed_window_out,
1453 PetscReal *mean_speed_window_prev_out,
1454 PetscReal *mean_speed_window_abs_out,
1455 PetscReal *mean_speed_window_rel_out,
1456 PetscReal *mean_speed_rms_window_out,
1457 PetscReal *mean_speed_rms_window_prev_out,
1458 PetscReal *mean_speed_rms_window_abs_out,
1459 PetscReal *mean_speed_rms_window_rel_out,
1460 PetscReal *mean_ke_window_out,
1461 PetscReal *mean_ke_window_prev_out,
1462 PetscReal *mean_ke_window_abs_out,
1463 PetscReal *mean_ke_window_rel_out,
1464 PetscReal *mean_ke_rms_window_out,
1465 PetscReal *mean_ke_rms_window_prev_out,
1466 PetscReal *mean_ke_rms_window_abs_out,
1467 PetscReal *mean_ke_rms_window_rel_out)
1468{
1469 PetscInt w = 0;
1470 PetscInt history_capacity = 0;
1471 PetscReal speed_sum = 0.0;
1472 PetscReal speed_sum_sq = 0.0;
1473 PetscReal speed_prev_sum = 0.0;
1474 PetscReal speed_prev_sum_sq = 0.0;
1475 PetscReal ke_sum = 0.0;
1476 PetscReal ke_sum_sq = 0.0;
1477 PetscReal ke_prev_sum = 0.0;
1478 PetscReal ke_prev_sum_sq = 0.0;
1479
1480 PetscFunctionBeginUser;
1481 if (!simCtx || !has_reference_out || !mean_speed_window_out || !mean_speed_window_prev_out ||
1482 !mean_speed_window_abs_out || !mean_speed_window_rel_out || !mean_speed_rms_window_out ||
1483 !mean_speed_rms_window_prev_out || !mean_speed_rms_window_abs_out || !mean_speed_rms_window_rel_out ||
1484 !mean_ke_window_out || !mean_ke_window_prev_out || !mean_ke_window_abs_out || !mean_ke_window_rel_out ||
1485 !mean_ke_rms_window_out || !mean_ke_rms_window_prev_out || !mean_ke_rms_window_abs_out ||
1486 !mean_ke_rms_window_rel_out) {
1487 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "ComputeStatisticalWindowMetrics received a NULL output pointer.");
1488 }
1489
1490 *has_reference_out = PETSC_FALSE;
1491 *mean_speed_window_out = 0.0;
1492 *mean_speed_window_prev_out = 0.0;
1493 *mean_speed_window_abs_out = 0.0;
1494 *mean_speed_window_rel_out = 0.0;
1495 *mean_speed_rms_window_out = 0.0;
1496 *mean_speed_rms_window_prev_out = 0.0;
1497 *mean_speed_rms_window_abs_out = 0.0;
1498 *mean_speed_rms_window_rel_out = 0.0;
1499 *mean_ke_window_out = 0.0;
1500 *mean_ke_window_prev_out = 0.0;
1501 *mean_ke_window_abs_out = 0.0;
1502 *mean_ke_window_rel_out = 0.0;
1503 *mean_ke_rms_window_out = 0.0;
1504 *mean_ke_rms_window_prev_out = 0.0;
1505 *mean_ke_rms_window_abs_out = 0.0;
1506 *mean_ke_rms_window_rel_out = 0.0;
1507
1508 w = simCtx->solutionConvergenceWindowSteps;
1509 history_capacity = 2 * w;
1510 if (w <= 0 || samples_available < w) PetscFunctionReturn(0);
1511
1512 for (PetscInt idx = 0; idx < w; ++idx) {
1513 PetscReal speed_value = SolutionConvergenceHistoryGet(simCtx->solutionConvergenceMeanSpeedHistory,
1514 history_capacity,
1515 samples_available,
1516 idx);
1517 PetscReal ke_value = SolutionConvergenceHistoryGet(simCtx->solutionConvergenceMeanKEHistory,
1518 history_capacity,
1519 samples_available,
1520 idx);
1521 speed_sum += speed_value;
1522 speed_sum_sq += speed_value * speed_value;
1523 ke_sum += ke_value;
1524 ke_sum_sq += ke_value * ke_value;
1525 }
1526
1527 *mean_speed_window_out = speed_sum / (PetscReal)w;
1528 *mean_speed_rms_window_out = PetscSqrtReal(PetscMax(0.0, speed_sum_sq / (PetscReal)w -
1529 (*mean_speed_window_out) * (*mean_speed_window_out)));
1530 *mean_ke_window_out = ke_sum / (PetscReal)w;
1531 *mean_ke_rms_window_out = PetscSqrtReal(PetscMax(0.0, ke_sum_sq / (PetscReal)w -
1532 (*mean_ke_window_out) * (*mean_ke_window_out)));
1533
1534 if (samples_available < 2 * w) PetscFunctionReturn(0);
1535
1536 for (PetscInt idx = w; idx < 2 * w; ++idx) {
1537 PetscReal speed_value = SolutionConvergenceHistoryGet(simCtx->solutionConvergenceMeanSpeedHistory,
1538 history_capacity,
1539 samples_available,
1540 idx);
1541 PetscReal ke_value = SolutionConvergenceHistoryGet(simCtx->solutionConvergenceMeanKEHistory,
1542 history_capacity,
1543 samples_available,
1544 idx);
1545 speed_prev_sum += speed_value;
1546 speed_prev_sum_sq += speed_value * speed_value;
1547 ke_prev_sum += ke_value;
1548 ke_prev_sum_sq += ke_value * ke_value;
1549 }
1550
1551 *has_reference_out = PETSC_TRUE;
1552 *mean_speed_window_prev_out = speed_prev_sum / (PetscReal)w;
1553 *mean_speed_window_abs_out = PetscAbsReal(*mean_speed_window_out - *mean_speed_window_prev_out);
1554 *mean_speed_window_rel_out = SolutionConvergenceSafeRelative(*mean_speed_window_abs_out, *mean_speed_window_out);
1555 *mean_speed_rms_window_prev_out = PetscSqrtReal(PetscMax(0.0, speed_prev_sum_sq / (PetscReal)w -
1556 (*mean_speed_window_prev_out) * (*mean_speed_window_prev_out)));
1557 *mean_speed_rms_window_abs_out = PetscAbsReal(*mean_speed_rms_window_out - *mean_speed_rms_window_prev_out);
1558 *mean_speed_rms_window_rel_out = SolutionConvergenceSafeRelative(*mean_speed_rms_window_abs_out, *mean_speed_rms_window_out);
1559
1560 *mean_ke_window_prev_out = ke_prev_sum / (PetscReal)w;
1561 *mean_ke_window_abs_out = PetscAbsReal(*mean_ke_window_out - *mean_ke_window_prev_out);
1562 *mean_ke_window_rel_out = SolutionConvergenceSafeRelative(*mean_ke_window_abs_out, *mean_ke_window_out);
1563 *mean_ke_rms_window_prev_out = PetscSqrtReal(PetscMax(0.0, ke_prev_sum_sq / (PetscReal)w -
1564 (*mean_ke_window_prev_out) * (*mean_ke_window_prev_out)));
1565 *mean_ke_rms_window_abs_out = PetscAbsReal(*mean_ke_rms_window_out - *mean_ke_rms_window_prev_out);
1566 *mean_ke_rms_window_rel_out = SolutionConvergenceSafeRelative(*mean_ke_rms_window_abs_out, *mean_ke_rms_window_out);
1567
1568 PetscFunctionReturn(0);
1569}
1570
1571/**
1572 * @brief Maps the internal solution-convergence mode enum to its log label.
1573 *
1574 * The logger writes a human-readable mode string into the
1575 * `solution_convergence.log` banner and `mode` column. This helper keeps the
1576 * formatting centralized so the file output stays consistent with the accepted
1577 * configuration names.
1578 *
1579 * @param[in] mode Internal solution-convergence mode selector.
1580 * @return Lowercase string label written to the log output.
1581 */
1582static const char *SolutionConvergenceModeToString(SolutionConvergenceMode mode)
1583{
1584 switch (mode) {
1585 case SOLUTION_CONVERGENCE_STEADY_DETERMINISTIC: return "steady_deterministic";
1586 case SOLUTION_CONVERGENCE_PERIODIC_DETERMINISTIC: return "periodic_deterministic";
1587 case SOLUTION_CONVERGENCE_STATISTICAL_STEADY: return "statistical_steady";
1588 case SOLUTION_CONVERGENCE_TRANSIENT: return "transient";
1589 default: return "unknown";
1590 }
1591}
1592
1593/**
1594 * @brief Implementation of \ref LOG_SOLUTION_CONVERGENCE().
1595 * @details Full API contract (arguments, ownership, side effects) is documented with
1596 * the header declaration in `include/logging.h`.
1597 * @see LOG_SOLUTION_CONVERGENCE()
1598 */
1599PetscErrorCode LOG_SOLUTION_CONVERGENCE(SimCtx *simCtx)
1600{
1601 PetscMPIInt rank = 0;
1602 PetscBool has_reference = PETSC_FALSE;
1603 PetscInt phase_step = -1;
1604 PetscInt samples_before = 0;
1605 PetscReal u_abs_l2 = 0.0, u_rel_l2 = 0.0, p_abs_l2 = 0.0, p_rel_l2 = 0.0;
1606 PetscReal mean_speed = 0.0, mean_speed_reference = 0.0, mean_speed_abs_drift = 0.0, mean_speed_rel_drift = 0.0;
1607 PetscReal mean_ke = 0.0, mean_ke_reference = 0.0, mean_ke_abs_drift = 0.0, mean_ke_rel_drift = 0.0;
1608 PetscReal mean_speed_window = 0.0, mean_speed_window_prev = 0.0, mean_speed_window_abs_drift = 0.0, mean_speed_window_rel_drift = 0.0;
1609 PetscReal mean_speed_rms_window = 0.0, mean_speed_rms_window_prev = 0.0, mean_speed_rms_window_abs_drift = 0.0, mean_speed_rms_window_rel_drift = 0.0;
1610 PetscReal mean_ke_window = 0.0, mean_ke_window_prev = 0.0, mean_ke_window_abs_drift = 0.0, mean_ke_window_rel_drift = 0.0;
1611 PetscReal mean_ke_rms_window = 0.0, mean_ke_rms_window_prev = 0.0, mean_ke_rms_window_abs_drift = 0.0, mean_ke_rms_window_rel_drift = 0.0;
1612
1613 PetscFunctionBeginUser;
1614 if (!simCtx) PetscFunctionReturn(0);
1615 if (simCtx->exec_mode != EXEC_MODE_SOLVER) PetscFunctionReturn(0);
1616
1617 samples_before = simCtx->solutionConvergenceSamplesRecorded;
1618
1619 switch (simCtx->solutionConvergenceMode) {
1620 case SOLUTION_CONVERGENCE_STEADY_DETERMINISTIC:
1621 case SOLUTION_CONVERGENCE_TRANSIENT:
1622 PetscCall(ComputeDeterministicSolutionMetrics(simCtx, PETSC_FALSE, -1, samples_before,
1623 &has_reference,
1624 &u_abs_l2, &u_rel_l2,
1625 &p_abs_l2, &p_rel_l2,
1626 &mean_speed, &mean_speed_reference,
1627 &mean_speed_abs_drift, &mean_speed_rel_drift,
1628 &mean_ke, &mean_ke_reference,
1629 &mean_ke_abs_drift, &mean_ke_rel_drift));
1630 break;
1631 case SOLUTION_CONVERGENCE_PERIODIC_DETERMINISTIC:
1632 phase_step = simCtx->solutionConvergencePeriodSteps > 0 ? (simCtx->step % simCtx->solutionConvergencePeriodSteps) : -1;
1633 PetscCall(ComputeDeterministicSolutionMetrics(simCtx, PETSC_TRUE, phase_step, samples_before,
1634 &has_reference,
1635 &u_abs_l2, &u_rel_l2,
1636 &p_abs_l2, &p_rel_l2,
1637 &mean_speed, &mean_speed_reference,
1638 &mean_speed_abs_drift, &mean_speed_rel_drift,
1639 &mean_ke, &mean_ke_reference,
1640 &mean_ke_abs_drift, &mean_ke_rel_drift));
1641 break;
1642 case SOLUTION_CONVERGENCE_STATISTICAL_STEADY:
1643 PetscCall(ComputeCurrentFlowObservables(simCtx, &mean_speed, &mean_ke));
1644 PetscCall(AppendStatisticalObservableSample(simCtx, samples_before, mean_speed, mean_ke));
1645 PetscCall(ComputeStatisticalWindowMetrics(simCtx, samples_before + 1,
1646 &has_reference,
1647 &mean_speed_window, &mean_speed_window_prev,
1648 &mean_speed_window_abs_drift, &mean_speed_window_rel_drift,
1649 &mean_speed_rms_window, &mean_speed_rms_window_prev,
1650 &mean_speed_rms_window_abs_drift, &mean_speed_rms_window_rel_drift,
1651 &mean_ke_window, &mean_ke_window_prev,
1652 &mean_ke_window_abs_drift, &mean_ke_window_rel_drift,
1653 &mean_ke_rms_window, &mean_ke_rms_window_prev,
1654 &mean_ke_rms_window_abs_drift, &mean_ke_rms_window_rel_drift));
1655 break;
1656 default:
1657 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown solution convergence mode %d.", (int)simCtx->solutionConvergenceMode);
1658 }
1659
1660 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
1661 if (rank == 0) {
1662 char log_path[PETSC_MAX_PATH_LEN + 32];
1663 FILE *f = NULL;
1664 const char *mode_str = SolutionConvergenceModeToString(simCtx->solutionConvergenceMode);
1665
1666 PetscCall(PetscSNPrintf(log_path, sizeof(log_path), "%s/solution_convergence.log", simCtx->log_dir));
1667 f = fopen(log_path, "a");
1668 if (!f) {
1669 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Cannot open solution convergence log: %s", log_path);
1670 }
1671
1672 if (ftell(f) == 0) {
1673 switch (simCtx->solutionConvergenceMode) {
1674 case SOLUTION_CONVERGENCE_STEADY_DETERMINISTIC:
1675 case SOLUTION_CONVERGENCE_TRANSIENT:
1676 fprintf(f, "==================== Solution Convergence Log [mode: %s] ====================\n", mode_str);
1677 /* 16 columns; header width = 314 chars */
1678 fprintf(f, "%-10s | %-18s | %-22s | %-3s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s\n",
1679 "step", "time", "mode", "ref",
1680 "u_abs_l2", "u_rel_l2", "p_abs_l2", "p_rel_l2",
1681 "mean_speed", "spd_ref", "spd_abs", "spd_rel",
1682 "mean_ke", "ke_ref", "ke_abs", "ke_rel");
1683 fprintf(f, "----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n");
1684 break;
1685 case SOLUTION_CONVERGENCE_PERIODIC_DETERMINISTIC:
1686 fprintf(f, "==================== Solution Convergence Log [mode: %s | period_steps: %d] ====================\n",
1687 mode_str, (int)simCtx->solutionConvergencePeriodSteps);
1688 /* 18 columns; header width = 330 chars */
1689 fprintf(f, "%-10s | %-18s | %-22s | %-3s | %-5s | %-5s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s\n",
1690 "step", "time", "mode", "ref", "ph", "per",
1691 "u_abs_l2", "u_rel_l2", "p_abs_l2", "p_rel_l2",
1692 "mean_speed", "spd_ref", "spd_abs", "spd_rel",
1693 "mean_ke", "ke_ref", "ke_abs", "ke_rel");
1694 fprintf(f, "----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n");
1695 break;
1696 case SOLUTION_CONVERGENCE_STATISTICAL_STEADY:
1697 fprintf(f, "==================== Solution Convergence Log [mode: %s | window_steps: %d] ====================\n",
1698 mode_str, (int)simCtx->solutionConvergenceWindowSteps);
1699 /* 21 columns; header width = 406 chars */
1700 fprintf(f, "%-10s | %-18s | %-22s | %-3s | %-5s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s\n",
1701 "step", "time", "mode", "ref", "win",
1702 "mean_speed", "mean_ke",
1703 "spd_win", "spd_win_prev", "spd_win_abs", "spd_win_rel",
1704 "spd_rms_win", "spd_rms_abs", "spd_rms_rel",
1705 "ke_win", "ke_win_prev", "ke_win_abs", "ke_win_rel",
1706 "ke_rms_win", "ke_rms_abs", "ke_rms_rel");
1707 fprintf(f, "------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n");
1708 break;
1709 default: break;
1710 }
1711 }
1712 if (simCtx->continueMode && simCtx->step == simCtx->StartStep + 1) {
1713 fprintf(f, "# Continuation from step %" PetscInt_FMT "\n", simCtx->StartStep);
1714 }
1715
1716 switch (simCtx->solutionConvergenceMode) {
1717 case SOLUTION_CONVERGENCE_STEADY_DETERMINISTIC:
1718 case SOLUTION_CONVERGENCE_TRANSIENT:
1719 fprintf(f,
1720 "%-10d | %-18.10e | %-22s | %-3d | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e\n",
1721 (int)simCtx->step, (double)simCtx->ti, mode_str, has_reference ? 1 : 0,
1722 (double)u_abs_l2, (double)u_rel_l2, (double)p_abs_l2, (double)p_rel_l2,
1723 (double)mean_speed, (double)mean_speed_reference,
1724 (double)mean_speed_abs_drift, (double)mean_speed_rel_drift,
1725 (double)mean_ke, (double)mean_ke_reference,
1726 (double)mean_ke_abs_drift, (double)mean_ke_rel_drift);
1727 break;
1728 case SOLUTION_CONVERGENCE_PERIODIC_DETERMINISTIC:
1729 fprintf(f,
1730 "%-10d | %-18.10e | %-22s | %-3d | %-5d | %-5d | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e\n",
1731 (int)simCtx->step, (double)simCtx->ti, mode_str, has_reference ? 1 : 0,
1732 (int)phase_step, (int)simCtx->solutionConvergencePeriodSteps,
1733 (double)u_abs_l2, (double)u_rel_l2, (double)p_abs_l2, (double)p_rel_l2,
1734 (double)mean_speed, (double)mean_speed_reference,
1735 (double)mean_speed_abs_drift, (double)mean_speed_rel_drift,
1736 (double)mean_ke, (double)mean_ke_reference,
1737 (double)mean_ke_abs_drift, (double)mean_ke_rel_drift);
1738 break;
1739 case SOLUTION_CONVERGENCE_STATISTICAL_STEADY:
1740 fprintf(f,
1741 "%-10d | %-18.10e | %-22s | %-3d | %-5d | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e | %-18.10e\n",
1742 (int)simCtx->step, (double)simCtx->ti, mode_str, has_reference ? 1 : 0,
1743 (int)simCtx->solutionConvergenceWindowSteps,
1744 (double)mean_speed, (double)mean_ke,
1745 (double)mean_speed_window, (double)mean_speed_window_prev,
1746 (double)mean_speed_window_abs_drift, (double)mean_speed_window_rel_drift,
1747 (double)mean_speed_rms_window,
1748 (double)mean_speed_rms_window_abs_drift, (double)mean_speed_rms_window_rel_drift,
1749 (double)mean_ke_window, (double)mean_ke_window_prev,
1750 (double)mean_ke_window_abs_drift, (double)mean_ke_window_rel_drift,
1751 (double)mean_ke_rms_window,
1752 (double)mean_ke_rms_window_abs_drift, (double)mean_ke_rms_window_rel_drift);
1753 break;
1754 default: break;
1755 }
1756 fclose(f);
1757 }
1758
1759 if (simCtx->solutionConvergenceMode == SOLUTION_CONVERGENCE_PERIODIC_DETERMINISTIC &&
1760 phase_step >= 0 && phase_step < simCtx->solutionConvergencePeriodSteps) {
1761 UserCtx *user = simCtx->usermg.mgctx[simCtx->usermg.mglevels - 1].user;
1762 for (PetscInt bi = 0; bi < simCtx->block_number; ++bi) {
1763 PetscCall(VecCopy(user[bi].Ucat, user[bi].solutionConvergencePeriodicUcatRef[phase_step]));
1764 PetscCall(VecCopy(user[bi].P, user[bi].solutionConvergencePeriodicPRef[phase_step]));
1765 }
1766 }
1767
1768 simCtx->solutionConvergenceSamplesRecorded = samples_before + 1;
1769
1770 PetscFunctionReturn(0);
1771}
1772
1773/**
1774 * @brief Logs continuity metrics for a single block to a file.
1775 *
1776 * This function should be called for each block, once per timestep. It opens a
1777 * central log file in append mode. To ensure the header is written only once,
1778 * it checks if it is processing block 0 on the simulation's start step.
1779 *
1780 * @param user A pointer to the UserCtx for the specific block whose metrics
1781 * are to be logged. The function accesses both global (SimCtx)
1782 * and local (user->...) data.
1783 * @return PetscErrorCode 0 on success.
1784 */
1785#undef __FUNCT__
1786#define __FUNCT__ "LOG_CONTINUITY_METRICS"
1787/**
1788 * @brief Implementation of \ref LOG_CONTINUITY_METRICS().
1789 * @details Full API contract (arguments, ownership, side effects) is documented with
1790 * the header declaration in `include/logging.h`.
1791 * @see LOG_CONTINUITY_METRICS()
1792 */
1793
1794PetscErrorCode LOG_CONTINUITY_METRICS(UserCtx *user)
1795{
1796 PetscErrorCode ierr;
1797 PetscMPIInt rank;
1798 SimCtx *simCtx = user->simCtx; // Get the shared SimCtx
1799 const PetscInt bi = user->_this; // Get this block's specific ID
1800 const PetscInt ti = simCtx->step; // Get the current timestep
1801
1802 PetscFunctionBeginUser;
1803 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
1804
1805 // Only rank 0 performs file I/O.
1806 if (!rank) {
1807 FILE *f;
1808 char filen[PETSC_MAX_PATH_LEN + 64];
1809 ierr = PetscSNPrintf(filen, sizeof(filen), "%s/Continuity_Metrics.log", simCtx->log_dir); CHKERRQ(ierr);
1810
1811 // Open the log file in append mode.
1812 f = fopen(filen, "a");
1813 if (!f) {
1814 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Cannot open log file: %s", filen);
1815 }
1816
1817 // Write a header only when the file is empty and it's the first block (bi=0).
1818 // Using ftell() instead of step comparison ensures correctness across continuations.
1819 if (ftell(f) == 0 && bi == 0) {
1820 PetscFPrintf(PETSC_COMM_SELF, f, "%-10s | %-6s | %-18s | %-30s | %-18s | %-18s | %-18s | %-18s\n",
1821 "Timestep", "Block", "Max Divergence", "Max Divergence Location ([k][j][i]=idx)", "Sum(RHS)","Total Flux In", "Total Flux Out", "Net Flux");
1822 PetscFPrintf(PETSC_COMM_SELF, f, "------------------------------------------------------------------------------------------------------------------------------------------\n");
1823 }
1824 if (simCtx->continueMode && ti == simCtx->StartStep + 1 && bi == 0) {
1825 PetscFPrintf(PETSC_COMM_SELF, f, "# Continuation from step %" PetscInt_FMT "\n", simCtx->StartStep);
1826 }
1827
1828 // Prepare the data strings and values for the current block.
1829 PetscReal net_flux = simCtx->FluxInSum - simCtx->FluxOutSum;
1830 char location_str[64];
1831 sprintf(location_str, "([%d][%d][%d] = %d)", (int)simCtx->MaxDivz, (int)simCtx->MaxDivy, (int)simCtx->MaxDivx, (int)simCtx->MaxDivFlatArg);
1832
1833 // Write the formatted line for the current block.
1834 PetscFPrintf(PETSC_COMM_SELF, f, "%-10d | %-6d | %-18.10e | %-39s | %-18.10e | %-18.10e | %-18.10e | %-18.10e\n",
1835 (int)ti,
1836 (int)bi,
1837 (double)simCtx->MaxDiv,
1838 location_str,
1839 (double)simCtx->summationRHS,
1840 (double)simCtx->FluxInSum,
1841 (double)simCtx->FluxOutSum,
1842 (double)net_flux);
1843
1844 fclose(f);
1845 }
1846
1847 PetscFunctionReturn(0);
1848}
1849
1850/**
1851 * @brief Implementation of \ref ParticleLocationStatusToString().
1852 * @details Full API contract (arguments, ownership, side effects) is documented with
1853 * the header declaration in `include/logging.h`.
1854 * @see ParticleLocationStatusToString()
1855 */
1856const char* ParticleLocationStatusToString(ParticleLocationStatus level)
1857{
1858 switch (level) {
1859 case NEEDS_LOCATION: return "NEEDS_LOCATION";
1860 case ACTIVE_AND_LOCATED: return "ACTIVE_AND_LOCATED";
1861 case MIGRATING_OUT: return "MIGRATING_OUT";
1862 case LOST: return "LOST";
1863 case UNINITIALIZED: return "UNINITIALIZED";
1864 default: return "UNKNOWN_LEVEL";
1865 }
1866}
1867
1868///////// Profiling System /////////
1869
1870// Data structure to hold profiling info for one function
1871typedef struct {
1872 const char *name;
1873 double total_time;
1874 double current_step_time;
1875 long long total_call_count;
1876 long long current_step_call_count;
1877 double start_time; // Timer for the current call
1878 PetscBool always_log;
1879} ProfiledFunction;
1880
1881// Global registry for all profiled functions
1882static ProfiledFunction *g_profiler_registry = NULL;
1883static PetscInt g_profiler_count = 0;
1884static PetscInt g_profiler_capacity = 0;
1885
1886// Internal helper to find a function in the registry or create it
1887/**
1888 * @brief Internal helper implementation: `_FindOrCreateEntry()`.
1889 * @details Local to this translation unit.
1890 */
1891static PetscErrorCode _FindOrCreateEntry(const char *func_name, PetscInt *idx)
1892{
1893 PetscFunctionBeginUser;
1894 // Search for existing entry
1895 for (PetscInt i = 0; i < g_profiler_count; ++i) {
1896 if (strcmp(g_profiler_registry[i].name, func_name) == 0) {
1897 *idx = i;
1898 PetscFunctionReturn(0);
1899 }
1900 }
1901
1902 // Not found, create a new entry
1903 if (g_profiler_count >= g_profiler_capacity) {
1904 PetscInt new_capacity = g_profiler_capacity == 0 ? 16 : g_profiler_capacity * 2;
1905 PetscErrorCode ierr = PetscRealloc(sizeof(ProfiledFunction) * new_capacity, &g_profiler_registry); CHKERRQ(ierr);
1906 g_profiler_capacity = new_capacity;
1907 }
1908
1909 *idx = g_profiler_count;
1910 g_profiler_registry[*idx].name = func_name;
1911 g_profiler_registry[*idx].total_time = 0.0;
1912 g_profiler_registry[*idx].current_step_time = 0.0;
1913 g_profiler_registry[*idx].total_call_count = 0;
1914 g_profiler_registry[*idx].current_step_call_count = 0;
1915 g_profiler_registry[*idx].start_time = 0.0;
1916 g_profiler_registry[*idx].always_log = PETSC_FALSE;
1917 g_profiler_count++;
1918
1919 PetscFunctionReturn(0);
1920}
1921
1922// --- Public API Implementation ---
1923/**
1924 * @brief Internal helper implementation: `ProfilingInitialize()`.
1925 * @details Local to this translation unit.
1926 */
1927PetscErrorCode ProfilingInitialize(SimCtx *simCtx)
1928{
1929 PetscFunctionBeginUser;
1930 if (!simCtx) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "SimCtx cannot be null for ProfilingInitialize");
1931
1932 // Iterate through the list of critical functions provided in SimCtx
1933 for (PetscInt i = 0; i < simCtx->nProfilingSelectedFuncs; ++i) {
1934 PetscInt idx;
1935 const char *func_name = simCtx->profilingSelectedFuncs[i];
1936 PetscErrorCode ierr = _FindOrCreateEntry(func_name, &idx); CHKERRQ(ierr);
1937 g_profiler_registry[idx].always_log = PETSC_TRUE;
1938
1939 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Marked '%s' as a critical function for profiling.\n", func_name);
1940 }
1941 PetscFunctionReturn(0);
1942}
1943
1944/**
1945 * @brief Implementation of \ref _ProfilingStart().
1946 * @details Full API contract (arguments, ownership, side effects) is documented with
1947 * the header declaration in `include/logging.h`.
1948 * @see _ProfilingStart()
1949 */
1950
1951void _ProfilingStart(const char *func_name)
1952{
1953 PetscInt idx;
1954 if (_FindOrCreateEntry(func_name, &idx) != 0) return; // Fail silently
1955 PetscTime(&g_profiler_registry[idx].start_time);
1956}
1957
1958/**
1959 * @brief Implementation of \ref _ProfilingEnd().
1960 * @details Full API contract (arguments, ownership, side effects) is documented with
1961 * the header declaration in `include/logging.h`.
1962 * @see _ProfilingEnd()
1963 */
1964
1965void _ProfilingEnd(const char *func_name)
1966{
1967 double end_time;
1968 PetscTime(&end_time);
1969
1970 PetscInt idx;
1971 if (_FindOrCreateEntry(func_name, &idx) != 0) return; // Fail silently
1972
1973 double elapsed = end_time - g_profiler_registry[idx].start_time;
1974 g_profiler_registry[idx].total_time += elapsed;
1975 g_profiler_registry[idx].current_step_time += elapsed;
1976 g_profiler_registry[idx].total_call_count++;
1977 g_profiler_registry[idx].current_step_call_count++;
1978}
1979
1980/**
1981 * @brief Implementation of \ref ProfilingResetTimestepCounters().
1982 * @details Full API contract (arguments, ownership, side effects) is documented with
1983 * the header declaration in `include/logging.h`.
1984 * @see ProfilingResetTimestepCounters()
1985 */
1986
1987PetscErrorCode ProfilingResetTimestepCounters(void)
1988{
1989 PetscFunctionBeginUser;
1990 for (PetscInt i = 0; i < g_profiler_count; ++i) {
1991 g_profiler_registry[i].current_step_time = 0.0;
1992 g_profiler_registry[i].current_step_call_count = 0;
1993 }
1994 PetscFunctionReturn(0);
1995}
1996
1997/**
1998 * @brief Implementation of \ref ProfilingLogTimestepSummary().
1999 * @details Full API contract (arguments, ownership, side effects) is documented with
2000 * the header declaration in `include/logging.h`.
2001 * @see ProfilingLogTimestepSummary()
2002 */
2003
2004PetscErrorCode ProfilingLogTimestepSummary(SimCtx *simCtx, PetscInt step)
2005{
2006 PetscBool should_write = PETSC_FALSE;
2007 FILE *f = NULL;
2008 char filen[(2 * PETSC_MAX_PATH_LEN) + 16];
2009
2010 PetscFunctionBeginUser;
2011 if (!simCtx) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "SimCtx cannot be null for ProfilingLogTimestepSummary");
2012
2013 if (strcmp(simCtx->profilingTimestepMode, "off") == 0) {
2014 for (PetscInt i = 0; i < g_profiler_count; ++i) {
2015 g_profiler_registry[i].current_step_time = 0.0;
2016 g_profiler_registry[i].current_step_call_count = 0;
2017 }
2018 PetscFunctionReturn(0);
2019 }
2020
2021 for (PetscInt i = 0; i < g_profiler_count; ++i) {
2022 if (g_profiler_registry[i].current_step_call_count <= 0) {
2023 continue;
2024 }
2025 if (strcmp(simCtx->profilingTimestepMode, "all") == 0 || g_profiler_registry[i].always_log) {
2026 should_write = PETSC_TRUE;
2027 break;
2028 }
2029 }
2030
2031 if (should_write && simCtx->rank == 0) {
2032 snprintf(filen, sizeof(filen), "%s/%s", simCtx->log_dir, simCtx->profilingTimestepFile);
2033 if (step == simCtx->StartStep + 1 && !simCtx->continueMode) {
2034 f = fopen(filen, "w");
2035 if (!f) {
2036 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Cannot open profiling timestep log file: %s", filen);
2037 }
2038 PetscFPrintf(PETSC_COMM_SELF, f, "step,function,calls,step_time_s\n");
2039 } else {
2040 f = fopen(filen, "a");
2041 if (!f) {
2042 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Cannot open profiling timestep log file: %s", filen);
2043 }
2044 if (step == simCtx->StartStep + 1 && ftell(f) == 0) {
2045 PetscFPrintf(PETSC_COMM_SELF, f, "step,function,calls,step_time_s\n");
2046 }
2047 }
2048 if (simCtx->continueMode && step == simCtx->StartStep + 1) {
2049 PetscFPrintf(PETSC_COMM_SELF, f, "# Continuation from step %" PetscInt_FMT "\n", simCtx->StartStep);
2050 }
2051
2052 for (PetscInt i = 0; i < g_profiler_count; ++i) {
2053 if (g_profiler_registry[i].current_step_call_count <= 0) {
2054 continue;
2055 }
2056 if (strcmp(simCtx->profilingTimestepMode, "all") == 0 || g_profiler_registry[i].always_log) {
2057 PetscFPrintf(
2058 PETSC_COMM_SELF,
2059 f,
2060 "%d,%s,%lld,%.6f\n",
2061 (int)step,
2062 g_profiler_registry[i].name,
2063 g_profiler_registry[i].current_step_call_count,
2064 g_profiler_registry[i].current_step_time
2065 );
2066 }
2067 }
2068 fclose(f);
2069 }
2070
2071 // Reset per-step counters for the next iteration
2072 for (PetscInt i = 0; i < g_profiler_count; ++i) {
2073 g_profiler_registry[i].current_step_time = 0.0;
2074 g_profiler_registry[i].current_step_call_count = 0;
2075 }
2076 PetscFunctionReturn(0);
2077}
2078
2079/**
2080 * @brief Implementation of \ref RuntimeMemoryLogSample().
2081 * @details Full API contract (arguments, ownership, side effects) is documented with
2082 * the header declaration in `include/logging.h`.
2083 * @see RuntimeMemoryLogSample()
2084 */
2085PetscErrorCode RuntimeMemoryLogSample(SimCtx *simCtx, PetscInt step, const char *event, const char *reason)
2086{
2087 PetscErrorCode ierr;
2088 PetscLogDouble process_current_bytes = 0.0;
2089 PetscLogDouble process_peak_bytes = 0.0;
2090 PetscLogDouble petsc_current_bytes = 0.0;
2091 PetscLogDouble petsc_peak_bytes = 0.0;
2092 PetscReal local_values[5];
2093 PetscReal global_values[5];
2094 PetscReal process_current_mb = 0.0;
2095 PetscReal process_peak_mb = 0.0;
2096 PetscReal petsc_current_mb = 0.0;
2097 PetscReal petsc_peak_mb = 0.0;
2098 PetscReal process_change_mb = 0.0;
2099 char path[(2 * PETSC_MAX_PATH_LEN) + 16];
2100 FILE *f = NULL;
2101
2102 PetscFunctionBeginUser;
2103 if (!simCtx) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "SimCtx cannot be null for RuntimeMemoryLogSample");
2104 if (!simCtx->runtimeMemoryLogEnabled) PetscFunctionReturn(0);
2105
2106 ierr = PetscMemoryGetCurrentUsage(&process_current_bytes); CHKERRQ(ierr);
2107 ierr = PetscMemoryGetMaximumUsage(&process_peak_bytes); CHKERRQ(ierr);
2108 ierr = PetscMallocGetCurrentUsage(&petsc_current_bytes); CHKERRQ(ierr);
2109 ierr = PetscMallocGetMaximumUsage(&petsc_peak_bytes); CHKERRQ(ierr);
2110
2111 process_current_mb = (PetscReal)(process_current_bytes / (1024.0 * 1024.0));
2112 process_peak_mb = (PetscReal)(process_peak_bytes / (1024.0 * 1024.0));
2113 petsc_current_mb = (PetscReal)(petsc_current_bytes / (1024.0 * 1024.0));
2114 petsc_peak_mb = (PetscReal)(petsc_peak_bytes / (1024.0 * 1024.0));
2115 if (simCtx->runtimeMemoryLogHasPrevious) {
2116 process_change_mb = process_current_mb - simCtx->runtimeMemoryLogPreviousProcessMB;
2117 }
2118
2119 local_values[0] = process_current_mb;
2120 local_values[1] = process_peak_mb;
2121 local_values[2] = petsc_current_mb;
2122 local_values[3] = petsc_peak_mb;
2123 local_values[4] = process_change_mb;
2124 ierr = MPI_Allreduce(local_values, global_values, 5, MPIU_REAL, MPI_MAX, PETSC_COMM_WORLD); CHKERRMPI(ierr);
2125
2126 simCtx->runtimeMemoryLogPreviousProcessMB = process_current_mb;
2127 simCtx->runtimeMemoryLogHasPrevious = PETSC_TRUE;
2128
2129 if (simCtx->rank == 0) {
2130 ierr = PetscSNPrintf(path, sizeof(path), "%s/%s", simCtx->log_dir, simCtx->runtimeMemoryLogFile); CHKERRQ(ierr);
2131 f = fopen(path, "a");
2132 if (!f) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Cannot open runtime memory log file: %s", path);
2133
2134 if (!simCtx->runtimeMemoryLogStarted) {
2135 fprintf(f, "# PICurv runtime memory log\n");
2136 if (simCtx->continueMode) {
2137 fprintf(f, "# Continuation from step %" PetscInt_FMT "\n", simCtx->StartStep);
2138 }
2139 fprintf(
2140 f,
2141 "%-8s %-10s %22s %20s %22s %28s %22s %-18s\n",
2142 "Step",
2143 "Event",
2144 "Process Current MB Max",
2145 "Process Peak MB Max",
2146 "PETSc Allocated MB Max",
2147 "PETSc Peak Allocated MB Max",
2148 "Process Change MB Max",
2149 "Reason"
2150 );
2151 simCtx->runtimeMemoryLogStarted = PETSC_TRUE;
2152 }
2153
2154 fprintf(
2155 f,
2156 "%-8" PetscInt_FMT " %-10s %22.3f %20.3f %22.3f %28.3f %22.3f %-18s\n",
2157 step,
2158 event ? event : "-",
2159 (double)global_values[0],
2160 (double)global_values[1],
2161 (double)global_values[2],
2162 (double)global_values[3],
2163 (double)global_values[4],
2164 (reason && reason[0]) ? reason : "-"
2165 );
2166 if ((event && (strcmp(event, "Shutdown") == 0 || strcmp(event, "Final") == 0))) {
2167 fflush(f);
2168 }
2169 fclose(f);
2170 }
2171
2172 PetscFunctionReturn(0);
2173}
2174
2175// Comparison function for qsort to sort by total_time in descending order
2176/**
2177 * @brief Internal helper implementation: `_CompareProfiledFunctions()`.
2178 * @details Local to this translation unit.
2179 */
2180static int _CompareProfiledFunctions(const void *a, const void *b)
2181{
2182 const ProfiledFunction *func_a = (const ProfiledFunction *)a;
2183 const ProfiledFunction *func_b = (const ProfiledFunction *)b;
2184
2185 if (func_a->total_time < func_b->total_time) return 1;
2186 if (func_a->total_time > func_b->total_time) return -1;
2187 return 0;
2188}
2189
2190/**
2191 * @brief Implementation of \ref ProfilingFinalize().
2192 * @details Full API contract (arguments, ownership, side effects) is documented with
2193 * the header declaration in `include/logging.h`.
2194 * @see ProfilingFinalize()
2195 */
2196PetscErrorCode ProfilingFinalize(SimCtx *simCtx)
2197{
2198 PetscErrorCode ierr;
2199 PetscInt rank = simCtx->rank;
2200 PetscFunctionBeginUser;
2201 if (!simCtx->profilingFinalSummary) PetscFunctionReturn(0);
2202 if (!rank) {
2203
2204 char exec_mode_modifier[32] = "Unknown";
2205 if(simCtx->exec_mode == EXEC_MODE_SOLVER) PetscCall(PetscStrncpy(exec_mode_modifier, "Solver", sizeof(exec_mode_modifier)));
2206 else if(simCtx->exec_mode == EXEC_MODE_POSTPROCESSOR) PetscCall(PetscStrncpy(exec_mode_modifier, "PostProcessor", sizeof(exec_mode_modifier)));
2207 //--- Step 0: Create a file viewer for log file
2208 FILE *f;
2209 char filen[PETSC_MAX_PATH_LEN + 128];
2210 ierr = PetscSNPrintf(filen, sizeof(filen), "%s/ProfilingSummary_%s.log",simCtx->log_dir,exec_mode_modifier); CHKERRQ(ierr);
2211
2212 // Open the log file: append with section label in continue mode, truncate otherwise.
2213 if (simCtx->continueMode) {
2214 f = fopen(filen, "a");
2215 if (!f) {
2216 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Cannot open log file: %s", filen);
2217 }
2218 fprintf(f, "\n=== Continuation from step %" PetscInt_FMT " ===\n", simCtx->StartStep);
2219 } else {
2220 f = fopen(filen, "w");
2221 if (!f) {
2222 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Cannot open log file: %s", filen);
2223 }
2224 }
2225
2226 // --- Step 1: Sort the data for readability ---
2227 qsort(g_profiler_registry, g_profiler_count, sizeof(ProfiledFunction), _CompareProfiledFunctions);
2228
2229 // --- Step 2: Dynamically determine the width for the function name column ---
2230 PetscInt max_name_len = strlen("Function"); // Start with the header's length
2231 for (PetscInt i = 0; i < g_profiler_count; ++i) {
2232 if (g_profiler_registry[i].total_call_count > 0) {
2233 PetscInt len = strlen(g_profiler_registry[i].name);
2234 if (len > max_name_len) {
2235 max_name_len = len;
2236 }
2237 }
2238 }
2239 // Add a little padding
2240 max_name_len += 2;
2241
2242 // --- Step 3: Define fixed widths for numeric columns for consistent alignment ---
2243 const int time_width = 18;
2244 const int count_width = 15;
2245 const int avg_width = 22;
2246
2247 // --- Step 4: Print the formatted table ---
2248 PetscFPrintf(PETSC_COMM_SELF, f, "=================================================================================================================\n");
2249 PetscFPrintf(PETSC_COMM_SELF, f, " FINAL PROFILING SUMMARY (Sorted by Total Time)\n");
2250 PetscFPrintf(PETSC_COMM_SELF, f, "=================================================================================================================\n");
2251
2252 // Header Row
2253 PetscFPrintf(PETSC_COMM_SELF, f, "%-*s | %-*s | %-*s | %-*s\n",
2254 max_name_len, "Function",
2255 time_width, "Total Time (s)",
2256 count_width, "Call Count",
2257 avg_width, "Avg. Time/Call (ms)");
2258
2259 // Separator Line (dynamically sized)
2260 for (int i = 0; i < max_name_len; i++) PetscFPrintf(PETSC_COMM_SELF, f, "-");
2261 PetscFPrintf(PETSC_COMM_SELF, f, "-|-");
2262 for (int i = 0; i < time_width; i++) PetscFPrintf(PETSC_COMM_SELF, f, "-");
2263 PetscFPrintf(PETSC_COMM_SELF, f, "-|-");
2264 for (int i = 0; i < count_width; i++) PetscFPrintf(PETSC_COMM_SELF, f, "-");
2265 PetscFPrintf(PETSC_COMM_SELF, f, "-|-");
2266 for (int i = 0; i < avg_width; i++) PetscFPrintf(PETSC_COMM_SELF, f, "-");
2267 PetscFPrintf(PETSC_COMM_SELF, f, "\n");
2268
2269 // Data Rows
2270 for (PetscInt i = 0; i < g_profiler_count; ++i) {
2271 if (g_profiler_registry[i].total_call_count > 0) {
2272 double avg_time_ms = (g_profiler_registry[i].total_time / g_profiler_registry[i].total_call_count) * 1000.0;
2273 PetscFPrintf(PETSC_COMM_SELF, f, "%-*s | %*.*f | %*lld | %*.*f\n",
2274 max_name_len, g_profiler_registry[i].name,
2275 time_width, 6, g_profiler_registry[i].total_time,
2276 count_width, g_profiler_registry[i].total_call_count,
2277 avg_width, 6, avg_time_ms);
2278 PetscFPrintf(PETSC_COMM_SELF, f, "------------------------------------------------------------------------------------------------------------------\n");
2279 }
2280 }
2281 PetscFPrintf(PETSC_COMM_SELF, f, "==================================================================================================================\n");
2282
2283 fclose(f);
2284 }
2285
2286 // --- Final Cleanup ---
2287 PetscFree(g_profiler_registry);
2288 g_profiler_registry = NULL;
2289 g_profiler_count = 0;
2290 g_profiler_capacity = 0;
2291 PetscFunctionReturn(0);
2292}
2293
2294/*================================================================================*
2295 * PROGRESS BAR UTILITY *
2296 *================================================================================*/
2297
2298/**
2299 * @brief Internal helper implementation: `PrintProgressBar()`.
2300 * @details Local to this translation unit.
2301 */
2302void PrintProgressBar(PetscInt step, PetscInt startStep, PetscInt totalSteps, PetscReal currentTime)
2303{
2304 if (totalSteps <= 0) return;
2305
2306 // --- Configuration ---
2307 const int barWidth = 50;
2308
2309 // --- Calculation ---
2310 // Calculate progress as a fraction from 0.0 to 1.0
2311 PetscReal progress = (PetscReal)(step - startStep + 1) / totalSteps;
2312 // Ensure progress doesn't exceed 1.0 due to floating point inaccuracies
2313 if (progress > 1.0) progress = 1.0;
2314
2315 int pos = (int)(barWidth * progress);
2316
2317 // --- Printing ---
2318 // Carriage return moves cursor to the beginning of the line
2319 PetscPrintf(PETSC_COMM_SELF, "\rProgress: [");
2320
2321 for (int i = 0; i < barWidth; ++i) {
2322 if (i < pos) {
2323 PetscPrintf(PETSC_COMM_SELF, "=");
2324 } else if (i == pos) {
2325 PetscPrintf(PETSC_COMM_SELF, ">");
2326 } else {
2327 PetscPrintf(PETSC_COMM_SELF, " ");
2328 }
2329 }
2330
2331 // Print percentage, step count, and current time
2332 PetscPrintf(PETSC_COMM_SELF, "] %3d%% (Step %" PetscInt_FMT "/%" PetscInt_FMT ", t=%.4f)",
2333 (int)(progress * 100.0),
2334 step + 1,
2335 startStep + totalSteps,
2336 currentTime);
2337
2338 // Flush the output buffer to ensure the bar is displayed immediately
2339 fflush(stdout);
2340}
2341
2342#undef __FUNCT__
2343#define __FUNCT__ "LOG_FIELD_MIN_MAX"
2344/**
2345 * @brief Implementation of \ref LOG_FIELD_MIN_MAX().
2346 * @details Full API contract is documented with the header declaration in `include/logging.h`.
2347 * @see LOG_FIELD_MIN_MAX()
2348 */
2349PetscErrorCode LOG_FIELD_MIN_MAX(UserCtx *user, const char *fieldName)
2350{
2351 PetscErrorCode ierr;
2352 PetscInt i, j, k;
2353 DMDALocalInfo info;
2354
2355 Vec fieldVec = NULL;
2356 DM dm = NULL;
2357 PetscInt dof;
2358 char data_layout[20];
2359
2360 PetscFunctionBeginUser;
2361
2362 // --- 1. Map string name to PETSc objects and determine data layout ---
2363 if (strcasecmp(fieldName, "Ucat") == 0) {
2364 fieldVec = user->Ucat; dm = user->fda; dof = 3; strcpy(data_layout, "Cell-Centered");
2365 } else if (strcasecmp(fieldName, "P") == 0) {
2366 fieldVec = user->P; dm = user->da; dof = 1; strcpy(data_layout, "Cell-Centered");
2367 } else if (strcasecmp(fieldName, "Diffusivity") == 0) {
2368 fieldVec = user->Diffusivity; dm = user->da; dof = 1; strcpy(data_layout, "Cell-Centered");
2369 } else if (strcasecmp(fieldName, "DiffusivityGradient") == 0) {
2370 fieldVec = user->DiffusivityGradient; dm = user->fda; dof = 3; strcpy(data_layout, "Cell-Centered");
2371 } else if (strcasecmp(fieldName, "Phi") == 0) {
2372 fieldVec = user->Phi; dm = user->da; dof = 1; strcpy(data_layout, "Cell-Centered");
2373 } else if (strcasecmp(fieldName, "Nvert") == 0) {
2374 fieldVec = user->Nvert; dm = user->da; dof = 1; strcpy(data_layout, "Cell-Centered");
2375 } else if (strcasecmp(fieldName, "Aj") == 0) {
2376 fieldVec = user->Aj; dm = user->da; dof = 1; strcpy(data_layout, "Cell-Centered");
2377 } else if (strcasecmp(fieldName, "Cent") == 0 || strcasecmp(fieldName, "Center-Coordinates") == 0) {
2378 fieldVec = user->Cent; dm = user->fda; dof = 3; strcpy(data_layout, "Cell-Centered");
2379 } else if (strcasecmp(fieldName, "Ucont") == 0) {
2380 fieldVec = user->lUcont; dm = user->fda; dof = 3; strcpy(data_layout, "Face-Centered");
2381 } else if (strcasecmp(fieldName, "Centx") == 0 || strcasecmp(fieldName, "X-Face-Centers") == 0) {
2382 fieldVec = user->Centx; dm = user->fda; dof = 3; strcpy(data_layout, "I-Face");
2383 } else if (strcasecmp(fieldName, "Centy") == 0 || strcasecmp(fieldName, "Y-Face-Centers") == 0) {
2384 fieldVec = user->Centy; dm = user->fda; dof = 3; strcpy(data_layout, "J-Face");
2385 } else if (strcasecmp(fieldName, "Centz") == 0 || strcasecmp(fieldName, "Z-Face-Centers") == 0) {
2386 fieldVec = user->Centz; dm = user->fda; dof = 3; strcpy(data_layout, "K-Face");
2387 } else if (strcasecmp(fieldName, "Csi") == 0 || strcasecmp(fieldName, "ICsi") == 0 ||
2388 strcasecmp(fieldName, "IEta") == 0 || strcasecmp(fieldName, "IZet") == 0) {
2389 fieldVec = strcasecmp(fieldName, "Csi") == 0 ? user->Csi :
2390 (strcasecmp(fieldName, "ICsi") == 0 ? user->ICsi :
2391 (strcasecmp(fieldName, "IEta") == 0 ? user->IEta : user->IZet));
2392 dm = user->fda; dof = 3; strcpy(data_layout, "I-Face");
2393 } else if (strcasecmp(fieldName, "Eta") == 0 || strcasecmp(fieldName, "JCsi") == 0 ||
2394 strcasecmp(fieldName, "JEta") == 0 || strcasecmp(fieldName, "JZet") == 0) {
2395 fieldVec = strcasecmp(fieldName, "Eta") == 0 ? user->Eta :
2396 (strcasecmp(fieldName, "JCsi") == 0 ? user->JCsi :
2397 (strcasecmp(fieldName, "JEta") == 0 ? user->JEta : user->JZet));
2398 dm = user->fda; dof = 3; strcpy(data_layout, "J-Face");
2399 } else if (strcasecmp(fieldName, "Zet") == 0 || strcasecmp(fieldName, "KCsi") == 0 ||
2400 strcasecmp(fieldName, "KEta") == 0 || strcasecmp(fieldName, "KZet") == 0) {
2401 fieldVec = strcasecmp(fieldName, "Zet") == 0 ? user->Zet :
2402 (strcasecmp(fieldName, "KCsi") == 0 ? user->KCsi :
2403 (strcasecmp(fieldName, "KEta") == 0 ? user->KEta : user->KZet));
2404 dm = user->fda; dof = 3; strcpy(data_layout, "K-Face");
2405 } else if (strcasecmp(fieldName, "IAj") == 0 || strcasecmp(fieldName, "JAj") == 0 ||
2406 strcasecmp(fieldName, "KAj") == 0) {
2407 fieldVec = strcasecmp(fieldName, "IAj") == 0 ? user->IAj :
2408 (strcasecmp(fieldName, "JAj") == 0 ? user->JAj : user->KAj);
2409 dm = user->da; dof = 1;
2410 strcpy(data_layout, strcasecmp(fieldName, "IAj") == 0 ? "I-Face" :
2411 (strcasecmp(fieldName, "JAj") == 0 ? "J-Face" : "K-Face"));
2412 } else if (strcasecmp(fieldName, "Coordinates") == 0) {
2413 ierr = DMGetCoordinates(user->da, &fieldVec); CHKERRQ(ierr);
2414 dm = user->fda; dof = 3; strcpy(data_layout, "Node-Centered");
2415 } else if (strcasecmp(fieldName,"Psi") == 0) {
2416 fieldVec = user->Psi; dm = user->da; dof = 1; strcpy(data_layout, "Cell-Centered"); // Assuming Psi is cell-centered
2417 } else {
2418 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Field %s not recognized.", fieldName);
2419 }
2420
2421 if (!fieldVec) {
2422 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vector for field '%s' is NULL.", fieldName);
2423 }
2424 if (!dm) {
2425 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "DM for field '%s' is NULL.", fieldName);
2426 }
2427
2428 ierr = DMDAGetLocalInfo(dm, &info); CHKERRQ(ierr);
2429
2430 // --- 2. Define Architecture-Aware Loop Bounds ---
2431 PetscInt i_start, i_end, j_start, j_end, k_start, k_end;
2432
2433 if (strcmp(data_layout, "Cell-Centered") == 0) {
2434 // For cell-centered data, the physical values are stored from index 1 to N-1.
2435 // We find the intersection of the rank's owned range [xs, xe) with the
2436 // physical data range [1, IM-1).
2437 i_start = PetscMax(info.xs, 1); i_end = PetscMin(info.xs + info.xm, user->IM);
2438 j_start = PetscMax(info.ys, 1); j_end = PetscMin(info.ys + info.ym, user->JM);
2439 k_start = PetscMax(info.zs, 1); k_end = PetscMin(info.zs + info.zm, user->KM);
2440 } else { // For Node- or Face-Centered data
2441 // The physical values are stored from index 0 to N-1.
2442 // We find the intersection of the rank's owned range [xs, xe) with the
2443 // physical data range [0, IM-1].
2444 i_start = PetscMax(info.xs, 0); i_end = PetscMin(info.xs + info.xm, user->IM);
2445 j_start = PetscMax(info.ys, 0); j_end = PetscMin(info.ys + info.ym, user->JM);
2446 k_start = PetscMax(info.zs, 0); k_end = PetscMin(info.zs + info.zm, user->KM);
2447 }
2448
2449 // --- 3. Barrier for clean, grouped output ---
2450 ierr = MPI_Barrier(PETSC_COMM_WORLD); CHKERRQ(ierr);
2451 if (user->simCtx->rank == 0) {
2452 PetscPrintf(PETSC_COMM_SELF, "\n--- Field Ranges: [%s] (Layout: %s) ---\n", fieldName, data_layout);
2453 }
2454
2455 // --- 4. Branch on DoF and perform calculation with correct bounds ---
2456 if (dof == 1) {
2457 PetscReal localMin = PETSC_MAX_REAL, localMax = PETSC_MIN_REAL;
2458 PetscReal globalMin, globalMax;
2459 const PetscScalar ***array;
2460
2461 ierr = DMDAVecGetArrayRead(dm, fieldVec, &array); CHKERRQ(ierr);
2462 for (k = k_start; k < k_end; k++) {
2463 for (j = j_start; j < j_end; j++) {
2464 for (i = i_start; i < i_end; i++) {
2465 localMin = PetscMin(localMin, array[k][j][i]);
2466 localMax = PetscMax(localMax, array[k][j][i]);
2467 }
2468 }
2469 }
2470 ierr = DMDAVecRestoreArrayRead(dm, fieldVec, &array); CHKERRQ(ierr);
2471
2472 ierr = MPI_Allreduce(&localMin, &globalMin, 1, MPIU_REAL, MPI_MIN, PETSC_COMM_WORLD); CHKERRQ(ierr);
2473 ierr = MPI_Allreduce(&localMax, &globalMax, 1, MPIU_REAL, MPI_MAX, PETSC_COMM_WORLD); CHKERRQ(ierr);
2474
2475 PetscSynchronizedPrintf(PETSC_COMM_WORLD, " [Rank %d] Local Range: [ %11.4e , %11.4e ]\n", user->simCtx->rank, localMin, localMax);
2476 ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT); CHKERRQ(ierr);
2477 if (user->simCtx->rank == 0) {
2478 PetscPrintf(PETSC_COMM_SELF, " Global Range: [ %11.4e , %11.4e ]\n", globalMin, globalMax);
2479 }
2480
2481 } else if (dof == 3) {
2482 Cmpnts localMin = {PETSC_MAX_REAL, PETSC_MAX_REAL, PETSC_MAX_REAL};
2483 Cmpnts localMax = {PETSC_MIN_REAL, PETSC_MIN_REAL, PETSC_MIN_REAL};
2484 Cmpnts globalMin, globalMax;
2485 const Cmpnts ***array;
2486
2487 ierr = DMDAVecGetArrayRead(dm, fieldVec, &array); CHKERRQ(ierr);
2488 for (k = k_start; k < k_end; k++) {
2489 for (j = j_start; j < j_end; j++) {
2490 for (i = i_start; i < i_end; i++) {
2491 localMin.x = PetscMin(localMin.x, array[k][j][i].x);
2492 localMin.y = PetscMin(localMin.y, array[k][j][i].y);
2493 localMin.z = PetscMin(localMin.z, array[k][j][i].z);
2494 localMax.x = PetscMax(localMax.x, array[k][j][i].x);
2495 localMax.y = PetscMax(localMax.y, array[k][j][i].y);
2496 localMax.z = PetscMax(localMax.z, array[k][j][i].z);
2497 }
2498 }
2499 }
2500 ierr = DMDAVecRestoreArrayRead(dm, fieldVec, &array); CHKERRQ(ierr);
2501
2502 ierr = MPI_Allreduce(&localMin, &globalMin, 3, MPIU_REAL, MPI_MIN, PETSC_COMM_WORLD); CHKERRQ(ierr);
2503 ierr = MPI_Allreduce(&localMax, &globalMax, 3, MPIU_REAL, MPI_MAX, PETSC_COMM_WORLD); CHKERRQ(ierr);
2504
2505 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, " [Rank %d] Local X-Range: [ %11.4e , %11.4e ]\n", user->simCtx->rank, localMin.x, localMax.x);
2506 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, " [Rank %d] Local Y-Range: [ %11.4e , %11.4e ]\n", user->simCtx->rank, localMin.y, localMax.y);
2507 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, " [Rank %d] Local Z-Range: [ %11.4e , %11.4e ]\n", user->simCtx->rank, localMin.z, localMax.z);
2508 ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT); CHKERRQ(ierr);
2509
2510 if (user->simCtx->rank == 0) {
2511 PetscPrintf(PETSC_COMM_SELF, " [Global] X-Range: [ %11.4e , %11.4e ]\n", globalMin.x, globalMax.x);
2512 PetscPrintf(PETSC_COMM_SELF, " [Global] Y-Range: [ %11.4e , %11.4e ]\n", globalMin.y, globalMax.y);
2513 PetscPrintf(PETSC_COMM_SELF, " [Global] Z-Range: [ %11.4e , %11.4e ]\n", globalMin.z, globalMax.z);
2514 }
2515
2516 } else {
2517 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "LogFieldStatistics only supports fields with 1 or 3 components, but field '%s' has %" PetscInt_FMT ".", fieldName, dof);
2518 }
2519
2520 // --- 5. Final barrier for clean output ordering ---
2521 ierr = MPI_Barrier(PETSC_COMM_WORLD); CHKERRQ(ierr);
2522 if (user->simCtx->rank == 0) {
2523 PetscPrintf(PETSC_COMM_SELF, "--------------------------------------------\n\n");
2524 }
2525
2526 PetscFunctionReturn(0);
2527}
2528
2529#undef __FUNCT__
2530#define __FUNCT__ "LOG_FIELD_ANATOMY"
2531/**
2532 * @brief Implementation of \ref LOG_FIELD_ANATOMY().
2533 * @details Full API contract is documented with the header declaration in `include/logging.h`.
2534 * @see LOG_FIELD_ANATOMY()
2535 */
2536PetscErrorCode LOG_FIELD_ANATOMY(UserCtx *user, const char *field_name, const char *stage_name)
2537{
2538 PetscErrorCode ierr;
2539 DMDALocalInfo info;
2540 PetscMPIInt rank;
2541
2542 Vec vec_local = NULL;
2543 DM dm = NULL;
2544 PetscInt dof = 0;
2545 char data_layout[20];
2546 char dominant_dir = '\0'; // 'x', 'y', 'z' for face-centered, 'm' for mixed (Ucont)
2547
2548 PetscFunctionBeginUser;
2549 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
2550
2551 // --- 1. Map string name to PETSc objects and determine data layout ---
2552 if (strcasecmp(field_name, "Ucat") == 0) {
2553 vec_local = user->lUcat; dm = user->fda; dof = 3; strcpy(data_layout, "Cell-Centered");
2554 } else if (strcasecmp(field_name, "P") == 0) {
2555 vec_local = user->lP; dm = user->da; dof = 1; strcpy(data_layout, "Cell-Centered");
2556 } else if (strcasecmp(field_name, "Diffusivity") == 0) {
2557 vec_local = user->lDiffusivity; dm = user->da; dof = 1; strcpy(data_layout, "Cell-Centered");
2558 } else if (strcasecmp(field_name, "DiffusivityGradient") == 0) {
2559 vec_local = user->lDiffusivityGradient; dm = user->fda; dof = 3; strcpy(data_layout, "Cell-Centered");
2560 } else if (strcasecmp(field_name, "Psi") == 0) {
2561 vec_local = user->lPsi; dm = user->da; dof = 1; strcpy(data_layout, "Cell-Centered");
2562 } else if (strcasecmp(field_name, "Center-Coordinates") == 0) {
2563 vec_local = user->lCent; dm = user->fda; dof = 3; strcpy(data_layout, "Cell-Centered");
2564 } else if (strcasecmp(field_name, "Ucont") == 0 ||
2565 strcasecmp(field_name, "Ucont_o") == 0 ||
2566 strcasecmp(field_name, "Ucont_rm1") == 0) {
2567 vec_local = strcasecmp(field_name, "Ucont") == 0 ? user->lUcont :
2568 (strcasecmp(field_name, "Ucont_o") == 0 ? user->lUcont_o : user->lUcont_rm1);
2569 dm = user->fda; dof = 3; strcpy(data_layout, "Face-Centered"); dominant_dir = 'm'; // Mixed
2570 } else if (strcasecmp(field_name, "Csi") == 0 || strcasecmp(field_name, "X-Face-Centers") == 0 ||
2571 strcasecmp(field_name, "ICsi") == 0 || strcasecmp(field_name, "IEta") == 0 ||
2572 strcasecmp(field_name, "IZet") == 0) {
2573 vec_local = strcasecmp(field_name, "Csi") == 0 ? user->lCsi :
2574 (strcasecmp(field_name, "X-Face-Centers") == 0 ? user->lCentx :
2575 (strcasecmp(field_name, "ICsi") == 0 ? user->lICsi :
2576 (strcasecmp(field_name, "IEta") == 0 ? user->lIEta : user->lIZet)));
2577 dm = user->fda; dof = 3; strcpy(data_layout, "Face-Centered"); dominant_dir = 'x';
2578 } else if (strcasecmp(field_name, "Eta") == 0 || strcasecmp(field_name, "Y-Face-Centers") == 0 ||
2579 strcasecmp(field_name, "JCsi") == 0 || strcasecmp(field_name, "JEta") == 0 ||
2580 strcasecmp(field_name, "JZet") == 0) {
2581 vec_local = strcasecmp(field_name, "Eta") == 0 ? user->lEta :
2582 (strcasecmp(field_name, "Y-Face-Centers") == 0 ? user->lCenty :
2583 (strcasecmp(field_name, "JCsi") == 0 ? user->lJCsi :
2584 (strcasecmp(field_name, "JEta") == 0 ? user->lJEta : user->lJZet)));
2585 dm = user->fda; dof = 3; strcpy(data_layout, "Face-Centered"); dominant_dir = 'y';
2586 } else if (strcasecmp(field_name, "Zet") == 0 || strcasecmp(field_name, "Z-Face-Centers") == 0 ||
2587 strcasecmp(field_name, "KCsi") == 0 || strcasecmp(field_name, "KEta") == 0 ||
2588 strcasecmp(field_name, "KZet") == 0) {
2589 vec_local = strcasecmp(field_name, "Zet") == 0 ? user->lZet :
2590 (strcasecmp(field_name, "Z-Face-Centers") == 0 ? user->lCentz :
2591 (strcasecmp(field_name, "KCsi") == 0 ? user->lKCsi :
2592 (strcasecmp(field_name, "KEta") == 0 ? user->lKEta : user->lKZet)));
2593 dm = user->fda; dof = 3; strcpy(data_layout, "Face-Centered"); dominant_dir = 'z';
2594 } else if (strcasecmp(field_name, "Coordinates") == 0) {
2595 ierr = DMGetCoordinatesLocal(user->da, &vec_local); CHKERRQ(ierr);
2596 dm = user->fda; dof = 3; strcpy(data_layout, "Node-Centered");
2597 } else if (strcasecmp(field_name, "CornerField")== 0){
2598 vec_local = user->lCellFieldAtCorner; strcpy(data_layout, "Node-Centered");
2599 PetscInt bs = 1;
2600 ierr = VecGetBlockSize(user->CellFieldAtCorner, &bs); CHKERRQ(ierr);
2601 dof = bs;
2602 if(dof == 1) dm = user->da;
2603 else dm = user->fda;
2604 } else {
2605 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Unknown field name for LOG_FIELD_ANATOMY: %s", field_name);
2606 }
2607
2608 // --- 2. Get Grid Info and Array Pointers ---
2609 ierr = DMDAGetLocalInfo(dm, &info); CHKERRQ(ierr);
2610
2611 ierr = PetscBarrier(NULL);
2612 PetscPrintf(PETSC_COMM_WORLD, "\n--- Field Anatomy Log: [%s] | Stage: [%s] | Layout: [%s] ---\n", field_name, stage_name, data_layout);
2613
2614 // Global physical dimensions (number of cells)
2615 PetscInt im_phys = user->IM;
2616 PetscInt jm_phys = user->JM;
2617 PetscInt km_phys = user->KM;
2618
2619 // Slice through the center of the local domain
2620 PetscInt i_mid = (PetscInt)(info.xs + 0.5 * info.xm) - 1;
2621 PetscInt j_mid = (PetscInt)(info.ys + 0.5 * info.ym) - 1;
2622 PetscInt k_mid = (PetscInt)(info.zs + 0.5 * info.zm) - 1;
2623
2624 // --- 3. Print Boundary Information based on Data Layout ---
2625
2626 // ======================================================================
2627 // === CASE 1: Cell-Centered Fields (Ucat, P) - USES SHIFTED INDEX ===
2628 // ======================================================================
2629 if (strcmp(data_layout, "Cell-Centered") == 0) {
2630 const void *l_arr;
2631 ierr = DMDAVecGetArrayRead(dm, vec_local, (void*)&l_arr); CHKERRQ(ierr);
2632
2633
2634 // --- I-Direction Boundaries ---
2635 if (info.xs == 0) { // Rank on -Xi boundary
2636 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, I-DIR]: Idx %2d (Ghost for Cell[k][j][0]) = ", rank, 0);
2637 if(dof==1) PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f)\n", ((const PetscReal***)l_arr)[k_mid][j_mid][0]);
2638 else PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f, %.5f, %.5f)\n", ((const Cmpnts***)l_arr)[k_mid][j_mid][0].x, ((const Cmpnts***)l_arr)[k_mid][j_mid][0].y, ((const Cmpnts***)l_arr)[k_mid][j_mid][0].z);
2639
2640 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, I-DIR]: Idx %2d (Value for Cell[k][j][0]) = ", rank, 1);
2641 if(dof==1) PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f)\n", ((const PetscReal***)l_arr)[k_mid][j_mid][1]);
2642 else PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f, %.5f, %.5f)\n", ((const Cmpnts***)l_arr)[k_mid][j_mid][1].x, ((const Cmpnts***)l_arr)[k_mid][j_mid][1].y, ((const Cmpnts***)l_arr)[k_mid][j_mid][1].z);
2643 }
2644 if (info.xs + info.xm == info.mx) { // Rank on +Xi boundary
2645 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, I-DIR]: Idx %2d (Value for Cell[k][j][%d]) = ", rank, im_phys - 1, im_phys - 2);
2646 if(dof==1) PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f)\n", ((const PetscReal***)l_arr)[k_mid][j_mid][im_phys - 1]);
2647 else PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f, %.5f, %.5f)\n", ((const Cmpnts***)l_arr)[k_mid][j_mid][im_phys - 1].x, ((const Cmpnts***)l_arr)[k_mid][j_mid][im_phys - 1].y, ((const Cmpnts***)l_arr)[k_mid][j_mid][im_phys - 1].z);
2648
2649 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, I-DIR]: Idx %2d (Ghost for Cell[k][j][%d]) = ", rank, im_phys, im_phys - 2);
2650 if(dof==1) PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f)\n", ((const PetscReal***)l_arr)[k_mid][j_mid][im_phys]);
2651 else PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f, %.5f, %.5f)\n", ((const Cmpnts***)l_arr)[k_mid][j_mid][im_phys].x, ((const Cmpnts***)l_arr)[k_mid][j_mid][im_phys].y, ((const Cmpnts***)l_arr)[k_mid][j_mid][im_phys].z);
2652 }
2653
2654 // --- J-Direction Boundaries ---
2655 if (info.ys == 0) { // Rank on -Eta boundary
2656 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, J-DIR]: Jdx %2d (Ghost for Cell[k][0][i]) = ", rank, 0);
2657 if(dof==1) PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f)\n", ((const PetscReal***)l_arr)[k_mid][0][i_mid]);
2658 else PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f, %.5f, %.5f)\n", ((const Cmpnts***)l_arr)[k_mid][0][i_mid].x, ((const Cmpnts***)l_arr)[k_mid][0][i_mid].y, ((const Cmpnts***)l_arr)[k_mid][0][i_mid].z);
2659
2660 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, J-DIR]: Jdx %2d (Value for Cell[k][0][i]) = ", rank, 1);
2661 if(dof==1) PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f)\n", ((const PetscReal***)l_arr)[k_mid][1][i_mid]);
2662 else PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f, %.5f, %.5f)\n", ((const Cmpnts***)l_arr)[k_mid][1][i_mid].x, ((const Cmpnts***)l_arr)[k_mid][1][i_mid].y, ((const Cmpnts***)l_arr)[k_mid][1][i_mid].z);
2663 }
2664
2665 if (info.ys + info.ym == info.my) { // Rank on +Eta boundary
2666 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, J-DIR]: Jdx %2d (Value for Cell[k][%d][i]) = ", rank, jm_phys - 1, jm_phys - 2);
2667 if(dof==1) PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f)\n", ((const PetscReal***)l_arr)[k_mid][jm_phys - 1][i_mid]);
2668 else PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f, %.5f, %.5f)\n", ((const Cmpnts***)l_arr)[k_mid][jm_phys - 1][i_mid].x, ((const Cmpnts***)l_arr)[k_mid][jm_phys - 1][i_mid].y, ((const Cmpnts***)l_arr)[k_mid][jm_phys - 1][i_mid].z);
2669
2670 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, J-DIR]: Jdx %2d (Ghost for Cell[k][%d][i]) = ", rank, jm_phys, jm_phys - 2);
2671 if(dof==1) PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f)\n", ((const PetscReal***)l_arr)[k_mid][jm_phys][i_mid]);
2672 else PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f, %.5f, %.5f)\n", ((const Cmpnts***)l_arr)[k_mid][jm_phys][i_mid].x, ((const Cmpnts***)l_arr)[k_mid][jm_phys][i_mid].y, ((const Cmpnts***)l_arr)[k_mid][jm_phys][i_mid].z);
2673 }
2674
2675 // --- K-Direction Boundaries ---
2676 if (info.zs == 0) { // Rank on -Zeta boundary
2677 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, K-DIR]: Kdx %2d (Ghost for Cell[0][j][i]) = ", rank, 0);
2678 if(dof==1) PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f)\n", ((const PetscReal***)l_arr)[0][j_mid][i_mid]);
2679 else PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f, %.5f, %.5f)\n", ((const Cmpnts***)l_arr)[0][j_mid][i_mid].x, ((const Cmpnts***)l_arr)[0][j_mid][i_mid].y, ((const Cmpnts***)l_arr)[0][j_mid][i_mid].z);
2680 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, K-DIR]: Kdx %2d (Value for Cell[0][j][i]) = ", rank, 1);
2681 if(dof==1) PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f)\n", ((const PetscReal***)l_arr)[1][j_mid][i_mid]);
2682 else PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f, %.5f, %.5f)\n", ((const Cmpnts***)l_arr)[1][j_mid][i_mid].x, ((const Cmpnts***)l_arr)[1][j_mid][i_mid].y, ((const Cmpnts***)l_arr)[1][j_mid][i_mid].z);
2683 }
2684 if (info.zs + info.zm == info.mz) { // Rank on +Zeta boundary
2685 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, K-DIR]: Kdx %2d (Value for Cell[%d][j][i]) = ", rank, km_phys - 1, km_phys - 2);
2686 if(dof==1) PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f)\n", ((const PetscReal***)l_arr)[km_phys - 1][j_mid][i_mid]);
2687 else PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f, %.5f, %.5f)\n", ((const Cmpnts***)l_arr)[km_phys - 1][j_mid][i_mid].x, ((const Cmpnts***)l_arr)[km_phys - 1][j_mid][i_mid].y, ((const Cmpnts***)l_arr)[km_phys - 1][j_mid][i_mid].z);
2688 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, K-DIR]: Kdx %2d (Ghost for Cell[%d][j][i]) = ", rank, km_phys, km_phys - 2);
2689 if(dof==1) PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f)\n", ((const PetscReal***)l_arr)[km_phys][j_mid][i_mid]);
2690 else PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f, %.5f, %.5f)\n", ((const Cmpnts***)l_arr)[km_phys][j_mid][i_mid].x, ((const Cmpnts***)l_arr)[km_phys][j_mid][i_mid].y, ((const Cmpnts***)l_arr)[km_phys][j_mid][i_mid].z);
2691 }
2692 ierr = DMDAVecRestoreArrayRead(dm, vec_local, (void*)&l_arr); CHKERRQ(ierr);
2693 }
2694 // ======================================================================
2695 // === CASE 2: Face-Centered Fields - NUANCED DIRECTIONAL LOGIC ===
2696 // ======================================================================
2697 else if (strcmp(data_layout, "Face-Centered") == 0) {
2698 const Cmpnts ***l_arr;
2699 ierr = DMDAVecGetArrayRead(dm, vec_local, (void*)&l_arr); CHKERRQ(ierr);
2700
2701 // --- I-Direction Boundaries ---
2702 if (info.xs == 0) { // Rank on -Xi boundary
2703 if (dominant_dir == 'x') { // Node-like in I-dir
2704 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, I-DIR]: Idx %2d (First Phys. X-Face) = (%.5f, %.5f, %.5f)\n", rank, 0, l_arr[k_mid][j_mid][0].x, l_arr[k_mid][j_mid][0].y, l_arr[k_mid][j_mid][0].z);
2705 } else if (dominant_dir == 'y' || dominant_dir == 'z') { // Cell-like in I-dir
2706 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, I-DIR]: Idx %2d (Ghost for Cell[k][j][0]) = (%.5f, %.5f, %.5f)\n", rank, 0, l_arr[k_mid][j_mid][0].x, l_arr[k_mid][j_mid][0].y, l_arr[k_mid][j_mid][0].z);
2707 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, I-DIR]: Idx %2d (Value for Cell[k][j][0]) = (%.5f, %.5f, %.5f)\n", rank, 1, l_arr[k_mid][j_mid][1].x, l_arr[k_mid][j_mid][1].y, l_arr[k_mid][j_mid][1].z);
2708 } else if (dominant_dir == 'm') { // Ucont: Mixed
2709 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, I-DIR]: u-comp @ Idx %2d (1st X-Face) = %.5f\n", rank, 0, l_arr[k_mid][j_mid][0].x);
2710 }
2711 }
2712 if (info.xs + info.xm == info.mx) { // Rank on +Xi boundary
2713 if (dominant_dir == 'x') { // Node-like in I-dir
2714 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, I-DIR]: Idx %2d (Last Phys. X-Face) = (%.5f, %.5f, %.5f)\n", rank, im_phys - 1, l_arr[k_mid][j_mid][im_phys - 1].x, l_arr[k_mid][j_mid][im_phys-1].y, l_arr[k_mid][j_mid][im_phys - 1].z);
2715 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, I-DIR]: Idx %2d (Ghost Location) = (%.5f, %.5f, %.5f)\n", rank, im_phys, l_arr[k_mid][j_mid][im_phys].x, l_arr[k_mid][j_mid][im_phys].y, l_arr[k_mid][j_mid][im_phys].z);
2716 } else if (dominant_dir == 'y' || dominant_dir == 'z') { // Cell-like in I-dir
2717 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, I-DIR]: Idx %2d (Value for Cell[k][j][%d]) = (%.5f, %.5f, %.5f)\n", rank, im_phys - 1, im_phys - 2, l_arr[k_mid][j_mid][im_phys - 1].x, l_arr[k_mid][j_mid][im_phys - 1].y, l_arr[k_mid][j_mid][im_phys-1].z);
2718 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, I-DIR]: Idx %2d (Ghost for Cell[k][j][%d]) = (%.5f, %.5f, %.5f)\n", rank, im_phys, im_phys - 2, l_arr[k_mid][j_mid][im_phys].x, l_arr[k_mid][j_mid][im_phys].y, l_arr[k_mid][j_mid][im_phys].z);
2719 } else if (dominant_dir == 'm') { // Ucont: Mixed
2720 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, I-DIR]: u-comp @ Idx %2d (Last X-Face) = %.5f\n", rank, im_phys - 1, l_arr[k_mid][j_mid][im_phys - 1].x);
2721 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, I-DIR]: u-comp @ Idx %2d (Ghost Location) = %.5f\n", rank, im_phys, l_arr[k_mid][j_mid][im_phys].x);
2722 }
2723 }
2724
2725 // --- J-Direction Boundaries ---
2726 if (info.ys == 0) { // Rank on -Eta boundary
2727 if (dominant_dir == 'y') { // Node-like in J-dir
2728 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, J-DIR]: Jdx %2d (First Phys. Y-Face) = (%.5f, %.5f, %.5f)\n", rank, 0, l_arr[k_mid][0][i_mid].x, l_arr[k_mid][0][i_mid].y, l_arr[k_mid][0][i_mid].z);
2729 } else if (dominant_dir == 'x' || dominant_dir == 'z') { // Cell-like in J-dir
2730 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, J-DIR]: Jdx %2d (Ghost for Cell[k][0][i]) = (%.5f, %.5f, %.5f)\n", rank, 0, l_arr[k_mid][0][i_mid].x, l_arr[k_mid][0][i_mid].y, l_arr[k_mid][0][i_mid].z);
2731 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, J-DIR]: Jdx %2d (Value for Cell[k][0][i]) = (%.5f, %.5f, %.5f)\n", rank, 1, l_arr[k_mid][1][i_mid].x, l_arr[k_mid][1][i_mid].y, l_arr[k_mid][1][i_mid].z);
2732 } else if (dominant_dir == 'm') { // Ucont: Mixed
2733 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, J-DIR]: v-comp @ Jdx %2d (1st Y-Face) = %.5f\n", rank, 0, l_arr[k_mid][0][i_mid].y);
2734 }
2735 }
2736 if (info.ys + info.ym == info.my) { // Rank on +Eta boundary
2737 if (dominant_dir == 'y') { // Node-like in J-dir
2738 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, J-DIR]: Jdx %2d (Last Phys. Y-Face) = (%.5f, %.5f, %.5f)\n", rank, jm_phys - 1, l_arr[k_mid][jm_phys - 1][i_mid].x, l_arr[k_mid][jm_phys - 1][i_mid].y, l_arr[k_mid][jm_phys - 1][i_mid].z);
2739 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, J-DIR]: Jdx %2d (Ghost Location) = (%.5f, %.5f, %.5f)\n", rank, jm_phys, l_arr[k_mid][jm_phys][i_mid].x, l_arr[k_mid][jm_phys][i_mid].y, l_arr[k_mid][jm_phys][i_mid].z);
2740 } else if (dominant_dir == 'x' || dominant_dir == 'z') { // Cell-like in J-dir
2741 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, J-DIR]: Jdx %2d (Value for Cell[k][%d][i]) = (%.5f, %.5f, %.5f)\n", rank, jm_phys-1, jm_phys-2, l_arr[k_mid][jm_phys - 1][i_mid].x, l_arr[k_mid][jm_phys - 1][i_mid].y, l_arr[k_mid][jm_phys - 1][i_mid].z);
2742 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, J-DIR]: Jdx %2d (Ghost for Cell[k][%d][i]) = (%.5f, %.5f, %.5f)\n", rank, jm_phys, jm_phys-2, l_arr[k_mid][jm_phys][i_mid].x, l_arr[k_mid][jm_phys][i_mid].y, l_arr[k_mid][jm_phys][i_mid].z);
2743 } else if (dominant_dir == 'm') { // Ucont: Mixed
2744 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, J-DIR]: v-comp @ Jdx %2d (Last Y-Face) = %.5f\n", rank, jm_phys - 1, l_arr[k_mid][jm_phys - 1][i_mid].y);
2745 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, J-DIR]: v-comp @ Jdx %2d (Ghost Location) = %.5f\n", rank, jm_phys, l_arr[k_mid][jm_phys][i_mid].y);
2746 }
2747 }
2748
2749 // --- K-Direction Boundaries ---
2750 if (info.zs == 0) { // Rank on -Zeta boundary
2751 if (dominant_dir == 'z') { // Node-like in K-dir
2752 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, K-DIR]: Kdx %2d (First Phys. Z-Face) = (%.5f, %.5f, %.5f)\n", rank, 0, l_arr[0][j_mid][i_mid].x, l_arr[0][j_mid][i_mid].y, l_arr[0][j_mid][i_mid].z);
2753 } else if (dominant_dir == 'x' || dominant_dir == 'y') { // Cell-like in K-dir
2754 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, K-DIR]: Kdx %2d (Ghost for Cell[0][j][i]) = (%.5f, %.5f, %.5f)\n", rank, 0, l_arr[0][j_mid][i_mid].x, l_arr[0][j_mid][i_mid].y, l_arr[0][j_mid][i_mid].z);
2755 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, K-DIR]: Kdx %2d (Value for Cell[0][j][i]) = (%.5f, %.5f, %.5f)\n", rank, 1, l_arr[1][j_mid][i_mid].x, l_arr[1][j_mid][i_mid].y, l_arr[1][j_mid][i_mid].z);
2756 } else if (dominant_dir == 'm') { // Ucont: Mixed
2757 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, K-DIR]: w-comp @ Idx %2d (1st Z-Face) = %.5f\n", rank, 0, l_arr[0][j_mid][i_mid].z);
2758 }
2759 }
2760 if (info.zs + info.zm == info.mz) { // Rank on +Zeta boundary
2761 if (dominant_dir == 'z') { // Node-like in K-dir
2762 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, K-DIR]: Idx %2d (Last Phys. Z-Face) = (%.5f, %.5f, %.5f)\n", rank, km_phys - 1, l_arr[km_phys - 1][j_mid][i_mid].x, l_arr[km_phys - 1][j_mid][i_mid].y, l_arr[km_phys - 1][j_mid][i_mid].z);
2763 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, K-DIR]: Idx %2d (Ghost Location) = (%.5f, %.5f, %.5f)\n", rank, km_phys, l_arr[km_phys][j_mid][i_mid].x, l_arr[km_phys][j_mid][i_mid].y, l_arr[km_phys][j_mid][i_mid].z);
2764 } else if (dominant_dir == 'x' || dominant_dir == 'y') { // Cell-like in K-dir
2765 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, K-DIR]: Idx %2d (Value for Cell[%d][j][i]) = (%.5f, %.5f, %.5f)\n", rank, km_phys-1, km_phys-2, l_arr[km_phys-1][j_mid][i_mid].x, l_arr[km_phys-1][j_mid][i_mid].y, l_arr[km_phys - 1][j_mid][i_mid].z);
2766 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, K-DIR]: Idx %2d (Ghost for Cell[%d][j][i]) = (%.5f, %.5f, %.5f)\n", rank, km_phys, km_phys-2, l_arr[km_phys][j_mid][i_mid].x, l_arr[km_phys][j_mid][i_mid].y, l_arr[km_phys][j_mid][i_mid].z);
2767 } else if (dominant_dir == 'm') { // Ucont: Mixed
2768 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, K-DIR]: w-comp @ Idx %2d (Last Z-Face) = %.5f\n", rank, km_phys - 1, l_arr[km_phys - 1][j_mid][i_mid].z);
2769 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, K-DIR]: w-comp @ Idx %2d (Ghost Loc.) = %.5f\n", rank, km_phys, l_arr[km_phys][j_mid][i_mid].z);
2770
2771 }
2772 }
2773 ierr = DMDAVecRestoreArrayRead(dm, vec_local, (void*)&l_arr); CHKERRQ(ierr);
2774 }
2775 // ======================================================================
2776 // === CASE 3: Node-Centered Fields - USES DIRECT INDEX ===
2777 // ======================================================================
2778 else if (strcmp(data_layout, "Node-Centered") == 0) {
2779 const Cmpnts ***l_arr;
2780 ierr = DMDAVecGetArrayRead(dm, vec_local, (void*)&l_arr); CHKERRQ(ierr);
2781
2782 // --- I-Direction Boundaries ---
2783 if (info.xs == 0) {
2784 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, I-DIR]: Idx %2d (First Phys. Node) = (%.5f, %.5f, %.5f)\n", rank, 0, l_arr[k_mid][j_mid][0].x, l_arr[k_mid][j_mid][0].y, l_arr[k_mid][j_mid][0].z);
2785 }
2786 if (info.xs + info.xm == info.mx) {
2787 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, I-DIR]: Idx %2d (Last Phys. Node) = (%.5f, %.5f, %.5f)\n", rank, im_phys - 1, l_arr[k_mid][j_mid][im_phys - 1].x, l_arr[k_mid][j_mid][im_phys - 1].y, l_arr[k_mid][j_mid][im_phys - 1].z);
2788 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, I-DIR]: Idx %2d (Unused/Ghost Loc) = (%.5f, %.5f, %.5f)\n", rank, im_phys, l_arr[k_mid][j_mid][im_phys].x, l_arr[k_mid][j_mid][im_phys].y, l_arr[k_mid][j_mid][im_phys].z);
2789 }
2790 // --- J-Direction Boundaries ---
2791 if (info.ys == 0) {
2792 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, J-DIR]: Jdx %2d (First Phys. Node) = (%.5f, %.5f, %.5f)\n", rank, 0, l_arr[k_mid][0][i_mid].x, l_arr[k_mid][0][i_mid].y, l_arr[k_mid][0][i_mid].z);
2793 }
2794 if (info.ys + info.ym == info.my) {
2795 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, J-DIR]: Jdx %2d (Last Phys. Node) = (%.5f, %.5f, %.5f)\n", rank, jm_phys - 1, l_arr[k_mid][jm_phys - 1][i_mid].x, l_arr[k_mid][jm_phys - 1][i_mid].y, l_arr[k_mid][jm_phys - 1][i_mid].z);
2796 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, J-DIR]: Jdx %2d (Unused/Ghost Loc) = (%.5f, %.5f, %.5f)\n", rank, jm_phys, l_arr[k_mid][jm_phys][i_mid].x, l_arr[k_mid][jm_phys][i_mid].y, l_arr[k_mid][jm_phys][i_mid].z);
2797 }
2798 // --- K-Direction Boundaries ---
2799 if (info.zs == 0) {
2800 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, K-DIR]: Kdx %2d (First Phys. Node) = (%.5f, %.5f, %.5f)\n", rank, 0, l_arr[0][j_mid][i_mid].x, l_arr[0][j_mid][i_mid].y, l_arr[0][j_mid][i_mid].z);
2801 }
2802 if(info.zs + info.zm == info.mz) {
2803 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, K-DIR]: Kdx %2d (Last Phys. Node) = (%.5f, %.5f, %.5f)\n", rank, km_phys - 1, l_arr[km_phys - 1][j_mid][i_mid].x, l_arr[km_phys - 1][j_mid][i_mid].y, l_arr[km_phys - 1][j_mid][i_mid].z);
2804 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, K-DIR]: Kdx %2d (Unused/Ghost Loc) = (%.5f, %.5f, %.5f)\n", rank, km_phys, l_arr[km_phys][j_mid][i_mid].x, l_arr[km_phys][j_mid][i_mid].y, l_arr[km_phys][j_mid][i_mid].z);
2805 }
2806 ierr = DMDAVecRestoreArrayRead(dm, vec_local, (void*)&l_arr); CHKERRQ(ierr);
2807 }
2808 else {
2809 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "LOG_FIELD_ANATOMY encountered an unknown data layout: %s", data_layout);
2810 }
2811
2812 ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT); CHKERRQ(ierr);
2813 ierr = PetscBarrier(NULL);
2814 PetscFunctionReturn(0);
2815}
2816
2817#undef __FUNCT__
2818#define __FUNCT__ "LOG_INTERPOLATION_ERROR"
2819/**
2820 * @brief Implementation of \ref LOG_INTERPOLATION_ERROR().
2821 * @details Full API contract (arguments, ownership, side effects) is documented with
2822 * the header declaration in `include/logging.h`.
2823 * @see LOG_INTERPOLATION_ERROR()
2824 */
2825PetscErrorCode LOG_INTERPOLATION_ERROR(UserCtx *user)
2826{
2827 SimCtx *simCtx = user->simCtx;
2828 PetscErrorCode ierr;
2829 DM swarm = user->swarm;
2830 Vec positionVec, analyticalvelocityVec, velocityVec, errorVec;
2831 PetscReal Interpolation_error = 0.0;
2832 PetscReal Maximum_Interpolation_error = 0.0;
2833 PetscReal AnalyticalSolution_magnitude = 0.0;
2834 PetscReal ErrorPercentage = 0.0;
2835
2836 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Creating global vectors.\n");
2837 ierr = DMSwarmCreateGlobalVectorFromField(swarm, "position", &positionVec); CHKERRQ(ierr);
2838 ierr = DMSwarmCreateGlobalVectorFromField(swarm, "velocity", &velocityVec); CHKERRQ(ierr);
2839
2840 ierr = VecDuplicate(positionVec, &analyticalvelocityVec); CHKERRQ(ierr);
2841 ierr = VecCopy(positionVec, analyticalvelocityVec); CHKERRQ(ierr);
2842
2843 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Computing analytical solution.\n");
2844 ierr = SetAnalyticalSolutionForParticles(analyticalvelocityVec, simCtx); CHKERRQ(ierr);
2845
2846 ierr = VecDuplicate(analyticalvelocityVec, &errorVec); CHKERRQ(ierr);
2847 ierr = VecCopy(analyticalvelocityVec, errorVec); CHKERRQ(ierr);
2848
2849 ierr = VecNorm(analyticalvelocityVec, NORM_2, &AnalyticalSolution_magnitude); CHKERRQ(ierr);
2850
2851 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Computing error.\n");
2852 ierr = VecAXPY(errorVec, -1.0, velocityVec); CHKERRQ(ierr);
2853 ierr = VecNorm(errorVec, NORM_2, &Interpolation_error); CHKERRQ(ierr);
2854 ierr = VecNorm(errorVec,NORM_INFINITY,&Maximum_Interpolation_error); CHKERRQ(ierr);
2855
2856 ErrorPercentage = (AnalyticalSolution_magnitude > 0) ?
2857 (Interpolation_error / AnalyticalSolution_magnitude * 100.0) : 0.0;
2858
2859 /* --- CSV output (always, rank 0 only) --- */
2860 if (simCtx->rank == 0) {
2861 char csv_path[PETSC_MAX_PATH_LEN + 32];
2862 ierr = PetscSNPrintf(csv_path, sizeof(csv_path), "%s/interpolation_error.csv", simCtx->log_dir); CHKERRQ(ierr);
2863 FILE *f = fopen(csv_path, "a");
2864 if (f) {
2865 if (ftell(f) == 0) {
2866 fprintf(f, "step,time,L2_error,Linf_error,L2_analytical,error_pct\n");
2867 }
2868 if (simCtx->continueMode && simCtx->step == simCtx->StartStep + 1) {
2869 fprintf(f, "# Continuation from step %" PetscInt_FMT "\n", simCtx->StartStep);
2870 }
2871 PetscReal t = (PetscReal)simCtx->ti * simCtx->dt;
2872 fprintf(f, "%d,%.6e,%.6e,%.6e,%.6e,%.4f\n",
2873 (int)simCtx->step, t,
2874 Interpolation_error, Maximum_Interpolation_error,
2875 AnalyticalSolution_magnitude, ErrorPercentage);
2876 fclose(f);
2877 }
2878 }
2879
2880 /* --- Console output (only at INFO level or above) --- */
2881 if (get_log_level() >= LOG_INFO) {
2882 LOG_ALLOW(GLOBAL, LOG_INFO, "Interpolation error (%%): %g\n", ErrorPercentage);
2883 PetscPrintf(PETSC_COMM_WORLD, "Interpolation error (%%): %g\n", ErrorPercentage);
2884 LOG_ALLOW(GLOBAL, LOG_INFO, "Maximum Interpolation error: %g\n", Maximum_Interpolation_error);
2885 PetscPrintf(PETSC_COMM_WORLD, "Maximum Interpolation error: %g\n", Maximum_Interpolation_error);
2886 }
2887
2888 ierr = VecDestroy(&analyticalvelocityVec); CHKERRQ(ierr);
2889 ierr = VecDestroy(&errorVec); CHKERRQ(ierr);
2890 ierr = DMSwarmDestroyGlobalVectorFromField(swarm, "position", &positionVec); CHKERRQ(ierr);
2891 ierr = DMSwarmDestroyGlobalVectorFromField(swarm, "velocity", &velocityVec); CHKERRQ(ierr);
2892
2893 return 0;
2894}
2895
2896#undef __FUNCT__
2897#define __FUNCT__ "LOG_SCATTER_METRICS"
2898/**
2899 * @brief Implementation of \ref LOG_SCATTER_METRICS().
2900 * @details Full API contract (arguments, ownership, side effects) is documented with
2901 * the header declaration in `include/logging.h`.
2902 * @see LOG_SCATTER_METRICS()
2903 */
2904PetscErrorCode LOG_SCATTER_METRICS(UserCtx *user)
2905{
2906 PetscErrorCode ierr;
2907 SimCtx *simCtx = NULL;
2908 DMDALocalInfo info;
2909 PetscInt xs, xe, ys, ye, zs, ze, mx, my, mz;
2910 PetscInt lxs, lxe, lys, lye, lzs, lze;
2911 Vec reference_vec = NULL;
2912 PetscReal ***psi = NULL;
2913 PetscReal ***psi_ref = NULL;
2914 PetscReal ***aj = NULL;
2915 PetscReal ***count = NULL;
2916 PetscReal *particle_psi = NULL;
2917 PetscInt nlocal = 0;
2918 PetscReal local_l1 = 0.0, global_l1 = 0.0;
2919 PetscReal local_l2_sq = 0.0, global_l2_sq = 0.0;
2920 PetscReal local_linf = 0.0, global_linf = 0.0;
2921 PetscReal local_ref_l2_sq = 0.0, global_ref_l2_sq = 0.0;
2922 PetscReal local_grid_integral = 0.0, global_grid_integral = 0.0;
2923 PetscReal local_domain_volume = 0.0, global_domain_volume = 0.0;
2924 PetscReal local_particle_sum = 0.0, global_particle_sum = 0.0;
2925 PetscInt64 local_particle_count = 0, global_particle_count = 0;
2926 PetscInt64 local_cell_count = 0, global_cell_count = 0;
2927 PetscInt64 local_occupied_count = 0, global_occupied_count = 0;
2928 PetscReal particle_integral = 0.0;
2929 PetscReal occupancy_fraction = 0.0;
2930 PetscReal mean_particles_per_occupied_cell = 0.0;
2931 PetscReal l2_error = 0.0;
2932 PetscReal relative_l2_error = 0.0;
2933
2934 PetscFunctionBeginUser;
2935 if (!user) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "UserCtx cannot be NULL.");
2936 simCtx = user->simCtx;
2937 if (!VerificationScalarOverrideActive(simCtx) || !user->swarm || !user->Psi || !user->ParticleCount) {
2938 PetscFunctionReturn(0);
2939 }
2940
2941 info = user->info;
2942 xs = info.xs; xe = info.xs + info.xm;
2943 ys = info.ys; ye = info.ys + info.ym;
2944 zs = info.zs; ze = info.zs + info.zm;
2945 mx = info.mx; my = info.my; mz = info.mz;
2946 lxs = (xs == 0) ? xs + 1 : xs; lxe = (xe == mx) ? xe - 1 : xe;
2947 lys = (ys == 0) ? ys + 1 : ys; lye = (ye == my) ? ye - 1 : ye;
2948 lzs = (zs == 0) ? zs + 1 : zs; lze = (ze == mz) ? ze - 1 : ze;
2949
2950 ierr = VecDuplicate(user->Psi, &reference_vec); CHKERRQ(ierr);
2951 ierr = SetAnalyticalScalarFieldAtCellCenters(user, reference_vec); CHKERRQ(ierr);
2952
2953 ierr = DMDAVecGetArrayRead(user->da, user->Psi, &psi); CHKERRQ(ierr);
2954 ierr = DMDAVecGetArrayRead(user->da, reference_vec, &psi_ref); CHKERRQ(ierr);
2955 ierr = DMDAVecGetArrayRead(user->da, user->Aj, &aj); CHKERRQ(ierr);
2956 ierr = DMDAVecGetArrayRead(user->da, user->ParticleCount, &count); CHKERRQ(ierr);
2957
2958 for (PetscInt k = lzs; k < lze; ++k) {
2959 for (PetscInt j = lys; j < lye; ++j) {
2960 for (PetscInt i = lxs; i < lxe; ++i) {
2961 const PetscReal cell_volume = (PetscAbsReal(aj[k][j][i]) > 1.0e-14) ? (1.0 / aj[k][j][i]) : 0.0;
2962 const PetscReal err = psi[k][j][i] - psi_ref[k][j][i];
2963 local_cell_count += 1;
2964 local_domain_volume += cell_volume;
2965 local_grid_integral += psi[k][j][i] * cell_volume;
2966 local_l1 += PetscAbsReal(err) * cell_volume;
2967 local_l2_sq += err * err * cell_volume;
2968 local_ref_l2_sq += psi_ref[k][j][i] * psi_ref[k][j][i] * cell_volume;
2969 local_linf = PetscMax(local_linf, PetscAbsReal(err));
2970 if (count[k][j][i] > 0.0) local_occupied_count += 1;
2971 }
2972 }
2973 }
2974
2975 ierr = DMDAVecRestoreArrayRead(user->da, user->ParticleCount, &count); CHKERRQ(ierr);
2976 ierr = DMDAVecRestoreArrayRead(user->da, user->Aj, &aj); CHKERRQ(ierr);
2977 ierr = DMDAVecRestoreArrayRead(user->da, reference_vec, &psi_ref); CHKERRQ(ierr);
2978 ierr = DMDAVecRestoreArrayRead(user->da, user->Psi, &psi); CHKERRQ(ierr);
2979 ierr = VecDestroy(&reference_vec); CHKERRQ(ierr);
2980
2981 ierr = DMSwarmGetLocalSize(user->swarm, &nlocal); CHKERRQ(ierr);
2982 local_particle_count = (PetscInt64)nlocal;
2983 if (nlocal > 0) {
2984 ierr = DMSwarmGetField(user->swarm, "Psi", NULL, NULL, (void **)&particle_psi); CHKERRQ(ierr);
2985 for (PetscInt p = 0; p < nlocal; ++p) local_particle_sum += particle_psi[p];
2986 ierr = DMSwarmRestoreField(user->swarm, "Psi", NULL, NULL, (void **)&particle_psi); CHKERRQ(ierr);
2987 }
2988
2989 ierr = MPI_Allreduce(&local_l1, &global_l1, 1, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD); CHKERRMPI(ierr);
2990 ierr = MPI_Allreduce(&local_l2_sq, &global_l2_sq, 1, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD); CHKERRMPI(ierr);
2991 ierr = MPI_Allreduce(&local_linf, &global_linf, 1, MPIU_REAL, MPI_MAX, PETSC_COMM_WORLD); CHKERRMPI(ierr);
2992 ierr = MPI_Allreduce(&local_ref_l2_sq, &global_ref_l2_sq, 1, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD); CHKERRMPI(ierr);
2993 ierr = MPI_Allreduce(&local_grid_integral, &global_grid_integral, 1, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD); CHKERRMPI(ierr);
2994 ierr = MPI_Allreduce(&local_domain_volume, &global_domain_volume, 1, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD); CHKERRMPI(ierr);
2995 ierr = MPI_Allreduce(&local_particle_sum, &global_particle_sum, 1, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD); CHKERRMPI(ierr);
2996 ierr = MPI_Allreduce(&local_particle_count, &global_particle_count, 1, MPIU_INT64, MPI_SUM, PETSC_COMM_WORLD); CHKERRMPI(ierr);
2997 ierr = MPI_Allreduce(&local_cell_count, &global_cell_count, 1, MPIU_INT64, MPI_SUM, PETSC_COMM_WORLD); CHKERRMPI(ierr);
2998 ierr = MPI_Allreduce(&local_occupied_count, &global_occupied_count, 1, MPIU_INT64, MPI_SUM, PETSC_COMM_WORLD); CHKERRMPI(ierr);
2999
3000 l2_error = PetscSqrtReal(global_l2_sq);
3001 relative_l2_error = (global_ref_l2_sq > 0.0) ? (l2_error / PetscSqrtReal(global_ref_l2_sq)) : 0.0;
3002 occupancy_fraction = (global_cell_count > 0) ? ((PetscReal)global_occupied_count / (PetscReal)global_cell_count) : 0.0;
3003 mean_particles_per_occupied_cell =
3004 (global_occupied_count > 0) ? ((PetscReal)global_particle_count / (PetscReal)global_occupied_count) : 0.0;
3005 particle_integral =
3006 (global_particle_count > 0) ? (global_domain_volume * global_particle_sum / (PetscReal)global_particle_count) : 0.0;
3007
3008 if (simCtx->rank == 0) {
3009 char csv_path[PETSC_MAX_PATH_LEN + 32];
3010 FILE *f = NULL;
3011 ierr = PetscSNPrintf(csv_path, sizeof(csv_path), "%s/scatter_metrics.csv", simCtx->log_dir); CHKERRQ(ierr);
3012 f = fopen(csv_path, "a");
3013 if (f) {
3014 if (ftell(f) == 0) {
3015 fprintf(f,
3016 "step,time,total_particles,total_cells,occupied_cells,occupancy_fraction,"
3017 "mean_particles_per_occupied_cell,particle_integral,grid_integral,"
3018 "conservation_error_abs,L1_error,L2_error,Linf_error,relative_L2_error\n");
3019 }
3020 if (simCtx->continueMode && simCtx->step == simCtx->StartStep + 1) {
3021 fprintf(f, "# Continuation from step %" PetscInt_FMT "\n", simCtx->StartStep);
3022 }
3023 fprintf(f, "%d,%.6e,%lld,%lld,%lld,%.6e,%.6e,%.6e,%.6e,%.6e,%.6e,%.6e,%.6e,%.6e\n",
3024 (int)simCtx->step,
3025 (double)simCtx->ti,
3026 (long long)global_particle_count,
3027 (long long)global_cell_count,
3028 (long long)global_occupied_count,
3029 (double)occupancy_fraction,
3030 (double)mean_particles_per_occupied_cell,
3031 (double)particle_integral,
3032 (double)global_grid_integral,
3033 (double)PetscAbsReal(global_grid_integral - particle_integral),
3034 (double)global_l1,
3035 (double)l2_error,
3036 (double)global_linf,
3037 (double)relative_l2_error);
3038 fclose(f);
3039 }
3040 }
3041
3042 if (get_log_level() >= LOG_INFO) {
3043 LOG_ALLOW(GLOBAL, LOG_INFO, "Scatter relative L2 error: %.6e\n", (double)relative_l2_error);
3044 LOG_ALLOW(GLOBAL, LOG_INFO, "Scatter occupancy fraction: %.6e\n", (double)occupancy_fraction);
3045 }
3046
3047 PetscFunctionReturn(0);
3048}
3049
3050#undef __FUNCT__
3051#define __FUNCT__ "ResetSearchMetrics"
3052/**
3053 * @brief Implementation of \ref ResetSearchMetrics().
3054 * @details Full API contract (arguments, ownership, side effects) is documented with
3055 * the header declaration in `include/logging.h`.
3056 * @see ResetSearchMetrics()
3057 */
3058PetscErrorCode ResetSearchMetrics(SimCtx *simCtx)
3059{
3060 PetscFunctionBeginUser;
3061 if (!simCtx) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "SimCtx cannot be NULL for ResetSearchMetrics.");
3062
3063 simCtx->searchMetrics.searchAttempts = 0;
3064 simCtx->searchMetrics.searchPopulation = 0;
3065 simCtx->searchMetrics.searchLocatedCount = 0;
3066 simCtx->searchMetrics.searchLostCount = 0;
3067 simCtx->searchMetrics.traversalStepsSum = 0;
3068 simCtx->searchMetrics.reSearchCount = 0;
3069 simCtx->searchMetrics.maxTraversalSteps = 0;
3070 simCtx->searchMetrics.maxTraversalFailCount = 0;
3071 simCtx->searchMetrics.tieBreakCount = 0;
3072 simCtx->searchMetrics.boundaryClampCount = 0;
3073 simCtx->searchMetrics.bboxGuessSuccessCount = 0;
3074 simCtx->searchMetrics.bboxGuessFallbackCount = 0;
3075 simCtx->searchMetrics.maxParticlePassDepth = 0;
3076 simCtx->searchMetrics.currentSettlementPass = 0;
3077
3078 PetscFunctionReturn(0);
3079}
3080
3081#undef __FUNCT__
3082#define __FUNCT__ "LOG_SEARCH_METRICS"
3083/**
3084 * @brief Implementation of \ref LOG_SEARCH_METRICS().
3085 * @details Full API contract (arguments, ownership, side effects) is documented with
3086 * the header declaration in `include/logging.h`.
3087 * @see LOG_SEARCH_METRICS()
3088 */
3089PetscErrorCode LOG_SEARCH_METRICS(UserCtx *user)
3090{
3091 PetscErrorCode ierr;
3092 SimCtx *simCtx = NULL;
3093 PetscInt totalParticles = 0;
3094 PetscReal local_metrics[SEARCH_METRIC_REDUCTION_LEN] = {0.0};
3095 PetscReal global_metrics[SEARCH_METRIC_REDUCTION_LEN] = {0.0};
3096 PetscReal meanTraversalSteps = 0.0;
3097 PetscReal searchFailureFraction = 0.0;
3098 PetscReal searchWorkIndex = 0.0;
3099 PetscReal reSearchFraction = 0.0;
3100 long long searchAttempts = 0;
3101 long long searchPopulation = 0;
3102 long long searchLocatedCount = 0;
3103 long long searchLostCount = 0;
3104 long long traversalStepsSum = 0;
3105 long long reSearchCount = 0;
3106 long long tieBreakCount = 0;
3107 long long boundaryClampCount = 0;
3108 long long bboxGuessSuccessCount = 0;
3109 long long bboxGuessFallbackCount = 0;
3110 long long maxTraversalFailCount = 0;
3111 long long maxTraversalSteps = 0;
3112 long long maxPassDepth = 0;
3113 MPI_Op reduction_op = MPI_OP_NULL;
3114
3115 PetscFunctionBeginUser;
3116 if (!user || !user->simCtx) {
3117 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "UserCtx and SimCtx are required for LOG_SEARCH_METRICS.");
3118 }
3119 simCtx = user->simCtx;
3120
3121 if (simCtx->np <= 0) {
3122 PetscFunctionReturn(0);
3123 }
3124
3125 ierr = DMSwarmGetSize(user->swarm, &totalParticles); CHKERRQ(ierr);
3126
3127 local_metrics[SEARCH_METRIC_SUM_SEARCH_ATTEMPTS] = (PetscReal)simCtx->searchMetrics.searchAttempts;
3128 local_metrics[SEARCH_METRIC_SUM_SEARCH_POPULATION] = (PetscReal)simCtx->searchMetrics.searchPopulation;
3129 local_metrics[SEARCH_METRIC_SUM_SEARCH_LOCATED] = (PetscReal)simCtx->searchMetrics.searchLocatedCount;
3130 local_metrics[SEARCH_METRIC_SUM_SEARCH_LOST] = (PetscReal)simCtx->searchMetrics.searchLostCount;
3131 local_metrics[SEARCH_METRIC_SUM_TRAVERSAL_STEPS] = (PetscReal)simCtx->searchMetrics.traversalStepsSum;
3132 local_metrics[SEARCH_METRIC_SUM_RESEARCH] = (PetscReal)simCtx->searchMetrics.reSearchCount;
3133 local_metrics[SEARCH_METRIC_SUM_TIE_BREAKS] = (PetscReal)simCtx->searchMetrics.tieBreakCount;
3134 local_metrics[SEARCH_METRIC_SUM_BOUNDARY_CLAMPS] = (PetscReal)simCtx->searchMetrics.boundaryClampCount;
3135 local_metrics[SEARCH_METRIC_SUM_BBOX_GUESS_SUCCESS] = (PetscReal)simCtx->searchMetrics.bboxGuessSuccessCount;
3136 local_metrics[SEARCH_METRIC_SUM_BBOX_GUESS_FALLBACK] = (PetscReal)simCtx->searchMetrics.bboxGuessFallbackCount;
3137 local_metrics[SEARCH_METRIC_SUM_MAX_TRAVERSAL_FAILS] = (PetscReal)simCtx->searchMetrics.maxTraversalFailCount;
3138 local_metrics[SEARCH_METRIC_MAX_TRAVERSAL_STEPS] = (PetscReal)simCtx->searchMetrics.maxTraversalSteps;
3139 local_metrics[SEARCH_METRIC_MAX_PASS_DEPTH] = (PetscReal)simCtx->searchMetrics.maxParticlePassDepth;
3140
3141 ierr = MPI_Op_create(SearchMetricsReduceOp, PETSC_TRUE, &reduction_op); CHKERRMPI(ierr);
3142 ierr = MPI_Allreduce(local_metrics, global_metrics, SEARCH_METRIC_REDUCTION_LEN, MPIU_REAL, reduction_op, PETSC_COMM_WORLD); CHKERRMPI(ierr);
3143 ierr = MPI_Op_free(&reduction_op); CHKERRMPI(ierr);
3144 reduction_op = MPI_OP_NULL;
3145
3146 searchAttempts = (long long)PetscFloorReal(global_metrics[SEARCH_METRIC_SUM_SEARCH_ATTEMPTS] + 0.5);
3147 searchPopulation = (long long)PetscFloorReal(global_metrics[SEARCH_METRIC_SUM_SEARCH_POPULATION] + 0.5);
3148 searchLocatedCount = (long long)PetscFloorReal(global_metrics[SEARCH_METRIC_SUM_SEARCH_LOCATED] + 0.5);
3149 searchLostCount = (long long)PetscFloorReal(global_metrics[SEARCH_METRIC_SUM_SEARCH_LOST] + 0.5);
3150 traversalStepsSum = (long long)PetscFloorReal(global_metrics[SEARCH_METRIC_SUM_TRAVERSAL_STEPS] + 0.5);
3151 reSearchCount = (long long)PetscFloorReal(global_metrics[SEARCH_METRIC_SUM_RESEARCH] + 0.5);
3152 tieBreakCount = (long long)PetscFloorReal(global_metrics[SEARCH_METRIC_SUM_TIE_BREAKS] + 0.5);
3153 boundaryClampCount = (long long)PetscFloorReal(global_metrics[SEARCH_METRIC_SUM_BOUNDARY_CLAMPS] + 0.5);
3154 bboxGuessSuccessCount = (long long)PetscFloorReal(global_metrics[SEARCH_METRIC_SUM_BBOX_GUESS_SUCCESS] + 0.5);
3155 bboxGuessFallbackCount = (long long)PetscFloorReal(global_metrics[SEARCH_METRIC_SUM_BBOX_GUESS_FALLBACK] + 0.5);
3156 maxTraversalFailCount = (long long)PetscFloorReal(global_metrics[SEARCH_METRIC_SUM_MAX_TRAVERSAL_FAILS] + 0.5);
3157 maxTraversalSteps = (long long)PetscFloorReal(global_metrics[SEARCH_METRIC_MAX_TRAVERSAL_STEPS] + 0.5);
3158 maxPassDepth = (long long)PetscFloorReal(global_metrics[SEARCH_METRIC_MAX_PASS_DEPTH] + 0.5);
3159
3160 if (searchAttempts > 0) {
3161 meanTraversalSteps = (PetscReal)traversalStepsSum / (PetscReal)searchAttempts;
3162 }
3163 if (searchPopulation > 0) {
3164 searchFailureFraction = (PetscReal)searchLostCount / (PetscReal)searchPopulation;
3165 searchWorkIndex = (PetscReal)traversalStepsSum / (PetscReal)searchPopulation;
3166 reSearchFraction = (PetscReal)reSearchCount / (PetscReal)searchPopulation;
3167 }
3168
3169 if (simCtx->rank == 0) {
3170 char csv_path[PETSC_MAX_PATH_LEN + 32];
3171 FILE *f = NULL;
3172
3173 ierr = PetscSNPrintf(csv_path, sizeof(csv_path), "%s/search_metrics.csv", simCtx->log_dir); CHKERRQ(ierr);
3174 f = fopen(csv_path, "a");
3175 if (!f) {
3176 LOG_ALLOW(GLOBAL, LOG_WARNING, "LOG_SEARCH_METRICS: could not open '%s' for writing.\n", csv_path);
3177 } else {
3178 if (ftell(f) == 0) {
3179 fprintf(f,
3180 "step,time,total_particles,lost,lost_cumulative,migrated,migration_passes,search_attempts,"
3181 "mean_traversal_steps,max_traversal_steps,tie_break_count,boundary_clamp_count,"
3182 "bbox_guess_success_count,bbox_guess_fallback_count,max_particle_pass_depth,load_imbalance,"
3183 "search_population,search_located_count,search_lost_count,traversal_steps_sum,re_search_count,"
3184 "max_traversal_fail_count,search_failure_fraction,search_work_index,re_search_fraction\n");
3185 }
3186 if (simCtx->continueMode && simCtx->step == simCtx->StartStep + 1) {
3187 fprintf(f, "# Continuation from step %" PetscInt_FMT "\n", simCtx->StartStep);
3188 }
3189 fprintf(f,
3190 "%d,%.6e,%d,%d,%d,%d,%d,%lld,%.6e,%lld,%lld,%lld,%lld,%lld,%lld,%.6e,%lld,%lld,%lld,%lld,%lld,%lld,%.6e,%.6e,%.6e\n",
3191 (int)simCtx->step,
3192 (double)simCtx->ti,
3193 (int)totalParticles,
3194 (int)simCtx->particlesLostLastStep,
3195 (int)simCtx->particlesLostCumulative,
3196 (int)simCtx->particlesMigratedLastStep,
3197 (int)simCtx->migrationPassesLastStep,
3198 searchAttempts,
3199 (double)meanTraversalSteps,
3200 maxTraversalSteps,
3201 tieBreakCount,
3202 boundaryClampCount,
3203 bboxGuessSuccessCount,
3204 bboxGuessFallbackCount,
3205 maxPassDepth,
3206 (double)simCtx->particleLoadImbalance,
3207 searchPopulation,
3208 searchLocatedCount,
3209 searchLostCount,
3210 traversalStepsSum,
3211 reSearchCount,
3212 maxTraversalFailCount,
3213 (double)searchFailureFraction,
3214 (double)searchWorkIndex,
3215 (double)reSearchFraction);
3216 fclose(f);
3217 }
3218 }
3219
3220 LOG_ALLOW(GLOBAL, LOG_DEBUG,
3221 "Search metrics: sff=%.3e swi=%.3e re_search=%.3e lost(step/total)=%d/%d migrated=%d passes=%d traversal(mean/max)=%.2f/%lld tie_breaks=%lld max_pass_depth=%lld\n",
3222 (double)searchFailureFraction,
3223 (double)searchWorkIndex,
3224 (double)reSearchFraction,
3225 (int)simCtx->particlesLostLastStep,
3226 (int)simCtx->particlesLostCumulative,
3227 (int)simCtx->particlesMigratedLastStep,
3228 (int)simCtx->migrationPassesLastStep,
3229 (double)meanTraversalSteps,
3230 maxTraversalSteps,
3231 tieBreakCount,
3232 maxPassDepth);
3233
3234 PetscFunctionReturn(0);
3235}
3236
3237#undef __FUNCT__
3238#define __FUNCT__ "CalculateAdvancedParticleMetrics"
3239/**
3240 * @brief Internal helper implementation: `CalculateAdvancedParticleMetrics()`.
3241 * @details Local to this translation unit.
3242 */
3243PetscErrorCode CalculateAdvancedParticleMetrics(UserCtx *user)
3244{
3245 PetscErrorCode ierr;
3246 SimCtx *simCtx = user->simCtx;
3247 PetscMPIInt size, rank;
3248
3249 PetscFunctionBeginUser;
3250 ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size); CHKERRQ(ierr);
3251 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
3252
3253 // --- 1. Particle Load Imbalance ---
3254 PetscInt nLocal, nGlobal, nLocalMax;
3255 ierr = DMSwarmGetLocalSize(user->swarm, &nLocal); CHKERRQ(ierr);
3256 ierr = DMSwarmGetSize(user->swarm, &nGlobal); CHKERRQ(ierr);
3257 ierr = MPI_Allreduce(&nLocal, &nLocalMax, 1, MPIU_INT, MPI_MAX, PETSC_COMM_WORLD); CHKERRQ(ierr);
3258
3259 PetscReal avg_per_rank = (size > 0) ? ((PetscReal)nGlobal / size) : 0.0;
3260 // Handle division by zero if there are no particles
3261 simCtx->particleLoadImbalance = (avg_per_rank > 1e-9) ? (nLocalMax / avg_per_rank) : 1.0;
3262
3263
3264 // --- 2. Number of Occupied Cells ---
3265 // This part requires access to the user->ParticleCount vector.
3266 PetscInt local_occupied_cells = 0;
3267 PetscInt global_occupied_cells;
3268 const PetscScalar *count_array;
3269 PetscInt vec_local_size;
3270
3271 ierr = VecGetLocalSize(user->ParticleCount, &vec_local_size); CHKERRQ(ierr);
3272 ierr = VecGetArrayRead(user->ParticleCount, &count_array); CHKERRQ(ierr);
3273
3274 for (PetscInt i = 0; i < vec_local_size; ++i) {
3275 if (count_array[i] > 0.5) { // Use 0.5 to be safe with floating point
3276 local_occupied_cells++;
3277 }
3278 }
3279 ierr = VecRestoreArrayRead(user->ParticleCount, &count_array); CHKERRQ(ierr);
3280
3281 ierr = MPI_Allreduce(&local_occupied_cells, &global_occupied_cells, 1, MPIU_INT, MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
3282 simCtx->occupiedCellCount = global_occupied_cells;
3283
3284 LOG_ALLOW_SYNC(GLOBAL, LOG_INFO, "[Rank %d] Advanced Metrics: Imbalance=%.2f, OccupiedCells=%d\n", rank, simCtx->particleLoadImbalance, simCtx->occupiedCellCount);
3285
3286 PetscFunctionReturn(0);
3287}
3288
3289#undef __FUNCT__
3290#define __FUNCT__ "LOG_PARTICLE_METRICS"
3291/**
3292 * @brief Implementation of \ref LOG_PARTICLE_METRICS().
3293 * @details Full API contract (arguments, ownership, side effects) is documented with
3294 * the header declaration in `include/logging.h`.
3295 * @see LOG_PARTICLE_METRICS()
3296 */
3297PetscErrorCode LOG_PARTICLE_METRICS(UserCtx *user, const char *stageName)
3298{
3299 PetscErrorCode ierr;
3300 PetscMPIInt rank;
3301 SimCtx *simCtx = user->simCtx;
3302 const char *stage_label = (stageName && stageName[0] != '\0') ? stageName : "N/A";
3303
3304 PetscFunctionBeginUser;
3305 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
3306
3307 PetscInt totalParticles;
3308 ierr = DMSwarmGetSize(user->swarm, &totalParticles); CHKERRQ(ierr);
3309
3310 if (!rank) {
3311 FILE *f;
3312 char filen[PETSC_MAX_PATH_LEN + 64];
3313 ierr = PetscSNPrintf(filen, sizeof(filen), "%s/Particle_Metrics.log", simCtx->log_dir); CHKERRQ(ierr);
3314 f = fopen(filen, "a");
3315 if (!f) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Cannot open particle log file: %s", filen);
3316
3317 if (ftell(f) == 0) {
3318 PetscFPrintf(PETSC_COMM_SELF, f, "%-18s | %-10s | %-12s | %-10s | %-10s | %-10s | %-15s | %-10s | %-10s\n",
3319 "Stage", "Timestep", "Total Ptls", "Lost", "Lost Total", "Migrated", "Occupied Cells", "Imbalance", "Mig Passes");
3320 PetscFPrintf(PETSC_COMM_SELF, f, "-------------------------------------------------------------------------------------------------------------------------------------------\n");
3321 }
3322 if (simCtx->continueMode && simCtx->step == simCtx->StartStep + 1) {
3323 PetscFPrintf(PETSC_COMM_SELF, f, "# Continuation from step %" PetscInt_FMT "\n", simCtx->StartStep);
3324 }
3325
3326 PetscFPrintf(PETSC_COMM_SELF, f, "%-18s | %-10d | %-12d | %-10d | %-10d | %-10d | %-15d | %-10.2f | %-10d\n",
3327 stage_label, (int)simCtx->step, (int)totalParticles, (int)simCtx->particlesLostLastStep,
3328 (int)simCtx->particlesLostCumulative, (int)simCtx->particlesMigratedLastStep, (int)simCtx->occupiedCellCount,
3329 (double)simCtx->particleLoadImbalance, (int)simCtx->migrationPassesLastStep);
3330 fclose(f);
3331 }
3332 PetscFunctionReturn(0);
3333}
SetAnalyticalScalarFieldAtCellCenters
PetscErrorCode SetAnalyticalScalarFieldAtCellCenters(UserCtx *user, Vec targetVec)
Writes the configured verification scalar profile at physical cell centers into a scalar Vec.
Definition AnalyticalSolutions.c:672
SetAnalyticalSolutionForParticles
PetscErrorCode SetAnalyticalSolutionForParticles(Vec tempVec, SimCtx *simCtx)
Applies the analytical solution to particle velocity vector.
Definition AnalyticalSolutions.c:541
SEARCH_METRIC_SUM_SEARCH_LOCATED
@ SEARCH_METRIC_SUM_SEARCH_LOCATED
Definition logging.c:33
SEARCH_METRIC_SUM_MAX_TRAVERSAL_FAILS
@ SEARCH_METRIC_SUM_MAX_TRAVERSAL_FAILS
Definition logging.c:41
SEARCH_METRIC_REDUCTION_LEN
@ SEARCH_METRIC_REDUCTION_LEN
Definition logging.c:44
SEARCH_METRIC_SUM_BBOX_GUESS_FALLBACK
@ SEARCH_METRIC_SUM_BBOX_GUESS_FALLBACK
Definition logging.c:40
SEARCH_METRIC_SUM_RESEARCH
@ SEARCH_METRIC_SUM_RESEARCH
Definition logging.c:36
SEARCH_METRIC_SUM_TRAVERSAL_STEPS
@ SEARCH_METRIC_SUM_TRAVERSAL_STEPS
Definition logging.c:35
SEARCH_METRIC_MAX_TRAVERSAL_STEPS
@ SEARCH_METRIC_MAX_TRAVERSAL_STEPS
Definition logging.c:42
SEARCH_METRIC_SUM_BOUNDARY_CLAMPS
@ SEARCH_METRIC_SUM_BOUNDARY_CLAMPS
Definition logging.c:38
SEARCH_METRIC_SUM_SEARCH_POPULATION
@ SEARCH_METRIC_SUM_SEARCH_POPULATION
Definition logging.c:32
SEARCH_METRIC_MAX_PASS_DEPTH
@ SEARCH_METRIC_MAX_PASS_DEPTH
Definition logging.c:43
SEARCH_METRIC_SUM_TIE_BREAKS
@ SEARCH_METRIC_SUM_TIE_BREAKS
Definition logging.c:37
SEARCH_METRIC_SUM_SEARCH_LOST
@ SEARCH_METRIC_SUM_SEARCH_LOST
Definition logging.c:34
SEARCH_METRIC_SUM_BBOX_GUESS_SUCCESS
@ SEARCH_METRIC_SUM_BBOX_GUESS_SUCCESS
Definition logging.c:39
SEARCH_METRIC_SUM_SEARCH_ATTEMPTS
@ SEARCH_METRIC_SUM_SEARCH_ATTEMPTS
Definition logging.c:31
set_allowed_functions
void set_allowed_functions(const char **functionList, int count)
Implementation of set_allowed_functions().
Definition logging.c:152
ProfiledFunction::always_log
PetscBool always_log
Definition logging.c:1878
AppendStatisticalObservableSample
static PetscErrorCode AppendStatisticalObservableSample(SimCtx *simCtx, PetscInt samples_before, PetscReal mean_speed, PetscReal mean_ke)
Appends one timestep's scalar observables to the statistical history.
Definition logging.c:1369
LOG_PARTICLE_METRICS
PetscErrorCode LOG_PARTICLE_METRICS(UserCtx *user, const char *stageName)
Implementation of LOG_PARTICLE_METRICS().
Definition logging.c:3297
BCHandlerTypeToString
const char * BCHandlerTypeToString(BCHandlerType handler_type)
Internal helper implementation: BCHandlerTypeToString().
Definition logging.c:792
is_function_allowed
PetscBool is_function_allowed(const char *functionName)
Implementation of is_function_allowed().
Definition logging.c:183
g_profiler_count
static PetscInt g_profiler_count
Definition logging.c:1883
DualMonitorDestroy
PetscErrorCode DualMonitorDestroy(void **ctx)
Implementation of DualMonitorDestroy().
Definition logging.c:830
TMP_BUF_SIZE
#define TMP_BUF_SIZE
Definition logging.c:7
gAllowedFunctions
static char ** gAllowedFunctions
Global/static array of function names allowed to log.
Definition logging.c:23
current_log_level
static LogLevel current_log_level
Static variable to cache the current logging level.
Definition logging.c:16
LOG_INTERPOLATION_ERROR
PetscErrorCode LOG_INTERPOLATION_ERROR(UserCtx *user)
Implementation of LOG_INTERPOLATION_ERROR().
Definition logging.c:2825
SolutionConvergenceHistoryGet
static PetscReal SolutionConvergenceHistoryGet(const PetscReal *history, PetscInt capacity, PetscInt samples_available, PetscInt offset_from_latest)
Reads one sample from the statistical ring buffer by age.
Definition logging.c:1334
ShouldEmitPeriodicParticleConsoleSnapshot
PetscBool ShouldEmitPeriodicParticleConsoleSnapshot(const SimCtx *simCtx, PetscInt completed_step)
Implementation of ShouldEmitPeriodicParticleConsoleSnapshot().
Definition logging.c:542
BCFaceToString
const char * BCFaceToString(BCFace face)
Implementation of BCFaceToString().
Definition logging.c:669
FreeAllowedFunctions
PetscErrorCode FreeAllowedFunctions(char **funcs, PetscInt n)
Internal helper implementation: FreeAllowedFunctions().
Definition logging.c:650
IsParticleConsoleSnapshotEnabled
PetscBool IsParticleConsoleSnapshotEnabled(const SimCtx *simCtx)
Implementation of IsParticleConsoleSnapshotEnabled().
Definition logging.c:525
print_log_level
PetscErrorCode print_log_level(void)
Internal helper implementation: print_log_level().
Definition logging.c:116
ProfiledFunction::total_call_count
long long total_call_count
Definition logging.c:1875
SolutionConvergenceSafeRelative
static PetscReal SolutionConvergenceSafeRelative(PetscReal numerator, PetscReal denominator)
Forms a guarded relative metric for solution-convergence logging.
Definition logging.c:948
EmitParticleConsoleSnapshot
PetscErrorCode EmitParticleConsoleSnapshot(UserCtx *user, SimCtx *simCtx, PetscInt step)
Implementation of EmitParticleConsoleSnapshot().
Definition logging.c:556
g_profiler_capacity
static PetscInt g_profiler_capacity
Definition logging.c:1884
ProfilingFinalize
PetscErrorCode ProfilingFinalize(SimCtx *simCtx)
Implementation of ProfilingFinalize().
Definition logging.c:2196
BuildRowFormatString
static void BuildRowFormatString(PetscMPIInt wRank, PetscInt wPID, PetscInt wCell, PetscInt wPos, PetscInt wVel, PetscInt wWt, char *fmtStr, size_t bufSize)
Definition logging.c:366
trim
static void trim(char *s)
Internal helper implementation: trim().
Definition logging.c:571
LoadAllowedFunctionsFromFile
PetscErrorCode LoadAllowedFunctionsFromFile(const char filename[], char ***funcsOut, PetscInt *nOut)
Implementation of LoadAllowedFunctionsFromFile().
Definition logging.c:596
SOLUTION_CONVERGENCE_REL_EPS
#define SOLUTION_CONVERGENCE_REL_EPS
Definition logging.c:917
LOG_FIELD_MIN_MAX
PetscErrorCode LOG_FIELD_MIN_MAX(UserCtx *user, const char *fieldName)
Implementation of LOG_FIELD_MIN_MAX().
Definition logging.c:2349
gNumAllowed
static int gNumAllowed
Number of entries in the gAllowedFunctions array.
Definition logging.c:28
ProfiledFunction::total_time
double total_time
Definition logging.c:1873
BuildHeaderString
static void BuildHeaderString(char *headerStr, size_t bufSize, PetscMPIInt wRank, PetscInt wPID, PetscInt wCell, PetscInt wPos, PetscInt wVel, PetscInt wWt)
Definition logging.c:379
PrintProgressBar
void PrintProgressBar(PetscInt step, PetscInt startStep, PetscInt totalSteps, PetscReal currentTime)
Internal helper implementation: PrintProgressBar().
Definition logging.c:2302
SolutionConvergenceModeToString
static const char * SolutionConvergenceModeToString(SolutionConvergenceMode mode)
Maps the internal solution-convergence mode enum to its log label.
Definition logging.c:1582
LOG_FIELD_ANATOMY
PetscErrorCode LOG_FIELD_ANATOMY(UserCtx *user, const char *field_name, const char *stage_name)
Implementation of LOG_FIELD_ANATOMY().
Definition logging.c:2536
_FindOrCreateEntry
static PetscErrorCode _FindOrCreateEntry(const char *func_name, PetscInt *idx)
Internal helper implementation: _FindOrCreateEntry().
Definition logging.c:1891
CellToStr
static void CellToStr(const PetscInt *cell, char *buf, size_t bufsize)
Definition logging.c:269
RuntimeMemoryLogSample
PetscErrorCode RuntimeMemoryLogSample(SimCtx *simCtx, PetscInt step, const char *event, const char *reason)
Implementation of RuntimeMemoryLogSample().
Definition logging.c:2085
get_log_level
LogLevel get_log_level()
Implementation of get_log_level().
Definition logging.c:84
ProfilingLogTimestepSummary
PetscErrorCode ProfilingLogTimestepSummary(SimCtx *simCtx, PetscInt step)
Implementation of ProfilingLogTimestepSummary().
Definition logging.c:2004
LOG_FACE_DISTANCES
PetscErrorCode LOG_FACE_DISTANCES(PetscReal *d)
Implementation of LOG_FACE_DISTANCES().
Definition logging.c:230
LOG_PARTICLE_FIELDS
PetscErrorCode LOG_PARTICLE_FIELDS(UserCtx *user, PetscInt printInterval)
Implementation of LOG_PARTICLE_FIELDS().
Definition logging.c:397
_ProfilingEnd
void _ProfilingEnd(const char *func_name)
Implementation of _ProfilingEnd().
Definition logging.c:1965
TripleRealToStr
static void TripleRealToStr(const PetscReal *arr, char *buf, size_t bufsize)
Definition logging.c:277
Int64ToStr
static void Int64ToStr(PetscInt64 value, char *buf, size_t bufsize)
Definition logging.c:261
BCTypeToString
const char * BCTypeToString(BCType type)
Implementation of BCTypeToString().
Definition logging.c:772
CalculateAdvancedParticleMetrics
PetscErrorCode CalculateAdvancedParticleMetrics(UserCtx *user)
Internal helper implementation: CalculateAdvancedParticleMetrics().
Definition logging.c:3243
ParticleLocationStatusToString
const char * ParticleLocationStatusToString(ParticleLocationStatus level)
Implementation of ParticleLocationStatusToString().
Definition logging.c:1856
LOG_SCATTER_METRICS
PetscErrorCode LOG_SCATTER_METRICS(UserCtx *user)
Implementation of LOG_SCATTER_METRICS().
Definition logging.c:2904
_CompareProfiledFunctions
static int _CompareProfiledFunctions(const void *a, const void *b)
Internal helper implementation: _CompareProfiledFunctions().
Definition logging.c:2180
LOG_SOLUTION_CONVERGENCE
PetscErrorCode LOG_SOLUTION_CONVERGENCE(SimCtx *simCtx)
Implementation of LOG_SOLUTION_CONVERGENCE().
Definition logging.c:1599
FlowDirectionToString
const char * FlowDirectionToString(FlowDirection fd)
Convert a FlowDirection enum value to its token string.
Definition logging.c:704
DualKSPMonitor
PetscErrorCode DualKSPMonitor(KSP ksp, PetscInt it, PetscReal rnorm, void *ctx)
Implementation of DualKSPMonitor().
Definition logging.c:869
LOG_CONTINUITY_METRICS
PetscErrorCode LOG_CONTINUITY_METRICS(UserCtx *user)
Implementation of LOG_CONTINUITY_METRICS().
Definition logging.c:1794
ProfiledFunction::current_step_time
double current_step_time
Definition logging.c:1874
ProfiledFunction::current_step_call_count
long long current_step_call_count
Definition logging.c:1876
ComputeStatisticalWindowMetrics
static PetscErrorCode ComputeStatisticalWindowMetrics(const SimCtx *simCtx, PetscInt samples_available, PetscBool *has_reference_out, PetscReal *mean_speed_window_out, PetscReal *mean_speed_window_prev_out, PetscReal *mean_speed_window_abs_out, PetscReal *mean_speed_window_rel_out, PetscReal *mean_speed_rms_window_out, PetscReal *mean_speed_rms_window_prev_out, PetscReal *mean_speed_rms_window_abs_out, PetscReal *mean_speed_rms_window_rel_out, PetscReal *mean_ke_window_out, PetscReal *mean_ke_window_prev_out, PetscReal *mean_ke_window_abs_out, PetscReal *mean_ke_window_rel_out, PetscReal *mean_ke_rms_window_out, PetscReal *mean_ke_rms_window_prev_out, PetscReal *mean_ke_rms_window_abs_out, PetscReal *mean_ke_rms_window_rel_out)
Computes adjacent-window drift metrics for statistical steady mode.
Definition logging.c:1449
LOG_SEARCH_METRICS
PetscErrorCode LOG_SEARCH_METRICS(UserCtx *user)
Implementation of LOG_SEARCH_METRICS().
Definition logging.c:3089
g_profiler_registry
static ProfiledFunction * g_profiler_registry
Definition logging.c:1882
InitialConditionModeToString
const char * InitialConditionModeToString(InitialConditionMode mode)
Implementation of InitialConditionModeToString().
Definition logging.c:687
ProfilingInitialize
PetscErrorCode ProfilingInitialize(SimCtx *simCtx)
Internal helper implementation: ProfilingInitialize().
Definition logging.c:1927
LESModelToString
const char * LESModelToString(LESModelType LESFlag)
Implementation of LESModelToString().
Definition logging.c:740
LOG_CELL_VERTICES
PetscErrorCode LOG_CELL_VERTICES(const Cell *cell, PetscMPIInt rank)
Implementation of LOG_CELL_VERTICES().
Definition logging.c:205
IntToStr
static void IntToStr(int value, char *buf, size_t bufsize)
Definition logging.c:253
ComputeDeterministicSolutionMetrics
static PetscErrorCode ComputeDeterministicSolutionMetrics(SimCtx *simCtx, PetscBool periodic_mode, PetscInt phase_step, PetscInt samples_before, PetscBool *has_reference_out, PetscReal *u_abs_l2_out, PetscReal *u_rel_l2_out, PetscReal *p_abs_l2_out, PetscReal *p_rel_l2_out, PetscReal *mean_speed_out, PetscReal *mean_speed_ref_out, PetscReal *mean_speed_abs_out, PetscReal *mean_speed_rel_out, PetscReal *mean_ke_out, PetscReal *mean_ke_ref_out, PetscReal *mean_ke_abs_out, PetscReal *mean_ke_rel_out)
Computes deterministic solution-drift metrics for the current step.
Definition logging.c:1091
ProfilingResetTimestepCounters
PetscErrorCode ProfilingResetTimestepCounters(void)
Implementation of ProfilingResetTimestepCounters().
Definition logging.c:1987
ResetSearchMetrics
PetscErrorCode ResetSearchMetrics(SimCtx *simCtx)
Implementation of ResetSearchMetrics().
Definition logging.c:3058
MomentumSolverTypeToString
const char * MomentumSolverTypeToString(MomentumSolverType SolverFlag)
Implementation of MomentumSolverTypeToString().
Definition logging.c:756
ComputeMaxColumnWidths
static PetscErrorCode ComputeMaxColumnWidths(PetscInt nParticles, const PetscMPIInt *ranks, const PetscInt64 *pids, const PetscInt *cellIDs, const PetscReal *positions, const PetscReal *velocities, const PetscReal *weights, int *wRank, int *wPID, int *wCell, int *wPos, int *wVel, int *wWt)
Definition logging.c:302
ComputeCurrentFlowObservables
static PetscErrorCode ComputeCurrentFlowObservables(SimCtx *simCtx, PetscReal *mean_speed_out, PetscReal *mean_ke_out)
Computes instantaneous global flow observables for statistical mode.
Definition logging.c:975
SOLUTION_CONVERGENCE_FLUID_THRESHOLD
#define SOLUTION_CONVERGENCE_FLUID_THRESHOLD
Definition logging.c:916
ProfiledFunction::start_time
double start_time
Definition logging.c:1877
ParticleInitializationToString
const char * ParticleInitializationToString(ParticleInitializationType ParticleInitialization)
Implementation of ParticleInitializationToString().
Definition logging.c:723
_ProfilingStart
void _ProfilingStart(const char *func_name)
Implementation of _ProfilingStart().
Definition logging.c:1951
SearchMetricsReduceOp
static void SearchMetricsReduceOp(void *invec, void *inoutvec, int *len, MPI_Datatype *datatype)
Internal reduction callback for packed search metrics.
Definition logging.c:51
ProfiledFunction::name
const char * name
Definition logging.c:1872
logging.h
Logging utilities and macros for PETSc-based applications.
LOG_ALLOW_SYNC
#define LOG_ALLOW_SYNC(scope, level, fmt,...)
Synchronized logging macro that checks both the log level and whether the calling function is in the ...
Definition logging.h:252
DualMonitorCtx::log_to_console
PetscBool log_to_console
Definition logging.h:57
LOCAL
#define LOCAL
Logging scope definitions for controlling message output.
Definition logging.h:44
GLOBAL
#define GLOBAL
Scope for global logging across all processes.
Definition logging.h:45
LOG_ALLOW
#define LOG_ALLOW(scope, level, fmt,...)
Logging macro that checks both the log level and whether the calling function is in the allowed-funct...
Definition logging.h:199
DualMonitorCtx::bnorm
PetscReal bnorm
Definition logging.h:58
DualMonitorCtx::step
PetscInt step
Definition logging.h:59
LOG
#define LOG(scope, level, fmt,...)
Logging macro for PETSc-based applications with scope control.
Definition logging.h:83
LogLevel
LogLevel
Enumeration of logging levels.
Definition logging.h:27
LOG_ERROR
@ LOG_ERROR
Critical errors that may halt the program.
Definition logging.h:28
LOG_TRACE
@ LOG_TRACE
Very fine-grained tracing information for in-depth debugging.
Definition logging.h:32
LOG_INFO
@ LOG_INFO
Informational messages about program execution.
Definition logging.h:30
LOG_WARNING
@ LOG_WARNING
Non-critical issues that warrant attention.
Definition logging.h:29
LOG_DEBUG
@ LOG_DEBUG
Detailed debugging information.
Definition logging.h:31
LOG_VERBOSE
@ LOG_VERBOSE
Extremely detailed logs, typically for development use only.
Definition logging.h:33
DualMonitorCtx::file_handle
FILE * file_handle
Definition logging.h:56
DualMonitorCtx::block_id
PetscInt block_id
Definition logging.h:60
DualMonitorCtx
Context for a dual-purpose KSP monitor.
Definition logging.h:55
LESModelType
LESModelType
Identifies the six logical faces of a structured computational block.
Definition variables.h:518
DYNAMIC_SMAGORINSKY
@ DYNAMIC_SMAGORINSKY
Definition variables.h:521
NO_LES_MODEL
@ NO_LES_MODEL
Definition variables.h:519
CONSTANT_SMAGORINSKY
@ CONSTANT_SMAGORINSKY
Definition variables.h:520
UserCtx::lDiffusivityGradient
Vec lDiffusivityGradient
Definition variables.h:908
UserCtx::lCent
Vec lCent
Definition variables.h:927
BCType
BCType
Defines the general mathematical/physical Category of a boundary.
Definition variables.h:281
INLET
@ INLET
Definition variables.h:288
INTERFACE
@ INTERFACE
Definition variables.h:283
FARFIELD
@ FARFIELD
Definition variables.h:289
OUTLET
@ OUTLET
Definition variables.h:287
PERIODIC
@ PERIODIC
Definition variables.h:290
WALL
@ WALL
Definition variables.h:284
SimCtx::continueMode
PetscBool continueMode
Definition variables.h:701
UserCtx::JCsi
Vec JCsi
Definition variables.h:931
UserCtx::KAj
Vec KAj
Definition variables.h:932
MGCtx::user
UserCtx * user
Definition variables.h:569
UserCtx::JEta
Vec JEta
Definition variables.h:931
UserCtx::Zet
Vec Zet
Definition variables.h:927
SimCtx::profilingFinalSummary
PetscBool profilingFinalSummary
Definition variables.h:837
SimCtx::rank
PetscMPIInt rank
Definition variables.h:687
SimCtx::profilingTimestepFile
char profilingTimestepFile[PETSC_MAX_PATH_LEN]
Definition variables.h:836
SearchMetricsState::searchLocatedCount
PetscInt64 searchLocatedCount
Definition variables.h:239
SimCtx::block_number
PetscInt block_number
Definition variables.h:768
UserCtx::lIEta
Vec lIEta
Definition variables.h:930
SearchMetricsState::searchLostCount
PetscInt64 searchLostCount
Definition variables.h:240
UserCtx::solutionConvergencePeriodicPRef
Vec * solutionConvergencePeriodicPRef
Definition variables.h:914
UserCtx::lIZet
Vec lIZet
Definition variables.h:930
UserCtx::Phi
Vec Phi
Definition variables.h:904
UserCtx::simCtx
SimCtx * simCtx
Back-pointer to the master simulation context.
Definition variables.h:879
ParticleInitializationType
ParticleInitializationType
Enumerator to identify the particle initialization strategy.
Definition variables.h:549
PARTICLE_INIT_SURFACE_RANDOM
@ PARTICLE_INIT_SURFACE_RANDOM
Random placement on the inlet face.
Definition variables.h:550
PARTICLE_INIT_SURFACE_EDGES
@ PARTICLE_INIT_SURFACE_EDGES
Deterministic placement at inlet face edges.
Definition variables.h:553
PARTICLE_INIT_POINT_SOURCE
@ PARTICLE_INIT_POINT_SOURCE
All particles at a fixed (psrc_x,psrc_y,psrc_z) — for validation.
Definition variables.h:552
PARTICLE_INIT_VOLUME
@ PARTICLE_INIT_VOLUME
Random volumetric distribution across the domain.
Definition variables.h:551
ParticleLocationStatus
ParticleLocationStatus
Defines the state of a particle with respect to its location and migration status during the iterativ...
Definition variables.h:135
LOST
@ LOST
Definition variables.h:139
NEEDS_LOCATION
@ NEEDS_LOCATION
Definition variables.h:136
ACTIVE_AND_LOCATED
@ ACTIVE_AND_LOCATED
Definition variables.h:137
UNINITIALIZED
@ UNINITIALIZED
Definition variables.h:140
MIGRATING_OUT
@ MIGRATING_OUT
Definition variables.h:138
SimCtx::solutionConvergenceMeanSpeedHistory
PetscReal * solutionConvergenceMeanSpeedHistory
Definition variables.h:753
SimCtx::FluxOutSum
PetscReal FluxOutSum
Definition variables.h:777
UserCtx::IZet
Vec IZet
Definition variables.h:930
UserCtx::Centz
Vec Centz
Definition variables.h:928
SimCtx::runtimeMemoryLogEnabled
PetscBool runtimeMemoryLogEnabled
Enable the rank-reduced runtime memory log.
Definition variables.h:852
UserCtx::IEta
Vec IEta
Definition variables.h:930
SearchMetricsState::boundaryClampCount
PetscInt64 boundaryClampCount
Definition variables.h:246
SimCtx::particlesLostLastStep
PetscInt particlesLostLastStep
Definition variables.h:803
UserCtx::KM
PetscInt KM
Definition variables.h:885
UserCtx::lZet
Vec lZet
Definition variables.h:927
SimCtx::usermg
UserMG usermg
Definition variables.h:821
UserCtx::Csi
Vec Csi
Definition variables.h:927
UserCtx::solutionConvergencePeriodicUcatRef
Vec * solutionConvergencePeriodicUcatRef
Definition variables.h:913
SearchMetricsState::traversalStepsSum
PetscInt64 traversalStepsSum
Definition variables.h:241
BCHandlerType
BCHandlerType
Defines the specific computational "strategy" for a boundary handler.
Definition variables.h:301
BC_HANDLER_INLET_PULSATILE_FLUX
@ BC_HANDLER_INLET_PULSATILE_FLUX
Definition variables.h:310
BC_HANDLER_PERIODIC_GEOMETRIC
@ BC_HANDLER_PERIODIC_GEOMETRIC
Definition variables.h:314
BC_HANDLER_INLET_PARABOLIC
@ BC_HANDLER_INLET_PARABOLIC
Definition variables.h:307
BC_HANDLER_INLET_CONSTANT_VELOCITY
@ BC_HANDLER_INLET_CONSTANT_VELOCITY
Definition variables.h:306
BC_HANDLER_PERIODIC_DRIVEN_INITIAL_FLUX
@ BC_HANDLER_PERIODIC_DRIVEN_INITIAL_FLUX
Definition variables.h:317
BC_HANDLER_INTERFACE_OVERSET
@ BC_HANDLER_INTERFACE_OVERSET
Definition variables.h:315
BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX
@ BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX
Definition variables.h:316
BC_HANDLER_WALL_MOVING
@ BC_HANDLER_WALL_MOVING
Definition variables.h:304
BC_HANDLER_INLET_PROFILE_FROM_FILE
@ BC_HANDLER_INLET_PROFILE_FROM_FILE
Definition variables.h:308
BC_HANDLER_WALL_NOSLIP
@ BC_HANDLER_WALL_NOSLIP
Definition variables.h:303
BC_HANDLER_OUTLET_CONSERVATION
@ BC_HANDLER_OUTLET_CONSERVATION
Definition variables.h:312
BC_HANDLER_FARFIELD_NONREFLECTING
@ BC_HANDLER_FARFIELD_NONREFLECTING
Definition variables.h:311
BC_HANDLER_OUTLET_PRESSURE
@ BC_HANDLER_OUTLET_PRESSURE
Definition variables.h:313
BC_HANDLER_SYMMETRY_PLANE
@ BC_HANDLER_SYMMETRY_PLANE
Definition variables.h:305
BC_HANDLER_UNDEFINED
@ BC_HANDLER_UNDEFINED
Definition variables.h:302
UserCtx::lUcont_rm1
Vec lUcont_rm1
Definition variables.h:912
SimCtx::solutionConvergenceSamplesRecorded
PetscInt solutionConvergenceSamplesRecorded
Definition variables.h:752
UserCtx::lCellFieldAtCorner
Vec lCellFieldAtCorner
Definition variables.h:915
UserCtx::_this
PetscInt _this
Definition variables.h:889
UserCtx::lKEta
Vec lKEta
Definition variables.h:932
SearchMetricsState::searchPopulation
PetscInt64 searchPopulation
Definition variables.h:238
SimCtx::solutionConvergenceMeanKEHistory
PetscReal * solutionConvergenceMeanKEHistory
Definition variables.h:754
SimCtx::dt
PetscReal dt
Definition variables.h:699
SimCtx::runtimeMemoryLogFile
char runtimeMemoryLogFile[PETSC_MAX_PATH_LEN]
File name written under log_dir.
Definition variables.h:853
SimCtx::runtimeMemoryLogStarted
PetscBool runtimeMemoryLogStarted
True after rank 0 writes the log header.
Definition variables.h:854
SimCtx::occupiedCellCount
PetscInt occupiedCellCount
Definition variables.h:807
SimCtx::profilingTimestepMode
char profilingTimestepMode[32]
Definition variables.h:835
SimCtx::k_periodic
PetscInt k_periodic
Definition variables.h:769
UserCtx::lPsi
Vec lPsi
Definition variables.h:953
SearchMetricsState::currentSettlementPass
PetscInt currentSettlementPass
Definition variables.h:250
SimCtx::np
PetscInt np
Definition variables.h:796
UserCtx::DiffusivityGradient
Vec DiffusivityGradient
Definition variables.h:908
UserCtx::lJCsi
Vec lJCsi
Definition variables.h:931
SimCtx::StartStep
PetscInt StartStep
Definition variables.h:694
MomentumSolverType
MomentumSolverType
Enumerator to identify the implemented momentum solver strategies.
Definition variables.h:532
MOMENTUM_SOLVER_DUALTIME_PICARD_JAMESON_RK
@ MOMENTUM_SOLVER_DUALTIME_PICARD_JAMESON_RK
Definition variables.h:534
MOMENTUM_SOLVER_EXPLICIT_RK
@ MOMENTUM_SOLVER_EXPLICIT_RK
Definition variables.h:533
MOMENTUM_SOLVER_NEWTON_KRYLOV
@ MOMENTUM_SOLVER_NEWTON_KRYLOV
Definition variables.h:535
SimCtx::solutionConvergencePeriodSteps
PetscInt solutionConvergencePeriodSteps
Definition variables.h:750
Cmpnts::x
PetscScalar x
Definition variables.h:101
UserCtx::JZet
Vec JZet
Definition variables.h:931
SearchMetricsState::reSearchCount
PetscInt64 reSearchCount
Definition variables.h:242
SimCtx::MaxDiv
PetscReal MaxDiv
Definition variables.h:828
UserCtx::Centx
Vec Centx
Definition variables.h:928
UserCtx::lUcont_o
Vec lUcont_o
Definition variables.h:911
SearchMetricsState::bboxGuessFallbackCount
PetscInt64 bboxGuessFallbackCount
Definition variables.h:248
UserCtx::Ucat_o
Vec Ucat_o
Definition variables.h:911
SimCtx::MaxDivx
PetscInt MaxDivx
Definition variables.h:829
SimCtx::MaxDivy
PetscInt MaxDivy
Definition variables.h:829
SearchMetricsState::bboxGuessSuccessCount
PetscInt64 bboxGuessSuccessCount
Definition variables.h:247
SimCtx::MaxDivz
PetscInt MaxDivz
Definition variables.h:829
UserCtx::lKZet
Vec lKZet
Definition variables.h:932
UserCtx::Eta
Vec Eta
Definition variables.h:927
SimCtx::log_dir
char log_dir[PETSC_MAX_PATH_LEN]
Definition variables.h:709
SimCtx::MaxDivFlatArg
PetscInt MaxDivFlatArg
Definition variables.h:829
SimCtx::FluxInSum
PetscReal FluxInSum
Definition variables.h:777
SearchMetricsState::maxParticlePassDepth
PetscInt64 maxParticlePassDepth
Definition variables.h:249
UserCtx::lJEta
Vec lJEta
Definition variables.h:931
UserCtx::lCsi
Vec lCsi
Definition variables.h:927
SearchMetricsState::maxTraversalSteps
PetscInt64 maxTraversalSteps
Definition variables.h:243
UserCtx::ICsi
Vec ICsi
Definition variables.h:930
Cmpnts::z
PetscScalar z
Definition variables.h:101
UserCtx::lKCsi
Vec lKCsi
Definition variables.h:932
UserCtx::Ucat
Vec Ucat
Definition variables.h:904
UserCtx::ParticleCount
Vec ParticleCount
Definition variables.h:952
UserCtx::CellFieldAtCorner
Vec CellFieldAtCorner
Definition variables.h:915
UserCtx::JM
PetscInt JM
Definition variables.h:885
UserCtx::lCenty
Vec lCenty
Definition variables.h:929
SimCtx::runtimeMemoryLogHasPrevious
PetscBool runtimeMemoryLogHasPrevious
True after the first process-memory sample.
Definition variables.h:855
UserMG::mglevels
PetscInt mglevels
Definition variables.h:576
SimCtx::profilingSelectedFuncs
char ** profilingSelectedFuncs
Definition variables.h:833
SimCtx::solutionConvergenceWindowSteps
PetscInt solutionConvergenceWindowSteps
Definition variables.h:751
UserCtx::lJZet
Vec lJZet
Definition variables.h:931
FlowDirection
FlowDirection
Primary flow direction for streamwise IC and Poiseuille modes.
Definition variables.h:270
FLOW_DIR_NEG_ZETA
@ FLOW_DIR_NEG_ZETA
Definition variables.h:276
FLOW_DIR_NEG_ETA
@ FLOW_DIR_NEG_ETA
Definition variables.h:274
FLOW_DIR_POS_ZETA
@ FLOW_DIR_POS_ZETA
Definition variables.h:275
FLOW_DIR_POS_XI
@ FLOW_DIR_POS_XI
Definition variables.h:271
FLOW_DIR_NEG_XI
@ FLOW_DIR_NEG_XI
Definition variables.h:272
FLOW_DIR_POS_ETA
@ FLOW_DIR_POS_ETA
Definition variables.h:273
SimCtx::particlesLostCumulative
PetscInt particlesLostCumulative
Definition variables.h:804
SimCtx::nProfilingSelectedFuncs
PetscInt nProfilingSelectedFuncs
Definition variables.h:834
UserCtx::IAj
Vec IAj
Definition variables.h:930
SimCtx::particlesMigratedLastStep
PetscInt particlesMigratedLastStep
Definition variables.h:806
UserCtx::JAj
Vec JAj
Definition variables.h:931
UserCtx::KEta
Vec KEta
Definition variables.h:932
InitialConditionMode
InitialConditionMode
Selects the algorithm used to populate a fresh Eulerian velocity field.
Definition variables.h:149
IC_MODE_CONSTANT_CARTESIAN
@ IC_MODE_CONSTANT_CARTESIAN
Definition variables.h:151
IC_MODE_POISEUILLE
@ IC_MODE_POISEUILLE
Definition variables.h:152
IC_MODE_CONSTANT_STREAMWISE
@ IC_MODE_CONSTANT_STREAMWISE
Definition variables.h:153
IC_MODE_FILE
@ IC_MODE_FILE
Definition variables.h:154
IC_MODE_ZERO
@ IC_MODE_ZERO
Definition variables.h:150
SimCtx::particleConsoleOutputFreq
PetscInt particleConsoleOutputFreq
Definition variables.h:697
UserCtx::lCentx
Vec lCentx
Definition variables.h:929
SimCtx::searchMetrics
SearchMetricsState searchMetrics
Definition variables.h:809
SimCtx::i_periodic
PetscInt i_periodic
Definition variables.h:769
UserCtx::lUcont
Vec lUcont
Definition variables.h:904
SimCtx::runtimeMemoryLogPreviousProcessMB
PetscReal runtimeMemoryLogPreviousProcessMB
Previous local process memory sample in MB.
Definition variables.h:856
SimCtx::step
PetscInt step
Definition variables.h:692
UserCtx::Diffusivity
Vec Diffusivity
Definition variables.h:907
UserCtx::lICsi
Vec lICsi
Definition variables.h:930
UserCtx::info
DMDALocalInfo info
Definition variables.h:883
UserCtx::lUcat
Vec lUcat
Definition variables.h:904
SimCtx::migrationPassesLastStep
PetscInt migrationPassesLastStep
Definition variables.h:805
Cmpnts::y
PetscScalar y
Definition variables.h:101
EXEC_MODE_SOLVER
@ EXEC_MODE_SOLVER
Definition variables.h:657
EXEC_MODE_POSTPROCESSOR
@ EXEC_MODE_POSTPROCESSOR
Definition variables.h:658
UserCtx::IM
PetscInt IM
Definition variables.h:885
TOP
@ TOP
Definition variables.h:145
FRONT
@ FRONT
Definition variables.h:145
BOTTOM
@ BOTTOM
Definition variables.h:145
BACK
@ BACK
Definition variables.h:145
LEFT
@ LEFT
Definition variables.h:145
RIGHT
@ RIGHT
Definition variables.h:145
UserCtx::lEta
Vec lEta
Definition variables.h:927
UserCtx::KZet
Vec KZet
Definition variables.h:932
UserCtx::Cent
Vec Cent
Definition variables.h:927
UserCtx::Nvert
Vec Nvert
Definition variables.h:904
UserCtx::KCsi
Vec KCsi
Definition variables.h:932
UserMG::mgctx
MGCtx * mgctx
Definition variables.h:579
SolutionConvergenceMode
SolutionConvergenceMode
Selects the runtime solution-convergence diagnostics mode.
Definition variables.h:541
SOLUTION_CONVERGENCE_TRANSIENT
@ SOLUTION_CONVERGENCE_TRANSIENT
Definition variables.h:545
SOLUTION_CONVERGENCE_PERIODIC_DETERMINISTIC
@ SOLUTION_CONVERGENCE_PERIODIC_DETERMINISTIC
Definition variables.h:543
SOLUTION_CONVERGENCE_STATISTICAL_STEADY
@ SOLUTION_CONVERGENCE_STATISTICAL_STEADY
Definition variables.h:544
SOLUTION_CONVERGENCE_STEADY_DETERMINISTIC
@ SOLUTION_CONVERGENCE_STEADY_DETERMINISTIC
Definition variables.h:542
UserCtx::lDiffusivity
Vec lDiffusivity
Definition variables.h:907
UserCtx::Centy
Vec Centy
Definition variables.h:928
SimCtx::solutionConvergenceMode
SolutionConvergenceMode solutionConvergenceMode
Definition variables.h:749
UserCtx::lCentz
Vec lCentz
Definition variables.h:929
SearchMetricsState::searchAttempts
PetscInt64 searchAttempts
Definition variables.h:237
SimCtx::exec_mode
ExecutionMode exec_mode
Definition variables.h:703
SearchMetricsState::tieBreakCount
PetscInt64 tieBreakCount
Definition variables.h:245
SimCtx::ti
PetscReal ti
Definition variables.h:693
SimCtx::summationRHS
PetscReal summationRHS
Definition variables.h:827
SearchMetricsState::maxTraversalFailCount
PetscInt64 maxTraversalFailCount
Definition variables.h:244
SimCtx::LoggingFrequency
PetscInt LoggingFrequency
Definition variables.h:826
Cell::vertices
Cmpnts vertices[8]
Coordinates of the eight vertices of the cell.
Definition variables.h:176
UserCtx::Psi
Vec Psi
Definition variables.h:953
SimCtx::particleLoadImbalance
PetscReal particleLoadImbalance
Definition variables.h:808
UserCtx::P_o
Vec P_o
Definition variables.h:911
BCFace
BCFace
Identifies the six logical faces of a structured computational block.
Definition variables.h:259
BC_FACE_NEG_X
@ BC_FACE_NEG_X
Definition variables.h:260
BC_FACE_POS_Z
@ BC_FACE_POS_Z
Definition variables.h:262
BC_FACE_POS_Y
@ BC_FACE_POS_Y
Definition variables.h:261
BC_FACE_NEG_Z
@ BC_FACE_NEG_Z
Definition variables.h:262
BC_FACE_POS_X
@ BC_FACE_POS_X
Definition variables.h:260
BC_FACE_NEG_Y
@ BC_FACE_NEG_Y
Definition variables.h:261
SimCtx::j_periodic
PetscInt j_periodic
Definition variables.h:769
Cell
Defines the vertices of a single hexahedral grid cell.
Definition variables.h:175
Cmpnts
A 3D point or vector with PetscScalar components.
Definition variables.h:100
SimCtx
The master context for the entire simulation.
Definition variables.h:684
UserCtx
User-defined context containing data specific to a single computational grid level.
Definition variables.h:876
verification_sources.h
VerificationScalarOverrideActive
PetscBool VerificationScalarOverrideActive(const SimCtx *simCtx)
Reports whether a verification-only scalar override is active.
Definition verification_sources.c:21
Generated by doxygen 1.9.8