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 FieldInitializationToString().
683 * @details Full API contract (arguments, ownership, side effects) is documented with
684 * the header declaration in `include/logging.h`.
685 * @see FieldInitializationToString()
686 */
687const char* FieldInitializationToString(PetscInt FieldInitialization)
688{
689 switch(FieldInitialization){
690 case 0: return "Zero";
691 case 1: return "Constant Normal velocity";
692 case 2: return "Poiseuille Normal velocity";
693 default: return "Unknown Field Initialization";
694 }
695}
696
697/**
698 * @brief Implementation of \ref ParticleInitializationToString().
699 * @details Full API contract (arguments, ownership, side effects) is documented with
700 * the header declaration in `include/logging.h`.
701 * @see ParticleInitializationToString()
702 */
703const char* ParticleInitializationToString(ParticleInitializationType ParticleInitialization)
704{
705 switch(ParticleInitialization){
706 case PARTICLE_INIT_SURFACE_RANDOM: return "Surface: Random";
707 case PARTICLE_INIT_VOLUME: return "Volume";
708 case PARTICLE_INIT_POINT_SOURCE: return "Point Source";
709 case PARTICLE_INIT_SURFACE_EDGES: return "Surface: At edges";
710 default: return "Unknown Particle Initialization";
711 }
712}
713
714/**
715 * @brief Implementation of \ref LESModelToString().
716 * @details Full API contract (arguments, ownership, side effects) is documented with
717 * the header declaration in `include/logging.h`.
718 * @see LESModelToString()
719 */
720const char* LESModelToString(LESModelType LESFlag)
721{
722 switch(LESFlag){
723 case NO_LES_MODEL: return "No LES";
724 case CONSTANT_SMAGORINSKY: return "Constant Smagorinsky";
725 case DYNAMIC_SMAGORINSKY: return "Dynamic Smagorinsky";
726 default: return "Unknown LES Flag";
727 }
728}
729
730/**
731 * @brief Implementation of \ref MomentumSolverTypeToString().
732 * @details Full API contract (arguments, ownership, side effects) is documented with
733 * the header declaration in `include/logging.h`.
734 * @see MomentumSolverTypeToString()
735 */
736const char* MomentumSolverTypeToString(MomentumSolverType SolverFlag)
737{
738 switch(SolverFlag){
739 case MOMENTUM_SOLVER_EXPLICIT_RK: return "Explicit 4 stage Runge-Kutta ";
740 case MOMENTUM_SOLVER_DUALTIME_PICARD_RK4: return "Dual Time Stepping with Picard Iterations and 4 stage Runge-Kutta Smoothing";
741 default: return "Unknown Momentum Solver Type";
742 }
743}
744
745/**
746 * @brief Implementation of \ref BCTypeToString().
747 * @details Full API contract (arguments, ownership, side effects) is documented with
748 * the header declaration in `include/logging.h`.
749 * @see BCTypeToString()
750 */
751const char* BCTypeToString(BCType type) {
752 switch (type) {
753 // case DIRICHLET: return "DIRICHLET";
754 // case NEUMANN: return "NEUMANN";
755 case WALL: return "WALL";
756 case INLET: return "INLET";
757 case OUTLET: return "OUTLET";
758 case FARFIELD: return "FARFIELD";
759 case PERIODIC: return "PERIODIC";
760 case INTERFACE: return "INTERFACE";
761
762 // case CUSTOM: return "CUSTOM";
763 default: return "Unknown BC Type";
764 }
765}
766
767/**
768 * @brief Internal helper implementation: `BCHandlerTypeToString()`.
769 * @details Local to this translation unit.
770 */
771const char* BCHandlerTypeToString(BCHandlerType handler_type) {
772 switch (handler_type) {
773 // Wall & Symmetry Handlers
774 case BC_HANDLER_WALL_NOSLIP: return "noslip";
775 case BC_HANDLER_WALL_MOVING: return "moving";
776 case BC_HANDLER_SYMMETRY_PLANE: return "symmetry_plane";
777
778 // Inlet Handlers
779 case BC_HANDLER_INLET_CONSTANT_VELOCITY: return "constant_velocity";
780 case BC_HANDLER_INLET_PULSATILE_FLUX: return "pulsatile_flux";
781 case BC_HANDLER_INLET_PARABOLIC: return "parabolic";
782 case BC_HANDLER_INLET_PROFILE_FROM_FILE: return "prescribed_flow";
783
784 // Outlet Handlers
785 case BC_HANDLER_OUTLET_CONSERVATION: return "conservation";
786 case BC_HANDLER_OUTLET_PRESSURE: return "pressure";
787
788 // Other Physical Handlers
789 case BC_HANDLER_FARFIELD_NONREFLECTING: return "nonreflecting";
790
791 // Multi-Block / Interface Handlers
792 case BC_HANDLER_PERIODIC_GEOMETRIC: return "geometric";
793 case BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX: return "constant flux";
794 case BC_HANDLER_PERIODIC_DRIVEN_INITIAL_FLUX: return "initial flux";
795 case BC_HANDLER_INTERFACE_OVERSET: return "overset";
796
797 // Default case
798 case BC_HANDLER_UNDEFINED:
799 default: return "UNKNOWN_HANDLER";
800 }
801}
802
803/**
804 * @brief Implementation of \ref DualMonitorDestroy().
805 * @details Full API contract (arguments, ownership, side effects) is documented with
806 * the header declaration in `include/logging.h`.
807 * @see DualMonitorDestroy()
808 */
809PetscErrorCode DualMonitorDestroy(void **ctx)
810{
811 DualMonitorCtx *monctx = (DualMonitorCtx*)*ctx;
812 PetscErrorCode ierr;
813 PetscMPIInt rank;
814
815 PetscFunctionBeginUser;
816 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank); CHKERRQ(ierr);
817 if(!rank && monctx->file_handle){
818 fclose(monctx->file_handle);
819 }
820
821 ierr = PetscFree(monctx); CHKERRQ(ierr);
822 *ctx = NULL;
823 PetscFunctionReturn(0);
824}
825
826/**
827 * @brief A custom KSP monitor that logs the true residual to a file and optionally to the console.
828 *
829 * This function replicates the behavior of KSPMonitorTrueResidualNorm by calculating
830 * the true residual norm ||b - Ax|| itself. It unconditionally logs to a file
831 * viewer and conditionally logs to the console based on a flag in the context.
832 *
833 * @param ksp The Krylov subspace context.
834 * @param it The current iteration number.
835 * @param rnorm The preconditioned residual norm (ignored, we compute our own).
836 * @param ctx A pointer to the DualMonitorCtx structure.
837 * @return PetscErrorCode 0 on success.
838 */
839#undef __FUNCT__
840#define __FUNCT__ "DualKSPMonitor"
841/**
842 * @brief Implementation of \ref DualKSPMonitor().
843 * @details Full API contract (arguments, ownership, side effects) is documented with
844 * the header declaration in `include/logging.h`.
845 * @see DualKSPMonitor()
846 */
847
848PetscErrorCode DualKSPMonitor(KSP ksp, PetscInt it, PetscReal rnorm, void *ctx)
849{
850 DualMonitorCtx *monctx = (DualMonitorCtx*)ctx;
851 PetscErrorCode ierr;
852 PetscReal trnorm, relnorm;
853 Vec r;
854 char norm_buf[256];
855 PetscMPIInt rank;
856
857 PetscFunctionBeginUser;
858 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank); CHKERRQ(ierr);
859
860 // 1. Calculate the true residual norm.
861 ierr = KSPBuildResidual(ksp, NULL, NULL, &r); CHKERRQ(ierr);
862 ierr = VecNorm(r, NORM_2, &trnorm); CHKERRQ(ierr);
863 ierr = VecDestroy(&r); CHKERRQ(ierr);
864
865 // 2. On the first iteration, compute and store the norm of the RHS vector `b`.
866 if (it == 0) {
867 Vec b;
868 ierr = KSPGetRhs(ksp, &b); CHKERRQ(ierr);
869 ierr = VecNorm(b, NORM_2, &monctx->bnorm); CHKERRQ(ierr);
870 }
871
872 if(!rank){
873 // 3. Compute the relative norm and format the output string.
874 if (monctx->bnorm > 1.e-15) {
875 relnorm = trnorm / monctx->bnorm;
876 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);
877 } else {
878 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);
879 }
880
881 // 4. Log to the file viewer (unconditionally).
882 if(monctx->file_handle){
883 ierr = PetscFPrintf(PETSC_COMM_SELF,monctx->file_handle,"%s\n", norm_buf); CHKERRQ(ierr);
884 }
885 // 5. Log to the console (conditionally).
886 if (monctx->log_to_console) {
887 PetscFPrintf(PETSC_COMM_SELF,stdout, "%s\n", norm_buf); CHKERRQ(ierr);
888 }
889
890 } //rank
891
892 PetscFunctionReturn(0);
893}
894
895#define SOLUTION_CONVERGENCE_FLUID_THRESHOLD 0.1
896#define SOLUTION_CONVERGENCE_REL_EPS 1.0e-30
897
909
914
915/**
916 * @brief Forms a guarded relative metric for solution-convergence logging.
917 *
918 * This helper centralizes the divide-by-nearly-zero protection used by the
919 * solution-convergence logger when turning an absolute drift into a relative
920 * one. The denominator is clamped away from zero so warmup rows, quiescent
921 * fields, and statistically small observables do not generate infinities.
922 *
923 * @param[in] numerator Absolute quantity or drift magnitude.
924 * @param[in] denominator Reference magnitude used for normalization.
925 * @return Guarded relative value `numerator / max(|denominator|, eps)`.
926 */
927static PetscReal SolutionConvergenceSafeRelative(PetscReal numerator, PetscReal denominator)
928{
929 return numerator / PetscMax(PetscAbsReal(denominator), SOLUTION_CONVERGENCE_REL_EPS);
930}
931
932/**
933 * @brief Computes instantaneous global flow observables for statistical mode.
934 *
935 * The statistical solution-convergence path does not compare full Eulerian
936 * fields. Instead, it tracks a compact history of global observables derived
937 * from the completed Eulerian state. This helper computes the current
938 * volume-weighted fluid-domain mean speed and mean kinetic energy from `Ucat`.
939 *
940 * Only physical fluid cells contribute:
941 * - solid/immersed cells are excluded using `Nvert`
942 * - cells with near-zero metric Jacobian are ignored to avoid invalid volume
943 * weights
944 *
945 * Each MPI rank accumulates local partial sums and the routine reduces them to
946 * one global pair of observables.
947 *
948 * @param[in] simCtx Simulation context owning the finest-level flow
949 * fields.
950 * @param[out] mean_speed_out Volume-weighted domain mean of `|u|`.
951 * @param[out] mean_ke_out Volume-weighted domain mean of `0.5 |u|^2`.
952 * @return PetscErrorCode 0 on success.
953 */
954static PetscErrorCode ComputeCurrentFlowObservables(SimCtx *simCtx, PetscReal *mean_speed_out, PetscReal *mean_ke_out)
955{
956 PetscReal local[3] = {0.0, 0.0, 0.0};
957 PetscReal global[3] = {0.0, 0.0, 0.0};
958 UserCtx *user = NULL;
959
960 PetscFunctionBeginUser;
961 if (!simCtx || !mean_speed_out || !mean_ke_out) {
962 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "ComputeCurrentFlowObservables received a NULL argument.");
963 }
964
965 user = simCtx->usermg.mgctx[simCtx->usermg.mglevels - 1].user;
966
967 for (PetscInt bi = 0; bi < simCtx->block_number; ++bi) {
968 const DMDALocalInfo info = user[bi].info;
969 Cmpnts ***ucat = NULL;
970 PetscReal ***aj = NULL;
971 PetscReal ***nvert = NULL;
972
973 PetscCall(DMDAVecGetArrayRead(user[bi].fda, user[bi].Ucat, &ucat));
974 PetscCall(DMDAVecGetArrayRead(user[bi].da, user[bi].Aj, &aj));
975 PetscCall(DMDAVecGetArrayRead(user[bi].da, user[bi].Nvert, &nvert));
976
977 for (PetscInt k = info.zs; k < info.zs + info.zm; ++k) {
978 for (PetscInt j = info.ys; j < info.ys + info.ym; ++j) {
979 for (PetscInt i = info.xs; i < info.xs + info.xm; ++i) {
980 PetscReal jac = aj[k][j][i];
981 PetscReal cell_volume = 0.0;
982 PetscReal speed = 0.0;
983 PetscReal ke = 0.0;
984
985 if (nvert[k][j][i] > SOLUTION_CONVERGENCE_FLUID_THRESHOLD) continue;
986 if (PetscAbsReal(jac) <= 1.0e-14) continue;
987
988 cell_volume = 1.0 / jac;
989 speed = PetscSqrtReal(ucat[k][j][i].x * ucat[k][j][i].x +
990 ucat[k][j][i].y * ucat[k][j][i].y +
991 ucat[k][j][i].z * ucat[k][j][i].z);
992 ke = 0.5 * speed * speed;
993
994 local[0] += cell_volume;
995 local[1] += speed * cell_volume;
996 local[2] += ke * cell_volume;
997 }
998 }
999 }
1000
1001 PetscCall(DMDAVecRestoreArrayRead(user[bi].da, user[bi].Nvert, &nvert));
1002 PetscCall(DMDAVecRestoreArrayRead(user[bi].da, user[bi].Aj, &aj));
1003 PetscCall(DMDAVecRestoreArrayRead(user[bi].fda, user[bi].Ucat, &ucat));
1004 }
1005
1006 PetscCallMPI(MPI_Allreduce(local, global, 3, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD));
1007
1008 if (global[0] <= 0.0) {
1009 *mean_speed_out = 0.0;
1010 *mean_ke_out = 0.0;
1011 } else {
1012 *mean_speed_out = global[1] / global[0];
1013 *mean_ke_out = global[2] / global[0];
1014 }
1015
1016 PetscFunctionReturn(0);
1017}
1018
1019/**
1020 * @brief Computes deterministic solution-drift metrics for the current step.
1021 *
1022 * This helper powers both `steady_deterministic` and
1023 * `periodic_deterministic` solution-convergence modes. It compares the
1024 * completed Eulerian state against either:
1025 * - the previous physical timestep (`Ucat_o`, `P_o`) for steady/transient use
1026 * - the stored phase-aligned snapshot ring for periodic use
1027 *
1028 * The routine performs two passes over the fluid cells:
1029 * 1. velocity/observable pass
1030 * - computes current and reference mean speed / mean KE
1031 * - computes absolute/relative velocity L2 drift
1032 * - accumulates pressure means needed for gauge removal
1033 * 2. pressure-only pass
1034 * - subtracts the volume-weighted mean pressure from both states
1035 * - computes gauge-invariant pressure L2 drift
1036 *
1037 * Warmup behavior is handled here. If no valid reference exists yet, all drift
1038 * outputs are left at zero and `has_reference_out` is set to `PETSC_FALSE`,
1039 * while current observables are still reported.
1040 *
1041 * @param[in] simCtx Simulation context owning the current state.
1042 * @param[in] periodic_mode `PETSC_TRUE` when comparing against
1043 * phase-aligned periodic storage.
1044 * @param[in] phase_step Active phase slot for periodic mode, or `-1`
1045 * when unused.
1046 * @param[in] samples_before Number of solution-convergence samples
1047 * already recorded before this timestep.
1048 * @param[out] has_reference_out Whether a valid comparison state existed.
1049 * @param[out] u_abs_l2_out Absolute L2 drift of Cartesian velocity.
1050 * @param[out] u_rel_l2_out Relative L2 drift of Cartesian velocity.
1051 * @param[out] p_abs_l2_out Gauge-invariant absolute L2 pressure drift.
1052 * @param[out] p_rel_l2_out Gauge-invariant relative L2 pressure drift.
1053 * @param[out] mean_speed_out Current volume-weighted mean speed.
1054 * @param[out] mean_speed_ref_out Reference volume-weighted mean speed.
1055 * @param[out] mean_speed_abs_out Absolute drift of mean speed.
1056 * @param[out] mean_speed_rel_out Relative drift of mean speed.
1057 * @param[out] mean_ke_out Current volume-weighted mean kinetic energy.
1058 * @param[out] mean_ke_ref_out Reference volume-weighted mean kinetic
1059 * energy.
1060 * @param[out] mean_ke_abs_out Absolute drift of mean kinetic energy.
1061 * @param[out] mean_ke_rel_out Relative drift of mean kinetic energy.
1062 * @return PetscErrorCode 0 on success.
1063 */
1064static PetscErrorCode ComputeDeterministicSolutionMetrics(SimCtx *simCtx,
1065 PetscBool periodic_mode,
1066 PetscInt phase_step,
1067 PetscInt samples_before,
1068 PetscBool *has_reference_out,
1069 PetscReal *u_abs_l2_out,
1070 PetscReal *u_rel_l2_out,
1071 PetscReal *p_abs_l2_out,
1072 PetscReal *p_rel_l2_out,
1073 PetscReal *mean_speed_out,
1074 PetscReal *mean_speed_ref_out,
1075 PetscReal *mean_speed_abs_out,
1076 PetscReal *mean_speed_rel_out,
1077 PetscReal *mean_ke_out,
1078 PetscReal *mean_ke_ref_out,
1079 PetscReal *mean_ke_abs_out,
1080 PetscReal *mean_ke_rel_out)
1081{
1082 SolutionConvergenceDeterministicPass1 local_pass1 = {0};
1083 SolutionConvergenceDeterministicPass1 global_pass1 = {0};
1084 SolutionConvergenceDeterministicPass2 local_pass2 = {0};
1085 SolutionConvergenceDeterministicPass2 global_pass2 = {0};
1086 PetscReal current_pressure_mean = 0.0;
1087 PetscReal reference_pressure_mean = 0.0;
1088 UserCtx *user = NULL;
1089
1090 PetscFunctionBeginUser;
1091 if (!simCtx || !has_reference_out || !u_abs_l2_out || !u_rel_l2_out || !p_abs_l2_out || !p_rel_l2_out ||
1092 !mean_speed_out || !mean_speed_ref_out || !mean_speed_abs_out || !mean_speed_rel_out ||
1093 !mean_ke_out || !mean_ke_ref_out || !mean_ke_abs_out || !mean_ke_rel_out) {
1094 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "ComputeDeterministicSolutionMetrics received a NULL output pointer.");
1095 }
1096
1097 *has_reference_out = periodic_mode
1098 ? (PetscBool)(simCtx->solutionConvergencePeriodSteps > 0 &&
1099 phase_step >= 0 &&
1100 phase_step < simCtx->solutionConvergencePeriodSteps &&
1101 samples_before >= simCtx->solutionConvergencePeriodSteps)
1102 : (PetscBool)(samples_before > 0);
1103
1104 *u_abs_l2_out = 0.0;
1105 *u_rel_l2_out = 0.0;
1106 *p_abs_l2_out = 0.0;
1107 *p_rel_l2_out = 0.0;
1108 *mean_speed_out = 0.0;
1109 *mean_speed_ref_out = 0.0;
1110 *mean_speed_abs_out = 0.0;
1111 *mean_speed_rel_out = 0.0;
1112 *mean_ke_out = 0.0;
1113 *mean_ke_ref_out = 0.0;
1114 *mean_ke_abs_out = 0.0;
1115 *mean_ke_rel_out = 0.0;
1116
1117 user = simCtx->usermg.mgctx[simCtx->usermg.mglevels - 1].user;
1118
1119 for (PetscInt bi = 0; bi < simCtx->block_number; ++bi) {
1120 const DMDALocalInfo info = user[bi].info;
1121 Cmpnts ***ucat = NULL;
1122 Cmpnts ***ucat_ref = NULL;
1123 PetscReal ***pressure = NULL;
1124 PetscReal ***pressure_ref = NULL;
1125 PetscReal ***aj = NULL;
1126 PetscReal ***nvert = NULL;
1127 Vec ucat_reference_vec = NULL;
1128 Vec pressure_reference_vec = NULL;
1129
1130 if (*has_reference_out) {
1131 if (periodic_mode) {
1132 ucat_reference_vec = user[bi].solutionConvergencePeriodicUcatRef[phase_step];
1133 pressure_reference_vec = user[bi].solutionConvergencePeriodicPRef[phase_step];
1134 } else {
1135 ucat_reference_vec = user[bi].Ucat_o;
1136 pressure_reference_vec = user[bi].P_o;
1137 }
1138 }
1139
1140 PetscCall(DMDAVecGetArrayRead(user[bi].fda, user[bi].Ucat, &ucat));
1141 PetscCall(DMDAVecGetArrayRead(user[bi].da, user[bi].P, &pressure));
1142 PetscCall(DMDAVecGetArrayRead(user[bi].da, user[bi].Aj, &aj));
1143 PetscCall(DMDAVecGetArrayRead(user[bi].da, user[bi].Nvert, &nvert));
1144 if (*has_reference_out) {
1145 PetscCall(DMDAVecGetArrayRead(user[bi].fda, ucat_reference_vec, &ucat_ref));
1146 PetscCall(DMDAVecGetArrayRead(user[bi].da, pressure_reference_vec, &pressure_ref));
1147 }
1148
1149 for (PetscInt k = info.zs; k < info.zs + info.zm; ++k) {
1150 for (PetscInt j = info.ys; j < info.ys + info.ym; ++j) {
1151 for (PetscInt i = info.xs; i < info.xs + info.xm; ++i) {
1152 PetscReal jac = aj[k][j][i];
1153 PetscReal cell_volume = 0.0;
1154 PetscReal speed = 0.0;
1155 PetscReal ke = 0.0;
1156
1157 if (nvert[k][j][i] > SOLUTION_CONVERGENCE_FLUID_THRESHOLD) continue;
1158 if (PetscAbsReal(jac) <= 1.0e-14) continue;
1159
1160 cell_volume = 1.0 / jac;
1161 speed = PetscSqrtReal(ucat[k][j][i].x * ucat[k][j][i].x +
1162 ucat[k][j][i].y * ucat[k][j][i].y +
1163 ucat[k][j][i].z * ucat[k][j][i].z);
1164 ke = 0.5 * speed * speed;
1165
1166 local_pass1.fluid_volume += cell_volume;
1167 local_pass1.current_speed_sum += speed * cell_volume;
1168 local_pass1.current_ke_sum += ke * cell_volume;
1169 local_pass1.current_u_norm_sq += (ucat[k][j][i].x * ucat[k][j][i].x +
1170 ucat[k][j][i].y * ucat[k][j][i].y +
1171 ucat[k][j][i].z * ucat[k][j][i].z) * cell_volume;
1172
1173 if (*has_reference_out) {
1174 PetscReal ref_speed = PetscSqrtReal(ucat_ref[k][j][i].x * ucat_ref[k][j][i].x +
1175 ucat_ref[k][j][i].y * ucat_ref[k][j][i].y +
1176 ucat_ref[k][j][i].z * ucat_ref[k][j][i].z);
1177 PetscReal ref_ke = 0.5 * ref_speed * ref_speed;
1178 PetscReal dux = ucat[k][j][i].x - ucat_ref[k][j][i].x;
1179 PetscReal duy = ucat[k][j][i].y - ucat_ref[k][j][i].y;
1180 PetscReal duz = ucat[k][j][i].z - ucat_ref[k][j][i].z;
1181
1182 local_pass1.reference_speed_sum += ref_speed * cell_volume;
1183 local_pass1.reference_ke_sum += ref_ke * cell_volume;
1184 local_pass1.delta_u_norm_sq += (dux * dux + duy * duy + duz * duz) * cell_volume;
1185 local_pass1.current_pressure_sum += pressure[k][j][i] * cell_volume;
1186 local_pass1.reference_pressure_sum += pressure_ref[k][j][i] * cell_volume;
1187 }
1188 }
1189 }
1190 }
1191
1192 if (*has_reference_out) {
1193 PetscCall(DMDAVecRestoreArrayRead(user[bi].da, pressure_reference_vec, &pressure_ref));
1194 PetscCall(DMDAVecRestoreArrayRead(user[bi].fda, ucat_reference_vec, &ucat_ref));
1195 }
1196 PetscCall(DMDAVecRestoreArrayRead(user[bi].da, user[bi].Nvert, &nvert));
1197 PetscCall(DMDAVecRestoreArrayRead(user[bi].da, user[bi].Aj, &aj));
1198 PetscCall(DMDAVecRestoreArrayRead(user[bi].da, user[bi].P, &pressure));
1199 PetscCall(DMDAVecRestoreArrayRead(user[bi].fda, user[bi].Ucat, &ucat));
1200 }
1201
1202 PetscCallMPI(MPI_Allreduce(&local_pass1, &global_pass1,
1203 sizeof(SolutionConvergenceDeterministicPass1) / sizeof(PetscReal),
1204 MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD));
1205
1206 if (global_pass1.fluid_volume <= 0.0) PetscFunctionReturn(0);
1207
1208 *mean_speed_out = global_pass1.current_speed_sum / global_pass1.fluid_volume;
1209 *mean_ke_out = global_pass1.current_ke_sum / global_pass1.fluid_volume;
1210
1211 if (!*has_reference_out) PetscFunctionReturn(0);
1212
1213 *mean_speed_ref_out = global_pass1.reference_speed_sum / global_pass1.fluid_volume;
1214 *mean_speed_abs_out = PetscAbsReal(*mean_speed_out - *mean_speed_ref_out);
1215 *mean_speed_rel_out = SolutionConvergenceSafeRelative(*mean_speed_abs_out, *mean_speed_out);
1216 *mean_ke_ref_out = global_pass1.reference_ke_sum / global_pass1.fluid_volume;
1217 *mean_ke_abs_out = PetscAbsReal(*mean_ke_out - *mean_ke_ref_out);
1218 *mean_ke_rel_out = SolutionConvergenceSafeRelative(*mean_ke_abs_out, *mean_ke_out);
1219
1220 current_pressure_mean = global_pass1.current_pressure_sum / global_pass1.fluid_volume;
1221 reference_pressure_mean = global_pass1.reference_pressure_sum / global_pass1.fluid_volume;
1222
1223 for (PetscInt bi = 0; bi < simCtx->block_number; ++bi) {
1224 const DMDALocalInfo info = user[bi].info;
1225 PetscReal ***pressure = NULL;
1226 PetscReal ***pressure_ref = NULL;
1227 PetscReal ***aj = NULL;
1228 PetscReal ***nvert = NULL;
1229 Vec pressure_reference_vec = periodic_mode ? user[bi].solutionConvergencePeriodicPRef[phase_step] : user[bi].P_o;
1230
1231 PetscCall(DMDAVecGetArrayRead(user[bi].da, user[bi].P, &pressure));
1232 PetscCall(DMDAVecGetArrayRead(user[bi].da, pressure_reference_vec, &pressure_ref));
1233 PetscCall(DMDAVecGetArrayRead(user[bi].da, user[bi].Aj, &aj));
1234 PetscCall(DMDAVecGetArrayRead(user[bi].da, user[bi].Nvert, &nvert));
1235
1236 for (PetscInt k = info.zs; k < info.zs + info.zm; ++k) {
1237 for (PetscInt j = info.ys; j < info.ys + info.ym; ++j) {
1238 for (PetscInt i = info.xs; i < info.xs + info.xm; ++i) {
1239 PetscReal jac = aj[k][j][i];
1240 PetscReal cell_volume = 0.0;
1241 PetscReal current_pressure = 0.0;
1242 PetscReal reference_pressure = 0.0;
1243 PetscReal delta_pressure = 0.0;
1244
1245 if (nvert[k][j][i] > SOLUTION_CONVERGENCE_FLUID_THRESHOLD) continue;
1246 if (PetscAbsReal(jac) <= 1.0e-14) continue;
1247
1248 cell_volume = 1.0 / jac;
1249 current_pressure = pressure[k][j][i] - current_pressure_mean;
1250 reference_pressure = pressure_ref[k][j][i] - reference_pressure_mean;
1251 delta_pressure = current_pressure - reference_pressure;
1252
1253 local_pass2.current_pressure_norm_sq += current_pressure * current_pressure * cell_volume;
1254 local_pass2.delta_pressure_norm_sq += delta_pressure * delta_pressure * cell_volume;
1255 }
1256 }
1257 }
1258
1259 PetscCall(DMDAVecRestoreArrayRead(user[bi].da, user[bi].Nvert, &nvert));
1260 PetscCall(DMDAVecRestoreArrayRead(user[bi].da, user[bi].Aj, &aj));
1261 PetscCall(DMDAVecRestoreArrayRead(user[bi].da, pressure_reference_vec, &pressure_ref));
1262 PetscCall(DMDAVecRestoreArrayRead(user[bi].da, user[bi].P, &pressure));
1263 }
1264
1265 PetscCallMPI(MPI_Allreduce(&local_pass2, &global_pass2,
1266 sizeof(SolutionConvergenceDeterministicPass2) / sizeof(PetscReal),
1267 MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD));
1268
1269 *u_abs_l2_out = PetscSqrtReal(global_pass1.delta_u_norm_sq);
1270 *u_rel_l2_out = SolutionConvergenceSafeRelative(*u_abs_l2_out, PetscSqrtReal(global_pass1.current_u_norm_sq));
1271 *p_abs_l2_out = PetscSqrtReal(global_pass2.delta_pressure_norm_sq);
1272 *p_rel_l2_out = SolutionConvergenceSafeRelative(*p_abs_l2_out, PetscSqrtReal(global_pass2.current_pressure_norm_sq));
1273
1274 PetscFunctionReturn(0);
1275}
1276
1277/**
1278 * @brief Reads one sample from the statistical ring buffer by age.
1279 *
1280 * The statistical logger stores scalar observables in a compact circular
1281 * buffer. This helper interprets the buffer using `samples_available` as the
1282 * logical end of the history and returns the entry `offset_from_latest` steps
1283 * back from the newest stored sample.
1284 *
1285 * Out-of-range requests return zero so warmup handling can remain simple and
1286 * deterministic.
1287 *
1288 * @param[in] history Ring-buffer storage array.
1289 * @param[in] capacity Total ring-buffer capacity.
1290 * @param[in] samples_available Number of logical samples available to read.
1291 * @param[in] offset_from_latest `0` means newest sample, `1` previous sample,
1292 * and so on.
1293 * @return Requested historical sample, or `0.0` if unavailable.
1294 */
1295static PetscReal SolutionConvergenceHistoryGet(const PetscReal *history,
1296 PetscInt capacity,
1297 PetscInt samples_available,
1298 PetscInt offset_from_latest)
1299{
1300 PetscInt count = 0;
1301 PetscInt index = 0;
1302
1303 if (!history || capacity <= 0 || samples_available <= 0 || offset_from_latest < 0) {
1304 return 0.0;
1305 }
1306
1307 count = PetscMin(samples_available, capacity);
1308 if (offset_from_latest >= count) return 0.0;
1309
1310 index = (samples_available - 1 - offset_from_latest) % capacity;
1311 if (index < 0) index += capacity;
1312 return history[index];
1313}
1314
1315/**
1316 * @brief Appends one timestep's scalar observables to the statistical history.
1317 *
1318 * Statistical solution-convergence compares adjacent windows of scalar
1319 * observables rather than full fields. This helper writes the current
1320 * `mean_speed` and `mean_ke` into the rolling history arrays using the current
1321 * sample count to choose the circular-buffer slot.
1322 *
1323 * @param[in,out] simCtx Simulation context owning the history arrays.
1324 * @param[in] samples_before Number of samples present before appending the
1325 * current timestep.
1326 * @param[in] mean_speed Current timestep mean-speed observable.
1327 * @param[in] mean_ke Current timestep mean-KE observable.
1328 * @return PetscErrorCode 0 on success.
1329 */
1330static PetscErrorCode AppendStatisticalObservableSample(SimCtx *simCtx,
1331 PetscInt samples_before,
1332 PetscReal mean_speed,
1333 PetscReal mean_ke)
1334{
1335 PetscInt history_capacity = 0;
1336 PetscInt slot = 0;
1337
1338 PetscFunctionBeginUser;
1339 if (!simCtx || !simCtx->solutionConvergenceMeanSpeedHistory || !simCtx->solutionConvergenceMeanKEHistory) {
1340 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Statistical solution-convergence history is not allocated.");
1341 }
1342
1343 history_capacity = 2 * simCtx->solutionConvergenceWindowSteps;
1344 if (history_capacity <= 0) {
1345 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Statistical solution-convergence history capacity must be positive.");
1346 }
1347
1348 slot = samples_before % history_capacity;
1349 simCtx->solutionConvergenceMeanSpeedHistory[slot] = mean_speed;
1350 simCtx->solutionConvergenceMeanKEHistory[slot] = mean_ke;
1351
1352 PetscFunctionReturn(0);
1353}
1354
1355/**
1356 * @brief Computes adjacent-window drift metrics for statistical steady mode.
1357 *
1358 * Once enough samples have been accumulated, this helper forms:
1359 * - the current window over the most recent `window_steps` samples
1360 * - the previous adjacent window over the preceding `window_steps` samples
1361 *
1362 * From those windows it computes means, RMS values, and absolute/relative
1363 * drift for both tracked observables (`mean_speed` and `mean_ke`). When the
1364 * history is still warming up:
1365 * - fewer than `window_steps` samples: no window metrics are available
1366 * - between `window_steps` and `2*window_steps - 1` samples: current-window
1367 * metrics are available, but no reference window exists yet
1368 *
1369 * @param[in] simCtx Simulation context owning the
1370 * statistical history.
1371 * @param[in] samples_available Number of samples available after
1372 * appending the current timestep.
1373 * @param[out] has_reference_out Whether both adjacent windows
1374 * exist and drift metrics are
1375 * meaningful.
1376 * @param[out] mean_speed_window_out Mean speed over the current
1377 * window.
1378 * @param[out] mean_speed_window_prev_out Mean speed over the previous
1379 * window.
1380 * @param[out] mean_speed_window_abs_out Absolute drift between current
1381 * and previous window means.
1382 * @param[out] mean_speed_window_rel_out Relative drift between current
1383 * and previous window means.
1384 * @param[out] mean_speed_rms_window_out RMS of mean-speed samples in the
1385 * current window.
1386 * @param[out] mean_speed_rms_window_prev_out RMS of mean-speed samples in the
1387 * previous window.
1388 * @param[out] mean_speed_rms_window_abs_out Absolute drift between window RMS
1389 * values.
1390 * @param[out] mean_speed_rms_window_rel_out Relative drift between window RMS
1391 * values.
1392 * @param[out] mean_ke_window_out Mean kinetic energy over the
1393 * current window.
1394 * @param[out] mean_ke_window_prev_out Mean kinetic energy over the
1395 * previous window.
1396 * @param[out] mean_ke_window_abs_out Absolute drift between current
1397 * and previous KE-window means.
1398 * @param[out] mean_ke_window_rel_out Relative drift between current
1399 * and previous KE-window means.
1400 * @param[out] mean_ke_rms_window_out RMS of mean-KE samples in the
1401 * current window.
1402 * @param[out] mean_ke_rms_window_prev_out RMS of mean-KE samples in the
1403 * previous window.
1404 * @param[out] mean_ke_rms_window_abs_out Absolute drift between KE-window
1405 * RMS values.
1406 * @param[out] mean_ke_rms_window_rel_out Relative drift between KE-window
1407 * RMS values.
1408 * @return PetscErrorCode 0 on success.
1409 */
1410static PetscErrorCode ComputeStatisticalWindowMetrics(const SimCtx *simCtx,
1411 PetscInt samples_available,
1412 PetscBool *has_reference_out,
1413 PetscReal *mean_speed_window_out,
1414 PetscReal *mean_speed_window_prev_out,
1415 PetscReal *mean_speed_window_abs_out,
1416 PetscReal *mean_speed_window_rel_out,
1417 PetscReal *mean_speed_rms_window_out,
1418 PetscReal *mean_speed_rms_window_prev_out,
1419 PetscReal *mean_speed_rms_window_abs_out,
1420 PetscReal *mean_speed_rms_window_rel_out,
1421 PetscReal *mean_ke_window_out,
1422 PetscReal *mean_ke_window_prev_out,
1423 PetscReal *mean_ke_window_abs_out,
1424 PetscReal *mean_ke_window_rel_out,
1425 PetscReal *mean_ke_rms_window_out,
1426 PetscReal *mean_ke_rms_window_prev_out,
1427 PetscReal *mean_ke_rms_window_abs_out,
1428 PetscReal *mean_ke_rms_window_rel_out)
1429{
1430 PetscInt w = 0;
1431 PetscInt history_capacity = 0;
1432 PetscReal speed_sum = 0.0;
1433 PetscReal speed_sum_sq = 0.0;
1434 PetscReal speed_prev_sum = 0.0;
1435 PetscReal speed_prev_sum_sq = 0.0;
1436 PetscReal ke_sum = 0.0;
1437 PetscReal ke_sum_sq = 0.0;
1438 PetscReal ke_prev_sum = 0.0;
1439 PetscReal ke_prev_sum_sq = 0.0;
1440
1441 PetscFunctionBeginUser;
1442 if (!simCtx || !has_reference_out || !mean_speed_window_out || !mean_speed_window_prev_out ||
1443 !mean_speed_window_abs_out || !mean_speed_window_rel_out || !mean_speed_rms_window_out ||
1444 !mean_speed_rms_window_prev_out || !mean_speed_rms_window_abs_out || !mean_speed_rms_window_rel_out ||
1445 !mean_ke_window_out || !mean_ke_window_prev_out || !mean_ke_window_abs_out || !mean_ke_window_rel_out ||
1446 !mean_ke_rms_window_out || !mean_ke_rms_window_prev_out || !mean_ke_rms_window_abs_out ||
1447 !mean_ke_rms_window_rel_out) {
1448 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "ComputeStatisticalWindowMetrics received a NULL output pointer.");
1449 }
1450
1451 *has_reference_out = PETSC_FALSE;
1452 *mean_speed_window_out = 0.0;
1453 *mean_speed_window_prev_out = 0.0;
1454 *mean_speed_window_abs_out = 0.0;
1455 *mean_speed_window_rel_out = 0.0;
1456 *mean_speed_rms_window_out = 0.0;
1457 *mean_speed_rms_window_prev_out = 0.0;
1458 *mean_speed_rms_window_abs_out = 0.0;
1459 *mean_speed_rms_window_rel_out = 0.0;
1460 *mean_ke_window_out = 0.0;
1461 *mean_ke_window_prev_out = 0.0;
1462 *mean_ke_window_abs_out = 0.0;
1463 *mean_ke_window_rel_out = 0.0;
1464 *mean_ke_rms_window_out = 0.0;
1465 *mean_ke_rms_window_prev_out = 0.0;
1466 *mean_ke_rms_window_abs_out = 0.0;
1467 *mean_ke_rms_window_rel_out = 0.0;
1468
1469 w = simCtx->solutionConvergenceWindowSteps;
1470 history_capacity = 2 * w;
1471 if (w <= 0 || samples_available < w) PetscFunctionReturn(0);
1472
1473 for (PetscInt idx = 0; idx < w; ++idx) {
1474 PetscReal speed_value = SolutionConvergenceHistoryGet(simCtx->solutionConvergenceMeanSpeedHistory,
1475 history_capacity,
1476 samples_available,
1477 idx);
1478 PetscReal ke_value = SolutionConvergenceHistoryGet(simCtx->solutionConvergenceMeanKEHistory,
1479 history_capacity,
1480 samples_available,
1481 idx);
1482 speed_sum += speed_value;
1483 speed_sum_sq += speed_value * speed_value;
1484 ke_sum += ke_value;
1485 ke_sum_sq += ke_value * ke_value;
1486 }
1487
1488 *mean_speed_window_out = speed_sum / (PetscReal)w;
1489 *mean_speed_rms_window_out = PetscSqrtReal(PetscMax(0.0, speed_sum_sq / (PetscReal)w -
1490 (*mean_speed_window_out) * (*mean_speed_window_out)));
1491 *mean_ke_window_out = ke_sum / (PetscReal)w;
1492 *mean_ke_rms_window_out = PetscSqrtReal(PetscMax(0.0, ke_sum_sq / (PetscReal)w -
1493 (*mean_ke_window_out) * (*mean_ke_window_out)));
1494
1495 if (samples_available < 2 * w) PetscFunctionReturn(0);
1496
1497 for (PetscInt idx = w; idx < 2 * w; ++idx) {
1498 PetscReal speed_value = SolutionConvergenceHistoryGet(simCtx->solutionConvergenceMeanSpeedHistory,
1499 history_capacity,
1500 samples_available,
1501 idx);
1502 PetscReal ke_value = SolutionConvergenceHistoryGet(simCtx->solutionConvergenceMeanKEHistory,
1503 history_capacity,
1504 samples_available,
1505 idx);
1506 speed_prev_sum += speed_value;
1507 speed_prev_sum_sq += speed_value * speed_value;
1508 ke_prev_sum += ke_value;
1509 ke_prev_sum_sq += ke_value * ke_value;
1510 }
1511
1512 *has_reference_out = PETSC_TRUE;
1513 *mean_speed_window_prev_out = speed_prev_sum / (PetscReal)w;
1514 *mean_speed_window_abs_out = PetscAbsReal(*mean_speed_window_out - *mean_speed_window_prev_out);
1515 *mean_speed_window_rel_out = SolutionConvergenceSafeRelative(*mean_speed_window_abs_out, *mean_speed_window_out);
1516 *mean_speed_rms_window_prev_out = PetscSqrtReal(PetscMax(0.0, speed_prev_sum_sq / (PetscReal)w -
1517 (*mean_speed_window_prev_out) * (*mean_speed_window_prev_out)));
1518 *mean_speed_rms_window_abs_out = PetscAbsReal(*mean_speed_rms_window_out - *mean_speed_rms_window_prev_out);
1519 *mean_speed_rms_window_rel_out = SolutionConvergenceSafeRelative(*mean_speed_rms_window_abs_out, *mean_speed_rms_window_out);
1520
1521 *mean_ke_window_prev_out = ke_prev_sum / (PetscReal)w;
1522 *mean_ke_window_abs_out = PetscAbsReal(*mean_ke_window_out - *mean_ke_window_prev_out);
1523 *mean_ke_window_rel_out = SolutionConvergenceSafeRelative(*mean_ke_window_abs_out, *mean_ke_window_out);
1524 *mean_ke_rms_window_prev_out = PetscSqrtReal(PetscMax(0.0, ke_prev_sum_sq / (PetscReal)w -
1525 (*mean_ke_window_prev_out) * (*mean_ke_window_prev_out)));
1526 *mean_ke_rms_window_abs_out = PetscAbsReal(*mean_ke_rms_window_out - *mean_ke_rms_window_prev_out);
1527 *mean_ke_rms_window_rel_out = SolutionConvergenceSafeRelative(*mean_ke_rms_window_abs_out, *mean_ke_rms_window_out);
1528
1529 PetscFunctionReturn(0);
1530}
1531
1532/**
1533 * @brief Maps the internal solution-convergence mode enum to its log label.
1534 *
1535 * The logger writes a human-readable mode string into the
1536 * `solution_convergence.log` banner and `mode` column. This helper keeps the
1537 * formatting centralized so the file output stays consistent with the accepted
1538 * configuration names.
1539 *
1540 * @param[in] mode Internal solution-convergence mode selector.
1541 * @return Lowercase string label written to the log output.
1542 */
1543static const char *SolutionConvergenceModeToString(SolutionConvergenceMode mode)
1544{
1545 switch (mode) {
1546 case SOLUTION_CONVERGENCE_STEADY_DETERMINISTIC: return "steady_deterministic";
1547 case SOLUTION_CONVERGENCE_PERIODIC_DETERMINISTIC: return "periodic_deterministic";
1548 case SOLUTION_CONVERGENCE_STATISTICAL_STEADY: return "statistical_steady";
1549 case SOLUTION_CONVERGENCE_TRANSIENT: return "transient";
1550 default: return "unknown";
1551 }
1552}
1553
1554/**
1555 * @brief Implementation of \ref LOG_SOLUTION_CONVERGENCE().
1556 * @details Full API contract (arguments, ownership, side effects) is documented with
1557 * the header declaration in `include/logging.h`.
1558 * @see LOG_SOLUTION_CONVERGENCE()
1559 */
1560PetscErrorCode LOG_SOLUTION_CONVERGENCE(SimCtx *simCtx)
1561{
1562 PetscMPIInt rank = 0;
1563 PetscBool has_reference = PETSC_FALSE;
1564 PetscInt phase_step = -1;
1565 PetscInt samples_before = 0;
1566 PetscReal u_abs_l2 = 0.0, u_rel_l2 = 0.0, p_abs_l2 = 0.0, p_rel_l2 = 0.0;
1567 PetscReal mean_speed = 0.0, mean_speed_reference = 0.0, mean_speed_abs_drift = 0.0, mean_speed_rel_drift = 0.0;
1568 PetscReal mean_ke = 0.0, mean_ke_reference = 0.0, mean_ke_abs_drift = 0.0, mean_ke_rel_drift = 0.0;
1569 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;
1570 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;
1571 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;
1572 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;
1573
1574 PetscFunctionBeginUser;
1575 if (!simCtx) PetscFunctionReturn(0);
1576 if (simCtx->exec_mode != EXEC_MODE_SOLVER) PetscFunctionReturn(0);
1577
1578 samples_before = simCtx->solutionConvergenceSamplesRecorded;
1579
1580 switch (simCtx->solutionConvergenceMode) {
1581 case SOLUTION_CONVERGENCE_STEADY_DETERMINISTIC:
1582 case SOLUTION_CONVERGENCE_TRANSIENT:
1583 PetscCall(ComputeDeterministicSolutionMetrics(simCtx, PETSC_FALSE, -1, samples_before,
1584 &has_reference,
1585 &u_abs_l2, &u_rel_l2,
1586 &p_abs_l2, &p_rel_l2,
1587 &mean_speed, &mean_speed_reference,
1588 &mean_speed_abs_drift, &mean_speed_rel_drift,
1589 &mean_ke, &mean_ke_reference,
1590 &mean_ke_abs_drift, &mean_ke_rel_drift));
1591 break;
1592 case SOLUTION_CONVERGENCE_PERIODIC_DETERMINISTIC:
1593 phase_step = simCtx->solutionConvergencePeriodSteps > 0 ? (simCtx->step % simCtx->solutionConvergencePeriodSteps) : -1;
1594 PetscCall(ComputeDeterministicSolutionMetrics(simCtx, PETSC_TRUE, phase_step, samples_before,
1595 &has_reference,
1596 &u_abs_l2, &u_rel_l2,
1597 &p_abs_l2, &p_rel_l2,
1598 &mean_speed, &mean_speed_reference,
1599 &mean_speed_abs_drift, &mean_speed_rel_drift,
1600 &mean_ke, &mean_ke_reference,
1601 &mean_ke_abs_drift, &mean_ke_rel_drift));
1602 break;
1603 case SOLUTION_CONVERGENCE_STATISTICAL_STEADY:
1604 PetscCall(ComputeCurrentFlowObservables(simCtx, &mean_speed, &mean_ke));
1605 PetscCall(AppendStatisticalObservableSample(simCtx, samples_before, mean_speed, mean_ke));
1606 PetscCall(ComputeStatisticalWindowMetrics(simCtx, samples_before + 1,
1607 &has_reference,
1608 &mean_speed_window, &mean_speed_window_prev,
1609 &mean_speed_window_abs_drift, &mean_speed_window_rel_drift,
1610 &mean_speed_rms_window, &mean_speed_rms_window_prev,
1611 &mean_speed_rms_window_abs_drift, &mean_speed_rms_window_rel_drift,
1612 &mean_ke_window, &mean_ke_window_prev,
1613 &mean_ke_window_abs_drift, &mean_ke_window_rel_drift,
1614 &mean_ke_rms_window, &mean_ke_rms_window_prev,
1615 &mean_ke_rms_window_abs_drift, &mean_ke_rms_window_rel_drift));
1616 break;
1617 default:
1618 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown solution convergence mode %d.", (int)simCtx->solutionConvergenceMode);
1619 }
1620
1621 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
1622 if (rank == 0) {
1623 char log_path[PETSC_MAX_PATH_LEN + 32];
1624 FILE *f = NULL;
1625 const char *mode_str = SolutionConvergenceModeToString(simCtx->solutionConvergenceMode);
1626
1627 PetscCall(PetscSNPrintf(log_path, sizeof(log_path), "%s/solution_convergence.log", simCtx->log_dir));
1628 f = fopen(log_path, "a");
1629 if (!f) {
1630 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Cannot open solution convergence log: %s", log_path);
1631 }
1632
1633 if (ftell(f) == 0) {
1634 switch (simCtx->solutionConvergenceMode) {
1635 case SOLUTION_CONVERGENCE_STEADY_DETERMINISTIC:
1636 case SOLUTION_CONVERGENCE_TRANSIENT:
1637 fprintf(f, "==================== Solution Convergence Log [mode: %s] ====================\n", mode_str);
1638 /* 16 columns; header width = 314 chars */
1639 fprintf(f, "%-10s | %-18s | %-22s | %-3s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s\n",
1640 "step", "time", "mode", "ref",
1641 "u_abs_l2", "u_rel_l2", "p_abs_l2", "p_rel_l2",
1642 "mean_speed", "spd_ref", "spd_abs", "spd_rel",
1643 "mean_ke", "ke_ref", "ke_abs", "ke_rel");
1644 fprintf(f, "----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n");
1645 break;
1646 case SOLUTION_CONVERGENCE_PERIODIC_DETERMINISTIC:
1647 fprintf(f, "==================== Solution Convergence Log [mode: %s | period_steps: %d] ====================\n",
1648 mode_str, (int)simCtx->solutionConvergencePeriodSteps);
1649 /* 18 columns; header width = 330 chars */
1650 fprintf(f, "%-10s | %-18s | %-22s | %-3s | %-5s | %-5s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s\n",
1651 "step", "time", "mode", "ref", "ph", "per",
1652 "u_abs_l2", "u_rel_l2", "p_abs_l2", "p_rel_l2",
1653 "mean_speed", "spd_ref", "spd_abs", "spd_rel",
1654 "mean_ke", "ke_ref", "ke_abs", "ke_rel");
1655 fprintf(f, "----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n");
1656 break;
1657 case SOLUTION_CONVERGENCE_STATISTICAL_STEADY:
1658 fprintf(f, "==================== Solution Convergence Log [mode: %s | window_steps: %d] ====================\n",
1659 mode_str, (int)simCtx->solutionConvergenceWindowSteps);
1660 /* 21 columns; header width = 406 chars */
1661 fprintf(f, "%-10s | %-18s | %-22s | %-3s | %-5s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s | %-18s\n",
1662 "step", "time", "mode", "ref", "win",
1663 "mean_speed", "mean_ke",
1664 "spd_win", "spd_win_prev", "spd_win_abs", "spd_win_rel",
1665 "spd_rms_win", "spd_rms_abs", "spd_rms_rel",
1666 "ke_win", "ke_win_prev", "ke_win_abs", "ke_win_rel",
1667 "ke_rms_win", "ke_rms_abs", "ke_rms_rel");
1668 fprintf(f, "------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n");
1669 break;
1670 default: break;
1671 }
1672 }
1673
1674 switch (simCtx->solutionConvergenceMode) {
1675 case SOLUTION_CONVERGENCE_STEADY_DETERMINISTIC:
1676 case SOLUTION_CONVERGENCE_TRANSIENT:
1677 fprintf(f,
1678 "%-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",
1679 (int)simCtx->step, (double)simCtx->ti, mode_str, has_reference ? 1 : 0,
1680 (double)u_abs_l2, (double)u_rel_l2, (double)p_abs_l2, (double)p_rel_l2,
1681 (double)mean_speed, (double)mean_speed_reference,
1682 (double)mean_speed_abs_drift, (double)mean_speed_rel_drift,
1683 (double)mean_ke, (double)mean_ke_reference,
1684 (double)mean_ke_abs_drift, (double)mean_ke_rel_drift);
1685 break;
1686 case SOLUTION_CONVERGENCE_PERIODIC_DETERMINISTIC:
1687 fprintf(f,
1688 "%-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",
1689 (int)simCtx->step, (double)simCtx->ti, mode_str, has_reference ? 1 : 0,
1690 (int)phase_step, (int)simCtx->solutionConvergencePeriodSteps,
1691 (double)u_abs_l2, (double)u_rel_l2, (double)p_abs_l2, (double)p_rel_l2,
1692 (double)mean_speed, (double)mean_speed_reference,
1693 (double)mean_speed_abs_drift, (double)mean_speed_rel_drift,
1694 (double)mean_ke, (double)mean_ke_reference,
1695 (double)mean_ke_abs_drift, (double)mean_ke_rel_drift);
1696 break;
1697 case SOLUTION_CONVERGENCE_STATISTICAL_STEADY:
1698 fprintf(f,
1699 "%-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",
1700 (int)simCtx->step, (double)simCtx->ti, mode_str, has_reference ? 1 : 0,
1701 (int)simCtx->solutionConvergenceWindowSteps,
1702 (double)mean_speed, (double)mean_ke,
1703 (double)mean_speed_window, (double)mean_speed_window_prev,
1704 (double)mean_speed_window_abs_drift, (double)mean_speed_window_rel_drift,
1705 (double)mean_speed_rms_window,
1706 (double)mean_speed_rms_window_abs_drift, (double)mean_speed_rms_window_rel_drift,
1707 (double)mean_ke_window, (double)mean_ke_window_prev,
1708 (double)mean_ke_window_abs_drift, (double)mean_ke_window_rel_drift,
1709 (double)mean_ke_rms_window,
1710 (double)mean_ke_rms_window_abs_drift, (double)mean_ke_rms_window_rel_drift);
1711 break;
1712 default: break;
1713 }
1714 fclose(f);
1715 }
1716
1717 if (simCtx->solutionConvergenceMode == SOLUTION_CONVERGENCE_PERIODIC_DETERMINISTIC &&
1718 phase_step >= 0 && phase_step < simCtx->solutionConvergencePeriodSteps) {
1719 UserCtx *user = simCtx->usermg.mgctx[simCtx->usermg.mglevels - 1].user;
1720 for (PetscInt bi = 0; bi < simCtx->block_number; ++bi) {
1721 PetscCall(VecCopy(user[bi].Ucat, user[bi].solutionConvergencePeriodicUcatRef[phase_step]));
1722 PetscCall(VecCopy(user[bi].P, user[bi].solutionConvergencePeriodicPRef[phase_step]));
1723 }
1724 }
1725
1726 simCtx->solutionConvergenceSamplesRecorded = samples_before + 1;
1727
1728 PetscFunctionReturn(0);
1729}
1730
1731/**
1732 * @brief Logs continuity metrics for a single block to a file.
1733 *
1734 * This function should be called for each block, once per timestep. It opens a
1735 * central log file in append mode. To ensure the header is written only once,
1736 * it checks if it is processing block 0 on the simulation's start step.
1737 *
1738 * @param user A pointer to the UserCtx for the specific block whose metrics
1739 * are to be logged. The function accesses both global (SimCtx)
1740 * and local (user->...) data.
1741 * @return PetscErrorCode 0 on success.
1742 */
1743#undef __FUNCT__
1744#define __FUNCT__ "LOG_CONTINUITY_METRICS"
1745/**
1746 * @brief Implementation of \ref LOG_CONTINUITY_METRICS().
1747 * @details Full API contract (arguments, ownership, side effects) is documented with
1748 * the header declaration in `include/logging.h`.
1749 * @see LOG_CONTINUITY_METRICS()
1750 */
1751
1752PetscErrorCode LOG_CONTINUITY_METRICS(UserCtx *user)
1753{
1754 PetscErrorCode ierr;
1755 PetscMPIInt rank;
1756 SimCtx *simCtx = user->simCtx; // Get the shared SimCtx
1757 const PetscInt bi = user->_this; // Get this block's specific ID
1758 const PetscInt ti = simCtx->step; // Get the current timestep
1759
1760 PetscFunctionBeginUser;
1761 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
1762
1763 // Only rank 0 performs file I/O.
1764 if (!rank) {
1765 FILE *f;
1766 char filen[PETSC_MAX_PATH_LEN + 64];
1767 ierr = PetscSNPrintf(filen, sizeof(filen), "%s/Continuity_Metrics.log", simCtx->log_dir); CHKERRQ(ierr);
1768
1769 // Open the log file in append mode.
1770 f = fopen(filen, "a");
1771 if (!f) {
1772 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Cannot open log file: %s", filen);
1773 }
1774
1775 // Write a header only when the file is empty and it's the first block (bi=0).
1776 // Using ftell() instead of step comparison ensures correctness across continuations.
1777 if (ftell(f) == 0 && bi == 0) {
1778 PetscFPrintf(PETSC_COMM_SELF, f, "%-10s | %-6s | %-18s | %-30s | %-18s | %-18s | %-18s | %-18s\n",
1779 "Timestep", "Block", "Max Divergence", "Max Divergence Location ([k][j][i]=idx)", "Sum(RHS)","Total Flux In", "Total Flux Out", "Net Flux");
1780 PetscFPrintf(PETSC_COMM_SELF, f, "------------------------------------------------------------------------------------------------------------------------------------------\n");
1781 }
1782
1783 // Prepare the data strings and values for the current block.
1784 PetscReal net_flux = simCtx->FluxInSum - simCtx->FluxOutSum;
1785 char location_str[64];
1786 sprintf(location_str, "([%d][%d][%d] = %d)", (int)simCtx->MaxDivz, (int)simCtx->MaxDivy, (int)simCtx->MaxDivx, (int)simCtx->MaxDivFlatArg);
1787
1788 // Write the formatted line for the current block.
1789 PetscFPrintf(PETSC_COMM_SELF, f, "%-10d | %-6d | %-18.10e | %-39s | %-18.10e | %-18.10e | %-18.10e | %-18.10e\n",
1790 (int)ti,
1791 (int)bi,
1792 (double)simCtx->MaxDiv,
1793 location_str,
1794 (double)simCtx->summationRHS,
1795 (double)simCtx->FluxInSum,
1796 (double)simCtx->FluxOutSum,
1797 (double)net_flux);
1798
1799 fclose(f);
1800 }
1801
1802 PetscFunctionReturn(0);
1803}
1804
1805/**
1806 * @brief Implementation of \ref ParticleLocationStatusToString().
1807 * @details Full API contract (arguments, ownership, side effects) is documented with
1808 * the header declaration in `include/logging.h`.
1809 * @see ParticleLocationStatusToString()
1810 */
1811const char* ParticleLocationStatusToString(ParticleLocationStatus level)
1812{
1813 switch (level) {
1814 case NEEDS_LOCATION: return "NEEDS_LOCATION";
1815 case ACTIVE_AND_LOCATED: return "ACTIVE_AND_LOCATED";
1816 case MIGRATING_OUT: return "MIGRATING_OUT";
1817 case LOST: return "LOST";
1818 case UNINITIALIZED: return "UNINITIALIZED";
1819 default: return "UNKNOWN_LEVEL";
1820 }
1821}
1822
1823///////// Profiling System /////////
1824
1825// Data structure to hold profiling info for one function
1826typedef struct {
1827 const char *name;
1828 double total_time;
1829 double current_step_time;
1830 long long total_call_count;
1831 long long current_step_call_count;
1832 double start_time; // Timer for the current call
1833 PetscBool always_log;
1834} ProfiledFunction;
1835
1836// Global registry for all profiled functions
1837static ProfiledFunction *g_profiler_registry = NULL;
1838static PetscInt g_profiler_count = 0;
1839static PetscInt g_profiler_capacity = 0;
1840
1841// Internal helper to find a function in the registry or create it
1842/**
1843 * @brief Internal helper implementation: `_FindOrCreateEntry()`.
1844 * @details Local to this translation unit.
1845 */
1846static PetscErrorCode _FindOrCreateEntry(const char *func_name, PetscInt *idx)
1847{
1848 PetscFunctionBeginUser;
1849 // Search for existing entry
1850 for (PetscInt i = 0; i < g_profiler_count; ++i) {
1851 if (strcmp(g_profiler_registry[i].name, func_name) == 0) {
1852 *idx = i;
1853 PetscFunctionReturn(0);
1854 }
1855 }
1856
1857 // Not found, create a new entry
1858 if (g_profiler_count >= g_profiler_capacity) {
1859 PetscInt new_capacity = g_profiler_capacity == 0 ? 16 : g_profiler_capacity * 2;
1860 PetscErrorCode ierr = PetscRealloc(sizeof(ProfiledFunction) * new_capacity, &g_profiler_registry); CHKERRQ(ierr);
1861 g_profiler_capacity = new_capacity;
1862 }
1863
1864 *idx = g_profiler_count;
1865 g_profiler_registry[*idx].name = func_name;
1866 g_profiler_registry[*idx].total_time = 0.0;
1867 g_profiler_registry[*idx].current_step_time = 0.0;
1868 g_profiler_registry[*idx].total_call_count = 0;
1869 g_profiler_registry[*idx].current_step_call_count = 0;
1870 g_profiler_registry[*idx].start_time = 0.0;
1871 g_profiler_registry[*idx].always_log = PETSC_FALSE;
1872 g_profiler_count++;
1873
1874 PetscFunctionReturn(0);
1875}
1876
1877// --- Public API Implementation ---
1878/**
1879 * @brief Internal helper implementation: `ProfilingInitialize()`.
1880 * @details Local to this translation unit.
1881 */
1882PetscErrorCode ProfilingInitialize(SimCtx *simCtx)
1883{
1884 PetscFunctionBeginUser;
1885 if (!simCtx) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "SimCtx cannot be null for ProfilingInitialize");
1886
1887 // Iterate through the list of critical functions provided in SimCtx
1888 for (PetscInt i = 0; i < simCtx->nProfilingSelectedFuncs; ++i) {
1889 PetscInt idx;
1890 const char *func_name = simCtx->profilingSelectedFuncs[i];
1891 PetscErrorCode ierr = _FindOrCreateEntry(func_name, &idx); CHKERRQ(ierr);
1892 g_profiler_registry[idx].always_log = PETSC_TRUE;
1893
1894 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Marked '%s' as a critical function for profiling.\n", func_name);
1895 }
1896 PetscFunctionReturn(0);
1897}
1898
1899/**
1900 * @brief Implementation of \ref _ProfilingStart().
1901 * @details Full API contract (arguments, ownership, side effects) is documented with
1902 * the header declaration in `include/logging.h`.
1903 * @see _ProfilingStart()
1904 */
1905
1906void _ProfilingStart(const char *func_name)
1907{
1908 PetscInt idx;
1909 if (_FindOrCreateEntry(func_name, &idx) != 0) return; // Fail silently
1910 PetscTime(&g_profiler_registry[idx].start_time);
1911}
1912
1913/**
1914 * @brief Implementation of \ref _ProfilingEnd().
1915 * @details Full API contract (arguments, ownership, side effects) is documented with
1916 * the header declaration in `include/logging.h`.
1917 * @see _ProfilingEnd()
1918 */
1919
1920void _ProfilingEnd(const char *func_name)
1921{
1922 double end_time;
1923 PetscTime(&end_time);
1924
1925 PetscInt idx;
1926 if (_FindOrCreateEntry(func_name, &idx) != 0) return; // Fail silently
1927
1928 double elapsed = end_time - g_profiler_registry[idx].start_time;
1929 g_profiler_registry[idx].total_time += elapsed;
1930 g_profiler_registry[idx].current_step_time += elapsed;
1931 g_profiler_registry[idx].total_call_count++;
1932 g_profiler_registry[idx].current_step_call_count++;
1933}
1934
1935/**
1936 * @brief Implementation of \ref ProfilingResetTimestepCounters().
1937 * @details Full API contract (arguments, ownership, side effects) is documented with
1938 * the header declaration in `include/logging.h`.
1939 * @see ProfilingResetTimestepCounters()
1940 */
1941
1942PetscErrorCode ProfilingResetTimestepCounters(void)
1943{
1944 PetscFunctionBeginUser;
1945 for (PetscInt i = 0; i < g_profiler_count; ++i) {
1946 g_profiler_registry[i].current_step_time = 0.0;
1947 g_profiler_registry[i].current_step_call_count = 0;
1948 }
1949 PetscFunctionReturn(0);
1950}
1951
1952/**
1953 * @brief Implementation of \ref ProfilingLogTimestepSummary().
1954 * @details Full API contract (arguments, ownership, side effects) is documented with
1955 * the header declaration in `include/logging.h`.
1956 * @see ProfilingLogTimestepSummary()
1957 */
1958
1959PetscErrorCode ProfilingLogTimestepSummary(SimCtx *simCtx, PetscInt step)
1960{
1961 PetscBool should_write = PETSC_FALSE;
1962 FILE *f = NULL;
1963 char filen[(2 * PETSC_MAX_PATH_LEN) + 16];
1964
1965 PetscFunctionBeginUser;
1966 if (!simCtx) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "SimCtx cannot be null for ProfilingLogTimestepSummary");
1967
1968 if (strcmp(simCtx->profilingTimestepMode, "off") == 0) {
1969 for (PetscInt i = 0; i < g_profiler_count; ++i) {
1970 g_profiler_registry[i].current_step_time = 0.0;
1971 g_profiler_registry[i].current_step_call_count = 0;
1972 }
1973 PetscFunctionReturn(0);
1974 }
1975
1976 for (PetscInt i = 0; i < g_profiler_count; ++i) {
1977 if (g_profiler_registry[i].current_step_call_count <= 0) {
1978 continue;
1979 }
1980 if (strcmp(simCtx->profilingTimestepMode, "all") == 0 || g_profiler_registry[i].always_log) {
1981 should_write = PETSC_TRUE;
1982 break;
1983 }
1984 }
1985
1986 if (should_write && simCtx->rank == 0) {
1987 snprintf(filen, sizeof(filen), "%s/%s", simCtx->log_dir, simCtx->profilingTimestepFile);
1988 if (step == simCtx->StartStep + 1 && !simCtx->continueMode) {
1989 f = fopen(filen, "w");
1990 if (!f) {
1991 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Cannot open profiling timestep log file: %s", filen);
1992 }
1993 PetscFPrintf(PETSC_COMM_SELF, f, "step,function,calls,step_time_s\n");
1994 } else {
1995 f = fopen(filen, "a");
1996 if (!f) {
1997 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Cannot open profiling timestep log file: %s", filen);
1998 }
1999 if (step == simCtx->StartStep + 1 && ftell(f) == 0) {
2000 PetscFPrintf(PETSC_COMM_SELF, f, "step,function,calls,step_time_s\n");
2001 }
2002 }
2003
2004 for (PetscInt i = 0; i < g_profiler_count; ++i) {
2005 if (g_profiler_registry[i].current_step_call_count <= 0) {
2006 continue;
2007 }
2008 if (strcmp(simCtx->profilingTimestepMode, "all") == 0 || g_profiler_registry[i].always_log) {
2009 PetscFPrintf(
2010 PETSC_COMM_SELF,
2011 f,
2012 "%d,%s,%lld,%.6f\n",
2013 (int)step,
2014 g_profiler_registry[i].name,
2015 g_profiler_registry[i].current_step_call_count,
2016 g_profiler_registry[i].current_step_time
2017 );
2018 }
2019 }
2020 fclose(f);
2021 }
2022
2023 // Reset per-step counters for the next iteration
2024 for (PetscInt i = 0; i < g_profiler_count; ++i) {
2025 g_profiler_registry[i].current_step_time = 0.0;
2026 g_profiler_registry[i].current_step_call_count = 0;
2027 }
2028 PetscFunctionReturn(0);
2029}
2030
2031/**
2032 * @brief Implementation of \ref RuntimeMemoryLogSample().
2033 * @details Full API contract (arguments, ownership, side effects) is documented with
2034 * the header declaration in `include/logging.h`.
2035 * @see RuntimeMemoryLogSample()
2036 */
2037PetscErrorCode RuntimeMemoryLogSample(SimCtx *simCtx, PetscInt step, const char *event, const char *reason)
2038{
2039 PetscErrorCode ierr;
2040 PetscLogDouble process_current_bytes = 0.0;
2041 PetscLogDouble process_peak_bytes = 0.0;
2042 PetscLogDouble petsc_current_bytes = 0.0;
2043 PetscLogDouble petsc_peak_bytes = 0.0;
2044 PetscReal local_values[5];
2045 PetscReal global_values[5];
2046 PetscReal process_current_mb = 0.0;
2047 PetscReal process_peak_mb = 0.0;
2048 PetscReal petsc_current_mb = 0.0;
2049 PetscReal petsc_peak_mb = 0.0;
2050 PetscReal process_change_mb = 0.0;
2051 char path[(2 * PETSC_MAX_PATH_LEN) + 16];
2052 FILE *f = NULL;
2053
2054 PetscFunctionBeginUser;
2055 if (!simCtx) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "SimCtx cannot be null for RuntimeMemoryLogSample");
2056 if (!simCtx->runtimeMemoryLogEnabled) PetscFunctionReturn(0);
2057
2058 ierr = PetscMemoryGetCurrentUsage(&process_current_bytes); CHKERRQ(ierr);
2059 ierr = PetscMemoryGetMaximumUsage(&process_peak_bytes); CHKERRQ(ierr);
2060 ierr = PetscMallocGetCurrentUsage(&petsc_current_bytes); CHKERRQ(ierr);
2061 ierr = PetscMallocGetMaximumUsage(&petsc_peak_bytes); CHKERRQ(ierr);
2062
2063 process_current_mb = (PetscReal)(process_current_bytes / (1024.0 * 1024.0));
2064 process_peak_mb = (PetscReal)(process_peak_bytes / (1024.0 * 1024.0));
2065 petsc_current_mb = (PetscReal)(petsc_current_bytes / (1024.0 * 1024.0));
2066 petsc_peak_mb = (PetscReal)(petsc_peak_bytes / (1024.0 * 1024.0));
2067 if (simCtx->runtimeMemoryLogHasPrevious) {
2068 process_change_mb = process_current_mb - simCtx->runtimeMemoryLogPreviousProcessMB;
2069 }
2070
2071 local_values[0] = process_current_mb;
2072 local_values[1] = process_peak_mb;
2073 local_values[2] = petsc_current_mb;
2074 local_values[3] = petsc_peak_mb;
2075 local_values[4] = process_change_mb;
2076 ierr = MPI_Allreduce(local_values, global_values, 5, MPIU_REAL, MPI_MAX, PETSC_COMM_WORLD); CHKERRMPI(ierr);
2077
2078 simCtx->runtimeMemoryLogPreviousProcessMB = process_current_mb;
2079 simCtx->runtimeMemoryLogHasPrevious = PETSC_TRUE;
2080
2081 if (simCtx->rank == 0) {
2082 ierr = PetscSNPrintf(path, sizeof(path), "%s/%s", simCtx->log_dir, simCtx->runtimeMemoryLogFile); CHKERRQ(ierr);
2083 f = fopen(path, "a");
2084 if (!f) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Cannot open runtime memory log file: %s", path);
2085
2086 if (!simCtx->runtimeMemoryLogStarted) {
2087 fprintf(f, "# PICurv runtime memory log\n");
2088 if (simCtx->continueMode) {
2089 fprintf(f, "# Continuation from step %" PetscInt_FMT "\n", simCtx->StartStep);
2090 }
2091 fprintf(
2092 f,
2093 "%-8s %-10s %22s %20s %22s %28s %22s %-18s\n",
2094 "Step",
2095 "Event",
2096 "Process Current MB Max",
2097 "Process Peak MB Max",
2098 "PETSc Allocated MB Max",
2099 "PETSc Peak Allocated MB Max",
2100 "Process Change MB Max",
2101 "Reason"
2102 );
2103 simCtx->runtimeMemoryLogStarted = PETSC_TRUE;
2104 }
2105
2106 fprintf(
2107 f,
2108 "%-8" PetscInt_FMT " %-10s %22.3f %20.3f %22.3f %28.3f %22.3f %-18s\n",
2109 step,
2110 event ? event : "-",
2111 (double)global_values[0],
2112 (double)global_values[1],
2113 (double)global_values[2],
2114 (double)global_values[3],
2115 (double)global_values[4],
2116 (reason && reason[0]) ? reason : "-"
2117 );
2118 if ((event && (strcmp(event, "Shutdown") == 0 || strcmp(event, "Final") == 0))) {
2119 fflush(f);
2120 }
2121 fclose(f);
2122 }
2123
2124 PetscFunctionReturn(0);
2125}
2126
2127// Comparison function for qsort to sort by total_time in descending order
2128/**
2129 * @brief Internal helper implementation: `_CompareProfiledFunctions()`.
2130 * @details Local to this translation unit.
2131 */
2132static int _CompareProfiledFunctions(const void *a, const void *b)
2133{
2134 const ProfiledFunction *func_a = (const ProfiledFunction *)a;
2135 const ProfiledFunction *func_b = (const ProfiledFunction *)b;
2136
2137 if (func_a->total_time < func_b->total_time) return 1;
2138 if (func_a->total_time > func_b->total_time) return -1;
2139 return 0;
2140}
2141
2142/**
2143 * @brief Implementation of \ref ProfilingFinalize().
2144 * @details Full API contract (arguments, ownership, side effects) is documented with
2145 * the header declaration in `include/logging.h`.
2146 * @see ProfilingFinalize()
2147 */
2148PetscErrorCode ProfilingFinalize(SimCtx *simCtx)
2149{
2150 PetscErrorCode ierr;
2151 PetscInt rank = simCtx->rank;
2152 PetscFunctionBeginUser;
2153 if (!simCtx->profilingFinalSummary) PetscFunctionReturn(0);
2154 if (!rank) {
2155
2156 char exec_mode_modifier[32] = "Unknown";
2157 if(simCtx->exec_mode == EXEC_MODE_SOLVER) PetscCall(PetscStrncpy(exec_mode_modifier, "Solver", sizeof(exec_mode_modifier)));
2158 else if(simCtx->exec_mode == EXEC_MODE_POSTPROCESSOR) PetscCall(PetscStrncpy(exec_mode_modifier, "PostProcessor", sizeof(exec_mode_modifier)));
2159 //--- Step 0: Create a file viewer for log file
2160 FILE *f;
2161 char filen[PETSC_MAX_PATH_LEN + 128];
2162 ierr = PetscSNPrintf(filen, sizeof(filen), "%s/ProfilingSummary_%s.log",simCtx->log_dir,exec_mode_modifier); CHKERRQ(ierr);
2163
2164 // Open the log file: append with section label in continue mode, truncate otherwise.
2165 if (simCtx->continueMode) {
2166 f = fopen(filen, "a");
2167 if (!f) {
2168 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Cannot open log file: %s", filen);
2169 }
2170 fprintf(f, "\n=== Continuation from step %" PetscInt_FMT " ===\n", simCtx->StartStep);
2171 } else {
2172 f = fopen(filen, "w");
2173 if (!f) {
2174 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Cannot open log file: %s", filen);
2175 }
2176 }
2177
2178 // --- Step 1: Sort the data for readability ---
2179 qsort(g_profiler_registry, g_profiler_count, sizeof(ProfiledFunction), _CompareProfiledFunctions);
2180
2181 // --- Step 2: Dynamically determine the width for the function name column ---
2182 PetscInt max_name_len = strlen("Function"); // Start with the header's length
2183 for (PetscInt i = 0; i < g_profiler_count; ++i) {
2184 if (g_profiler_registry[i].total_call_count > 0) {
2185 PetscInt len = strlen(g_profiler_registry[i].name);
2186 if (len > max_name_len) {
2187 max_name_len = len;
2188 }
2189 }
2190 }
2191 // Add a little padding
2192 max_name_len += 2;
2193
2194 // --- Step 3: Define fixed widths for numeric columns for consistent alignment ---
2195 const int time_width = 18;
2196 const int count_width = 15;
2197 const int avg_width = 22;
2198
2199 // --- Step 4: Print the formatted table ---
2200 PetscFPrintf(PETSC_COMM_SELF, f, "=================================================================================================================\n");
2201 PetscFPrintf(PETSC_COMM_SELF, f, " FINAL PROFILING SUMMARY (Sorted by Total Time)\n");
2202 PetscFPrintf(PETSC_COMM_SELF, f, "=================================================================================================================\n");
2203
2204 // Header Row
2205 PetscFPrintf(PETSC_COMM_SELF, f, "%-*s | %-*s | %-*s | %-*s\n",
2206 max_name_len, "Function",
2207 time_width, "Total Time (s)",
2208 count_width, "Call Count",
2209 avg_width, "Avg. Time/Call (ms)");
2210
2211 // Separator Line (dynamically sized)
2212 for (int i = 0; i < max_name_len; i++) PetscFPrintf(PETSC_COMM_SELF, f, "-");
2213 PetscFPrintf(PETSC_COMM_SELF, f, "-|-");
2214 for (int i = 0; i < time_width; i++) PetscFPrintf(PETSC_COMM_SELF, f, "-");
2215 PetscFPrintf(PETSC_COMM_SELF, f, "-|-");
2216 for (int i = 0; i < count_width; i++) PetscFPrintf(PETSC_COMM_SELF, f, "-");
2217 PetscFPrintf(PETSC_COMM_SELF, f, "-|-");
2218 for (int i = 0; i < avg_width; i++) PetscFPrintf(PETSC_COMM_SELF, f, "-");
2219 PetscFPrintf(PETSC_COMM_SELF, f, "\n");
2220
2221 // Data Rows
2222 for (PetscInt i = 0; i < g_profiler_count; ++i) {
2223 if (g_profiler_registry[i].total_call_count > 0) {
2224 double avg_time_ms = (g_profiler_registry[i].total_time / g_profiler_registry[i].total_call_count) * 1000.0;
2225 PetscFPrintf(PETSC_COMM_SELF, f, "%-*s | %*.*f | %*lld | %*.*f\n",
2226 max_name_len, g_profiler_registry[i].name,
2227 time_width, 6, g_profiler_registry[i].total_time,
2228 count_width, g_profiler_registry[i].total_call_count,
2229 avg_width, 6, avg_time_ms);
2230 PetscFPrintf(PETSC_COMM_SELF, f, "------------------------------------------------------------------------------------------------------------------\n");
2231 }
2232 }
2233 PetscFPrintf(PETSC_COMM_SELF, f, "==================================================================================================================\n");
2234
2235 fclose(f);
2236 }
2237
2238 // --- Final Cleanup ---
2239 PetscFree(g_profiler_registry);
2240 g_profiler_registry = NULL;
2241 g_profiler_count = 0;
2242 g_profiler_capacity = 0;
2243 PetscFunctionReturn(0);
2244}
2245
2246/*================================================================================*
2247 * PROGRESS BAR UTILITY *
2248 *================================================================================*/
2249
2250/**
2251 * @brief Internal helper implementation: `PrintProgressBar()`.
2252 * @details Local to this translation unit.
2253 */
2254void PrintProgressBar(PetscInt step, PetscInt startStep, PetscInt totalSteps, PetscReal currentTime)
2255{
2256 if (totalSteps <= 0) return;
2257
2258 // --- Configuration ---
2259 const int barWidth = 50;
2260
2261 // --- Calculation ---
2262 // Calculate progress as a fraction from 0.0 to 1.0
2263 PetscReal progress = (PetscReal)(step - startStep + 1) / totalSteps;
2264 // Ensure progress doesn't exceed 1.0 due to floating point inaccuracies
2265 if (progress > 1.0) progress = 1.0;
2266
2267 int pos = (int)(barWidth * progress);
2268
2269 // --- Printing ---
2270 // Carriage return moves cursor to the beginning of the line
2271 PetscPrintf(PETSC_COMM_SELF, "\rProgress: [");
2272
2273 for (int i = 0; i < barWidth; ++i) {
2274 if (i < pos) {
2275 PetscPrintf(PETSC_COMM_SELF, "=");
2276 } else if (i == pos) {
2277 PetscPrintf(PETSC_COMM_SELF, ">");
2278 } else {
2279 PetscPrintf(PETSC_COMM_SELF, " ");
2280 }
2281 }
2282
2283 // Print percentage, step count, and current time
2284 PetscPrintf(PETSC_COMM_SELF, "] %3d%% (Step %" PetscInt_FMT "/%" PetscInt_FMT ", t=%.4f)",
2285 (int)(progress * 100.0),
2286 step + 1,
2287 startStep + totalSteps,
2288 currentTime);
2289
2290 // Flush the output buffer to ensure the bar is displayed immediately
2291 fflush(stdout);
2292}
2293
2294#undef __FUNCT__
2295#define __FUNCT__ "LOG_FIELD_MIN_MAX"
2296/**
2297 * @brief Implementation of \ref LOG_FIELD_MIN_MAX().
2298 * @details Full API contract is documented with the header declaration in `include/logging.h`.
2299 * @see LOG_FIELD_MIN_MAX()
2300 */
2301PetscErrorCode LOG_FIELD_MIN_MAX(UserCtx *user, const char *fieldName)
2302{
2303 PetscErrorCode ierr;
2304 PetscInt i, j, k;
2305 DMDALocalInfo info;
2306
2307 Vec fieldVec = NULL;
2308 DM dm = NULL;
2309 PetscInt dof;
2310 char data_layout[20];
2311
2312 PetscFunctionBeginUser;
2313
2314 // --- 1. Map string name to PETSc objects and determine data layout ---
2315 if (strcasecmp(fieldName, "Ucat") == 0) {
2316 fieldVec = user->Ucat; dm = user->fda; dof = 3; strcpy(data_layout, "Cell-Centered");
2317 } else if (strcasecmp(fieldName, "P") == 0) {
2318 fieldVec = user->P; dm = user->da; dof = 1; strcpy(data_layout, "Cell-Centered");
2319 } else if (strcasecmp(fieldName, "Diffusivity") == 0) {
2320 fieldVec = user->Diffusivity; dm = user->da; dof = 1; strcpy(data_layout, "Cell-Centered");
2321 } else if (strcasecmp(fieldName, "DiffusivityGradient") == 0) {
2322 fieldVec = user->DiffusivityGradient; dm = user->fda; dof = 3; strcpy(data_layout, "Cell-Centered");
2323 } else if (strcasecmp(fieldName, "Ucont") == 0) {
2324 fieldVec = user->lUcont; dm = user->fda; dof = 3; strcpy(data_layout, "Face-Centered");
2325 } else if (strcasecmp(fieldName, "Coordinates") == 0) {
2326 ierr = DMGetCoordinates(user->da, &fieldVec); CHKERRQ(ierr);
2327 dm = user->fda; dof = 3; strcpy(data_layout, "Node-Centered");
2328 } else if (strcasecmp(fieldName,"Psi") == 0) {
2329 fieldVec = user->Psi; dm = user->da; dof = 1; strcpy(data_layout, "Cell-Centered"); // Assuming Psi is cell-centered
2330 } else {
2331 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Field %s not recognized.", fieldName);
2332 }
2333
2334 if (!fieldVec) {
2335 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vector for field '%s' is NULL.", fieldName);
2336 }
2337 if (!dm) {
2338 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "DM for field '%s' is NULL.", fieldName);
2339 }
2340
2341 ierr = DMDAGetLocalInfo(dm, &info); CHKERRQ(ierr);
2342
2343 // --- 2. Define Architecture-Aware Loop Bounds ---
2344 PetscInt i_start, i_end, j_start, j_end, k_start, k_end;
2345
2346 if (strcmp(data_layout, "Cell-Centered") == 0) {
2347 // For cell-centered data, the physical values are stored from index 1 to N-1.
2348 // We find the intersection of the rank's owned range [xs, xe) with the
2349 // physical data range [1, IM-1).
2350 i_start = PetscMax(info.xs, 1); i_end = PetscMin(info.xs + info.xm, user->IM);
2351 j_start = PetscMax(info.ys, 1); j_end = PetscMin(info.ys + info.ym, user->JM);
2352 k_start = PetscMax(info.zs, 1); k_end = PetscMin(info.zs + info.zm, user->KM);
2353 } else { // For Node- or Face-Centered data
2354 // The physical values are stored from index 0 to N-1.
2355 // We find the intersection of the rank's owned range [xs, xe) with the
2356 // physical data range [0, IM-1].
2357 i_start = PetscMax(info.xs, 0); i_end = PetscMin(info.xs + info.xm, user->IM);
2358 j_start = PetscMax(info.ys, 0); j_end = PetscMin(info.ys + info.ym, user->JM);
2359 k_start = PetscMax(info.zs, 0); k_end = PetscMin(info.zs + info.zm, user->KM);
2360 }
2361
2362 // --- 3. Barrier for clean, grouped output ---
2363 ierr = MPI_Barrier(PETSC_COMM_WORLD); CHKERRQ(ierr);
2364 if (user->simCtx->rank == 0) {
2365 PetscPrintf(PETSC_COMM_SELF, "\n--- Field Ranges: [%s] (Layout: %s) ---\n", fieldName, data_layout);
2366 }
2367
2368 // --- 4. Branch on DoF and perform calculation with correct bounds ---
2369 if (dof == 1) {
2370 PetscReal localMin = PETSC_MAX_REAL, localMax = PETSC_MIN_REAL;
2371 PetscReal globalMin, globalMax;
2372 const PetscScalar ***array;
2373
2374 ierr = DMDAVecGetArrayRead(dm, fieldVec, &array); CHKERRQ(ierr);
2375 for (k = k_start; k < k_end; k++) {
2376 for (j = j_start; j < j_end; j++) {
2377 for (i = i_start; i < i_end; i++) {
2378 localMin = PetscMin(localMin, array[k][j][i]);
2379 localMax = PetscMax(localMax, array[k][j][i]);
2380 }
2381 }
2382 }
2383 ierr = DMDAVecRestoreArrayRead(dm, fieldVec, &array); CHKERRQ(ierr);
2384
2385 ierr = MPI_Allreduce(&localMin, &globalMin, 1, MPIU_REAL, MPI_MIN, PETSC_COMM_WORLD); CHKERRQ(ierr);
2386 ierr = MPI_Allreduce(&localMax, &globalMax, 1, MPIU_REAL, MPI_MAX, PETSC_COMM_WORLD); CHKERRQ(ierr);
2387
2388 PetscSynchronizedPrintf(PETSC_COMM_WORLD, " [Rank %d] Local Range: [ %11.4e , %11.4e ]\n", user->simCtx->rank, localMin, localMax);
2389 ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT); CHKERRQ(ierr);
2390 if (user->simCtx->rank == 0) {
2391 PetscPrintf(PETSC_COMM_SELF, " Global Range: [ %11.4e , %11.4e ]\n", globalMin, globalMax);
2392 }
2393
2394 } else if (dof == 3) {
2395 Cmpnts localMin = {PETSC_MAX_REAL, PETSC_MAX_REAL, PETSC_MAX_REAL};
2396 Cmpnts localMax = {PETSC_MIN_REAL, PETSC_MIN_REAL, PETSC_MIN_REAL};
2397 Cmpnts globalMin, globalMax;
2398 const Cmpnts ***array;
2399
2400 ierr = DMDAVecGetArrayRead(dm, fieldVec, &array); CHKERRQ(ierr);
2401 for (k = k_start; k < k_end; k++) {
2402 for (j = j_start; j < j_end; j++) {
2403 for (i = i_start; i < i_end; i++) {
2404 localMin.x = PetscMin(localMin.x, array[k][j][i].x);
2405 localMin.y = PetscMin(localMin.y, array[k][j][i].y);
2406 localMin.z = PetscMin(localMin.z, array[k][j][i].z);
2407 localMax.x = PetscMax(localMax.x, array[k][j][i].x);
2408 localMax.y = PetscMax(localMax.y, array[k][j][i].y);
2409 localMax.z = PetscMax(localMax.z, array[k][j][i].z);
2410 }
2411 }
2412 }
2413 ierr = DMDAVecRestoreArrayRead(dm, fieldVec, &array); CHKERRQ(ierr);
2414
2415 ierr = MPI_Allreduce(&localMin, &globalMin, 3, MPIU_REAL, MPI_MIN, PETSC_COMM_WORLD); CHKERRQ(ierr);
2416 ierr = MPI_Allreduce(&localMax, &globalMax, 3, MPIU_REAL, MPI_MAX, PETSC_COMM_WORLD); CHKERRQ(ierr);
2417
2418 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, " [Rank %d] Local X-Range: [ %11.4e , %11.4e ]\n", user->simCtx->rank, localMin.x, localMax.x);
2419 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, " [Rank %d] Local Y-Range: [ %11.4e , %11.4e ]\n", user->simCtx->rank, localMin.y, localMax.y);
2420 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, " [Rank %d] Local Z-Range: [ %11.4e , %11.4e ]\n", user->simCtx->rank, localMin.z, localMax.z);
2421 ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT); CHKERRQ(ierr);
2422
2423 if (user->simCtx->rank == 0) {
2424 PetscPrintf(PETSC_COMM_SELF, " [Global] X-Range: [ %11.4e , %11.4e ]\n", globalMin.x, globalMax.x);
2425 PetscPrintf(PETSC_COMM_SELF, " [Global] Y-Range: [ %11.4e , %11.4e ]\n", globalMin.y, globalMax.y);
2426 PetscPrintf(PETSC_COMM_SELF, " [Global] Z-Range: [ %11.4e , %11.4e ]\n", globalMin.z, globalMax.z);
2427 }
2428
2429 } else {
2430 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);
2431 }
2432
2433 // --- 5. Final barrier for clean output ordering ---
2434 ierr = MPI_Barrier(PETSC_COMM_WORLD); CHKERRQ(ierr);
2435 if (user->simCtx->rank == 0) {
2436 PetscPrintf(PETSC_COMM_SELF, "--------------------------------------------\n\n");
2437 }
2438
2439 PetscFunctionReturn(0);
2440}
2441
2442#undef __FUNCT__
2443#define __FUNCT__ "LOG_FIELD_ANATOMY"
2444/**
2445 * @brief Implementation of \ref LOG_FIELD_ANATOMY().
2446 * @details Full API contract is documented with the header declaration in `include/logging.h`.
2447 * @see LOG_FIELD_ANATOMY()
2448 */
2449PetscErrorCode LOG_FIELD_ANATOMY(UserCtx *user, const char *field_name, const char *stage_name)
2450{
2451 PetscErrorCode ierr;
2452 DMDALocalInfo info;
2453 PetscMPIInt rank;
2454
2455 Vec vec_local = NULL;
2456 DM dm = NULL;
2457 PetscInt dof = 0;
2458 char data_layout[20];
2459 char dominant_dir = '\0'; // 'x', 'y', 'z' for face-centered, 'm' for mixed (Ucont)
2460
2461 PetscFunctionBeginUser;
2462 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
2463
2464 // --- 1. Map string name to PETSc objects and determine data layout ---
2465 if (strcasecmp(field_name, "Ucat") == 0) {
2466 vec_local = user->lUcat; dm = user->fda; dof = 3; strcpy(data_layout, "Cell-Centered");
2467 } else if (strcasecmp(field_name, "P") == 0) {
2468 vec_local = user->lP; dm = user->da; dof = 1; strcpy(data_layout, "Cell-Centered");
2469 } else if (strcasecmp(field_name, "Diffusivity") == 0) {
2470 vec_local = user->lDiffusivity; dm = user->da; dof = 1; strcpy(data_layout, "Cell-Centered");
2471 } else if (strcasecmp(field_name, "DiffusivityGradient") == 0) {
2472 vec_local = user->lDiffusivityGradient; dm = user->fda; dof = 3; strcpy(data_layout, "Cell-Centered");
2473 } else if (strcasecmp(field_name, "Psi") == 0) {
2474 vec_local = user->lPsi; dm = user->da; dof = 1; strcpy(data_layout, "Cell-Centered");
2475 } else if (strcasecmp(field_name, "Center-Coordinates") == 0) {
2476 vec_local = user->lCent; dm = user->fda; dof = 3; strcpy(data_layout, "Cell-Centered");
2477 } else if (strcasecmp(field_name, "Ucont") == 0) {
2478 vec_local = user->lUcont; dm = user->fda; dof = 3; strcpy(data_layout, "Face-Centered"); dominant_dir = 'm'; // Mixed
2479 } else if (strcasecmp(field_name, "Csi") == 0 || strcasecmp(field_name, "X-Face-Centers") == 0) {
2480 vec_local = (strcasecmp(field_name, "Csi") == 0) ? user->lCsi : user->Centx;
2481 dm = user->fda; dof = 3; strcpy(data_layout, "Face-Centered"); dominant_dir = 'x';
2482 } else if (strcasecmp(field_name, "Eta") == 0 || strcasecmp(field_name, "Y-Face-Centers") == 0) {
2483 vec_local = (strcasecmp(field_name, "Eta") == 0) ? user->lEta : user->Centy;
2484 dm = user->fda; dof = 3; strcpy(data_layout, "Face-Centered"); dominant_dir = 'y';
2485 } else if (strcasecmp(field_name, "Zet") == 0 || strcasecmp(field_name, "Z-Face-Centers") == 0) {
2486 vec_local = (strcasecmp(field_name, "Zet") == 0) ? user->lZet : user->Centz;
2487 dm = user->fda; dof = 3; strcpy(data_layout, "Face-Centered"); dominant_dir = 'z';
2488 } else if (strcasecmp(field_name, "Coordinates") == 0) {
2489 ierr = DMGetCoordinatesLocal(user->da, &vec_local); CHKERRQ(ierr);
2490 dm = user->fda; dof = 3; strcpy(data_layout, "Node-Centered");
2491 } else if (strcasecmp(field_name, "CornerField")== 0){
2492 vec_local = user->lCellFieldAtCorner; strcpy(data_layout, "Node-Centered");
2493 PetscInt bs = 1;
2494 ierr = VecGetBlockSize(user->CellFieldAtCorner, &bs); CHKERRQ(ierr);
2495 dof = bs;
2496 if(dof == 1) dm = user->da;
2497 else dm = user->fda;
2498 } else {
2499 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Unknown field name for LOG_FIELD_ANATOMY: %s", field_name);
2500 }
2501
2502 // --- 2. Get Grid Info and Array Pointers ---
2503 ierr = DMDAGetLocalInfo(dm, &info); CHKERRQ(ierr);
2504
2505 ierr = PetscBarrier(NULL);
2506 PetscPrintf(PETSC_COMM_WORLD, "\n--- Field Anatomy Log: [%s] | Stage: [%s] | Layout: [%s] ---\n", field_name, stage_name, data_layout);
2507
2508 // Global physical dimensions (number of cells)
2509 PetscInt im_phys = user->IM;
2510 PetscInt jm_phys = user->JM;
2511 PetscInt km_phys = user->KM;
2512
2513 // Slice through the center of the local domain
2514 PetscInt i_mid = (PetscInt)(info.xs + 0.5 * info.xm) - 1;
2515 PetscInt j_mid = (PetscInt)(info.ys + 0.5 * info.ym) - 1;
2516 PetscInt k_mid = (PetscInt)(info.zs + 0.5 * info.zm) - 1;
2517
2518 // --- 3. Print Boundary Information based on Data Layout ---
2519
2520 // ======================================================================
2521 // === CASE 1: Cell-Centered Fields (Ucat, P) - USES SHIFTED INDEX ===
2522 // ======================================================================
2523 if (strcmp(data_layout, "Cell-Centered") == 0) {
2524 const void *l_arr;
2525 ierr = DMDAVecGetArrayRead(dm, vec_local, (void*)&l_arr); CHKERRQ(ierr);
2526
2527
2528 // --- I-Direction Boundaries ---
2529 if (info.xs == 0) { // Rank on -Xi boundary
2530 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, I-DIR]: Idx %2d (Ghost for Cell[k][j][0]) = ", rank, 0);
2531 if(dof==1) PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f)\n", ((const PetscReal***)l_arr)[k_mid][j_mid][0]);
2532 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);
2533
2534 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, I-DIR]: Idx %2d (Value for Cell[k][j][0]) = ", rank, 1);
2535 if(dof==1) PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f)\n", ((const PetscReal***)l_arr)[k_mid][j_mid][1]);
2536 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);
2537 }
2538 if (info.xs + info.xm == info.mx) { // Rank on +Xi boundary
2539 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, I-DIR]: Idx %2d (Value for Cell[k][j][%d]) = ", rank, im_phys - 1, im_phys - 2);
2540 if(dof==1) PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f)\n", ((const PetscReal***)l_arr)[k_mid][j_mid][im_phys - 1]);
2541 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);
2542
2543 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, I-DIR]: Idx %2d (Ghost for Cell[k][j][%d]) = ", rank, im_phys, im_phys - 2);
2544 if(dof==1) PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f)\n", ((const PetscReal***)l_arr)[k_mid][j_mid][im_phys]);
2545 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);
2546 }
2547
2548 // --- J-Direction Boundaries ---
2549 if (info.ys == 0) { // Rank on -Eta boundary
2550 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, J-DIR]: Jdx %2d (Ghost for Cell[k][0][i]) = ", rank, 0);
2551 if(dof==1) PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f)\n", ((const PetscReal***)l_arr)[k_mid][0][i_mid]);
2552 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);
2553
2554 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, J-DIR]: Jdx %2d (Value for Cell[k][0][i]) = ", rank, 1);
2555 if(dof==1) PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f)\n", ((const PetscReal***)l_arr)[k_mid][1][i_mid]);
2556 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);
2557 }
2558
2559 if (info.ys + info.ym == info.my) { // Rank on +Eta boundary
2560 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, J-DIR]: Jdx %2d (Value for Cell[k][%d][i]) = ", rank, jm_phys - 1, jm_phys - 2);
2561 if(dof==1) PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f)\n", ((const PetscReal***)l_arr)[k_mid][jm_phys - 1][i_mid]);
2562 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);
2563
2564 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, J-DIR]: Jdx %2d (Ghost for Cell[k][%d][i]) = ", rank, jm_phys, jm_phys - 2);
2565 if(dof==1) PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f)\n", ((const PetscReal***)l_arr)[k_mid][jm_phys][i_mid]);
2566 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);
2567 }
2568
2569 // --- K-Direction Boundaries ---
2570 if (info.zs == 0) { // Rank on -Zeta boundary
2571 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, K-DIR]: Kdx %2d (Ghost for Cell[0][j][i]) = ", rank, 0);
2572 if(dof==1) PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f)\n", ((const PetscReal***)l_arr)[0][j_mid][i_mid]);
2573 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);
2574 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, K-DIR]: Kdx %2d (Value for Cell[0][j][i]) = ", rank, 1);
2575 if(dof==1) PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f)\n", ((const PetscReal***)l_arr)[1][j_mid][i_mid]);
2576 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);
2577 }
2578 if (info.zs + info.zm == info.mz) { // Rank on +Zeta boundary
2579 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, K-DIR]: Kdx %2d (Value for Cell[%d][j][i]) = ", rank, km_phys - 1, km_phys - 2);
2580 if(dof==1) PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f)\n", ((const PetscReal***)l_arr)[km_phys - 1][j_mid][i_mid]);
2581 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);
2582 PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[Rank %d, K-DIR]: Kdx %2d (Ghost for Cell[%d][j][i]) = ", rank, km_phys, km_phys - 2);
2583 if(dof==1) PetscSynchronizedPrintf(PETSC_COMM_WORLD, "(%.5f)\n", ((const PetscReal***)l_arr)[km_phys][j_mid][i_mid]);
2584 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);
2585 }
2586 ierr = DMDAVecRestoreArrayRead(dm, vec_local, (void*)&l_arr); CHKERRQ(ierr);
2587 }
2588 // ======================================================================
2589 // === CASE 2: Face-Centered Fields - NUANCED DIRECTIONAL LOGIC ===
2590 // ======================================================================
2591 else if (strcmp(data_layout, "Face-Centered") == 0) {
2592 const Cmpnts ***l_arr;
2593 ierr = DMDAVecGetArrayRead(dm, vec_local, (void*)&l_arr); CHKERRQ(ierr);
2594
2595 // --- I-Direction Boundaries ---
2596 if (info.xs == 0) { // Rank on -Xi boundary
2597 if (dominant_dir == 'x') { // Node-like in I-dir
2598 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);
2599 } else if (dominant_dir == 'y' || dominant_dir == 'z') { // Cell-like in I-dir
2600 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);
2601 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);
2602 } else if (dominant_dir == 'm') { // Ucont: Mixed
2603 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);
2604 }
2605 }
2606 if (info.xs + info.xm == info.mx) { // Rank on +Xi boundary
2607 if (dominant_dir == 'x') { // Node-like in I-dir
2608 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);
2609 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);
2610 } else if (dominant_dir == 'y' || dominant_dir == 'z') { // Cell-like in I-dir
2611 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);
2612 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);
2613 } else if (dominant_dir == 'm') { // Ucont: Mixed
2614 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);
2615 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);
2616 }
2617 }
2618
2619 // --- J-Direction Boundaries ---
2620 if (info.ys == 0) { // Rank on -Eta boundary
2621 if (dominant_dir == 'y') { // Node-like in J-dir
2622 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);
2623 } else if (dominant_dir == 'x' || dominant_dir == 'z') { // Cell-like in J-dir
2624 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);
2625 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);
2626 } else if (dominant_dir == 'm') { // Ucont: Mixed
2627 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);
2628 }
2629 }
2630 if (info.ys + info.ym == info.my) { // Rank on +Eta boundary
2631 if (dominant_dir == 'y') { // Node-like in J-dir
2632 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);
2633 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);
2634 } else if (dominant_dir == 'x' || dominant_dir == 'z') { // Cell-like in J-dir
2635 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);
2636 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);
2637 } else if (dominant_dir == 'm') { // Ucont: Mixed
2638 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);
2639 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);
2640 }
2641 }
2642
2643 // --- K-Direction Boundaries ---
2644 if (info.zs == 0) { // Rank on -Zeta boundary
2645 if (dominant_dir == 'z') { // Node-like in K-dir
2646 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);
2647 } else if (dominant_dir == 'x' || dominant_dir == 'y') { // Cell-like in K-dir
2648 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);
2649 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);
2650 } else if (dominant_dir == 'm') { // Ucont: Mixed
2651 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);
2652 }
2653 }
2654 if (info.zs + info.zm == info.mz) { // Rank on +Zeta boundary
2655 if (dominant_dir == 'z') { // Node-like in K-dir
2656 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);
2657 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);
2658 } else if (dominant_dir == 'x' || dominant_dir == 'y') { // Cell-like in K-dir
2659 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);
2660 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);
2661 } else if (dominant_dir == 'm') { // Ucont: Mixed
2662 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);
2663 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);
2664
2665 }
2666 }
2667 ierr = DMDAVecRestoreArrayRead(dm, vec_local, (void*)&l_arr); CHKERRQ(ierr);
2668 }
2669 // ======================================================================
2670 // === CASE 3: Node-Centered Fields - USES DIRECT INDEX ===
2671 // ======================================================================
2672 else if (strcmp(data_layout, "Node-Centered") == 0) {
2673 const Cmpnts ***l_arr;
2674 ierr = DMDAVecGetArrayRead(dm, vec_local, (void*)&l_arr); CHKERRQ(ierr);
2675
2676 // --- I-Direction Boundaries ---
2677 if (info.xs == 0) {
2678 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);
2679 }
2680 if (info.xs + info.xm == info.mx) {
2681 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);
2682 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);
2683 }
2684 // --- J-Direction Boundaries ---
2685 if (info.ys == 0) {
2686 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);
2687 }
2688 if (info.ys + info.ym == info.my) {
2689 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);
2690 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);
2691 }
2692 // --- K-Direction Boundaries ---
2693 if (info.zs == 0) {
2694 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);
2695 }
2696 if(info.zs + info.zm == info.mz) {
2697 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);
2698 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);
2699 }
2700 ierr = DMDAVecRestoreArrayRead(dm, vec_local, (void*)&l_arr); CHKERRQ(ierr);
2701 }
2702 else {
2703 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "LOG_FIELD_ANATOMY encountered an unknown data layout: %s", data_layout);
2704 }
2705
2706 ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT); CHKERRQ(ierr);
2707 ierr = PetscBarrier(NULL);
2708 PetscFunctionReturn(0);
2709}
2710
2711#undef __FUNCT__
2712#define __FUNCT__ "LOG_INTERPOLATION_ERROR"
2713/**
2714 * @brief Implementation of \ref LOG_INTERPOLATION_ERROR().
2715 * @details Full API contract (arguments, ownership, side effects) is documented with
2716 * the header declaration in `include/logging.h`.
2717 * @see LOG_INTERPOLATION_ERROR()
2718 */
2719PetscErrorCode LOG_INTERPOLATION_ERROR(UserCtx *user)
2720{
2721 SimCtx *simCtx = user->simCtx;
2722 PetscErrorCode ierr;
2723 DM swarm = user->swarm;
2724 Vec positionVec, analyticalvelocityVec, velocityVec, errorVec;
2725 PetscReal Interpolation_error = 0.0;
2726 PetscReal Maximum_Interpolation_error = 0.0;
2727 PetscReal AnalyticalSolution_magnitude = 0.0;
2728 PetscReal ErrorPercentage = 0.0;
2729
2730 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Creating global vectors.\n");
2731 ierr = DMSwarmCreateGlobalVectorFromField(swarm, "position", &positionVec); CHKERRQ(ierr);
2732 ierr = DMSwarmCreateGlobalVectorFromField(swarm, "velocity", &velocityVec); CHKERRQ(ierr);
2733
2734 ierr = VecDuplicate(positionVec, &analyticalvelocityVec); CHKERRQ(ierr);
2735 ierr = VecCopy(positionVec, analyticalvelocityVec); CHKERRQ(ierr);
2736
2737 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Computing analytical solution.\n");
2738 ierr = SetAnalyticalSolutionForParticles(analyticalvelocityVec, simCtx); CHKERRQ(ierr);
2739
2740 ierr = VecDuplicate(analyticalvelocityVec, &errorVec); CHKERRQ(ierr);
2741 ierr = VecCopy(analyticalvelocityVec, errorVec); CHKERRQ(ierr);
2742
2743 ierr = VecNorm(analyticalvelocityVec, NORM_2, &AnalyticalSolution_magnitude); CHKERRQ(ierr);
2744
2745 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Computing error.\n");
2746 ierr = VecAXPY(errorVec, -1.0, velocityVec); CHKERRQ(ierr);
2747 ierr = VecNorm(errorVec, NORM_2, &Interpolation_error); CHKERRQ(ierr);
2748 ierr = VecNorm(errorVec,NORM_INFINITY,&Maximum_Interpolation_error); CHKERRQ(ierr);
2749
2750 ErrorPercentage = (AnalyticalSolution_magnitude > 0) ?
2751 (Interpolation_error / AnalyticalSolution_magnitude * 100.0) : 0.0;
2752
2753 /* --- CSV output (always, rank 0 only) --- */
2754 if (simCtx->rank == 0) {
2755 char csv_path[PETSC_MAX_PATH_LEN + 32];
2756 ierr = PetscSNPrintf(csv_path, sizeof(csv_path), "%s/interpolation_error.csv", simCtx->log_dir); CHKERRQ(ierr);
2757 FILE *f = fopen(csv_path, "a");
2758 if (f) {
2759 if (ftell(f) == 0) {
2760 fprintf(f, "step,time,L2_error,Linf_error,L2_analytical,error_pct\n");
2761 }
2762 PetscReal t = (PetscReal)simCtx->ti * simCtx->dt;
2763 fprintf(f, "%d,%.6e,%.6e,%.6e,%.6e,%.4f\n",
2764 (int)simCtx->step, t,
2765 Interpolation_error, Maximum_Interpolation_error,
2766 AnalyticalSolution_magnitude, ErrorPercentage);
2767 fclose(f);
2768 }
2769 }
2770
2771 /* --- Console output (only at INFO level or above) --- */
2772 if (get_log_level() >= LOG_INFO) {
2773 LOG_ALLOW(GLOBAL, LOG_INFO, "Interpolation error (%%): %g\n", ErrorPercentage);
2774 PetscPrintf(PETSC_COMM_WORLD, "Interpolation error (%%): %g\n", ErrorPercentage);
2775 LOG_ALLOW(GLOBAL, LOG_INFO, "Maximum Interpolation error: %g\n", Maximum_Interpolation_error);
2776 PetscPrintf(PETSC_COMM_WORLD, "Maximum Interpolation error: %g\n", Maximum_Interpolation_error);
2777 }
2778
2779 ierr = VecDestroy(&analyticalvelocityVec); CHKERRQ(ierr);
2780 ierr = VecDestroy(&errorVec); CHKERRQ(ierr);
2781 ierr = DMSwarmDestroyGlobalVectorFromField(swarm, "position", &positionVec); CHKERRQ(ierr);
2782 ierr = DMSwarmDestroyGlobalVectorFromField(swarm, "velocity", &velocityVec); CHKERRQ(ierr);
2783
2784 return 0;
2785}
2786
2787#undef __FUNCT__
2788#define __FUNCT__ "LOG_SCATTER_METRICS"
2789/**
2790 * @brief Implementation of \ref LOG_SCATTER_METRICS().
2791 * @details Full API contract (arguments, ownership, side effects) is documented with
2792 * the header declaration in `include/logging.h`.
2793 * @see LOG_SCATTER_METRICS()
2794 */
2795PetscErrorCode LOG_SCATTER_METRICS(UserCtx *user)
2796{
2797 PetscErrorCode ierr;
2798 SimCtx *simCtx = NULL;
2799 DMDALocalInfo info;
2800 PetscInt xs, xe, ys, ye, zs, ze, mx, my, mz;
2801 PetscInt lxs, lxe, lys, lye, lzs, lze;
2802 Vec reference_vec = NULL;
2803 PetscReal ***psi = NULL;
2804 PetscReal ***psi_ref = NULL;
2805 PetscReal ***aj = NULL;
2806 PetscReal ***count = NULL;
2807 PetscReal *particle_psi = NULL;
2808 PetscInt nlocal = 0;
2809 PetscReal local_l1 = 0.0, global_l1 = 0.0;
2810 PetscReal local_l2_sq = 0.0, global_l2_sq = 0.0;
2811 PetscReal local_linf = 0.0, global_linf = 0.0;
2812 PetscReal local_ref_l2_sq = 0.0, global_ref_l2_sq = 0.0;
2813 PetscReal local_grid_integral = 0.0, global_grid_integral = 0.0;
2814 PetscReal local_domain_volume = 0.0, global_domain_volume = 0.0;
2815 PetscReal local_particle_sum = 0.0, global_particle_sum = 0.0;
2816 PetscInt64 local_particle_count = 0, global_particle_count = 0;
2817 PetscInt64 local_cell_count = 0, global_cell_count = 0;
2818 PetscInt64 local_occupied_count = 0, global_occupied_count = 0;
2819 PetscReal particle_integral = 0.0;
2820 PetscReal occupancy_fraction = 0.0;
2821 PetscReal mean_particles_per_occupied_cell = 0.0;
2822 PetscReal l2_error = 0.0;
2823 PetscReal relative_l2_error = 0.0;
2824
2825 PetscFunctionBeginUser;
2826 if (!user) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "UserCtx cannot be NULL.");
2827 simCtx = user->simCtx;
2828 if (!VerificationScalarOverrideActive(simCtx) || !user->swarm || !user->Psi || !user->ParticleCount) {
2829 PetscFunctionReturn(0);
2830 }
2831
2832 info = user->info;
2833 xs = info.xs; xe = info.xs + info.xm;
2834 ys = info.ys; ye = info.ys + info.ym;
2835 zs = info.zs; ze = info.zs + info.zm;
2836 mx = info.mx; my = info.my; mz = info.mz;
2837 lxs = (xs == 0) ? xs + 1 : xs; lxe = (xe == mx) ? xe - 1 : xe;
2838 lys = (ys == 0) ? ys + 1 : ys; lye = (ye == my) ? ye - 1 : ye;
2839 lzs = (zs == 0) ? zs + 1 : zs; lze = (ze == mz) ? ze - 1 : ze;
2840
2841 ierr = VecDuplicate(user->Psi, &reference_vec); CHKERRQ(ierr);
2842 ierr = SetAnalyticalScalarFieldAtCellCenters(user, reference_vec); CHKERRQ(ierr);
2843
2844 ierr = DMDAVecGetArrayRead(user->da, user->Psi, &psi); CHKERRQ(ierr);
2845 ierr = DMDAVecGetArrayRead(user->da, reference_vec, &psi_ref); CHKERRQ(ierr);
2846 ierr = DMDAVecGetArrayRead(user->da, user->Aj, &aj); CHKERRQ(ierr);
2847 ierr = DMDAVecGetArrayRead(user->da, user->ParticleCount, &count); CHKERRQ(ierr);
2848
2849 for (PetscInt k = lzs; k < lze; ++k) {
2850 for (PetscInt j = lys; j < lye; ++j) {
2851 for (PetscInt i = lxs; i < lxe; ++i) {
2852 const PetscReal cell_volume = (PetscAbsReal(aj[k][j][i]) > 1.0e-14) ? (1.0 / aj[k][j][i]) : 0.0;
2853 const PetscReal err = psi[k][j][i] - psi_ref[k][j][i];
2854 local_cell_count += 1;
2855 local_domain_volume += cell_volume;
2856 local_grid_integral += psi[k][j][i] * cell_volume;
2857 local_l1 += PetscAbsReal(err) * cell_volume;
2858 local_l2_sq += err * err * cell_volume;
2859 local_ref_l2_sq += psi_ref[k][j][i] * psi_ref[k][j][i] * cell_volume;
2860 local_linf = PetscMax(local_linf, PetscAbsReal(err));
2861 if (count[k][j][i] > 0.0) local_occupied_count += 1;
2862 }
2863 }
2864 }
2865
2866 ierr = DMDAVecRestoreArrayRead(user->da, user->ParticleCount, &count); CHKERRQ(ierr);
2867 ierr = DMDAVecRestoreArrayRead(user->da, user->Aj, &aj); CHKERRQ(ierr);
2868 ierr = DMDAVecRestoreArrayRead(user->da, reference_vec, &psi_ref); CHKERRQ(ierr);
2869 ierr = DMDAVecRestoreArrayRead(user->da, user->Psi, &psi); CHKERRQ(ierr);
2870 ierr = VecDestroy(&reference_vec); CHKERRQ(ierr);
2871
2872 ierr = DMSwarmGetLocalSize(user->swarm, &nlocal); CHKERRQ(ierr);
2873 local_particle_count = (PetscInt64)nlocal;
2874 if (nlocal > 0) {
2875 ierr = DMSwarmGetField(user->swarm, "Psi", NULL, NULL, (void **)&particle_psi); CHKERRQ(ierr);
2876 for (PetscInt p = 0; p < nlocal; ++p) local_particle_sum += particle_psi[p];
2877 ierr = DMSwarmRestoreField(user->swarm, "Psi", NULL, NULL, (void **)&particle_psi); CHKERRQ(ierr);
2878 }
2879
2880 ierr = MPI_Allreduce(&local_l1, &global_l1, 1, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD); CHKERRMPI(ierr);
2881 ierr = MPI_Allreduce(&local_l2_sq, &global_l2_sq, 1, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD); CHKERRMPI(ierr);
2882 ierr = MPI_Allreduce(&local_linf, &global_linf, 1, MPIU_REAL, MPI_MAX, PETSC_COMM_WORLD); CHKERRMPI(ierr);
2883 ierr = MPI_Allreduce(&local_ref_l2_sq, &global_ref_l2_sq, 1, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD); CHKERRMPI(ierr);
2884 ierr = MPI_Allreduce(&local_grid_integral, &global_grid_integral, 1, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD); CHKERRMPI(ierr);
2885 ierr = MPI_Allreduce(&local_domain_volume, &global_domain_volume, 1, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD); CHKERRMPI(ierr);
2886 ierr = MPI_Allreduce(&local_particle_sum, &global_particle_sum, 1, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD); CHKERRMPI(ierr);
2887 ierr = MPI_Allreduce(&local_particle_count, &global_particle_count, 1, MPIU_INT64, MPI_SUM, PETSC_COMM_WORLD); CHKERRMPI(ierr);
2888 ierr = MPI_Allreduce(&local_cell_count, &global_cell_count, 1, MPIU_INT64, MPI_SUM, PETSC_COMM_WORLD); CHKERRMPI(ierr);
2889 ierr = MPI_Allreduce(&local_occupied_count, &global_occupied_count, 1, MPIU_INT64, MPI_SUM, PETSC_COMM_WORLD); CHKERRMPI(ierr);
2890
2891 l2_error = PetscSqrtReal(global_l2_sq);
2892 relative_l2_error = (global_ref_l2_sq > 0.0) ? (l2_error / PetscSqrtReal(global_ref_l2_sq)) : 0.0;
2893 occupancy_fraction = (global_cell_count > 0) ? ((PetscReal)global_occupied_count / (PetscReal)global_cell_count) : 0.0;
2894 mean_particles_per_occupied_cell =
2895 (global_occupied_count > 0) ? ((PetscReal)global_particle_count / (PetscReal)global_occupied_count) : 0.0;
2896 particle_integral =
2897 (global_particle_count > 0) ? (global_domain_volume * global_particle_sum / (PetscReal)global_particle_count) : 0.0;
2898
2899 if (simCtx->rank == 0) {
2900 char csv_path[PETSC_MAX_PATH_LEN + 32];
2901 FILE *f = NULL;
2902 ierr = PetscSNPrintf(csv_path, sizeof(csv_path), "%s/scatter_metrics.csv", simCtx->log_dir); CHKERRQ(ierr);
2903 f = fopen(csv_path, "a");
2904 if (f) {
2905 if (ftell(f) == 0) {
2906 fprintf(f,
2907 "step,time,total_particles,total_cells,occupied_cells,occupancy_fraction,"
2908 "mean_particles_per_occupied_cell,particle_integral,grid_integral,"
2909 "conservation_error_abs,L1_error,L2_error,Linf_error,relative_L2_error\n");
2910 }
2911 fprintf(f, "%d,%.6e,%lld,%lld,%lld,%.6e,%.6e,%.6e,%.6e,%.6e,%.6e,%.6e,%.6e,%.6e\n",
2912 (int)simCtx->step,
2913 (double)simCtx->ti,
2914 (long long)global_particle_count,
2915 (long long)global_cell_count,
2916 (long long)global_occupied_count,
2917 (double)occupancy_fraction,
2918 (double)mean_particles_per_occupied_cell,
2919 (double)particle_integral,
2920 (double)global_grid_integral,
2921 (double)PetscAbsReal(global_grid_integral - particle_integral),
2922 (double)global_l1,
2923 (double)l2_error,
2924 (double)global_linf,
2925 (double)relative_l2_error);
2926 fclose(f);
2927 }
2928 }
2929
2930 if (get_log_level() >= LOG_INFO) {
2931 LOG_ALLOW(GLOBAL, LOG_INFO, "Scatter relative L2 error: %.6e\n", (double)relative_l2_error);
2932 LOG_ALLOW(GLOBAL, LOG_INFO, "Scatter occupancy fraction: %.6e\n", (double)occupancy_fraction);
2933 }
2934
2935 PetscFunctionReturn(0);
2936}
2937
2938#undef __FUNCT__
2939#define __FUNCT__ "ResetSearchMetrics"
2940/**
2941 * @brief Implementation of \ref ResetSearchMetrics().
2942 * @details Full API contract (arguments, ownership, side effects) is documented with
2943 * the header declaration in `include/logging.h`.
2944 * @see ResetSearchMetrics()
2945 */
2946PetscErrorCode ResetSearchMetrics(SimCtx *simCtx)
2947{
2948 PetscFunctionBeginUser;
2949 if (!simCtx) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "SimCtx cannot be NULL for ResetSearchMetrics.");
2950
2951 simCtx->searchMetrics.searchAttempts = 0;
2952 simCtx->searchMetrics.searchPopulation = 0;
2953 simCtx->searchMetrics.searchLocatedCount = 0;
2954 simCtx->searchMetrics.searchLostCount = 0;
2955 simCtx->searchMetrics.traversalStepsSum = 0;
2956 simCtx->searchMetrics.reSearchCount = 0;
2957 simCtx->searchMetrics.maxTraversalSteps = 0;
2958 simCtx->searchMetrics.maxTraversalFailCount = 0;
2959 simCtx->searchMetrics.tieBreakCount = 0;
2960 simCtx->searchMetrics.boundaryClampCount = 0;
2961 simCtx->searchMetrics.bboxGuessSuccessCount = 0;
2962 simCtx->searchMetrics.bboxGuessFallbackCount = 0;
2963 simCtx->searchMetrics.maxParticlePassDepth = 0;
2964 simCtx->searchMetrics.currentSettlementPass = 0;
2965
2966 PetscFunctionReturn(0);
2967}
2968
2969#undef __FUNCT__
2970#define __FUNCT__ "LOG_SEARCH_METRICS"
2971/**
2972 * @brief Implementation of \ref LOG_SEARCH_METRICS().
2973 * @details Full API contract (arguments, ownership, side effects) is documented with
2974 * the header declaration in `include/logging.h`.
2975 * @see LOG_SEARCH_METRICS()
2976 */
2977PetscErrorCode LOG_SEARCH_METRICS(UserCtx *user)
2978{
2979 PetscErrorCode ierr;
2980 SimCtx *simCtx = NULL;
2981 PetscInt totalParticles = 0;
2982 PetscReal local_metrics[SEARCH_METRIC_REDUCTION_LEN] = {0.0};
2983 PetscReal global_metrics[SEARCH_METRIC_REDUCTION_LEN] = {0.0};
2984 PetscReal meanTraversalSteps = 0.0;
2985 PetscReal searchFailureFraction = 0.0;
2986 PetscReal searchWorkIndex = 0.0;
2987 PetscReal reSearchFraction = 0.0;
2988 long long searchAttempts = 0;
2989 long long searchPopulation = 0;
2990 long long searchLocatedCount = 0;
2991 long long searchLostCount = 0;
2992 long long traversalStepsSum = 0;
2993 long long reSearchCount = 0;
2994 long long tieBreakCount = 0;
2995 long long boundaryClampCount = 0;
2996 long long bboxGuessSuccessCount = 0;
2997 long long bboxGuessFallbackCount = 0;
2998 long long maxTraversalFailCount = 0;
2999 long long maxTraversalSteps = 0;
3000 long long maxPassDepth = 0;
3001 MPI_Op reduction_op = MPI_OP_NULL;
3002
3003 PetscFunctionBeginUser;
3004 if (!user || !user->simCtx) {
3005 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "UserCtx and SimCtx are required for LOG_SEARCH_METRICS.");
3006 }
3007 simCtx = user->simCtx;
3008
3009 if (simCtx->np <= 0) {
3010 PetscFunctionReturn(0);
3011 }
3012
3013 ierr = DMSwarmGetSize(user->swarm, &totalParticles); CHKERRQ(ierr);
3014
3015 local_metrics[SEARCH_METRIC_SUM_SEARCH_ATTEMPTS] = (PetscReal)simCtx->searchMetrics.searchAttempts;
3016 local_metrics[SEARCH_METRIC_SUM_SEARCH_POPULATION] = (PetscReal)simCtx->searchMetrics.searchPopulation;
3017 local_metrics[SEARCH_METRIC_SUM_SEARCH_LOCATED] = (PetscReal)simCtx->searchMetrics.searchLocatedCount;
3018 local_metrics[SEARCH_METRIC_SUM_SEARCH_LOST] = (PetscReal)simCtx->searchMetrics.searchLostCount;
3019 local_metrics[SEARCH_METRIC_SUM_TRAVERSAL_STEPS] = (PetscReal)simCtx->searchMetrics.traversalStepsSum;
3020 local_metrics[SEARCH_METRIC_SUM_RESEARCH] = (PetscReal)simCtx->searchMetrics.reSearchCount;
3021 local_metrics[SEARCH_METRIC_SUM_TIE_BREAKS] = (PetscReal)simCtx->searchMetrics.tieBreakCount;
3022 local_metrics[SEARCH_METRIC_SUM_BOUNDARY_CLAMPS] = (PetscReal)simCtx->searchMetrics.boundaryClampCount;
3023 local_metrics[SEARCH_METRIC_SUM_BBOX_GUESS_SUCCESS] = (PetscReal)simCtx->searchMetrics.bboxGuessSuccessCount;
3024 local_metrics[SEARCH_METRIC_SUM_BBOX_GUESS_FALLBACK] = (PetscReal)simCtx->searchMetrics.bboxGuessFallbackCount;
3025 local_metrics[SEARCH_METRIC_SUM_MAX_TRAVERSAL_FAILS] = (PetscReal)simCtx->searchMetrics.maxTraversalFailCount;
3026 local_metrics[SEARCH_METRIC_MAX_TRAVERSAL_STEPS] = (PetscReal)simCtx->searchMetrics.maxTraversalSteps;
3027 local_metrics[SEARCH_METRIC_MAX_PASS_DEPTH] = (PetscReal)simCtx->searchMetrics.maxParticlePassDepth;
3028
3029 ierr = MPI_Op_create(SearchMetricsReduceOp, PETSC_TRUE, &reduction_op); CHKERRMPI(ierr);
3030 ierr = MPI_Allreduce(local_metrics, global_metrics, SEARCH_METRIC_REDUCTION_LEN, MPIU_REAL, reduction_op, PETSC_COMM_WORLD); CHKERRMPI(ierr);
3031 ierr = MPI_Op_free(&reduction_op); CHKERRMPI(ierr);
3032 reduction_op = MPI_OP_NULL;
3033
3034 searchAttempts = (long long)PetscFloorReal(global_metrics[SEARCH_METRIC_SUM_SEARCH_ATTEMPTS] + 0.5);
3035 searchPopulation = (long long)PetscFloorReal(global_metrics[SEARCH_METRIC_SUM_SEARCH_POPULATION] + 0.5);
3036 searchLocatedCount = (long long)PetscFloorReal(global_metrics[SEARCH_METRIC_SUM_SEARCH_LOCATED] + 0.5);
3037 searchLostCount = (long long)PetscFloorReal(global_metrics[SEARCH_METRIC_SUM_SEARCH_LOST] + 0.5);
3038 traversalStepsSum = (long long)PetscFloorReal(global_metrics[SEARCH_METRIC_SUM_TRAVERSAL_STEPS] + 0.5);
3039 reSearchCount = (long long)PetscFloorReal(global_metrics[SEARCH_METRIC_SUM_RESEARCH] + 0.5);
3040 tieBreakCount = (long long)PetscFloorReal(global_metrics[SEARCH_METRIC_SUM_TIE_BREAKS] + 0.5);
3041 boundaryClampCount = (long long)PetscFloorReal(global_metrics[SEARCH_METRIC_SUM_BOUNDARY_CLAMPS] + 0.5);
3042 bboxGuessSuccessCount = (long long)PetscFloorReal(global_metrics[SEARCH_METRIC_SUM_BBOX_GUESS_SUCCESS] + 0.5);
3043 bboxGuessFallbackCount = (long long)PetscFloorReal(global_metrics[SEARCH_METRIC_SUM_BBOX_GUESS_FALLBACK] + 0.5);
3044 maxTraversalFailCount = (long long)PetscFloorReal(global_metrics[SEARCH_METRIC_SUM_MAX_TRAVERSAL_FAILS] + 0.5);
3045 maxTraversalSteps = (long long)PetscFloorReal(global_metrics[SEARCH_METRIC_MAX_TRAVERSAL_STEPS] + 0.5);
3046 maxPassDepth = (long long)PetscFloorReal(global_metrics[SEARCH_METRIC_MAX_PASS_DEPTH] + 0.5);
3047
3048 if (searchAttempts > 0) {
3049 meanTraversalSteps = (PetscReal)traversalStepsSum / (PetscReal)searchAttempts;
3050 }
3051 if (searchPopulation > 0) {
3052 searchFailureFraction = (PetscReal)searchLostCount / (PetscReal)searchPopulation;
3053 searchWorkIndex = (PetscReal)traversalStepsSum / (PetscReal)searchPopulation;
3054 reSearchFraction = (PetscReal)reSearchCount / (PetscReal)searchPopulation;
3055 }
3056
3057 if (simCtx->rank == 0) {
3058 char csv_path[PETSC_MAX_PATH_LEN + 32];
3059 FILE *f = NULL;
3060
3061 ierr = PetscSNPrintf(csv_path, sizeof(csv_path), "%s/search_metrics.csv", simCtx->log_dir); CHKERRQ(ierr);
3062 f = fopen(csv_path, "a");
3063 if (!f) {
3064 LOG_ALLOW(GLOBAL, LOG_WARNING, "LOG_SEARCH_METRICS: could not open '%s' for writing.\n", csv_path);
3065 } else {
3066 if (ftell(f) == 0) {
3067 fprintf(f,
3068 "step,time,total_particles,lost,lost_cumulative,migrated,migration_passes,search_attempts,"
3069 "mean_traversal_steps,max_traversal_steps,tie_break_count,boundary_clamp_count,"
3070 "bbox_guess_success_count,bbox_guess_fallback_count,max_particle_pass_depth,load_imbalance,"
3071 "search_population,search_located_count,search_lost_count,traversal_steps_sum,re_search_count,"
3072 "max_traversal_fail_count,search_failure_fraction,search_work_index,re_search_fraction\n");
3073 }
3074 fprintf(f,
3075 "%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",
3076 (int)simCtx->step,
3077 (double)simCtx->ti,
3078 (int)totalParticles,
3079 (int)simCtx->particlesLostLastStep,
3080 (int)simCtx->particlesLostCumulative,
3081 (int)simCtx->particlesMigratedLastStep,
3082 (int)simCtx->migrationPassesLastStep,
3083 searchAttempts,
3084 (double)meanTraversalSteps,
3085 maxTraversalSteps,
3086 tieBreakCount,
3087 boundaryClampCount,
3088 bboxGuessSuccessCount,
3089 bboxGuessFallbackCount,
3090 maxPassDepth,
3091 (double)simCtx->particleLoadImbalance,
3092 searchPopulation,
3093 searchLocatedCount,
3094 searchLostCount,
3095 traversalStepsSum,
3096 reSearchCount,
3097 maxTraversalFailCount,
3098 (double)searchFailureFraction,
3099 (double)searchWorkIndex,
3100 (double)reSearchFraction);
3101 fclose(f);
3102 }
3103 }
3104
3105 LOG_ALLOW(GLOBAL, LOG_DEBUG,
3106 "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",
3107 (double)searchFailureFraction,
3108 (double)searchWorkIndex,
3109 (double)reSearchFraction,
3110 (int)simCtx->particlesLostLastStep,
3111 (int)simCtx->particlesLostCumulative,
3112 (int)simCtx->particlesMigratedLastStep,
3113 (int)simCtx->migrationPassesLastStep,
3114 (double)meanTraversalSteps,
3115 maxTraversalSteps,
3116 tieBreakCount,
3117 maxPassDepth);
3118
3119 PetscFunctionReturn(0);
3120}
3121
3122#undef __FUNCT__
3123#define __FUNCT__ "CalculateAdvancedParticleMetrics"
3124/**
3125 * @brief Internal helper implementation: `CalculateAdvancedParticleMetrics()`.
3126 * @details Local to this translation unit.
3127 */
3128PetscErrorCode CalculateAdvancedParticleMetrics(UserCtx *user)
3129{
3130 PetscErrorCode ierr;
3131 SimCtx *simCtx = user->simCtx;
3132 PetscMPIInt size, rank;
3133
3134 PetscFunctionBeginUser;
3135 ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size); CHKERRQ(ierr);
3136 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
3137
3138 // --- 1. Particle Load Imbalance ---
3139 PetscInt nLocal, nGlobal, nLocalMax;
3140 ierr = DMSwarmGetLocalSize(user->swarm, &nLocal); CHKERRQ(ierr);
3141 ierr = DMSwarmGetSize(user->swarm, &nGlobal); CHKERRQ(ierr);
3142 ierr = MPI_Allreduce(&nLocal, &nLocalMax, 1, MPIU_INT, MPI_MAX, PETSC_COMM_WORLD); CHKERRQ(ierr);
3143
3144 PetscReal avg_per_rank = (size > 0) ? ((PetscReal)nGlobal / size) : 0.0;
3145 // Handle division by zero if there are no particles
3146 simCtx->particleLoadImbalance = (avg_per_rank > 1e-9) ? (nLocalMax / avg_per_rank) : 1.0;
3147
3148
3149 // --- 2. Number of Occupied Cells ---
3150 // This part requires access to the user->ParticleCount vector.
3151 PetscInt local_occupied_cells = 0;
3152 PetscInt global_occupied_cells;
3153 const PetscScalar *count_array;
3154 PetscInt vec_local_size;
3155
3156 ierr = VecGetLocalSize(user->ParticleCount, &vec_local_size); CHKERRQ(ierr);
3157 ierr = VecGetArrayRead(user->ParticleCount, &count_array); CHKERRQ(ierr);
3158
3159 for (PetscInt i = 0; i < vec_local_size; ++i) {
3160 if (count_array[i] > 0.5) { // Use 0.5 to be safe with floating point
3161 local_occupied_cells++;
3162 }
3163 }
3164 ierr = VecRestoreArrayRead(user->ParticleCount, &count_array); CHKERRQ(ierr);
3165
3166 ierr = MPI_Allreduce(&local_occupied_cells, &global_occupied_cells, 1, MPIU_INT, MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
3167 simCtx->occupiedCellCount = global_occupied_cells;
3168
3169 LOG_ALLOW_SYNC(GLOBAL, LOG_INFO, "[Rank %d] Advanced Metrics: Imbalance=%.2f, OccupiedCells=%d\n", rank, simCtx->particleLoadImbalance, simCtx->occupiedCellCount);
3170
3171 PetscFunctionReturn(0);
3172}
3173
3174#undef __FUNCT__
3175#define __FUNCT__ "LOG_PARTICLE_METRICS"
3176/**
3177 * @brief Implementation of \ref LOG_PARTICLE_METRICS().
3178 * @details Full API contract (arguments, ownership, side effects) is documented with
3179 * the header declaration in `include/logging.h`.
3180 * @see LOG_PARTICLE_METRICS()
3181 */
3182PetscErrorCode LOG_PARTICLE_METRICS(UserCtx *user, const char *stageName)
3183{
3184 PetscErrorCode ierr;
3185 PetscMPIInt rank;
3186 SimCtx *simCtx = user->simCtx;
3187 const char *stage_label = (stageName && stageName[0] != '\0') ? stageName : "N/A";
3188
3189 PetscFunctionBeginUser;
3190 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
3191
3192 PetscInt totalParticles;
3193 ierr = DMSwarmGetSize(user->swarm, &totalParticles); CHKERRQ(ierr);
3194
3195 if (!rank) {
3196 FILE *f;
3197 char filen[PETSC_MAX_PATH_LEN + 64];
3198 ierr = PetscSNPrintf(filen, sizeof(filen), "%s/Particle_Metrics.log", simCtx->log_dir); CHKERRQ(ierr);
3199 f = fopen(filen, "a");
3200 if (!f) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Cannot open particle log file: %s", filen);
3201
3202 if (ftell(f) == 0) {
3203 PetscFPrintf(PETSC_COMM_SELF, f, "%-18s | %-10s | %-12s | %-10s | %-10s | %-10s | %-15s | %-10s | %-10s\n",
3204 "Stage", "Timestep", "Total Ptls", "Lost", "Lost Total", "Migrated", "Occupied Cells", "Imbalance", "Mig Passes");
3205 PetscFPrintf(PETSC_COMM_SELF, f, "-------------------------------------------------------------------------------------------------------------------------------------------\n");
3206 }
3207
3208 PetscFPrintf(PETSC_COMM_SELF, f, "%-18s | %-10d | %-12d | %-10d | %-10d | %-10d | %-15d | %-10.2f | %-10d\n",
3209 stage_label, (int)simCtx->step, (int)totalParticles, (int)simCtx->particlesLostLastStep,
3210 (int)simCtx->particlesLostCumulative, (int)simCtx->particlesMigratedLastStep, (int)simCtx->occupiedCellCount,
3211 (double)simCtx->particleLoadImbalance, (int)simCtx->migrationPassesLastStep);
3212 fclose(f);
3213 }
3214 PetscFunctionReturn(0);
3215}
SetAnalyticalScalarFieldAtCellCenters
PetscErrorCode SetAnalyticalScalarFieldAtCellCenters(UserCtx *user, Vec targetVec)
Writes the configured verification scalar profile at physical cell centers into a scalar Vec.
Definition AnalyticalSolutions.c:695
SetAnalyticalSolutionForParticles
PetscErrorCode SetAnalyticalSolutionForParticles(Vec tempVec, SimCtx *simCtx)
Applies the analytical solution to particle velocity vector.
Definition AnalyticalSolutions.c:564
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:1833
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:1330
LOG_PARTICLE_METRICS
PetscErrorCode LOG_PARTICLE_METRICS(UserCtx *user, const char *stageName)
Implementation of LOG_PARTICLE_METRICS().
Definition logging.c:3182
BCHandlerTypeToString
const char * BCHandlerTypeToString(BCHandlerType handler_type)
Internal helper implementation: BCHandlerTypeToString().
Definition logging.c:771
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:1838
FieldInitializationToString
const char * FieldInitializationToString(PetscInt FieldInitialization)
Implementation of FieldInitializationToString().
Definition logging.c:687
DualMonitorDestroy
PetscErrorCode DualMonitorDestroy(void **ctx)
Implementation of DualMonitorDestroy().
Definition logging.c:809
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:2719
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:1295
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:1830
SolutionConvergenceSafeRelative
static PetscReal SolutionConvergenceSafeRelative(PetscReal numerator, PetscReal denominator)
Forms a guarded relative metric for solution-convergence logging.
Definition logging.c:927
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:1839
ProfilingFinalize
PetscErrorCode ProfilingFinalize(SimCtx *simCtx)
Implementation of ProfilingFinalize().
Definition logging.c:2148
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:896
LOG_FIELD_MIN_MAX
PetscErrorCode LOG_FIELD_MIN_MAX(UserCtx *user, const char *fieldName)
Implementation of LOG_FIELD_MIN_MAX().
Definition logging.c:2301
gNumAllowed
static int gNumAllowed
Number of entries in the gAllowedFunctions array.
Definition logging.c:28
ProfiledFunction::total_time
double total_time
Definition logging.c:1828
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:2254
SolutionConvergenceModeToString
static const char * SolutionConvergenceModeToString(SolutionConvergenceMode mode)
Maps the internal solution-convergence mode enum to its log label.
Definition logging.c:1543
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:2449
_FindOrCreateEntry
static PetscErrorCode _FindOrCreateEntry(const char *func_name, PetscInt *idx)
Internal helper implementation: _FindOrCreateEntry().
Definition logging.c:1846
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:2037
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:1959
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:1920
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:751
CalculateAdvancedParticleMetrics
PetscErrorCode CalculateAdvancedParticleMetrics(UserCtx *user)
Internal helper implementation: CalculateAdvancedParticleMetrics().
Definition logging.c:3128
ParticleLocationStatusToString
const char * ParticleLocationStatusToString(ParticleLocationStatus level)
Implementation of ParticleLocationStatusToString().
Definition logging.c:1811
LOG_SCATTER_METRICS
PetscErrorCode LOG_SCATTER_METRICS(UserCtx *user)
Implementation of LOG_SCATTER_METRICS().
Definition logging.c:2795
_CompareProfiledFunctions
static int _CompareProfiledFunctions(const void *a, const void *b)
Internal helper implementation: _CompareProfiledFunctions().
Definition logging.c:2132
LOG_SOLUTION_CONVERGENCE
PetscErrorCode LOG_SOLUTION_CONVERGENCE(SimCtx *simCtx)
Implementation of LOG_SOLUTION_CONVERGENCE().
Definition logging.c:1560
DualKSPMonitor
PetscErrorCode DualKSPMonitor(KSP ksp, PetscInt it, PetscReal rnorm, void *ctx)
Implementation of DualKSPMonitor().
Definition logging.c:848
LOG_CONTINUITY_METRICS
PetscErrorCode LOG_CONTINUITY_METRICS(UserCtx *user)
Implementation of LOG_CONTINUITY_METRICS().
Definition logging.c:1752
ProfiledFunction::current_step_time
double current_step_time
Definition logging.c:1829
ProfiledFunction::current_step_call_count
long long current_step_call_count
Definition logging.c:1831
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:1410
LOG_SEARCH_METRICS
PetscErrorCode LOG_SEARCH_METRICS(UserCtx *user)
Implementation of LOG_SEARCH_METRICS().
Definition logging.c:2977
g_profiler_registry
static ProfiledFunction * g_profiler_registry
Definition logging.c:1837
ProfilingInitialize
PetscErrorCode ProfilingInitialize(SimCtx *simCtx)
Internal helper implementation: ProfilingInitialize().
Definition logging.c:1882
LESModelToString
const char * LESModelToString(LESModelType LESFlag)
Implementation of LESModelToString().
Definition logging.c:720
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:1064
ProfilingResetTimestepCounters
PetscErrorCode ProfilingResetTimestepCounters(void)
Implementation of ProfilingResetTimestepCounters().
Definition logging.c:1942
ResetSearchMetrics
PetscErrorCode ResetSearchMetrics(SimCtx *simCtx)
Implementation of ResetSearchMetrics().
Definition logging.c:2946
MomentumSolverTypeToString
const char * MomentumSolverTypeToString(MomentumSolverType SolverFlag)
Implementation of MomentumSolverTypeToString().
Definition logging.c:736
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:954
SOLUTION_CONVERGENCE_FLUID_THRESHOLD
#define SOLUTION_CONVERGENCE_FLUID_THRESHOLD
Definition logging.c:895
ProfiledFunction::start_time
double start_time
Definition logging.c:1832
ParticleInitializationToString
const char * ParticleInitializationToString(ParticleInitializationType ParticleInitialization)
Implementation of ParticleInitializationToString().
Definition logging.c:703
_ProfilingStart
void _ProfilingStart(const char *func_name)
Implementation of _ProfilingStart().
Definition logging.c:1906
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:1827
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:488
DYNAMIC_SMAGORINSKY
@ DYNAMIC_SMAGORINSKY
Definition variables.h:491
NO_LES_MODEL
@ NO_LES_MODEL
Definition variables.h:489
CONSTANT_SMAGORINSKY
@ CONSTANT_SMAGORINSKY
Definition variables.h:490
UserCtx::lDiffusivityGradient
Vec lDiffusivityGradient
Definition variables.h:860
UserCtx::lCent
Vec lCent
Definition variables.h:879
BCType
BCType
Defines the general mathematical/physical Category of a boundary.
Definition variables.h:251
INLET
@ INLET
Definition variables.h:258
INTERFACE
@ INTERFACE
Definition variables.h:253
FARFIELD
@ FARFIELD
Definition variables.h:259
OUTLET
@ OUTLET
Definition variables.h:257
PERIODIC
@ PERIODIC
Definition variables.h:260
WALL
@ WALL
Definition variables.h:254
SimCtx::continueMode
PetscBool continueMode
Definition variables.h:668
MGCtx::user
UserCtx * user
Definition variables.h:536
SimCtx::profilingFinalSummary
PetscBool profilingFinalSummary
Definition variables.h:794
SimCtx::rank
PetscMPIInt rank
Definition variables.h:654
SimCtx::profilingTimestepFile
char profilingTimestepFile[PETSC_MAX_PATH_LEN]
Definition variables.h:793
SearchMetricsState::searchLocatedCount
PetscInt64 searchLocatedCount
Definition variables.h:224
SimCtx::block_number
PetscInt block_number
Definition variables.h:726
SearchMetricsState::searchLostCount
PetscInt64 searchLostCount
Definition variables.h:225
UserCtx::solutionConvergencePeriodicPRef
Vec * solutionConvergencePeriodicPRef
Definition variables.h:866
UserCtx::simCtx
SimCtx * simCtx
Back-pointer to the master simulation context.
Definition variables.h:833
ParticleInitializationType
ParticleInitializationType
Enumerator to identify the particle initialization strategy.
Definition variables.h:516
PARTICLE_INIT_SURFACE_RANDOM
@ PARTICLE_INIT_SURFACE_RANDOM
Random placement on the inlet face.
Definition variables.h:517
PARTICLE_INIT_SURFACE_EDGES
@ PARTICLE_INIT_SURFACE_EDGES
Deterministic placement at inlet face edges.
Definition variables.h:520
PARTICLE_INIT_POINT_SOURCE
@ PARTICLE_INIT_POINT_SOURCE
All particles at a fixed (psrc_x,psrc_y,psrc_z) — for validation.
Definition variables.h:519
PARTICLE_INIT_VOLUME
@ PARTICLE_INIT_VOLUME
Random volumetric distribution across the domain.
Definition variables.h:518
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:711
SimCtx::FluxOutSum
PetscReal FluxOutSum
Definition variables.h:735
UserCtx::Centz
Vec Centz
Definition variables.h:880
SimCtx::runtimeMemoryLogEnabled
PetscBool runtimeMemoryLogEnabled
Enable the rank-reduced runtime memory log.
Definition variables.h:809
SearchMetricsState::boundaryClampCount
PetscInt64 boundaryClampCount
Definition variables.h:231
SimCtx::particlesLostLastStep
PetscInt particlesLostLastStep
Definition variables.h:760
UserCtx::KM
PetscInt KM
Definition variables.h:839
UserCtx::lZet
Vec lZet
Definition variables.h:879
SimCtx::usermg
UserMG usermg
Definition variables.h:778
UserCtx::solutionConvergencePeriodicUcatRef
Vec * solutionConvergencePeriodicUcatRef
Definition variables.h:865
SearchMetricsState::traversalStepsSum
PetscInt64 traversalStepsSum
Definition variables.h:226
BCHandlerType
BCHandlerType
Defines the specific computational "strategy" for a boundary handler.
Definition variables.h:271
BC_HANDLER_INLET_PULSATILE_FLUX
@ BC_HANDLER_INLET_PULSATILE_FLUX
Definition variables.h:280
BC_HANDLER_PERIODIC_GEOMETRIC
@ BC_HANDLER_PERIODIC_GEOMETRIC
Definition variables.h:284
BC_HANDLER_INLET_PARABOLIC
@ BC_HANDLER_INLET_PARABOLIC
Definition variables.h:277
BC_HANDLER_INLET_CONSTANT_VELOCITY
@ BC_HANDLER_INLET_CONSTANT_VELOCITY
Definition variables.h:276
BC_HANDLER_PERIODIC_DRIVEN_INITIAL_FLUX
@ BC_HANDLER_PERIODIC_DRIVEN_INITIAL_FLUX
Definition variables.h:287
BC_HANDLER_INTERFACE_OVERSET
@ BC_HANDLER_INTERFACE_OVERSET
Definition variables.h:285
BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX
@ BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX
Definition variables.h:286
BC_HANDLER_WALL_MOVING
@ BC_HANDLER_WALL_MOVING
Definition variables.h:274
BC_HANDLER_INLET_PROFILE_FROM_FILE
@ BC_HANDLER_INLET_PROFILE_FROM_FILE
Definition variables.h:278
BC_HANDLER_WALL_NOSLIP
@ BC_HANDLER_WALL_NOSLIP
Definition variables.h:273
BC_HANDLER_OUTLET_CONSERVATION
@ BC_HANDLER_OUTLET_CONSERVATION
Definition variables.h:282
BC_HANDLER_FARFIELD_NONREFLECTING
@ BC_HANDLER_FARFIELD_NONREFLECTING
Definition variables.h:281
BC_HANDLER_OUTLET_PRESSURE
@ BC_HANDLER_OUTLET_PRESSURE
Definition variables.h:283
BC_HANDLER_SYMMETRY_PLANE
@ BC_HANDLER_SYMMETRY_PLANE
Definition variables.h:275
BC_HANDLER_UNDEFINED
@ BC_HANDLER_UNDEFINED
Definition variables.h:272
SimCtx::solutionConvergenceSamplesRecorded
PetscInt solutionConvergenceSamplesRecorded
Definition variables.h:710
UserCtx::lCellFieldAtCorner
Vec lCellFieldAtCorner
Definition variables.h:867
UserCtx::_this
PetscInt _this
Definition variables.h:843
SearchMetricsState::searchPopulation
PetscInt64 searchPopulation
Definition variables.h:223
SimCtx::solutionConvergenceMeanKEHistory
PetscReal * solutionConvergenceMeanKEHistory
Definition variables.h:712
SimCtx::dt
PetscReal dt
Definition variables.h:666
SimCtx::runtimeMemoryLogFile
char runtimeMemoryLogFile[PETSC_MAX_PATH_LEN]
File name written under log_dir.
Definition variables.h:810
SimCtx::runtimeMemoryLogStarted
PetscBool runtimeMemoryLogStarted
True after rank 0 writes the log header.
Definition variables.h:811
SimCtx::occupiedCellCount
PetscInt occupiedCellCount
Definition variables.h:764
SimCtx::profilingTimestepMode
char profilingTimestepMode[32]
Definition variables.h:792
UserCtx::lPsi
Vec lPsi
Definition variables.h:904
SearchMetricsState::currentSettlementPass
PetscInt currentSettlementPass
Definition variables.h:235
SimCtx::np
PetscInt np
Definition variables.h:753
UserCtx::DiffusivityGradient
Vec DiffusivityGradient
Definition variables.h:860
SimCtx::StartStep
PetscInt StartStep
Definition variables.h:661
MomentumSolverType
MomentumSolverType
Enumerator to identify the implemented momentum solver strategies.
Definition variables.h:502
MOMENTUM_SOLVER_EXPLICIT_RK
@ MOMENTUM_SOLVER_EXPLICIT_RK
Definition variables.h:503
MOMENTUM_SOLVER_DUALTIME_PICARD_RK4
@ MOMENTUM_SOLVER_DUALTIME_PICARD_RK4
Definition variables.h:504
SimCtx::solutionConvergencePeriodSteps
PetscInt solutionConvergencePeriodSteps
Definition variables.h:708
Cmpnts::x
PetscScalar x
Definition variables.h:101
SearchMetricsState::reSearchCount
PetscInt64 reSearchCount
Definition variables.h:227
SimCtx::MaxDiv
PetscReal MaxDiv
Definition variables.h:785
UserCtx::Centx
Vec Centx
Definition variables.h:880
SearchMetricsState::bboxGuessFallbackCount
PetscInt64 bboxGuessFallbackCount
Definition variables.h:233
UserCtx::Ucat_o
Vec Ucat_o
Definition variables.h:863
SimCtx::MaxDivx
PetscInt MaxDivx
Definition variables.h:786
SimCtx::MaxDivy
PetscInt MaxDivy
Definition variables.h:786
SearchMetricsState::bboxGuessSuccessCount
PetscInt64 bboxGuessSuccessCount
Definition variables.h:232
SimCtx::MaxDivz
PetscInt MaxDivz
Definition variables.h:786
SimCtx::log_dir
char log_dir[PETSC_MAX_PATH_LEN]
Definition variables.h:676
SimCtx::MaxDivFlatArg
PetscInt MaxDivFlatArg
Definition variables.h:786
SimCtx::FluxInSum
PetscReal FluxInSum
Definition variables.h:735
SearchMetricsState::maxParticlePassDepth
PetscInt64 maxParticlePassDepth
Definition variables.h:234
UserCtx::lCsi
Vec lCsi
Definition variables.h:879
SearchMetricsState::maxTraversalSteps
PetscInt64 maxTraversalSteps
Definition variables.h:228
Cmpnts::z
PetscScalar z
Definition variables.h:101
UserCtx::Ucat
Vec Ucat
Definition variables.h:856
UserCtx::ParticleCount
Vec ParticleCount
Definition variables.h:903
UserCtx::CellFieldAtCorner
Vec CellFieldAtCorner
Definition variables.h:867
UserCtx::JM
PetscInt JM
Definition variables.h:839
SimCtx::runtimeMemoryLogHasPrevious
PetscBool runtimeMemoryLogHasPrevious
True after the first process-memory sample.
Definition variables.h:812
UserMG::mglevels
PetscInt mglevels
Definition variables.h:543
SimCtx::profilingSelectedFuncs
char ** profilingSelectedFuncs
Definition variables.h:790
SimCtx::solutionConvergenceWindowSteps
PetscInt solutionConvergenceWindowSteps
Definition variables.h:709
SimCtx::particlesLostCumulative
PetscInt particlesLostCumulative
Definition variables.h:761
SimCtx::nProfilingSelectedFuncs
PetscInt nProfilingSelectedFuncs
Definition variables.h:791
SimCtx::particlesMigratedLastStep
PetscInt particlesMigratedLastStep
Definition variables.h:763
SimCtx::particleConsoleOutputFreq
PetscInt particleConsoleOutputFreq
Definition variables.h:664
SimCtx::searchMetrics
SearchMetricsState searchMetrics
Definition variables.h:766
UserCtx::lUcont
Vec lUcont
Definition variables.h:856
SimCtx::runtimeMemoryLogPreviousProcessMB
PetscReal runtimeMemoryLogPreviousProcessMB
Previous local process memory sample in MB.
Definition variables.h:813
SimCtx::step
PetscInt step
Definition variables.h:659
UserCtx::Diffusivity
Vec Diffusivity
Definition variables.h:859
UserCtx::info
DMDALocalInfo info
Definition variables.h:837
UserCtx::lUcat
Vec lUcat
Definition variables.h:856
SimCtx::migrationPassesLastStep
PetscInt migrationPassesLastStep
Definition variables.h:762
Cmpnts::y
PetscScalar y
Definition variables.h:101
EXEC_MODE_SOLVER
@ EXEC_MODE_SOLVER
Definition variables.h:624
EXEC_MODE_POSTPROCESSOR
@ EXEC_MODE_POSTPROCESSOR
Definition variables.h:625
UserCtx::IM
PetscInt IM
Definition variables.h:839
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:879
UserMG::mgctx
MGCtx * mgctx
Definition variables.h:546
SolutionConvergenceMode
SolutionConvergenceMode
Selects the runtime solution-convergence diagnostics mode.
Definition variables.h:508
SOLUTION_CONVERGENCE_TRANSIENT
@ SOLUTION_CONVERGENCE_TRANSIENT
Definition variables.h:512
SOLUTION_CONVERGENCE_PERIODIC_DETERMINISTIC
@ SOLUTION_CONVERGENCE_PERIODIC_DETERMINISTIC
Definition variables.h:510
SOLUTION_CONVERGENCE_STATISTICAL_STEADY
@ SOLUTION_CONVERGENCE_STATISTICAL_STEADY
Definition variables.h:511
SOLUTION_CONVERGENCE_STEADY_DETERMINISTIC
@ SOLUTION_CONVERGENCE_STEADY_DETERMINISTIC
Definition variables.h:509
UserCtx::lDiffusivity
Vec lDiffusivity
Definition variables.h:859
UserCtx::Centy
Vec Centy
Definition variables.h:880
SimCtx::solutionConvergenceMode
SolutionConvergenceMode solutionConvergenceMode
Definition variables.h:707
SearchMetricsState::searchAttempts
PetscInt64 searchAttempts
Definition variables.h:222
SimCtx::exec_mode
ExecutionMode exec_mode
Definition variables.h:670
SearchMetricsState::tieBreakCount
PetscInt64 tieBreakCount
Definition variables.h:230
SimCtx::ti
PetscReal ti
Definition variables.h:660
SimCtx::summationRHS
PetscReal summationRHS
Definition variables.h:784
SearchMetricsState::maxTraversalFailCount
PetscInt64 maxTraversalFailCount
Definition variables.h:229
SimCtx::LoggingFrequency
PetscInt LoggingFrequency
Definition variables.h:783
Cell::vertices
Cmpnts vertices[8]
Coordinates of the eight vertices of the cell.
Definition variables.h:161
UserCtx::Psi
Vec Psi
Definition variables.h:904
SimCtx::particleLoadImbalance
PetscReal particleLoadImbalance
Definition variables.h:765
UserCtx::P_o
Vec P_o
Definition variables.h:863
BCFace
BCFace
Identifies the six logical faces of a structured computational block.
Definition variables.h:244
BC_FACE_NEG_X
@ BC_FACE_NEG_X
Definition variables.h:245
BC_FACE_POS_Z
@ BC_FACE_POS_Z
Definition variables.h:247
BC_FACE_POS_Y
@ BC_FACE_POS_Y
Definition variables.h:246
BC_FACE_NEG_Z
@ BC_FACE_NEG_Z
Definition variables.h:247
BC_FACE_POS_X
@ BC_FACE_POS_X
Definition variables.h:245
BC_FACE_NEG_Y
@ BC_FACE_NEG_Y
Definition variables.h:246
Cell
Defines the vertices of a single hexahedral grid cell.
Definition variables.h:160
Cmpnts
A 3D point or vector with PetscScalar components.
Definition variables.h:100
SimCtx
The master context for the entire simulation.
Definition variables.h:651
UserCtx
User-defined context containing data specific to a single computational grid level.
Definition variables.h:830
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