PICurv 0.1.0
A Parallel Particle-In-Cell Solver for Curvilinear LES
Loading...
Searching...
No Matches
setup.c
Go to the documentation of this file.
1/**
2 * @file setup.c // Setup code for running any simulation
3 * @brief Test program for DMSwarm interpolation using the fdf-curvIB method.
4 * Provides the setup to start any simulation with DMSwarm and DMDAs.
5 **/
6
7#include <ctype.h>
8#include <errno.h>
9
10#include "setup.h"
11#include "statistics_config.h"
12#include "statistics_accumulator.h"
13
14/**
15 * @brief Implementation of \ref RuntimeWalltimeGuardParsePositiveSeconds().
16 * @details Full API contract (arguments, ownership, side effects) is documented with
17 * the header declaration in `include/setup.h`.
18 * @see RuntimeWalltimeGuardParsePositiveSeconds()
19 */
20PetscBool RuntimeWalltimeGuardParsePositiveSeconds(const char *text, PetscReal *seconds_out)
21{
22 char *endptr = NULL;
23 double parsed_value;
24
25 if (seconds_out) *seconds_out = 0.0;
26 if (!text || text[0] == '\0') return PETSC_FALSE;
27
28 errno = 0;
29 parsed_value = strtod(text, &endptr);
30 if (endptr == text || errno == ERANGE || !isfinite(parsed_value) || parsed_value <= 0.0) {
31 return PETSC_FALSE;
32 }
33
34 while (*endptr != '\0' && isspace((unsigned char)*endptr)) {
35 endptr++;
36 }
37 if (*endptr != '\0') return PETSC_FALSE;
38
39 if (seconds_out) *seconds_out = (PetscReal)parsed_value;
40 return PETSC_TRUE;
41}
42
43/**
44 * @brief Implementation of \ref InitializeSolutionConvergenceState().
45 * @details Full API contract (arguments, ownership, side effects) is documented with
46 * the header declaration in `include/setup.h`.
47 * @see InitializeSolutionConvergenceState()
48 */
49PetscErrorCode InitializeSolutionConvergenceState(SimCtx *simCtx)
50{
51 UserCtx *user = NULL;
52 PetscInt history_capacity = 0;
53
54 PetscFunctionBeginUser;
55 if (!simCtx) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "SimCtx cannot be NULL.");
56 if (simCtx->exec_mode != EXEC_MODE_SOLVER) PetscFunctionReturn(0);
57 if (!simCtx->solutionConvergenceEnabled) PetscFunctionReturn(0);
58 if (!simCtx->usermg.mgctx) {
59 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE,
60 "Multigrid hierarchy must exist before initializing solution convergence storage.");
61 }
62
63 user = simCtx->usermg.mgctx[simCtx->usermg.mglevels - 1].user;
64 if (!user) {
65 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE,
66 "Finest-level UserCtx must exist before initializing solution convergence storage.");
67 }
68
69 simCtx->solutionConvergenceSamplesRecorded = 0;
70
71 if (simCtx->solutionConvergenceMode == SOLUTION_CONVERGENCE_PERIODIC_DETERMINISTIC) {
72 for (PetscInt bi = 0; bi < simCtx->block_number; ++bi) {
73 PetscCall(PetscCalloc1((size_t)simCtx->solutionConvergencePeriodSteps,
74 &user[bi].solutionConvergencePeriodicUcatRef));
75 PetscCall(PetscCalloc1((size_t)simCtx->solutionConvergencePeriodSteps,
76 &user[bi].solutionConvergencePeriodicPRef));
77 for (PetscInt phase = 0; phase < simCtx->solutionConvergencePeriodSteps; ++phase) {
78 PetscCall(VecDuplicate(user[bi].Ucat, &user[bi].solutionConvergencePeriodicUcatRef[phase]));
79 PetscCall(VecSet(user[bi].solutionConvergencePeriodicUcatRef[phase], 0.0));
80 PetscCall(VecDuplicate(user[bi].P, &user[bi].solutionConvergencePeriodicPRef[phase]));
81 PetscCall(VecSet(user[bi].solutionConvergencePeriodicPRef[phase], 0.0));
82 }
83 }
84 }
85
86 if (simCtx->solutionConvergenceMode == SOLUTION_CONVERGENCE_STATISTICAL_STEADY) {
87 history_capacity = 2 * simCtx->solutionConvergenceWindowSteps;
88 PetscCall(PetscCalloc1((size_t)history_capacity, &simCtx->solutionConvergenceMeanSpeedHistory));
89 PetscCall(PetscCalloc1((size_t)history_capacity, &simCtx->solutionConvergenceMeanKEHistory));
90 }
91
92 PetscFunctionReturn(0);
93}
94
95/**
96 * @brief Implementation of \ref DestroySolutionConvergenceState().
97 * @details Full API contract (arguments, ownership, side effects) is documented with
98 * the header declaration in `include/setup.h`.
99 * @see DestroySolutionConvergenceState()
100 */
101PetscErrorCode DestroySolutionConvergenceState(SimCtx *simCtx)
102{
103 UserCtx *user = NULL;
104
105 PetscFunctionBeginUser;
106 if (!simCtx) PetscFunctionReturn(0);
107 if (simCtx->fieldStatisticsWindows) {
108 PetscCall(PetscFree(simCtx->fieldStatisticsWindows));
109 simCtx->fieldStatisticsWindows = NULL;
110 simCtx->fieldStatisticsWindowCount = 0;
111 }
112 if (!simCtx->usermg.mgctx) PetscFunctionReturn(0);
113
114 user = simCtx->usermg.mgctx[simCtx->usermg.mglevels - 1].user;
115 if (user) {
116 for (PetscInt bi = 0; bi < simCtx->block_number; ++bi) {
117 if (user[bi].solutionConvergencePeriodicUcatRef) {
118 for (PetscInt phase = 0; phase < simCtx->solutionConvergencePeriodSteps; ++phase) {
119 if (user[bi].solutionConvergencePeriodicUcatRef[phase]) {
120 PetscCall(VecDestroy(&user[bi].solutionConvergencePeriodicUcatRef[phase]));
121 }
122 }
123 PetscCall(PetscFree(user[bi].solutionConvergencePeriodicUcatRef));
124 user[bi].solutionConvergencePeriodicUcatRef = NULL;
125 }
126 if (user[bi].solutionConvergencePeriodicPRef) {
127 for (PetscInt phase = 0; phase < simCtx->solutionConvergencePeriodSteps; ++phase) {
128 if (user[bi].solutionConvergencePeriodicPRef[phase]) {
129 PetscCall(VecDestroy(&user[bi].solutionConvergencePeriodicPRef[phase]));
130 }
131 }
132 PetscCall(PetscFree(user[bi].solutionConvergencePeriodicPRef));
133 user[bi].solutionConvergencePeriodicPRef = NULL;
134 }
135 }
136 }
137
138 if (simCtx->solutionConvergenceMeanSpeedHistory) {
139 PetscCall(PetscFree(simCtx->solutionConvergenceMeanSpeedHistory));
140 simCtx->solutionConvergenceMeanSpeedHistory = NULL;
141 }
142 if (simCtx->solutionConvergenceMeanKEHistory) {
143 PetscCall(PetscFree(simCtx->solutionConvergenceMeanKEHistory));
144 simCtx->solutionConvergenceMeanKEHistory = NULL;
145 }
146 simCtx->solutionConvergenceSamplesRecorded = 0;
147
148 PetscFunctionReturn(0);
149}
150
151#undef __FUNCT__
152#define __FUNCT__ "CreateSimulationContext"
153
154/**
155 * @brief Implementation of \ref CreateSimulationContext().
156 * @details Full API contract (arguments, ownership, side effects) is documented with
157 * the header declaration in `include/setup.h`.
158 * @see CreateSimulationContext()
159 */
160PetscErrorCode CreateSimulationContext(int argc, char **argv, SimCtx **p_simCtx)
161{
162 PetscErrorCode ierr;
163 (void)argc;
164 (void)argv;
165 SimCtx *simCtx;
166 char control_filename[PETSC_MAX_PATH_LEN] = ""; // Temporary placeholder for control file name.
167 PetscBool control_flg; // Temporary placeholder for control file tag existence check flag.
168 PetscBool particle_console_output_freq_flg = PETSC_FALSE;
169
170 PetscFunctionBeginUser;
171
172 PROFILE_FUNCTION_BEGIN;
173
174 // === 1. Allocate the Context Struct and Set ALL Defaults ==================
175 ierr = PetscNew(p_simCtx); CHKERRQ(ierr);
176 simCtx = *p_simCtx;
177
178 // --- Group 1: Parallelism & MPI ---
179 simCtx->rank = 0; simCtx->size = 1;
180
181 // --- Group 2: Simulation Control, Time, and I/O ---
182 simCtx->step = 0; simCtx->ti = 0.0; simCtx->StartStep = 0; simCtx->StepsToRun = 10;
183 simCtx->tiout = 10; simCtx->particleConsoleOutputFreq = simCtx->tiout;
184 simCtx->StartTime = 0.0; simCtx->dt = 0.001;
185 simCtx->OnlySetup = PETSC_FALSE;
186 simCtx->continueMode = PETSC_FALSE;
187 simCtx->logviewer = NULL;
188 strcpy(simCtx->eulerianSource,"solve");
189 strcpy(simCtx->restart_dir,"restart");
190 strcpy(simCtx->output_dir,"output");
191 strcpy(simCtx->log_dir,"logs");
192 simCtx->_io_context_buffer[0] = '\0';
193 simCtx->current_io_directory = NULL;
194 simCtx->checkpointGeometrySHA256[0] = '\0';
195 simCtx->checkpointGeometryHashReady = PETSC_FALSE;
196 simCtx->restartHistoryAvailable = PETSC_FALSE;
197
198 // --- Group 3: High-Level Physics & Model Selection Flags ---
199 simCtx->immersed = 0; simCtx->movefsi = 0; simCtx->rotatefsi = 0;
200 simCtx->sediment = 0; simCtx->rheology = 0; simCtx->invicid = 0;
201 simCtx->TwoD = 0; simCtx->thin = 0; simCtx->moveframe = 0;
202 simCtx->rotateframe = 0; simCtx->blank = 0;
203 simCtx->dgf_x = 0; simCtx->dgf_y = 1; simCtx->dgf_z = 0;
204 simCtx->dgf_ax = 1; simCtx->dgf_ay = 0; simCtx->dgf_az = 0;
205 strcpy(simCtx->AnalyticalSolutionType,"TGV3D");
206
207 // --- Group 4: Specific Simulation Case Flags --- (DEPRICATED)
208 simCtx->cop=0; simCtx->fish=0; simCtx->fish_c=0; simCtx->fishcyl=0;
209 simCtx->eel=0; simCtx->pizza=0; simCtx->turbine=0; simCtx->Pipe=0;
210 simCtx->wing=0; simCtx->hydro=0; simCtx->MHV=0; simCtx->LV=0;
211 simCtx->channelz = 0;
212
213 // --- Group 5: Solver & Numerics Parameters ---
214 simCtx->mom_solver_type = MOMENTUM_SOLVER_DUALTIME_PICARD_JAMESON_RK; simCtx->mom_max_pseudo_steps = 50;
215 simCtx->mom_dt_jameson_residual_norm_noise_allowance_factor = 1.1; // raised from 1.05; less aggressive rejection
216 simCtx->mom_atol = 1e-7; simCtx->mom_rtol = 1e-4;
217 /* Residual-based convergence is the default. Leaving both at 0.0 selected the
218 * update-only branch, where |dU| <= mom_atol can pass purely because dtau collapsed
219 * (|dU| ~ dtau*|R|), converging on a state that does not satisfy the equations.
220 * Shipped configs all set these explicitly; these defaults protect a minimal
221 * user-written config. Set both non-positive to opt back out deliberately. */
222 simCtx->mom_resid_atol = 1e-8; simCtx->mom_resid_rtol = 1e-3;
223 simCtx->imp_stol = 1.e-8;
224 simCtx->mglevels = 3; simCtx->mg_MAX_IT = 30; simCtx->mg_idx = 1;
225 simCtx->mg_preItr = 1; simCtx->mg_poItr = 1;
226 simCtx->poisson = 0; simCtx->poisson_tol = 5.e-9;
227 simCtx->STRONG_COUPLING = 0;simCtx->central=0;
228 /* pseudo_cfl and its bounds are now dimensionless Courant numbers: CFL = dtau * lambda_max,
229 where lambda_max is the global spectral radius computed at each physical timestep.
230 Stable range for 4-stage Jameson RK: ~0–2.83. Initial 0.5 gives a comfortable margin. */
231 simCtx->ren = 100.0; simCtx->pseudo_cfl = 0.5;
232 simCtx->max_pseudo_cfl = 2.0; simCtx->min_pseudo_cfl = 0.001;
233 simCtx->pseudo_cfl_reduction_factor = 0.75;
234 simCtx->pseudo_cfl_growth_factor = 1.1; // raised from 1.0; controller can now increase CFL
235 simCtx->no_pseudo_cfl_backtrack = PETSC_FALSE;
236 simCtx->mom_ratio_ema_alpha = 0.3; /* moderate smoothing; set to 1.0 to recover original raw-ratio behavior */
237 simCtx->mom_last_converged = PETSC_TRUE;
238 simCtx->mom_last_lambda_max = 0.0; /* populated after first momentum solve */
239 simCtx->mom_nk_monitor_history = PETSC_FALSE;
240 simCtx->ps_ksp_pic_monitor_true_residual = PETSC_FALSE;
241 simCtx->cdisx = 0.0; simCtx->cdisy = 0.0; simCtx->cdisz = 0.0;
242 simCtx->initialConditionMode = IC_MODE_ZERO;
243 simCtx->initialConditionField = IC_FIELD_UCAT;
244 strcpy(simCtx->initialConditionDirectory, "config/initial_condition");
245 simCtx->InitialConstantContra.x = 0.0;
246 simCtx->InitialConstantContra.y = 0.0;
247 simCtx->InitialConstantContra.z = 0.0;
248 simCtx->flowDirection = FLOW_DIR_UNSET;
249 simCtx->icVelocityPhysical = 0.0;
250 simCtx->AnalyticalUniformVelocity.x = 0.0;
251 simCtx->AnalyticalUniformVelocity.y = 0.0;
252 simCtx->AnalyticalUniformVelocity.z = 0.0;
253 simCtx->fieldStatisticsEnabled = PETSC_FALSE;
254 simCtx->fieldStatisticsWindowCount = 0;
255 simCtx->fieldStatisticsWindows = NULL;
256 simCtx->statisticsConsoleOutputFreq = 0;
257 simCtx->fieldStatisticsContinue = PETSC_FALSE;
258 simCtx->solutionConvergenceEnabled = PETSC_TRUE;
259 simCtx->solutionConvergenceMode = SOLUTION_CONVERGENCE_STEADY_DETERMINISTIC;
260 simCtx->solutionConvergencePeriodSteps = 0;
261 simCtx->solutionConvergenceWindowSteps = 0;
262 simCtx->solutionConvergenceSamplesRecorded = 0;
263 simCtx->solutionConvergenceMeanSpeedHistory = NULL;
264 simCtx->solutionConvergenceMeanKEHistory = NULL;
265 simCtx->verificationDiffusivity.enabled = PETSC_FALSE;
266 strcpy(simCtx->verificationDiffusivity.mode, "");
267 strcpy(simCtx->verificationDiffusivity.profile, "");
268 simCtx->verificationDiffusivity.gamma0 = 0.0;
269 simCtx->verificationDiffusivity.slope_x = 0.0;
270
271 // --- Group 6: Physical & Geometric Parameters ---
272 simCtx->NumberOfBodies = 1; simCtx->Flux_in = 1.0; simCtx->angle = 0.0;
273 simCtx->max_angle = -54. * 3.1415926 / 180.;
274 simCtx->CMx_c=0.0; simCtx->CMy_c=0.0; simCtx->CMz_c=0.0;
275 simCtx->wall_roughness_height = 1e-16;
276 simCtx->schmidt_number = 1.0; simCtx->Turbulent_schmidt_number = 0.7;
277
278 // --- Group 7: Grid, Domain, and Boundary Condition Settings ---
279 simCtx->block_number = 1; simCtx->inletprofile = 1;
280 simCtx->grid1d = 0; simCtx->Ogrid = 0;
281 simCtx->i_periodic = 0; simCtx->j_periodic = 0; simCtx->k_periodic = 0;
282 simCtx->blkpbc = 10; simCtx->pseudo_periodic = 0;
283 strcpy(simCtx->grid_file, "config/grid.run");
284 simCtx->generate_grid = PETSC_FALSE;
285 simCtx->da_procs_x = PETSC_DECIDE;
286 simCtx->da_procs_y = PETSC_DECIDE;
287 simCtx->da_procs_z = PETSC_DECIDE;
288 simCtx->grid_rotation_angle = 0.0;
289 simCtx->Croty = 0.0; simCtx->Crotz = 0.0;
290 simCtx->num_bcs_files = 1;
291 ierr = PetscMalloc1(1, &simCtx->bcs_files); CHKERRQ(ierr);
292 ierr = PetscStrallocpy("config/bcs.run", &simCtx->bcs_files[0]); CHKERRQ(ierr);
293 simCtx->FluxInSum = 0.0; simCtx->FluxOutSum = 0.0; simCtx->Fluxsum = 0.0;
294 simCtx->drivingForceMagnitude = 0.0, simCtx->forceScalingFactor = 1.8;
295 simCtx->drivingForceStep = -1;
296 simCtx->targetVolumetricFlux = 0.0;
297 simCtx->drivenFluxTargetLatched = PETSC_FALSE;
298 simCtx->bulkVelocityCorrection = 0.0;
299 simCtx->boundaryVelocityCorrection = 0.0;
300 simCtx->AreaInSum = 0.0; simCtx->AreaOutSum = 0.0;
301 simCtx->U_bc = 0.0; simCtx->ccc = 0;
302 simCtx->ratio = 0.0;
303
304
305 // --- Group 8: Turbulence Modeling (LES/RANS) ---
306 simCtx->les = NO_LES_MODEL; simCtx->rans = 0;
307 simCtx->wallfunction = 0; simCtx->mixed = 0; simCtx->clark = 0;
308 simCtx->dynamic_freq = 1; simCtx->max_cs = 0.5;
309 simCtx->Const_CS = 0.03;
310 simCtx->testfilter_ik = 0; simCtx->testfilter_1d = 0;
311 simCtx->i_homo_filter = 0; simCtx->j_homo_filter = 0; simCtx->k_homo_filter = 0;
312
313 // --- Group 9: Particle / DMSwarm Data & Settings ---
314 simCtx->np = 0; simCtx->readFields = PETSC_FALSE;
315 simCtx->dm_swarm = NULL; simCtx->bboxlist = NULL;
316 simCtx->ParticleInitialization = PARTICLE_INIT_SURFACE_RANDOM;
317 simCtx->interpolationMethod = INTERP_TRILINEAR;
318 strcpy(simCtx->particleRestartMode,"load");
319 simCtx->particlesLostLastStep = 0;
320 simCtx->particlesLostCumulative = 0;
321 simCtx->particlesMigratedLastStep = 0;
322 simCtx->occupiedCellCount = 0;
323 simCtx->particleLoadImbalance = 0.0;
324 simCtx->migrationPassesLastStep = 0;
325 simCtx->searchMetrics.searchAttempts = 0;
326 simCtx->searchMetrics.searchPopulation = 0;
327 simCtx->searchMetrics.searchLocatedCount = 0;
328 simCtx->searchMetrics.searchLostCount = 0;
329 simCtx->searchMetrics.traversalStepsSum = 0;
330 simCtx->searchMetrics.reSearchCount = 0;
331 simCtx->searchMetrics.maxTraversalSteps = 0;
332 simCtx->searchMetrics.maxTraversalFailCount = 0;
333 simCtx->searchMetrics.tieBreakCount = 0;
334 simCtx->searchMetrics.boundaryClampCount = 0;
335 simCtx->searchMetrics.bboxGuessSuccessCount = 0;
336 simCtx->searchMetrics.bboxGuessFallbackCount = 0;
337 simCtx->searchMetrics.maxParticlePassDepth = 0;
338 simCtx->searchMetrics.currentSettlementPass = 0;
339 simCtx->BrownianMotionRNG = NULL;
340 simCtx->C_IEM = 2.0;
341
342 // --- Group 10: Immersed Boundary & FSI Data Object Pointers ---
343 simCtx->ibm = NULL; simCtx->ibmv = NULL; simCtx->fsi = NULL;
344 simCtx->rstart_fsi = PETSC_FALSE; simCtx->duplicate = 0;
345
346 // --- Group 11: Logging and Custom Configuration ---
347 strcpy(simCtx->allowedFile, "config/whitelist.run");
348 simCtx->useCfg = PETSC_FALSE;
349 simCtx->allowedFuncs = NULL;
350 simCtx->nAllowed = 0;
351 simCtx->LoggingFrequency = 10;
352 simCtx->summationRHS = 0.0;
353 simCtx->MaxDiv = 0.0;
354 simCtx->MaxDivFlatArg = 0; simCtx->MaxDivx = 0; simCtx->MaxDivy = 0; simCtx->MaxDivz = 0;
355 strcpy(simCtx->profilingSelectedFuncsFile, "config/profile.run");
356 simCtx->useProfilingSelectedFuncsCfg = PETSC_FALSE;
357 simCtx->profilingSelectedFuncs = NULL;
358 simCtx->nProfilingSelectedFuncs = 0;
359 strcpy(simCtx->profilingTimestepMode, "selected");
360 strcpy(simCtx->profilingTimestepFile, "Profiling_Timestep_Summary.csv");
361 simCtx->profilingFinalSummary = PETSC_TRUE;
362 simCtx->walltimeGuardEnabled = PETSC_FALSE;
363 simCtx->walltimeGuardActive = PETSC_FALSE;
364 simCtx->walltimeGuardWarmupSteps = 10;
365 simCtx->walltimeGuardMultiplier = 2.0;
366 simCtx->walltimeGuardMinSeconds = 60.0;
367 simCtx->walltimeGuardEstimatorAlpha = 0.35;
368 simCtx->walltimeGuardJobStartEpochSeconds = 0.0;
369 simCtx->walltimeGuardLimitSeconds = 0.0;
370 simCtx->walltimeGuardCompletedSteps = 0;
371 simCtx->walltimeGuardWarmupTotalSeconds = 0.0;
372 simCtx->walltimeGuardWarmupAverageSeconds = 0.0;
373 simCtx->walltimeGuardHasEWMA = PETSC_FALSE;
374 simCtx->walltimeGuardEWMASeconds = 0.0;
375 simCtx->walltimeGuardLatestStepSeconds = 0.0;
376 simCtx->runtimeMemoryLogEnabled = PETSC_TRUE;
377 strcpy(simCtx->runtimeMemoryLogFile, "Runtime_Memory.log");
378 simCtx->runtimeMemoryLogStarted = PETSC_FALSE;
379 simCtx->runtimeMemoryLogHasPrevious = PETSC_FALSE;
380 simCtx->runtimeMemoryLogPreviousProcessMB = 0.0;
381 // --- Group 11: Post-Processing Information ---
382 strcpy(simCtx->PostprocessingControlFile, "config/post.run");
383 ierr = PetscNew(&simCtx->pps); CHKERRQ(ierr);
384
385 // === 2. Get MPI Info and Handle Config File =============================
386 // -- Group 1: Parallelism & MPI Information
387 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &simCtx->rank); CHKERRQ(ierr);
388 ierr = MPI_Comm_size(PETSC_COMM_WORLD, &simCtx->size); CHKERRQ(ierr);
389
390 // First, check if the -control_file argument was provided by the user/script.
391 ierr = PetscOptionsGetString(NULL, NULL, "-control_file", control_filename, sizeof(control_filename), &control_flg); CHKERRQ(ierr);
392
393 // If the flag is NOT present or the filename is empty, abort with a helpful error.
394 if (!control_flg || strlen(control_filename) == 0) {
395 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
396 "\n\n*** MANDATORY ARGUMENT MISSING ***\n"
397 "The -control_file argument was not provided.\n"
398 "This program must be launched with a configuration file.\n"
399 "Example: mpiexec -n 4 ./simulator -control_file /path/to/your/config.control\n"
400 "This is typically handled automatically by the 'picurv' script.\n");
401 }
402
403 // At this point, we have a valid filename. Attempt to load it.
404 LOG(GLOBAL, LOG_INFO, "Loading mandatory configuration from: %s\n", control_filename);
405 ierr = PetscOptionsInsertFile(PETSC_COMM_WORLD, NULL, control_filename, PETSC_FALSE);
406 if (ierr == PETSC_ERR_FILE_OPEN) {
407 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_FILE_OPEN, "The specified control file was not found or could not be opened: %s", control_filename);
408 }
409 CHKERRQ(ierr);
410
411 {
412 PetscBool legacy_averaging = PETSC_FALSE;
413 ierr = PetscOptionsHasName(NULL, NULL, "-averaging", &legacy_averaging); CHKERRQ(ierr);
414 PetscCheck(!legacy_averaging, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
415 "Legacy -averaging was removed. Use instantaneous output and offline "
416 "postprocessing until the replacement field-statistics pipeline is available.");
417 }
418
419 // === 3. A Configure Logging System ========================================
420 // This logic determines the logging configuration and STORES it in simCtx for
421 // later reference and cleanup.
422 ierr = PetscOptionsGetString(NULL, NULL, "-whitelist_config_file", simCtx->allowedFile, PETSC_MAX_PATH_LEN, &simCtx->useCfg); CHKERRQ(ierr);
423
424 if (simCtx->useCfg) {
425 ierr = LoadAllowedFunctionsFromFile(simCtx->allowedFile, &simCtx->allowedFuncs, &simCtx->nAllowed);
426 if (ierr) {
427 // Use direct PetscPrintf as logging system isn't fully active yet.
428 PetscPrintf(PETSC_COMM_SELF, "[%s] WARNING: Failed to load allowed functions from '%s'. Falling back to default list.\n", __func__, simCtx->allowedFile);
429 simCtx->useCfg = PETSC_FALSE; // Mark as failed.
430 ierr = 0; // Clear the error to allow fallback.
431 } else if (simCtx->nAllowed == 0) {
432 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
433 "Whitelist config file '%s' is empty. Omit -whitelist_config_file to use the default allow-list, or list at least one function.",
434 simCtx->allowedFile);
435 }
436 }
437 if (!simCtx->useCfg) {
438 // Fallback to default logging functions if no file was used or if loading failed.
439 simCtx->nAllowed = 2;
440 ierr = PetscMalloc1(simCtx->nAllowed, &simCtx->allowedFuncs); CHKERRQ(ierr);
441 ierr = PetscStrallocpy("main", &simCtx->allowedFuncs[0]); CHKERRQ(ierr);
442 ierr = PetscStrallocpy("CreateSimulationContext", &simCtx->allowedFuncs[1]); CHKERRQ(ierr);
443 }
444
445 // Activate the configuration by passing it to the logging module's setup function.
446 set_allowed_functions((const char**)simCtx->allowedFuncs, (size_t)simCtx->nAllowed);
447
448 // Now that the logger is configured, we can use it.
449 LOG_ALLOW_SYNC(LOCAL, LOG_INFO, "Context created. Initializing on rank %d of %d.\n", simCtx->rank, simCtx->size);
450 print_log_level(); // This will now correctly reflect the LOG_LEVEL environment variable.
451
452 // === 3.B Configure Profiling System ========================================
453 ierr = PetscOptionsGetString(NULL, NULL, "-profiling_timestep_mode", simCtx->profilingTimestepMode, sizeof(simCtx->profilingTimestepMode), NULL); CHKERRQ(ierr);
454 ierr = PetscOptionsGetString(NULL, NULL, "-profiling_timestep_file", simCtx->profilingTimestepFile, PETSC_MAX_PATH_LEN, NULL); CHKERRQ(ierr);
455 ierr = PetscOptionsGetBool(NULL, NULL, "-profiling_final_summary", &simCtx->profilingFinalSummary, NULL); CHKERRQ(ierr);
456 if (strcmp(simCtx->profilingTimestepMode, "off") != 0 &&
457 strcmp(simCtx->profilingTimestepMode, "selected") != 0 &&
458 strcmp(simCtx->profilingTimestepMode, "all") != 0) {
459 PetscPrintf(PETSC_COMM_SELF, "[%s] WARNING: Unknown profiling timestep mode '%s'. Falling back to 'selected'.\n", __func__, simCtx->profilingTimestepMode);
460 strcpy(simCtx->profilingTimestepMode, "selected");
461 }
462
463 if (strcmp(simCtx->profilingTimestepMode, "selected") == 0) {
464 ierr = PetscOptionsGetString(NULL, NULL, "-profile_config_file", simCtx->profilingSelectedFuncsFile, PETSC_MAX_PATH_LEN, &simCtx->useProfilingSelectedFuncsCfg); CHKERRQ(ierr);
465 if (simCtx->useProfilingSelectedFuncsCfg) {
466 ierr = LoadAllowedFunctionsFromFile(simCtx->profilingSelectedFuncsFile, &simCtx->profilingSelectedFuncs, &simCtx->nProfilingSelectedFuncs);
467 if (ierr) {
468 PetscPrintf(PETSC_COMM_SELF, "[%s] WARNING: Failed to load selected profiling functions from '%s'. Falling back to default list.\n", __func__, simCtx->profilingSelectedFuncsFile);
469 simCtx->useProfilingSelectedFuncsCfg = PETSC_FALSE;
470 ierr = 0;
471 }
472 }
473 if (!simCtx->useProfilingSelectedFuncsCfg) {
474 // Fallback to a hardcoded default list if no file was provided or loading failed.
475 simCtx->nProfilingSelectedFuncs = 4;
476 ierr = PetscMalloc1(simCtx->nProfilingSelectedFuncs, &simCtx->profilingSelectedFuncs); CHKERRQ(ierr);
477 ierr = PetscStrallocpy("FlowSolver", &simCtx->profilingSelectedFuncs[0]); CHKERRQ(ierr);
478 ierr = PetscStrallocpy("AdvanceSimulation", &simCtx->profilingSelectedFuncs[1]); CHKERRQ(ierr);
479 ierr = PetscStrallocpy("LocateAllParticlesInGrid", &simCtx->profilingSelectedFuncs[2]); CHKERRQ(ierr);
480 ierr = PetscStrallocpy("InterpolateAllFieldsToSwarm", &simCtx->profilingSelectedFuncs[3]); CHKERRQ(ierr);
481 }
482 }
483
484 // Initialize the profiling system with the current updated simulation context.
485 ierr = ProfilingInitialize(simCtx); CHKERRQ(ierr);
486
487 // === 4. Parse All Command Line Options ==================================
488 LOG_ALLOW(GLOBAL, LOG_INFO, "Parsing command-line options...\n");
489
490 // --- Group 2
491 LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 2: Simulation Control,Time and I/O.\n");
492 // Read the physical time to start from.
493 // The default is already 0.0, so this will only be non-zero if the user provides it.
494 ierr = PetscOptionsGetInt(NULL, NULL, "-start_step", &simCtx->StartStep, NULL); CHKERRQ(ierr);
495 ierr = PetscOptionsGetInt(NULL,NULL, "-totalsteps", &simCtx->StepsToRun, NULL); CHKERRQ(ierr);
496 ierr = PetscOptionsGetBool(NULL, NULL, "-only_setup", &simCtx->OnlySetup, NULL); CHKERRQ(ierr);
497 ierr = PetscOptionsGetBool(NULL, NULL, "-continue_mode", &simCtx->continueMode, NULL); CHKERRQ(ierr);
498 /* Resuming in the same run directory continues accumulated statistics by
499 * default; branching with --restart-from must opt in, because the branch may
500 * follow a different physical trajectory than the samples already collected. */
501 simCtx->fieldStatisticsContinue = simCtx->continueMode;
502 ierr = PetscOptionsGetReal(NULL, NULL, "-dt", &simCtx->dt, NULL); CHKERRQ(ierr);
503 ierr = PetscOptionsGetInt(NULL, NULL, "-tio", &simCtx->tiout, NULL); CHKERRQ(ierr);
504 ierr = PetscOptionsGetInt(NULL, NULL, "-particle_console_output_freq", &simCtx->particleConsoleOutputFreq, &particle_console_output_freq_flg); CHKERRQ(ierr);
505 if (!particle_console_output_freq_flg) {
506 simCtx->particleConsoleOutputFreq = simCtx->tiout;
507 }
508 ierr = PetscOptionsGetString(NULL,NULL,"-euler_field_source",simCtx->eulerianSource,sizeof(simCtx->eulerianSource),NULL);CHKERRQ(ierr);
509 ierr = PetscOptionsGetString(NULL,NULL,"-output_dir",simCtx->output_dir,sizeof(simCtx->output_dir),NULL);CHKERRQ(ierr);
510 ierr = PetscOptionsGetString(NULL,NULL,"-restart_dir",simCtx->restart_dir,sizeof(simCtx->restart_dir),NULL);CHKERRQ(ierr);
511 ierr = PetscOptionsGetString(NULL,NULL,"-log_dir",simCtx->log_dir,sizeof(simCtx->log_dir),NULL);CHKERRQ(ierr);
512 ierr = PetscOptionsGetBool(NULL, NULL, "-walltime_guard_enabled", &simCtx->walltimeGuardEnabled, NULL); CHKERRQ(ierr);
513 ierr = PetscOptionsGetInt(NULL, NULL, "-walltime_guard_warmup_steps", &simCtx->walltimeGuardWarmupSteps, NULL); CHKERRQ(ierr);
514 ierr = PetscOptionsGetReal(NULL, NULL, "-walltime_guard_multiplier", &simCtx->walltimeGuardMultiplier, NULL); CHKERRQ(ierr);
515 ierr = PetscOptionsGetBool(NULL, NULL, "-runtime_memory_log_enabled", &simCtx->runtimeMemoryLogEnabled, NULL); CHKERRQ(ierr);
516 ierr = PetscOptionsGetString(NULL, NULL, "-runtime_memory_log_file", simCtx->runtimeMemoryLogFile, PETSC_MAX_PATH_LEN, NULL); CHKERRQ(ierr);
517 ierr = PetscOptionsGetReal(NULL, NULL, "-walltime_guard_min_seconds", &simCtx->walltimeGuardMinSeconds, NULL); CHKERRQ(ierr);
518 ierr = PetscOptionsGetReal(NULL, NULL, "-walltime_guard_estimator_alpha", &simCtx->walltimeGuardEstimatorAlpha, NULL); CHKERRQ(ierr);
519
520 if (simCtx->walltimeGuardWarmupSteps <= 0) {
521 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Invalid value for -walltime_guard_warmup_steps: %d. Must be > 0.", simCtx->walltimeGuardWarmupSteps);
522 }
523 if (simCtx->walltimeGuardMultiplier <= 0.0 || simCtx->walltimeGuardMultiplier > 5.0) {
524 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Invalid value for -walltime_guard_multiplier: %.6f. Must be in (0, 5].", (double)simCtx->walltimeGuardMultiplier);
525 }
526 if (simCtx->walltimeGuardMinSeconds <= 0.0) {
527 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Invalid value for -walltime_guard_min_seconds: %.6f. Must be > 0.", (double)simCtx->walltimeGuardMinSeconds);
528 }
529 if (simCtx->walltimeGuardEstimatorAlpha <= 0.0 || simCtx->walltimeGuardEstimatorAlpha > 1.0) {
530 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Invalid value for -walltime_guard_estimator_alpha: %.6f. Must be in (0, 1].", (double)simCtx->walltimeGuardEstimatorAlpha);
531 }
532
533 if(strcmp(simCtx->eulerianSource,"solve")!= 0 && strcmp(simCtx->eulerianSource,"load") != 0 && strcmp(simCtx->eulerianSource,"analytical")!=0){
534 SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Invalid value for -euler_field_source. Must be 'load','analytical' or 'solve'. You provided '%s'.",simCtx->eulerianSource);
535 }
536 if (simCtx->walltimeGuardEnabled) {
537 const char *job_start_env = getenv("PICURV_JOB_START_EPOCH");
538 const char *limit_env = getenv("PICURV_WALLTIME_LIMIT_SECONDS");
539 PetscBool job_start_ok = RuntimeWalltimeGuardParsePositiveSeconds(job_start_env, &simCtx->walltimeGuardJobStartEpochSeconds);
540 PetscBool limit_ok = RuntimeWalltimeGuardParsePositiveSeconds(limit_env, &simCtx->walltimeGuardLimitSeconds);
541
542 if (!job_start_ok || !limit_ok) {
543 simCtx->walltimeGuardActive = PETSC_FALSE;
544 simCtx->walltimeGuardJobStartEpochSeconds = 0.0;
545 simCtx->walltimeGuardLimitSeconds = 0.0;
546 LOG_ALLOW(
547 GLOBAL,
548 LOG_WARNING,
549 "Runtime walltime guard enabled but %s/%s are missing or invalid. Falling back to external shutdown signals only.\n",
550 "PICURV_JOB_START_EPOCH",
551 "PICURV_WALLTIME_LIMIT_SECONDS"
552 );
553 } else {
554 simCtx->walltimeGuardActive = PETSC_TRUE;
555 }
556 }
557
558 // --- Group 3
559 LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 3: High-Level Physics & Model Selection Flags\n");
560 ierr = PetscOptionsGetInt(NULL, NULL, "-imm", &simCtx->immersed, NULL); CHKERRQ(ierr);
561 ierr = PetscOptionsGetInt(NULL, NULL, "-fsi", &simCtx->movefsi, NULL); CHKERRQ(ierr);
562 ierr = PetscOptionsGetInt(NULL, NULL, "-rfsi", &simCtx->rotatefsi, NULL); CHKERRQ(ierr);
563 ierr = PetscOptionsGetInt(NULL, NULL, "-sediment", &simCtx->sediment, NULL); CHKERRQ(ierr);
564 ierr = PetscOptionsGetInt(NULL, NULL, "-rheology", &simCtx->rheology, NULL); CHKERRQ(ierr);
565 ierr = PetscOptionsGetInt(NULL, NULL, "-inv", &simCtx->invicid, NULL); CHKERRQ(ierr);
566 ierr = PetscOptionsGetInt(NULL, NULL, "-TwoD", &simCtx->TwoD, NULL); CHKERRQ(ierr);
567 ierr = PetscOptionsGetInt(NULL, NULL, "-thin", &simCtx->thin, NULL); CHKERRQ(ierr);
568 ierr = PetscOptionsGetInt(NULL, NULL, "-mframe", &simCtx->moveframe, NULL); CHKERRQ(ierr);
569 ierr = PetscOptionsGetInt(NULL, NULL, "-rframe", &simCtx->rotateframe, NULL); CHKERRQ(ierr);
570 ierr = PetscOptionsGetInt(NULL, NULL, "-blk", &simCtx->blank, NULL); CHKERRQ(ierr);
571 ierr = PetscOptionsGetInt(NULL, NULL, "-dgf_z", &simCtx->dgf_z, NULL); CHKERRQ(ierr);
572 ierr = PetscOptionsGetInt(NULL, NULL, "-dgf_y", &simCtx->dgf_y, NULL); CHKERRQ(ierr);
573 ierr = PetscOptionsGetInt(NULL, NULL, "-dgf_x", &simCtx->dgf_x, NULL); CHKERRQ(ierr);
574 ierr = PetscOptionsGetInt(NULL, NULL, "-dgf_az", &simCtx->dgf_az, NULL); CHKERRQ(ierr);
575 ierr = PetscOptionsGetInt(NULL, NULL, "-dgf_ay", &simCtx->dgf_ay, NULL); CHKERRQ(ierr);
576 ierr = PetscOptionsGetInt(NULL, NULL, "-dgf_ax", &simCtx->dgf_ax, NULL); CHKERRQ(ierr);
577 ierr = PetscOptionsGetString(NULL,NULL,"-analytical_type",simCtx->AnalyticalSolutionType,sizeof(simCtx->AnalyticalSolutionType),NULL);CHKERRQ(ierr);
578
579 // --- Group 4
580 LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 4: Specific Simulation Case Flags \n");
581 ierr = PetscOptionsGetInt(NULL, NULL, "-cop", &simCtx->cop, NULL); CHKERRQ(ierr);
582 ierr = PetscOptionsGetInt(NULL, NULL, "-fish", &simCtx->fish, NULL); CHKERRQ(ierr);
583 ierr = PetscOptionsGetInt(NULL, NULL, "-pizza", &simCtx->pizza, NULL); CHKERRQ(ierr);
584 ierr = PetscOptionsGetInt(NULL, NULL, "-turbine", &simCtx->turbine, NULL); CHKERRQ(ierr);
585 ierr = PetscOptionsGetInt(NULL, NULL, "-fishcyl", &simCtx->fishcyl, NULL); CHKERRQ(ierr);
586 ierr = PetscOptionsGetInt(NULL, NULL, "-eel", &simCtx->eel, NULL); CHKERRQ(ierr);
587 ierr = PetscOptionsGetInt(NULL, NULL, "-cstart", &simCtx->fish_c, NULL); CHKERRQ(ierr);
588 ierr = PetscOptionsGetInt(NULL, NULL, "-wing", &simCtx->wing, NULL); CHKERRQ(ierr);
589 ierr = PetscOptionsGetInt(NULL, NULL, "-mhv", &simCtx->MHV, NULL); CHKERRQ(ierr);
590 ierr = PetscOptionsGetInt(NULL, NULL, "-hydro", &simCtx->hydro, NULL); CHKERRQ(ierr);
591 ierr = PetscOptionsGetInt(NULL, NULL, "-lv", &simCtx->LV, NULL); CHKERRQ(ierr);
592 ierr = PetscOptionsGetInt(NULL, NULL, "-Pipe", &simCtx->Pipe, NULL); CHKERRQ(ierr);
593 ierr = PetscOptionsGetInt(NULL, NULL, "-Turbulent_Channel_z", &simCtx->channelz, NULL); CHKERRQ(ierr);
594 ierr = PetscOptionsGetReal(NULL,NULL,"-driven_flow_initial_force",&simCtx->drivingForceMagnitude,NULL);CHKERRQ(ierr);
595 ierr = PetscOptionsGetReal(NULL,NULL,"-driven_flow_scaling_factor",&simCtx->forceScalingFactor,NULL);CHKERRQ(ierr);
596 // --- Group 5
597 LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 5: Solver & Numerics Parameters \n");
598 char mom_solver_type_char[PETSC_MAX_PATH_LEN];
599 char solution_convergence_mode_char[PETSC_MAX_PATH_LEN];
600 PetscBool mom_solver_type_flg = PETSC_FALSE;
601 PetscBool solution_convergence_mode_flg = PETSC_FALSE;
602 ierr = PetscOptionsGetString(NULL, NULL, "-mom_solver_type", mom_solver_type_char, sizeof(mom_solver_type_char), &mom_solver_type_flg); CHKERRQ(ierr);
603 ierr = PetscOptionsGetInt(NULL, NULL, "-mom_max_pseudo_steps", &simCtx->mom_max_pseudo_steps, NULL); CHKERRQ(ierr);
604 ierr = PetscOptionsGetReal(NULL, NULL, "-mom_atol", &simCtx->mom_atol, NULL); CHKERRQ(ierr);
605 ierr = PetscOptionsGetReal(NULL, NULL, "-mom_rtol", &simCtx->mom_rtol, NULL); CHKERRQ(ierr);
606 ierr = PetscOptionsGetReal(NULL, NULL, "-mom_resid_atol", &simCtx->mom_resid_atol, NULL); CHKERRQ(ierr);
607 ierr = PetscOptionsGetReal(NULL, NULL, "-mom_resid_rtol", &simCtx->mom_resid_rtol, NULL); CHKERRQ(ierr);
608 ierr = PetscOptionsGetReal(NULL, NULL, "-imp_stol", &simCtx->imp_stol, NULL); CHKERRQ(ierr);
609 ierr = PetscOptionsGetInt(NULL, NULL, "-central", &simCtx->central, NULL); CHKERRQ(ierr);
610 ierr = PetscOptionsGetString(NULL, NULL, "-solution_convergence_mode",
611 solution_convergence_mode_char, sizeof(solution_convergence_mode_char),
612 &solution_convergence_mode_flg); CHKERRQ(ierr);
613 ierr = PetscOptionsGetBool(NULL, NULL, "-solution_convergence_enabled", &simCtx->solutionConvergenceEnabled, NULL); CHKERRQ(ierr);
614 ierr = PetscOptionsGetInt(NULL, NULL, "-solution_convergence_period_steps", &simCtx->solutionConvergencePeriodSteps, NULL); CHKERRQ(ierr);
615 ierr = PetscOptionsGetInt(NULL, NULL, "-solution_convergence_window_steps", &simCtx->solutionConvergenceWindowSteps, NULL); CHKERRQ(ierr);
616 /* Field statistics resolve here rather than inline, because a variable-arity
617 * window list needs its own parse site; see statistics_config.c. It must
618 * precede CreateAndInitializeAllVectors, which sizes the per-window
619 * accumulators from the window count resolved here. */
620 ierr = ParseFieldStatisticsConfig(simCtx); CHKERRQ(ierr);
621
622 // Keep parser acceptance aligned with the enum and FlowSolver dispatch.
623 if (mom_solver_type_flg) {
624 if(strcmp(mom_solver_type_char, "DUALTIME_PICARD_JAMESON_RK") == 0 ||
625 strcmp(mom_solver_type_char, "DUALTIME_PICARD_RK4") == 0) {
626 simCtx->mom_solver_type = MOMENTUM_SOLVER_DUALTIME_PICARD_JAMESON_RK;
627 } else if (strcmp(mom_solver_type_char, "EXPLICIT_RK") == 0) {
628 simCtx->mom_solver_type = MOMENTUM_SOLVER_EXPLICIT_RK;
629 } else if (strcmp(mom_solver_type_char, "newton_krylov") == 0) {
630 simCtx->mom_solver_type = MOMENTUM_SOLVER_NEWTON_KRYLOV;
631 } else {
632 LOG(GLOBAL, LOG_ERROR, "Invalid value for -mom_solver_type: '%s'. Valid options are: 'DUALTIME_PICARD_JAMESON_RK', 'EXPLICIT_RK', 'newton_krylov'.\n", mom_solver_type_char);
633 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Invalid value for -mom_solver_type: '%s'.", mom_solver_type_char);
634 }
635 }
636
637 if (solution_convergence_mode_flg) {
638 if (strcmp(solution_convergence_mode_char, "STEADY_DETERMINISTIC") == 0) {
639 simCtx->solutionConvergenceMode = SOLUTION_CONVERGENCE_STEADY_DETERMINISTIC;
640 } else if (strcmp(solution_convergence_mode_char, "PERIODIC_DETERMINISTIC") == 0) {
641 simCtx->solutionConvergenceMode = SOLUTION_CONVERGENCE_PERIODIC_DETERMINISTIC;
642 } else if (strcmp(solution_convergence_mode_char, "STATISTICAL_STEADY") == 0) {
643 simCtx->solutionConvergenceMode = SOLUTION_CONVERGENCE_STATISTICAL_STEADY;
644 } else if (strcmp(solution_convergence_mode_char, "TRANSIENT") == 0) {
645 simCtx->solutionConvergenceMode = SOLUTION_CONVERGENCE_TRANSIENT;
646 } else {
647 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
648 "Invalid value for -solution_convergence_mode: '%s'.", solution_convergence_mode_char);
649 }
650 }
651
652 if (simCtx->solutionConvergenceMode == SOLUTION_CONVERGENCE_PERIODIC_DETERMINISTIC &&
653 simCtx->solutionConvergencePeriodSteps <= 0) {
654 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
655 "solution convergence mode PERIODIC_DETERMINISTIC requires -solution_convergence_period_steps > 0.");
656 }
657 if (simCtx->solutionConvergenceMode == SOLUTION_CONVERGENCE_STATISTICAL_STEADY &&
658 simCtx->solutionConvergenceWindowSteps <= 0) {
659 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
660 "solution convergence mode STATISTICAL_STEADY requires -solution_convergence_window_steps > 0.");
661 }
662
663 // --- Multigrid Options ---
664 ierr = PetscOptionsGetInt(NULL, NULL, "-mg_level", &simCtx->mglevels, NULL); CHKERRQ(ierr);
665 ierr = PetscOptionsGetInt(NULL, NULL, "-mg_max_it", &simCtx->mg_MAX_IT, NULL); CHKERRQ(ierr);
666 ierr = PetscOptionsGetInt(NULL, NULL, "-mg_idx", &simCtx->mg_idx, NULL); CHKERRQ(ierr);
667 ierr = PetscOptionsGetInt(NULL, NULL, "-mg_pre_it", &simCtx->mg_preItr, NULL); CHKERRQ(ierr);
668 ierr = PetscOptionsGetInt(NULL, NULL, "-mg_post_it", &simCtx->mg_poItr, NULL); CHKERRQ(ierr);
669
670 // --- Other Solver Options ---
671 ierr = PetscOptionsGetInt(NULL, NULL, "-poisson", &simCtx->poisson, NULL); CHKERRQ(ierr);
672 ierr = PetscOptionsGetReal(NULL, NULL, "-poisson_tol", &simCtx->poisson_tol, NULL); CHKERRQ(ierr);
673 ierr = PetscOptionsGetInt(NULL, NULL, "-str", &simCtx->STRONG_COUPLING, NULL); CHKERRQ(ierr);
674 ierr = PetscOptionsGetReal(NULL, NULL, "-ren", &simCtx->ren, NULL); CHKERRQ(ierr);
675 ierr = PetscOptionsGetReal(NULL, NULL, "-pseudo_cfl", &simCtx->pseudo_cfl, NULL); CHKERRQ(ierr);
676 ierr = PetscOptionsGetReal(NULL, NULL, "-max_pseudo_cfl", &simCtx->max_pseudo_cfl, NULL); CHKERRQ(ierr);
677 ierr = PetscOptionsGetReal(NULL, NULL, "-min_pseudo_cfl", &simCtx->min_pseudo_cfl, NULL); CHKERRQ(ierr);
678 ierr = PetscOptionsGetReal(NULL, NULL, "-pseudo_cfl_reduction_factor", &simCtx->pseudo_cfl_reduction_factor, NULL); CHKERRQ(ierr);
679 ierr = PetscOptionsGetReal(NULL, NULL, "-pseudo_cfl_growth_factor", &simCtx->pseudo_cfl_growth_factor, NULL); CHKERRQ(ierr);
680 // Read the deprecated RK4 spelling first so the canonical Jameson option wins if both are present.
681 ierr = PetscOptionsGetReal(NULL,NULL, "-mom_dt_rk4_residual_norm_noise_allowance_factor",&simCtx->mom_dt_jameson_residual_norm_noise_allowance_factor,NULL);CHKERRQ(ierr);
682 ierr = PetscOptionsGetReal(NULL,NULL, "-mom_dt_jameson_residual_norm_noise_allowance_factor",&simCtx->mom_dt_jameson_residual_norm_noise_allowance_factor,NULL);CHKERRQ(ierr);
683 ierr = PetscOptionsGetBool(NULL, NULL, "-no_pseudo_cfl_backtrack", &simCtx->no_pseudo_cfl_backtrack, NULL); CHKERRQ(ierr);
684 ierr = PetscOptionsGetReal(NULL, NULL, "-mom_ratio_ema_alpha", &simCtx->mom_ratio_ema_alpha, NULL); CHKERRQ(ierr);
685 if (simCtx->min_pseudo_cfl <= 0.0 ||
686 simCtx->pseudo_cfl < simCtx->min_pseudo_cfl ||
687 simCtx->pseudo_cfl > simCtx->max_pseudo_cfl) {
688 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE,
689 "Pseudo-CFL controls require 0 < minimum <= initial <= maximum.");
690 }
691 if (simCtx->pseudo_cfl_growth_factor < 1.0 ||
692 simCtx->pseudo_cfl_reduction_factor <= 0.0 ||
693 simCtx->pseudo_cfl_reduction_factor >= 1.0 ||
694 simCtx->mom_dt_jameson_residual_norm_noise_allowance_factor < 1.0) {
695 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE,
696 "Pseudo-CFL controls require growth_factor >= 1, 0 < reduction_factor < 1, and noise allowance >= 1.");
697 }
698 if (simCtx->mom_ratio_ema_alpha < 0.0 || simCtx->mom_ratio_ema_alpha > 1.0) {
699 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE,
700 "-mom_ratio_ema_alpha must be in [0, 1].");
701 }
702 ierr = PetscOptionsHasName(NULL, NULL, "-ps_ksp_pic_monitor_true_residual", &simCtx->ps_ksp_pic_monitor_true_residual); CHKERRQ(ierr);
703 ierr = PetscOptionsGetBool(NULL, NULL, "-mom_nk_pic_monitor", &simCtx->mom_nk_monitor_history, NULL); CHKERRQ(ierr);
704 {
705 PetscInt ic_mode = (PetscInt)simCtx->initialConditionMode;
706 PetscInt ic_field = (PetscInt)simCtx->initialConditionField;
707 ierr = PetscOptionsGetInt(NULL, NULL, "-finit", &ic_mode, NULL); CHKERRQ(ierr);
708 ierr = PetscOptionsGetInt(NULL, NULL, "-ic_field", &ic_field, NULL); CHKERRQ(ierr);
709 simCtx->initialConditionMode = (InitialConditionMode)ic_mode;
710 simCtx->initialConditionField = (InitialConditionField)ic_field;
711 }
712 ierr = PetscOptionsGetString(NULL, NULL, "-ic_dir", simCtx->initialConditionDirectory,
713 sizeof(simCtx->initialConditionDirectory), NULL); CHKERRQ(ierr);
714 if (simCtx->initialConditionMode < IC_MODE_ZERO || simCtx->initialConditionMode > IC_MODE_FILE) {
715 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE,
716 "Invalid value for -finit. Expected an initial-condition mode in [0,4], got %d.",
717 simCtx->initialConditionMode);
718 }
719 if (simCtx->initialConditionField < IC_FIELD_UCAT || simCtx->initialConditionField > IC_FIELD_UCONT) {
720 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE,
721 "Invalid value for -ic_field. Expected 0 (Ucat) or 1 (Ucont), got %d.",
722 simCtx->initialConditionField);
723 }
724 ierr = PetscOptionsGetReal(NULL, NULL, "-ucont_x", &simCtx->InitialConstantContra.x, NULL); CHKERRQ(ierr);
725 ierr = PetscOptionsGetReal(NULL, NULL, "-ucont_y", &simCtx->InitialConstantContra.y, NULL); CHKERRQ(ierr);
726 ierr = PetscOptionsGetReal(NULL, NULL, "-ucont_z", &simCtx->InitialConstantContra.z, NULL); CHKERRQ(ierr);
727 {
728 PetscInt fd_int = (PetscInt)FLOW_DIR_UNSET;
729 PetscBool fd_set = PETSC_FALSE;
730 ierr = PetscOptionsGetInt(NULL, NULL, "-flow_direction", &fd_int, &fd_set); CHKERRQ(ierr);
731 if (fd_set) simCtx->flowDirection = (FlowDirection)fd_int;
732 }
733 ierr = PetscOptionsGetReal(NULL, NULL, "-ic_velocity_physical", &simCtx->icVelocityPhysical, NULL); CHKERRQ(ierr);
734 ierr = PetscOptionsGetReal(NULL, NULL, "-analytical_uniform_u", &simCtx->AnalyticalUniformVelocity.x, NULL); CHKERRQ(ierr);
735 ierr = PetscOptionsGetReal(NULL, NULL, "-analytical_uniform_v", &simCtx->AnalyticalUniformVelocity.y, NULL); CHKERRQ(ierr);
736 ierr = PetscOptionsGetReal(NULL, NULL, "-analytical_uniform_w", &simCtx->AnalyticalUniformVelocity.z, NULL); CHKERRQ(ierr);
737 PetscBool verification_scalar_value_set = PETSC_FALSE;
738 PetscBool verification_scalar_phi0_set = PETSC_FALSE;
739 PetscBool verification_scalar_slope_x_set = PETSC_FALSE;
740 PetscBool verification_scalar_amplitude_set = PETSC_FALSE;
741 PetscBool verification_scalar_kx_set = PETSC_FALSE;
742 PetscBool verification_scalar_ky_set = PETSC_FALSE;
743 PetscBool verification_scalar_kz_set = PETSC_FALSE;
744 ierr = PetscOptionsGetString(NULL, NULL, "-verification_diffusivity_mode",
745 simCtx->verificationDiffusivity.mode,
746 sizeof(simCtx->verificationDiffusivity.mode), NULL); CHKERRQ(ierr);
747 ierr = PetscOptionsGetString(NULL, NULL, "-verification_diffusivity_profile",
748 simCtx->verificationDiffusivity.profile,
749 sizeof(simCtx->verificationDiffusivity.profile), NULL); CHKERRQ(ierr);
750 ierr = PetscOptionsGetReal(NULL, NULL, "-verification_diffusivity_gamma0",
751 &simCtx->verificationDiffusivity.gamma0, NULL); CHKERRQ(ierr);
752 ierr = PetscOptionsGetReal(NULL, NULL, "-verification_diffusivity_slope_x",
753 &simCtx->verificationDiffusivity.slope_x, NULL); CHKERRQ(ierr);
754 ierr = PetscOptionsGetString(NULL, NULL, "-verification_scalar_mode",
755 simCtx->verificationScalar.mode,
756 sizeof(simCtx->verificationScalar.mode), NULL); CHKERRQ(ierr);
757 ierr = PetscOptionsGetString(NULL, NULL, "-verification_scalar_profile",
758 simCtx->verificationScalar.profile,
759 sizeof(simCtx->verificationScalar.profile), NULL); CHKERRQ(ierr);
760 ierr = PetscOptionsGetReal(NULL, NULL, "-verification_scalar_value",
761 &simCtx->verificationScalar.value, &verification_scalar_value_set); CHKERRQ(ierr);
762 ierr = PetscOptionsGetReal(NULL, NULL, "-verification_scalar_phi0",
763 &simCtx->verificationScalar.phi0, &verification_scalar_phi0_set); CHKERRQ(ierr);
764 ierr = PetscOptionsGetReal(NULL, NULL, "-verification_scalar_slope_x",
765 &simCtx->verificationScalar.slope_x, &verification_scalar_slope_x_set); CHKERRQ(ierr);
766 ierr = PetscOptionsGetReal(NULL, NULL, "-verification_scalar_amplitude",
767 &simCtx->verificationScalar.amplitude, &verification_scalar_amplitude_set); CHKERRQ(ierr);
768 ierr = PetscOptionsGetReal(NULL, NULL, "-verification_scalar_kx",
769 &simCtx->verificationScalar.kx, &verification_scalar_kx_set); CHKERRQ(ierr);
770 ierr = PetscOptionsGetReal(NULL, NULL, "-verification_scalar_ky",
771 &simCtx->verificationScalar.ky, &verification_scalar_ky_set); CHKERRQ(ierr);
772 ierr = PetscOptionsGetReal(NULL, NULL, "-verification_scalar_kz",
773 &simCtx->verificationScalar.kz, &verification_scalar_kz_set); CHKERRQ(ierr);
774 simCtx->verificationDiffusivity.enabled =
775 (PetscBool)(simCtx->verificationDiffusivity.mode[0] != '\0' ||
776 simCtx->verificationDiffusivity.profile[0] != '\0');
777 simCtx->verificationScalar.enabled =
778 (PetscBool)(simCtx->verificationScalar.mode[0] != '\0' ||
779 simCtx->verificationScalar.profile[0] != '\0');
780 if (simCtx->verificationDiffusivity.enabled) {
781 if (strcmp(simCtx->eulerianSource, "analytical") != 0) {
782 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE,
783 "verification diffusivity overrides require -euler_field_source \"analytical\".");
784 }
785 if (strcmp(simCtx->verificationDiffusivity.mode, "analytical") != 0) {
786 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
787 "Unsupported -verification_diffusivity_mode '%s'. Only 'analytical' is supported.",
788 simCtx->verificationDiffusivity.mode);
789 }
790 if (strcmp(simCtx->verificationDiffusivity.profile, "LINEAR_X") != 0) {
791 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
792 "Unsupported -verification_diffusivity_profile '%s'. Only 'LINEAR_X' is supported.",
793 simCtx->verificationDiffusivity.profile);
794 }
795 }
796 if (simCtx->verificationScalar.enabled) {
797 if (strcmp(simCtx->eulerianSource, "analytical") != 0) {
798 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE,
799 "verification scalar overrides require -euler_field_source \"analytical\".");
800 }
801 if (strcmp(simCtx->verificationScalar.mode, "analytical") != 0) {
802 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
803 "Unsupported -verification_scalar_mode '%s'. Only 'analytical' is supported.",
804 simCtx->verificationScalar.mode);
805 }
806 if (strcmp(simCtx->verificationScalar.profile, "CONSTANT") == 0) {
807 if (!verification_scalar_value_set) {
808 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
809 "verification scalar profile CONSTANT requires -verification_scalar_value.");
810 }
811 } else if (strcmp(simCtx->verificationScalar.profile, "LINEAR_X") == 0) {
812 if (!verification_scalar_phi0_set || !verification_scalar_slope_x_set) {
813 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
814 "verification scalar profile LINEAR_X requires -verification_scalar_phi0 and -verification_scalar_slope_x.");
815 }
816 } else if (strcmp(simCtx->verificationScalar.profile, "SIN_PRODUCT") == 0) {
817 if (!verification_scalar_amplitude_set || !verification_scalar_kx_set ||
818 !verification_scalar_ky_set || !verification_scalar_kz_set) {
819 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
820 "verification scalar profile SIN_PRODUCT requires -verification_scalar_amplitude, -verification_scalar_kx, -verification_scalar_ky, and -verification_scalar_kz.");
821 }
822 } else {
823 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
824 "Unsupported -verification_scalar_profile '%s'. Supported profiles: CONSTANT, LINEAR_X, SIN_PRODUCT.",
825 simCtx->verificationScalar.profile);
826 }
827 }
828 // NOTE: cdisx,cdisy,cdisz haven't been parsed, add if necessary.
829
830 // --- Group 6
831 LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 6: Physical & Geometric Parameters \n");
832 ierr = PetscOptionsGetReal(NULL,NULL,"-schmidt_number",&simCtx->schmidt_number,NULL);CHKERRQ(ierr);
833 ierr = PetscOptionsGetReal(NULL,NULL,"-turb_schmidt_number",&simCtx->Turbulent_schmidt_number,NULL);CHKERRQ(ierr);
834 ierr = PetscOptionsGetInt(NULL, NULL, "-no_of_bodies", &simCtx->NumberOfBodies, NULL); CHKERRQ(ierr);
835 ierr = PetscOptionsGetReal(NULL,NULL,"-wall_roughness",&simCtx->wall_roughness_height,NULL);CHKERRQ(ierr);
836 // NOTE: angle is not parsed in the original code, it set programmatically. We will follow that.
837 // NOTE: max_angle is calculated based on other flags (like MHV) in the legacy code.
838 // We will defer that logic to a later setup stage and not parse them directly.
839 // The Scaling Information is calculated here
840 ierr = ParseScalingInformation(simCtx); CHKERRQ(ierr);
841
842 // --- Group 7
843 LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 7: Grid, Domain, and Boundary Condition Settings \n");
844 ierr = PetscOptionsGetInt(NULL, NULL, "-nblk", &simCtx->block_number, NULL); CHKERRQ(ierr); // This is also a modern option
845 ierr = PetscOptionsGetInt(NULL, NULL, "-inlet", &simCtx->inletprofile, NULL); CHKERRQ(ierr);
846 ierr = PetscOptionsGetInt(NULL, NULL, "-Ogrid", &simCtx->Ogrid, NULL); CHKERRQ(ierr);
847 // NOTE: channelz was not parsed, likely set programmatically. We will omit its parsing call.
848 ierr = PetscOptionsGetInt(NULL, NULL, "-grid1d", &simCtx->grid1d, NULL); CHKERRQ(ierr);
849 ierr = PetscOptionsGetBool(NULL, NULL, "-grid", &simCtx->generate_grid, NULL); CHKERRQ(ierr);
850 ierr = PetscOptionsGetString(NULL, NULL, "-grid_file", simCtx->grid_file, PETSC_MAX_PATH_LEN, NULL); CHKERRQ(ierr);
851 ierr = PetscOptionsGetInt(NULL, NULL, "-da_processors_x", &simCtx->da_procs_x, NULL); CHKERRQ(ierr);
852 ierr = PetscOptionsGetInt(NULL, NULL, "-da_processors_y", &simCtx->da_procs_y, NULL); CHKERRQ(ierr);
853 ierr = PetscOptionsGetInt(NULL, NULL, "-da_processors_z", &simCtx->da_procs_z, NULL); CHKERRQ(ierr);
854 ierr = PetscOptionsGetInt(NULL, NULL, "-pbc_domain", &simCtx->blkpbc, NULL); CHKERRQ(ierr);
855 // NOTE: pseudo_periodic was not parsed. We will omit its parsing call.
856 ierr = PetscOptionsGetReal(NULL, NULL, "-grid_rotation_angle", &simCtx->grid_rotation_angle, NULL); CHKERRQ(ierr);
857 ierr = PetscOptionsGetReal(NULL, NULL, "-Croty", &simCtx->Croty, NULL); CHKERRQ(ierr);
858 ierr = PetscOptionsGetReal(NULL, NULL, "-Crotz", &simCtx->Crotz, NULL); CHKERRQ(ierr);
859 PetscBool bcs_flg;
860 char file_list_str[PETSC_MAX_PATH_LEN * 10]; // Buffer for comma-separated list
861
862 ierr = PetscOptionsGetString(NULL, NULL, "-bcs_files", file_list_str, sizeof(file_list_str), &bcs_flg); CHKERRQ(ierr);
863 ierr = PetscOptionsGetReal(NULL, NULL, "-U_bc", &simCtx->U_bc, NULL); CHKERRQ(ierr);
864
865 if (bcs_flg) {
866 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Found -bcs_files option, overriding default.\n");
867
868 // A. Clean up the default memory we allocated in Phase 1.
869 ierr = PetscFree(simCtx->bcs_files[0]); CHKERRQ(ierr);
870 ierr = PetscFree(simCtx->bcs_files); CHKERRQ(ierr);
871 simCtx->num_bcs_files = 0;
872 simCtx->bcs_files = NULL;
873
874 // B. Parse the user-provided comma-separated list.
875 char *token;
876 char *str_copy;
877 ierr = PetscStrallocpy(file_list_str, &str_copy); CHKERRQ(ierr);
878
879 // First pass: count the number of files.
880 token = strtok(str_copy, ",");
881 while (token) {
882 simCtx->num_bcs_files++;
883 token = strtok(NULL, ",");
884 }
885 ierr = PetscFree(str_copy); CHKERRQ(ierr);
886
887 // Second pass: allocate memory and store the filenames.
888 ierr = PetscMalloc1(simCtx->num_bcs_files, &simCtx->bcs_files); CHKERRQ(ierr);
889 ierr = PetscStrallocpy(file_list_str, &str_copy); CHKERRQ(ierr);
890 token = strtok(str_copy, ",");
891 for (PetscInt i = 0; i < simCtx->num_bcs_files; i++) {
892 ierr = PetscStrallocpy(token, &simCtx->bcs_files[i]); CHKERRQ(ierr);
893 token = strtok(NULL, ",");
894 }
895 ierr = PetscFree(str_copy); CHKERRQ(ierr);
896 }
897
898
899 // --- Group 8
900 LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 8: Turbulence Modeling (LES/RANS) \n");
901 PetscInt temp_les_model;
902 ierr = PetscOptionsGetInt(NULL, NULL, "-les", &temp_les_model, NULL); CHKERRQ(ierr);
903 simCtx->les = (LESModelType)temp_les_model;
904 ierr = PetscOptionsGetInt(NULL, NULL, "-rans", &simCtx->rans, NULL); CHKERRQ(ierr);
905 ierr = PetscOptionsGetInt(NULL, NULL, "-wallfunction", &simCtx->wallfunction, NULL); CHKERRQ(ierr);
906 ierr = PetscOptionsGetInt(NULL, NULL, "-mixed", &simCtx->mixed, NULL); CHKERRQ(ierr);
907 ierr = PetscOptionsGetInt(NULL, NULL, "-clark", &simCtx->clark, NULL); CHKERRQ(ierr);
908 ierr = PetscOptionsGetInt(NULL, NULL, "-dynamic_freq", &simCtx->dynamic_freq, NULL); CHKERRQ(ierr);
909 ierr = PetscOptionsGetReal(NULL, NULL, "-max_cs", &simCtx->max_cs, NULL); CHKERRQ(ierr);
910 ierr = PetscOptionsGetReal(NULL, NULL, "-const_cs", &simCtx->Const_CS, NULL); CHKERRQ(ierr);
911 ierr = PetscOptionsGetInt(NULL, NULL, "-testfilter_ik", &simCtx->testfilter_ik, NULL); CHKERRQ(ierr);
912 ierr = PetscOptionsGetInt(NULL, NULL, "-testfilter_1d", &simCtx->testfilter_1d, NULL); CHKERRQ(ierr);
913 ierr = PetscOptionsGetInt(NULL, NULL, "-i_homo_filter", &simCtx->i_homo_filter, NULL); CHKERRQ(ierr);
914 ierr = PetscOptionsGetInt(NULL, NULL, "-j_homo_filter", &simCtx->j_homo_filter, NULL); CHKERRQ(ierr);
915 ierr = PetscOptionsGetInt(NULL, NULL, "-k_homo_filter", &simCtx->k_homo_filter, NULL); CHKERRQ(ierr);
916
917 // --- Group 9
918 LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 9: Particle / DMSwarm Data & Settings \n");
919 ierr = PetscOptionsGetInt(NULL, NULL, "-numParticles", &simCtx->np, NULL); CHKERRQ(ierr);
920 ierr = PetscOptionsGetBool(NULL, NULL, "-read_fields", &simCtx->readFields, NULL); CHKERRQ(ierr);
921 PetscInt temp_pinit = (PetscInt)PARTICLE_INIT_SURFACE_RANDOM;
922 ierr = PetscOptionsGetInt(NULL, NULL, "-pinit", &temp_pinit, NULL); CHKERRQ(ierr);
923 simCtx->ParticleInitialization = (ParticleInitializationType)temp_pinit;
924 PetscInt temp_interp = (PetscInt)INTERP_TRILINEAR;
925 ierr = PetscOptionsGetInt(NULL, NULL, "-interpolation_method", &temp_interp, NULL); CHKERRQ(ierr);
926 simCtx->interpolationMethod = (InterpolationMethod)temp_interp;
927 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Interpolation method: %s\n",
928 simCtx->interpolationMethod == INTERP_TRILINEAR ? "Trilinear (direct cell-center)" : "CornerAveraged (legacy)");
929 ierr = PetscOptionsGetReal(NULL, NULL, "-psrc_x", &simCtx->psrc_x, NULL); CHKERRQ(ierr);
930 ierr = PetscOptionsGetReal(NULL, NULL, "-psrc_y", &simCtx->psrc_y, NULL); CHKERRQ(ierr);
931 ierr = PetscOptionsGetReal(NULL, NULL, "-psrc_z", &simCtx->psrc_z, NULL); CHKERRQ(ierr);
932 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Particle initialization mode: %s. Point source: (%.6f, %.6f, %.6f)\n",
933 ParticleInitializationToString(simCtx->ParticleInitialization),
934 simCtx->psrc_x, simCtx->psrc_y, simCtx->psrc_z);
935 ierr = PetscOptionsGetString(NULL,NULL,"-particle_restart_mode",simCtx->particleRestartMode,sizeof(simCtx->particleRestartMode),NULL); CHKERRQ(ierr);
936 // Validation for Particle Restart Mode
937 if (strcmp(simCtx->particleRestartMode, "load") != 0 && strcmp(simCtx->particleRestartMode, "init") != 0) {
938 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Invalid value for -particle_restart_mode. Must be 'load' or 'init'. You provided '%s'.", simCtx->particleRestartMode);
939 }
940 ierr = InitializeBrownianRNG(simCtx); CHKERRQ(ierr);
941 // --- Group 10
942 LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 10: Immersed Boundary & FSI Data Object Pointers \n");
943 ierr = PetscOptionsGetBool(NULL, NULL, "-rs_fsi", &simCtx->rstart_fsi, NULL); CHKERRQ(ierr);
944 ierr = PetscOptionsGetInt(NULL, NULL, "-duplicate", &simCtx->duplicate, NULL); CHKERRQ(ierr);
945
946 // --- Group 11
947 LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 11: Top-Level Managers & Custom Configuration \n");
948 ierr = PetscOptionsGetInt(NULL, NULL, "-logfreq", &simCtx->LoggingFrequency, NULL); CHKERRQ(ierr);
949
950 if (simCtx->num_bcs_files != simCtx->block_number) {
951 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_INCOMP, "Number of BC files (%d) does not match number of blocks (%d). Use -bcs_files \"file1.dat,file2.dat,...\".", simCtx->num_bcs_files, simCtx->block_number);
952 }
953
954 // --- Group 12
955 LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 12: Post-Processing Information.\n");
956 // This logic determines the Post Processing configuration and STORES it in simCtx for later reference and cleanup.
957 ierr = PetscOptionsGetString(NULL,NULL,"-postprocessing_config_file",simCtx->PostprocessingControlFile,PETSC_MAX_PATH_LEN,NULL); CHKERRQ(ierr);
958 /* Parse post settings for both solver and post-processor binaries using the single pre-allocated pps object. */
959 ierr = ParsePostProcessingSettings(simCtx);
960
961 // === 5. Dependent Parameter Calculations ================================
962 // Some parameters depend on others, so we calculate them here.
963 simCtx->StartTime = (PetscReal)simCtx->StartStep*simCtx->dt;
964 simCtx->ti = simCtx->StartTime;
965 simCtx->step = simCtx->StartStep;
966
967 // === 5. Log Summary and Finalize Setup ==================================
968 LOG_ALLOW(GLOBAL, LOG_DEBUG, "-- Console Output Functions [Total : %d] : --\n", simCtx->nAllowed);
969 for (PetscInt i = 0; i < simCtx->nAllowed; ++i) {
970 LOG_ALLOW(GLOBAL, LOG_DEBUG, " [%2d] «%s»\n", i, simCtx->allowedFuncs[i]);
971 }
972
973 LOG_ALLOW(GLOBAL, LOG_INFO, "Configuration complete. Key parameters:\n");
974 LOG_ALLOW(GLOBAL, LOG_INFO, " - Run mode: %s\n", simCtx->OnlySetup ? "SETUP ONLY" : "Full Simulation");
975 LOG_ALLOW(GLOBAL, LOG_INFO, " - Time steps: %d (from %d to %d)\n", simCtx->StepsToRun, simCtx->StartStep, simCtx->StartStep + simCtx->StepsToRun);
976 LOG_ALLOW(GLOBAL, LOG_INFO, " - Time step size (dt): %g\n", simCtx->dt);
977 if (simCtx->tiout > 0) {
978 LOG_ALLOW(GLOBAL, LOG_INFO, " - Field/restart output cadence: every %d step(s)\n", simCtx->tiout);
979 } else {
980 LOG_ALLOW(GLOBAL, LOG_INFO, " - Field/restart output cadence: DISABLED\n");
981 }
982 LOG_ALLOW(GLOBAL, LOG_INFO, " - Immersed Boundary: %s\n", simCtx->immersed ? "ENABLED" : "DISABLED");
983 LOG_ALLOW(GLOBAL, LOG_INFO, " - Particles: %d\n", simCtx->np);
984 if (simCtx->np > 0) {
985 if (simCtx->particleConsoleOutputFreq > 0) {
986 LOG_ALLOW(GLOBAL, LOG_INFO, " - Particle console cadence: every %d step(s)\n", simCtx->particleConsoleOutputFreq);
987 } else {
988 LOG_ALLOW(GLOBAL, LOG_INFO, " - Particle console cadence: DISABLED\n");
989 }
990 LOG_ALLOW(GLOBAL, LOG_INFO, " - Particle console row subsampling: every %d particle(s)\n", simCtx->LoggingFrequency);
991 }
992 if (simCtx->StartStep > 0 && simCtx->np > 0) {
993 LOG_ALLOW(GLOBAL, LOG_INFO, " - Particle Restart Mode: %s\n", simCtx->particleRestartMode);
994 }
995
996 // --- Initialize PETSc's internal performance logging stage ---
997 ierr = PetscLogDefaultBegin(); CHKERRQ(ierr); // REDUNDANT but safe.
998 ierr = PetscMemorySetGetMaximumUsage(); CHKERRQ(ierr);
999
1000 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Finished CreateSimulationContext successfully on rank %d.\n", simCtx->rank);
1001
1002 PROFILE_FUNCTION_END;
1003 PetscFunctionReturn(0);
1004}
1005
1006#undef __FUNCT__
1007#define __FUNCT__ "PetscMkdirRecursive"
1008/**
1009 * @brief Create a directory path recursively using PETSc-compatible error handling.
1010 */
1011static PetscErrorCode PetscMkdirRecursive(const char *path)
1012{
1013 PetscErrorCode ierr;
1014 char tmp_path[PETSC_MAX_PATH_LEN];
1015 char *p = NULL;
1016 size_t len;
1017 PetscBool exists;
1018
1019 PetscFunctionBeginUser;
1020
1021 // Create a mutable copy of the path
1022 len = strlen(path);
1023 if (len >= sizeof(tmp_path)) {
1024 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Path is too long to process: %s", path);
1025 }
1026 strcpy(tmp_path, path);
1027
1028 // If the path ends with a separator, remove it
1029 if (tmp_path[len - 1] == '/') {
1030 tmp_path[len - 1] = 0;
1031 }
1032
1033 // Iterate through the path, creating each directory level
1034 for (p = tmp_path + 1; *p; p++) {
1035 if (*p == '/') {
1036 *p = 0; // Temporarily terminate the string
1037
1038 // Check if this directory level exists
1039 ierr = PetscTestDirectory(tmp_path, 'r', &exists); CHKERRQ(ierr);
1040 if (!exists) {
1041 ierr = PetscMkdir(tmp_path); CHKERRQ(ierr);
1042 }
1043
1044 *p = '/'; // Restore the separator
1045 }
1046 }
1047
1048 // Create the final, full directory path
1049 ierr = PetscTestDirectory(tmp_path, 'r', &exists); CHKERRQ(ierr);
1050 if (!exists) {
1051 ierr = PetscMkdir(tmp_path); CHKERRQ(ierr);
1052 }
1053
1054 PetscFunctionReturn(0);
1055}
1056
1057#undef __FUNCT__
1058#define __FUNCT__ "SetupSimulationEnvironment"
1059/**
1060 * @brief Internal helper implementation: `SetupSimulationEnvironment()`.
1061 * @details Local to this translation unit.
1062 */
1063PetscErrorCode SetupSimulationEnvironment(SimCtx *simCtx)
1064{
1065 PetscErrorCode ierr;
1066 PetscMPIInt rank;
1067 PetscBool exists;
1068
1069 PetscFunctionBeginUser;
1070 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
1071
1072 LOG_ALLOW(GLOBAL, LOG_INFO, "--- Setting up simulation environment ---\n");
1073
1074 /* =====================================================================
1075 * Phase 1: Check for all required and optional INPUT files.
1076 * ===================================================================== */
1077 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Phase 1: Verifying input files...\n");
1078
1079 // --- Mandatory Inputs ---
1080 if (!simCtx->generate_grid) {
1081 ierr = VerifyPathExistence(simCtx->grid_file, PETSC_FALSE, PETSC_FALSE, "Grid file", &exists); CHKERRQ(ierr);
1082 }
1083 for (PetscInt i = 0; i < simCtx->num_bcs_files; i++) {
1084 char desc[128];
1085 ierr = PetscSNPrintf(desc, sizeof(desc), "BCS file #%d", i + 1); CHKERRQ(ierr);
1086 ierr = VerifyPathExistence(simCtx->bcs_files[i], PETSC_FALSE, PETSC_FALSE, desc, &exists); CHKERRQ(ierr);
1087 }
1088
1089 // --- Optional Inputs (these produce warnings if missing) ---
1090 if (simCtx->useCfg) {
1091 ierr = VerifyPathExistence(simCtx->allowedFile, PETSC_FALSE, PETSC_TRUE, "Whitelist config file", &exists); CHKERRQ(ierr);
1092 }
1093 if (simCtx->useProfilingSelectedFuncsCfg) {
1094 ierr = VerifyPathExistence(simCtx->profilingSelectedFuncsFile, PETSC_FALSE, PETSC_TRUE, "Profiling config file", &exists); CHKERRQ(ierr);
1095 }
1096 if (simCtx->exec_mode == EXEC_MODE_POSTPROCESSOR) {
1097 ierr = VerifyPathExistence(simCtx->PostprocessingControlFile, PETSC_FALSE, PETSC_TRUE, "Post-processing control file", &exists); CHKERRQ(ierr);
1098 }
1099
1100
1101 /* =====================================================================
1102 * Phase 2: Validate directories specific to the execution mode.
1103 * ===================================================================== */
1104 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Phase 2: Verifying execution mode directories...\n");
1105 // The data source directory must exist if we intend to load any data from it.
1106 // This is true if:
1107 // 1. We are restarting from a previous time step (StartStep > 0), which implies
1108 // loading Eulerian fields and/or particle fields.
1109 // 2. We are starting from t=0 but are explicitly told to load the initial
1110 // Eulerian fields from a file (eulerianSource == "load").
1111 if (simCtx->StartStep > 0 || strcmp(simCtx->eulerianSource,"load")== 0){ // If this is a restart run
1112 ierr = VerifyPathExistence(simCtx->restart_dir, PETSC_TRUE, PETSC_FALSE, "Restart source directory", &exists); CHKERRQ(ierr);
1113 }
1114 if (simCtx->StartStep == 0 && strcmp(simCtx->eulerianSource, "solve") == 0 &&
1115 simCtx->initialConditionMode == IC_MODE_FILE) {
1116 ierr = VerifyPathExistence(simCtx->initialConditionDirectory, PETSC_TRUE, PETSC_FALSE,
1117 "Initial-condition source directory", &exists); CHKERRQ(ierr);
1118 }
1119 if (simCtx->exec_mode == EXEC_MODE_POSTPROCESSOR) {
1120 ierr = VerifyPathExistence(simCtx->pps->source_dir, PETSC_TRUE, PETSC_FALSE, "Post-processing source directory", &exists); CHKERRQ(ierr);
1121 }
1122
1123 /* =====================================================================
1124 * Phase 3: Create and prepare all OUTPUT directories.
1125 * ===================================================================== */
1126 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Phase 3: Preparing output directories...\n");
1127
1128 if (rank == 0){
1129 if(simCtx->exec_mode == EXEC_MODE_SOLVER){
1130 // --- Prepare Log Directory ---
1131 if (!simCtx->continueMode) {
1132 // Only wipe logs on fresh runs; continue mode appends to existing logs.
1133 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Creating/cleaning log directory: %s\n", simCtx->log_dir);
1134 ierr = PetscRMTree(simCtx->log_dir); // Wipes the directory and its contents
1135 if (ierr) { /* Ignore file-not-found error, but fail on others */
1136 PetscError(PETSC_COMM_SELF, __LINE__, __FUNCT__, __FILE__, ierr, PETSC_ERROR_INITIAL, "Could not remove existing log directory '%s'. Check permissions.", simCtx->log_dir);
1137 }
1138 ierr = PetscMkdir(simCtx->log_dir); CHKERRQ(ierr);
1139 } else {
1140 // In continue mode, ensure log directory exists but don't wipe it.
1141 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Continue mode: preserving existing log directory: %s\n", simCtx->log_dir);
1142 ierr = PetscMkdir(simCtx->log_dir); CHKERRQ(ierr);
1143 }
1144
1145 // --- Prepare Output Directory ---
1146 // The checkpoint coordinator creates its fixed internal hierarchy
1147 // transactionally when a checkpoint is due.
1148 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Verifying main output directory: %s\n", simCtx->output_dir);
1149 ierr = PetscTestDirectory(simCtx->output_dir, 'r', &exists); CHKERRQ(ierr);
1150 if (!exists) {
1151 LOG_ALLOW(GLOBAL, LOG_INFO, "Output directory not found. Creating: %s\n", simCtx->output_dir);
1152 ierr = PetscMkdir(simCtx->output_dir); CHKERRQ(ierr);
1153 }
1154 } else if(simCtx->exec_mode == EXEC_MODE_POSTPROCESSOR){
1155 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Preparing post-processing output directories ...\n");
1156
1157 PostProcessParams *pps = simCtx->pps;
1158 char path_buffer[PETSC_MAX_PATH_LEN];
1159
1160 const char *last_slash_euler = strrchr(pps->output_prefix, '/');
1161 if(last_slash_euler){
1162 size_t dir_len = last_slash_euler - pps->output_prefix;
1163 if(dir_len > 0){
1164 if(dir_len >= sizeof(path_buffer)) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Post-processing output prefix path is too long.");
1165 strncpy(path_buffer, pps->output_prefix, dir_len);
1166 path_buffer[dir_len] = '\0';
1167
1168 ierr = PetscTestDirectory(path_buffer, 'r', &exists); CHKERRQ(ierr);
1169 if (!exists){
1170 LOG_ALLOW(GLOBAL, LOG_INFO, "Creating post-processing Eulerian output directory: %s\n", path_buffer);
1171 ierr = PetscMkdirRecursive(path_buffer); CHKERRQ(ierr);
1172 }
1173 }
1174 }
1175
1176 // Particle output directory
1177 if(pps->outputParticles){
1178 const char *last_slash_particle = strrchr(pps->particle_output_prefix, '/');
1179 if(last_slash_particle){
1180 size_t dir_len = last_slash_particle - pps->particle_output_prefix;
1181 if(dir_len > 0){
1182 if(dir_len > sizeof(path_buffer)) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Post-processing particle output prefix path is too long.");
1183 strncpy(path_buffer, pps->particle_output_prefix, dir_len);
1184 path_buffer[dir_len] = '\0';
1185
1186 ierr = PetscTestDirectory(path_buffer, 'r', &exists); CHKERRQ(ierr);
1187
1188 if (!exists){
1189 LOG_ALLOW(GLOBAL, LOG_INFO, "Creating post-processing Particle output directory: %s\n", path_buffer);
1190 ierr = PetscMkdirRecursive(path_buffer); CHKERRQ(ierr);
1191 }
1192 }
1193 }
1194 }
1195
1196 // Statistics output directory
1197 if(pps->statistics_pipeline[0] != '\0'){
1198 const char *last_slash_stats = strrchr(pps->statistics_output_prefix, '/');
1199 if(last_slash_stats){
1200 size_t dir_len = last_slash_stats - pps->statistics_output_prefix;
1201 if(dir_len > 0){
1202 if(dir_len >= sizeof(path_buffer)) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Post-processing statistics output prefix path is too long.");
1203 strncpy(path_buffer, pps->statistics_output_prefix, dir_len);
1204 path_buffer[dir_len] = '\0';
1205
1206 ierr = PetscTestDirectory(path_buffer, 'r', &exists); CHKERRQ(ierr);
1207 if (!exists){
1208 LOG_ALLOW(GLOBAL, LOG_INFO, "Creating post-processing Statistics output directory: %s\n", path_buffer);
1209 ierr = PetscMkdirRecursive(path_buffer); CHKERRQ(ierr);
1210 }
1211 }
1212 }
1213 }
1214 }
1215 }
1216
1217 // Synchronize all processes before proceeding
1218 ierr = MPI_Barrier(PETSC_COMM_WORLD); CHKERRMPI(ierr);
1219
1220 LOG_ALLOW(GLOBAL, LOG_INFO, "--- Environment setup complete ---\n");
1221
1222 PetscFunctionReturn(0);
1223}
1224
1225#undef __FUNCT__
1226#define __FUNCT__ "AllocateContextHeirarchy"
1227/**
1228 * @brief Allocate the user-context objects required by every multigrid level.
1229 */
1230static PetscErrorCode AllocateContextHierarchy(SimCtx *simCtx)
1231{
1232 PetscErrorCode ierr;
1233 UserMG *usermg = &simCtx->usermg;
1234 MGCtx *mgctx;
1235 PetscInt nblk = simCtx->block_number;
1236 PetscBool found;
1237 PetscFunctionBeginUser;
1238 PROFILE_FUNCTION_BEGIN;
1239
1240 LOG_ALLOW(GLOBAL, LOG_INFO, "Allocating context hierarchy for %d levels and %d blocks...\n", simCtx->mglevels, nblk);
1241
1242 // Store the number of levels in the UserMG struct itself
1243 usermg->mglevels = simCtx->mglevels;
1244
1245 // --- 1. Allocate the array of MGCtx structs ---
1246 ierr = PetscMalloc(usermg->mglevels * sizeof(MGCtx), &usermg->mgctx); CHKERRQ(ierr);
1247 // Zero-initialize to ensure all pointers (especially packer) are NULL
1248 ierr = PetscMemzero(usermg->mgctx, usermg->mglevels * sizeof(MGCtx)); CHKERRQ(ierr);
1249 mgctx = usermg->mgctx;
1250 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Allocated MGCtx array of size %d.\n", simCtx->rank, usermg->mglevels);
1251
1252 // --- 2. Parse semi-coarsening options (logic from MG_Initial) ---
1253 // These flags determine if a dimension is coarsened in the multigrid hierarchy.
1254 PetscInt *isc, *jsc, *ksc;
1255 ierr = PetscMalloc3(nblk, &isc, nblk, &jsc, nblk, &ksc); CHKERRQ(ierr);
1256 // Set defaults to FALSE (full coarsening)
1257 for (PetscInt i = 0; i < nblk; ++i) {
1258 isc[i] = 0; jsc[i] = 0; ksc[i] = 0;
1259 }
1260
1261// Use a temporary variable for the 'count' argument to the parsing function.
1262 // This protects the original 'nblk' which is needed for the loop bounds.
1263 PetscInt n_opts_found = nblk;
1264 ierr = PetscOptionsGetIntArray(NULL, NULL, "-mg_i_semi", isc, &n_opts_found, &found); CHKERRQ(ierr);
1265
1266 n_opts_found = nblk; // Reset the temp variable before the next call
1267 ierr = PetscOptionsGetIntArray(NULL, NULL, "-mg_j_semi", jsc, &n_opts_found, &found); CHKERRQ(ierr);
1268
1269 n_opts_found = nblk; // Reset the temp variable before the next call
1270 ierr = PetscOptionsGetIntArray(NULL, NULL, "-mg_k_semi", ksc, &n_opts_found, &found); CHKERRQ(ierr);
1271
1272 // --- 3. Loop over levels and blocks to allocate UserCtx arrays ---
1273 for (PetscInt level = 0; level < simCtx->mglevels; level++) {
1274
1275 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Setting up MG Level %d...\n", simCtx->rank, level);
1276 // Allocate the array of UserCtx structs for this level
1277 ierr = PetscMalloc(nblk * sizeof(UserCtx), &mgctx[level].user); CHKERRQ(ierr);
1278 // It's good practice to zero out the memory to avoid uninitialized values
1279 ierr = PetscMemzero(mgctx[level].user, nblk * sizeof(UserCtx)); CHKERRQ(ierr);
1280 mgctx[level].thislevel = level;
1281
1282 for (PetscInt bi = 0; bi < nblk; bi++) {
1283 UserCtx *currentUser = &mgctx[level].user[bi];
1284 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Initializing UserCtx for Level %d, Block %d.\n", simCtx->rank, level, bi);
1285
1286 // --- CRITICAL STEP: Set the back-pointer to the master context ---
1287 currentUser->simCtx = simCtx;
1288
1289 // Initialize other per-context values
1290 currentUser->thislevel = level;
1291 currentUser->_this = bi; //
1292 currentUser->mglevels = usermg->mglevels;
1293
1294 // Assign semi-coarsening flags
1295 currentUser->isc = isc[bi];
1296 currentUser->jsc = jsc[bi];
1297 currentUser->ksc = ksc[bi];
1298
1299 // Link to finer/coarser contexts for multigrid operations
1300 if (level > 0) {
1301 currentUser->user_c = &mgctx[level-1].user[bi];
1302 LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG, "Rank %d: -> Linked to coarser context (user_c).\n", simCtx->rank);
1303 }
1304 if (level < usermg->mglevels - 1) {
1305 currentUser->user_f = &mgctx[level+1].user[bi];
1306 LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG, "Rank %d: -> Linked to finer context (user_f).\n", simCtx->rank);
1307 }
1308 }
1309 }
1310
1311 // Log a summary of the parsed flags on each rank.
1312 if (get_log_level() >= LOG_DEBUG && nblk > 0) {
1313 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Final semi-coarsening configuration view:\n", simCtx->rank);
1314 for (PetscInt bi = 0; bi < nblk; ++bi) {
1315 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Block %d: i-semi=%d, j-semi=%d, k-semi=%d\n", simCtx->rank, bi, isc[bi], jsc[bi], ksc[bi]);
1316 }
1317 }
1318
1319 // Clean up temporary arrays
1320 ierr = PetscFree3(isc, jsc, ksc); CHKERRQ(ierr);
1321
1322 LOG_ALLOW(GLOBAL, LOG_INFO, "Context hierarchy allocation complete.\n");
1323 PROFILE_FUNCTION_END;
1324 PetscFunctionReturn(0);
1325}
1326
1327#undef __FUNCT__
1328#define __FUNCT__ "SetupSolverParameters"
1329/**
1330 * @brief Configure solver tolerances, operators, and time-integration controls from the context.
1331 */
1332static PetscErrorCode SetupSolverParameters(SimCtx *simCtx){
1333
1334 PetscFunctionBeginUser;
1335 PROFILE_FUNCTION_BEGIN;
1336
1337 LOG_ALLOW(GLOBAL,LOG_INFO, " -- Setting up solver parameters -- .\n");
1338
1339 UserMG *usermg = &simCtx->usermg;
1340 MGCtx *mgctx = usermg->mgctx;
1341 PetscInt nblk = simCtx->block_number;
1342
1343 for (PetscInt level = usermg->mglevels-1; level >=0; level--) {
1344 for (PetscInt bi = 0; bi < nblk; bi++) {
1345 UserCtx *user = &mgctx[level].user[bi];
1346 LOG_ALLOW_SYNC(LOCAL, LOG_DEBUG, "Rank %d: Setting up parameters for level %d, block %d\n", simCtx->rank, level, bi);
1347
1348 user->assignedA = PETSC_FALSE;
1349 user->multinullspace = PETSC_FALSE;
1350 }
1351 }
1352 PROFILE_FUNCTION_END;
1353 PetscFunctionReturn(0);
1354}
1355
1356#undef __FUNCT__
1357#define __FUNCT__ "SetupGridAndSolvers"
1358/**
1359 * @brief Implementation of \ref SetupGridAndSolvers().
1360 * @details Full API contract (arguments, ownership, side effects) is documented with
1361 * the header declaration in `include/setup.h`.
1362 * @see SetupGridAndSolvers()
1363 */
1364PetscErrorCode SetupGridAndSolvers(SimCtx *simCtx)
1365{
1366 PetscErrorCode ierr;
1367 PetscFunctionBeginUser;
1368
1369 PROFILE_FUNCTION_BEGIN;
1370
1371 LOG_ALLOW(GLOBAL, LOG_INFO, "--- Starting Grid and Solvers Setup ---\n");
1372
1373 // Phase 1: Allocate the UserMG and UserCtx hierarchy
1374 ierr = AllocateContextHierarchy(simCtx); CHKERRQ(ierr);
1375
1376 ierr = DefineAllGridDimensions(simCtx); CHKERRQ(ierr);
1377 ierr = InitializeAllGridDMs(simCtx); CHKERRQ(ierr);
1378 ierr = AssignAllGridCoordinates(simCtx); CHKERRQ(ierr);
1379 ierr = CreateAndInitializeAllVectors(simCtx); CHKERRQ(ierr);
1380 ierr = SetupSolverParameters(simCtx); CHKERRQ(ierr);
1381
1382 // NOTE: CalculateAllGridMetrics is now called inside SetupBoundaryConditions (not here) to ensure:
1383 // 1. Boundary condition configuration data (boundary_faces) is available for periodic BC corrections
1384 // 2. Computed metrics are available for inlet/outlet area calculations
1385 // This resolves the circular dependency between BC setup and metric calculations.
1386
1387 LOG_ALLOW(GLOBAL, LOG_INFO, "--- Grid and Solvers Setup Complete ---\n");
1388
1389 PROFILE_FUNCTION_END;
1390 PetscFunctionReturn(0);
1391}
1392
1393
1394#undef __FUNCT__
1395#define __FUNCT__ "CreateAndInitializeAllVectors"
1396/**
1397 * @brief Internal helper implementation: `CreateAndInitializeAllVectors()`.
1398 * @details Local to this translation unit.
1399 */
1400PetscErrorCode CreateAndInitializeAllVectors(SimCtx *simCtx)
1401{
1402 PetscErrorCode ierr;
1403 UserMG *usermg = &simCtx->usermg;
1404 MGCtx *mgctx = usermg->mgctx;
1405 PetscInt nblk = simCtx->block_number;
1406
1407 PetscFunctionBeginUser;
1408
1409 PROFILE_FUNCTION_BEGIN;
1410
1411 LOG_ALLOW(GLOBAL, LOG_INFO, "Creating and initializing all simulation vectors...\n");
1412
1413 for (PetscInt level = usermg->mglevels-1; level >=0; level--) {
1414 for (PetscInt bi = 0; bi < nblk; bi++) {
1415 UserCtx *user = &mgctx[level].user[bi];
1416
1417 if(!user->da || !user->fda) {
1418 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "DMs not properly initialized in UserCtx before vector creation.");
1419 }
1420
1421 LOG_ALLOW_SYNC(LOCAL, LOG_DEBUG, "Rank %d: Creating vectors for level %d, block %d\n", simCtx->rank, level, bi);
1422
1423 // --- Group A: Primary Flow Fields (Global and Local) ---
1424 // These are the core solution variables.
1425 ierr = DMCreateGlobalVector(user->fda, &user->Ucont); CHKERRQ(ierr); ierr = VecSet(user->Ucont, 0.0); CHKERRQ(ierr);
1426 ierr = DMCreateGlobalVector(user->fda, &user->Ucat); CHKERRQ(ierr); ierr = VecSet(user->Ucat, 0.0); CHKERRQ(ierr);
1427 ierr = DMCreateGlobalVector(user->da, &user->P); CHKERRQ(ierr); ierr = VecSet(user->P, 0.0); CHKERRQ(ierr);
1428 ierr = DMCreateGlobalVector(user->da, &user->Nvert); CHKERRQ(ierr); ierr = VecSet(user->Nvert, 0.0); CHKERRQ(ierr);
1429
1430 ierr = DMCreateLocalVector(user->fda, &user->lUcont); CHKERRQ(ierr); ierr = VecSet(user->lUcont, 0.0); CHKERRQ(ierr);
1431 ierr = DMCreateLocalVector(user->fda, &user->lUcat); CHKERRQ(ierr); ierr = VecSet(user->lUcat, 0.0); CHKERRQ(ierr);
1432 ierr = DMCreateLocalVector(user->da, &user->lP); CHKERRQ(ierr); ierr = VecSet(user->lP, 0.0); CHKERRQ(ierr);
1433 ierr = DMCreateLocalVector(user->da, &user->lNvert); CHKERRQ(ierr); ierr = VecSet(user->lNvert, 0.0); CHKERRQ(ierr);
1434
1435 // -- Group A2: Derived Flow Fields (Global and Local) ---
1436 ierr = VecDuplicate(user->P,&user->Diffusivity); CHKERRQ(ierr); ierr = VecSet(user->Diffusivity, 0.0); CHKERRQ(ierr);
1437 ierr = VecDuplicate(user->lP,&user->lDiffusivity); CHKERRQ(ierr); ierr = VecSet(user->lDiffusivity, 0.0); CHKERRQ(ierr);
1438 ierr = VecDuplicate(user->Ucat,&user->DiffusivityGradient); CHKERRQ(ierr); ierr = VecSet(user->DiffusivityGradient, 0.0); CHKERRQ(ierr);
1439 ierr = VecDuplicate(user->lUcat,&user->lDiffusivityGradient); CHKERRQ(ierr); ierr = VecSet(user->lDiffusivityGradient, 0.0); CHKERRQ(ierr);
1440
1441 // -- Group B: Solver Work Vectors (Global and Local) ---
1442 ierr = VecDuplicate(user->P, &user->Phi); CHKERRQ(ierr); ierr = VecSet(user->Phi, 0.0); CHKERRQ(ierr);
1443 ierr = VecDuplicate(user->lP, &user->lPhi); CHKERRQ(ierr); ierr = VecSet(user->lPhi, 0.0); CHKERRQ(ierr);
1444
1445 // --- Group C: Time-Stepping & Workspace Fields (Finest Level Only) ---
1446 if (level == usermg->mglevels - 1) {
1447 ierr = VecDuplicate(user->Ucont, &user->Ucont_o); CHKERRQ(ierr); ierr = VecSet(user->Ucont_o, 0.0); CHKERRQ(ierr);
1448 ierr = VecDuplicate(user->Ucont, &user->Ucont_rm1); CHKERRQ(ierr); ierr = VecSet(user->Ucont_rm1, 0.0); CHKERRQ(ierr);
1449 ierr = VecDuplicate(user->Ucat, &user->Ucat_o); CHKERRQ(ierr); ierr = VecSet(user->Ucat_o, 0.0); CHKERRQ(ierr);
1450 ierr = VecDuplicate(user->P, &user->P_o); CHKERRQ(ierr); ierr = VecSet(user->P_o, 0.0); CHKERRQ(ierr);
1451 ierr = VecDuplicate(user->lUcont, &user->lUcont_o); CHKERRQ(ierr); ierr = VecSet(user->lUcont_o, 0.0); CHKERRQ(ierr);
1452 ierr = VecDuplicate(user->lUcont, &user->lUcont_rm1); CHKERRQ(ierr); ierr = VecSet(user->lUcont_rm1, 0.0); CHKERRQ(ierr);
1453 ierr = DMCreateLocalVector(user->da, &user->lNvert_o); CHKERRQ(ierr); ierr = VecSet(user->lNvert_o, 0.0); CHKERRQ(ierr);
1454 ierr = VecDuplicate(user->Nvert, &user->Nvert_o); CHKERRQ(ierr); ierr = VecSet(user->Nvert_o, 0.0); CHKERRQ(ierr);
1455 }
1456
1457 // --- Group D: Grid Metrics (Face-Centered) ---
1458 ierr = DMCreateGlobalVector(user->fda, &user->Csi); CHKERRQ(ierr); ierr = VecSet(user->Csi, 0.0); CHKERRQ(ierr);
1459 ierr = VecDuplicate(user->Csi, &user->Eta); CHKERRQ(ierr); ierr = VecSet(user->Eta, 0.0); CHKERRQ(ierr);
1460 ierr = VecDuplicate(user->Csi, &user->Zet); CHKERRQ(ierr); ierr = VecSet(user->Zet, 0.0); CHKERRQ(ierr);
1461 ierr = DMCreateGlobalVector(user->da, &user->Aj); CHKERRQ(ierr); ierr = VecSet(user->Aj, 0.0); CHKERRQ(ierr);
1462
1463 ierr = DMCreateLocalVector(user->fda, &user->lCsi); CHKERRQ(ierr); ierr = VecSet(user->lCsi, 0.0); CHKERRQ(ierr);
1464 ierr = VecDuplicate(user->lCsi, &user->lEta); CHKERRQ(ierr); ierr = VecSet(user->lEta, 0.0); CHKERRQ(ierr);
1465 ierr = VecDuplicate(user->lCsi, &user->lZet); CHKERRQ(ierr); ierr = VecSet(user->lZet, 0.0); CHKERRQ(ierr);
1466 ierr = DMCreateLocalVector(user->da, &user->lAj); CHKERRQ(ierr); ierr = VecSet(user->lAj, 0.0); CHKERRQ(ierr);
1467
1468
1469 // --- Group E: Grid Metrics (Face-Centered) ---
1470 // Vector metrics are duplicated from Csi (DOF=3, fda-based)
1471 ierr = VecDuplicate(user->Csi, &user->ICsi); CHKERRQ(ierr); ierr = VecSet(user->ICsi, 0.0); CHKERRQ(ierr);
1472 ierr = VecDuplicate(user->Csi, &user->IEta); CHKERRQ(ierr); ierr = VecSet(user->IEta, 0.0); CHKERRQ(ierr);
1473 ierr = VecDuplicate(user->Csi, &user->IZet); CHKERRQ(ierr); ierr = VecSet(user->IZet, 0.0); CHKERRQ(ierr);
1474 ierr = VecDuplicate(user->Csi, &user->JCsi); CHKERRQ(ierr); ierr = VecSet(user->JCsi, 0.0); CHKERRQ(ierr);
1475 ierr = VecDuplicate(user->Csi, &user->JEta); CHKERRQ(ierr); ierr = VecSet(user->JEta, 0.0); CHKERRQ(ierr);
1476 ierr = VecDuplicate(user->Csi, &user->JZet); CHKERRQ(ierr); ierr = VecSet(user->JZet, 0.0); CHKERRQ(ierr);
1477 ierr = VecDuplicate(user->Csi, &user->KCsi); CHKERRQ(ierr); ierr = VecSet(user->KCsi, 0.0); CHKERRQ(ierr);
1478 ierr = VecDuplicate(user->Csi, &user->KEta); CHKERRQ(ierr); ierr = VecSet(user->KEta, 0.0); CHKERRQ(ierr);
1479 ierr = VecDuplicate(user->Csi, &user->KZet); CHKERRQ(ierr); ierr = VecSet(user->KZet, 0.0); CHKERRQ(ierr);
1480 // Scalar metrics are duplicated from Aj (DOF=1, da-based)
1481 ierr = VecDuplicate(user->Aj, &user->IAj); CHKERRQ(ierr); ierr = VecSet(user->IAj, 0.0); CHKERRQ(ierr);
1482 ierr = VecDuplicate(user->Aj, &user->JAj); CHKERRQ(ierr); ierr = VecSet(user->JAj, 0.0); CHKERRQ(ierr);
1483 ierr = VecDuplicate(user->Aj, &user->KAj); CHKERRQ(ierr); ierr = VecSet(user->KAj, 0.0); CHKERRQ(ierr);
1484
1485 ierr = VecDuplicate(user->lCsi, &user->lICsi); CHKERRQ(ierr); ierr = VecSet(user->lICsi, 0.0); CHKERRQ(ierr);
1486 ierr = VecDuplicate(user->lCsi, &user->lIEta); CHKERRQ(ierr); ierr = VecSet(user->lIEta, 0.0); CHKERRQ(ierr);
1487 ierr = VecDuplicate(user->lCsi, &user->lIZet); CHKERRQ(ierr); ierr = VecSet(user->lIZet, 0.0); CHKERRQ(ierr);
1488 ierr = VecDuplicate(user->lCsi, &user->lJCsi); CHKERRQ(ierr); ierr = VecSet(user->lJCsi, 0.0); CHKERRQ(ierr);
1489 ierr = VecDuplicate(user->lCsi, &user->lJEta); CHKERRQ(ierr); ierr = VecSet(user->lJEta, 0.0); CHKERRQ(ierr);
1490 ierr = VecDuplicate(user->lCsi, &user->lJZet); CHKERRQ(ierr); ierr = VecSet(user->lJZet, 0.0); CHKERRQ(ierr);
1491 ierr = VecDuplicate(user->lCsi, &user->lKCsi); CHKERRQ(ierr); ierr = VecSet(user->lKCsi, 0.0); CHKERRQ(ierr);
1492 ierr = VecDuplicate(user->lCsi, &user->lKEta); CHKERRQ(ierr); ierr = VecSet(user->lKEta, 0.0); CHKERRQ(ierr);
1493 ierr = VecDuplicate(user->lCsi, &user->lKZet); CHKERRQ(ierr); ierr = VecSet(user->lKZet, 0.0); CHKERRQ(ierr);
1494
1495 ierr = VecDuplicate(user->lAj, &user->lIAj); CHKERRQ(ierr); ierr = VecSet(user->lIAj, 0.0); CHKERRQ(ierr);
1496 ierr = VecDuplicate(user->lAj, &user->lJAj); CHKERRQ(ierr); ierr = VecSet(user->lJAj, 0.0); CHKERRQ(ierr);
1497 ierr = VecDuplicate(user->lAj, &user->lKAj); CHKERRQ(ierr); ierr = VecSet(user->lKAj, 0.0); CHKERRQ(ierr);
1498
1499 // --- Group F: Cell/Face Center Coordinates and Grid Spacing ---
1500 ierr = DMCreateGlobalVector(user->fda, &user->Cent); CHKERRQ(ierr); ierr = VecSet(user->Cent, 0.0); CHKERRQ(ierr);
1501 ierr = DMCreateLocalVector(user->fda, &user->lCent); CHKERRQ(ierr); ierr = VecSet(user->lCent, 0.0); CHKERRQ(ierr);
1502
1503 ierr = VecDuplicate(user->Cent, &user->GridSpace); CHKERRQ(ierr); ierr = VecSet(user->GridSpace, 0.0); CHKERRQ(ierr);
1504 ierr = VecDuplicate(user->lCent, &user->lGridSpace); CHKERRQ(ierr); ierr = VecSet(user->lGridSpace, 0.0); CHKERRQ(ierr);
1505
1506 ierr = VecDuplicate(user->Cent, &user->Centx); CHKERRQ(ierr); ierr = VecSet(user->Centx, 0.0); CHKERRQ(ierr);
1507 ierr = VecDuplicate(user->Cent, &user->Centy); CHKERRQ(ierr); ierr = VecSet(user->Centy, 0.0); CHKERRQ(ierr);
1508 ierr = VecDuplicate(user->Cent, &user->Centz); CHKERRQ(ierr); ierr = VecSet(user->Centz, 0.0); CHKERRQ(ierr);
1509 ierr = VecDuplicate(user->lCent, &user->lCentx); CHKERRQ(ierr); ierr = VecSet(user->lCentx, 0.0); CHKERRQ(ierr);
1510 ierr = VecDuplicate(user->lCent, &user->lCenty); CHKERRQ(ierr); ierr = VecSet(user->lCenty, 0.0); CHKERRQ(ierr);
1511 ierr = VecDuplicate(user->lCent, &user->lCentz); CHKERRQ(ierr); ierr = VecSet(user->lCentz, 0.0); CHKERRQ(ierr);
1512
1513 if(level == usermg->mglevels -1){
1514 // --- Group G: Turbulence Models (Finest Level Only) ---
1515 if (simCtx->les || simCtx->rans) {
1516 ierr = DMCreateGlobalVector(user->da, &user->Nu_t); CHKERRQ(ierr); ierr = VecSet(user->Nu_t, 0.0); CHKERRQ(ierr);
1517 ierr = DMCreateLocalVector(user->da, &user->lNu_t); CHKERRQ(ierr); ierr = VecSet(user->lNu_t, 0.0); CHKERRQ(ierr);
1518 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Turbulence viscosity (Nu_t) vectors created for LES/RANS model.\n");
1519 if(simCtx->les){
1520 ierr = DMCreateGlobalVector(user->da,&user->CS); CHKERRQ(ierr); ierr = VecSet(user->CS,0.0); CHKERRQ(ierr);
1521 ierr = DMCreateLocalVector(user->da,&user->lCs); CHKERRQ(ierr); ierr = VecSet(user->lCs,0.0); CHKERRQ(ierr);
1522 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Smagorinsky constant (CS) vectors created for LES model.\n");
1523 }
1524
1525 if(simCtx->wallfunction){
1526 ierr = DMCreateLocalVector(user->fda,&user->lFriction_Velocity); CHKERRQ(ierr); ierr = VecSet(user->lFriction_Velocity,0.0);
1527 }
1528 // Add K_Omega etc. here as needed
1529
1530 // Note: Add any other vectors from the legacy MG_Initial here as needed.
1531 // For example: Rhs, Forcing, turbulence Vecs (K_Omega, Nu_t)...
1532
1533 }
1534 // --- Group H: Particle Methods
1535 if(simCtx->np>0){
1536 ierr = DMCreateGlobalVector(user->da,&user->ParticleCount); CHKERRQ(ierr); ierr = VecSet(user->ParticleCount,0.0); CHKERRQ(ierr);
1537 ierr = DMCreateLocalVector(user->da,&user->lParticleCount); CHKERRQ(ierr); ierr = VecSet(user->lParticleCount,0.0); CHKERRQ(ierr);
1538 // Scalar field to hold particle scalar property (e.g., temperature, concentration)
1539 ierr = DMCreateGlobalVector(user->da,&user->Psi); CHKERRQ(ierr); ierr = VecSet(user->Psi,0.0); CHKERRQ(ierr);
1540 ierr = DMCreateLocalVector(user->da,&user->lPsi); CHKERRQ(ierr); ierr = VecSet(user->lPsi,0.0); CHKERRQ(ierr);
1541 LOG_ALLOW(GLOBAL,LOG_DEBUG,"ParticleCount & Scalar(Psi) created for %d particles.\n",simCtx->np);
1542 }
1543 }
1544 // --- Group I: Boundary Condition vectors ---
1545 ierr = DMCreateGlobalVector(user->fda, &user->Bcs.Ubcs); CHKERRQ(ierr);
1546 ierr = VecSet(user->Bcs.Ubcs, 0.0); CHKERRQ(ierr);
1547 ierr = DMCreateGlobalVector(user->fda, &user->Bcs.Uch); CHKERRQ(ierr);
1548 ierr = VecSet(user->Bcs.Uch, 0.0); CHKERRQ(ierr);
1549
1550 // --- Group J: Coarsest-Level Immersed-Boundary Workspace ---
1551 // Sized from the coarsest-level DMDA, which is fixed once the DMs exist,
1552 // so this is allocated once here rather than on every Poisson solve.
1553 // FullyBlocked() guards every read with its own Blocked[] flag, so the
1554 // buffer carries no state between solves.
1555 if (level == 0 && simCtx->immersed) {
1556 ierr = PetscCalloc1((size_t)(user->info.mx * user->info.my * 2), &user->KSKE); CHKERRQ(ierr);
1557 LOG_ALLOW(LOCAL, LOG_DEBUG, "Coarsest-level KSKE workspace allocated for immersed boundaries.\n");
1558 }
1559
1560 // --- Group L: Field-Statistics Accumulators (Finest Level Only) ---
1561 // Config-counted, like the convergence-state reference fields: the window
1562 // count is resolved before this factory runs, and each accumulator is
1563 // duplicated from a vector created above so it inherits DM and layout.
1564 if (level == usermg->mglevels - 1 && FieldStatisticsIsActive(simCtx)) {
1565 ierr = PetscCalloc1((size_t)simCtx->fieldStatisticsWindowCount,
1566 &user->fieldStatisticsStorage); CHKERRQ(ierr);
1567 for (PetscInt w = 0; w < simCtx->fieldStatisticsWindowCount; ++w) {
1568 ierr = PicurvWindowStorageCreate(user,
1569 &simCtx->fieldStatisticsWindows[w].definition,
1570 &user->fieldStatisticsStorage[w]); CHKERRQ(ierr);
1571 }
1572 LOG_ALLOW(LOCAL, LOG_INFO,
1573 "Allocated accumulators for %d statistics window(s).\n",
1574 simCtx->fieldStatisticsWindowCount);
1575 }
1576
1577 // --- Group K: Corner-Staging Workspace (Finest Level Only) ---
1578 // One pair per block size. Previously these were created lazily inside
1579 // the interpolation routine and rebuilt whenever the block size changed;
1580 // allocating both here removes that churn and the cached-size check.
1581 if (level == usermg->mglevels - 1) {
1582 ierr = DMCreateGlobalVector(user->da, &user->CellScalarAtCorner); CHKERRQ(ierr);
1583 ierr = VecSet(user->CellScalarAtCorner, 0.0); CHKERRQ(ierr);
1584 ierr = DMCreateLocalVector(user->da, &user->lCellScalarAtCorner); CHKERRQ(ierr);
1585 ierr = VecSet(user->lCellScalarAtCorner, 0.0); CHKERRQ(ierr);
1586 ierr = DMCreateGlobalVector(user->fda, &user->CellVectorAtCorner); CHKERRQ(ierr);
1587 ierr = VecSet(user->CellVectorAtCorner, 0.0); CHKERRQ(ierr);
1588 ierr = DMCreateLocalVector(user->fda, &user->lCellVectorAtCorner); CHKERRQ(ierr);
1589 ierr = VecSet(user->lCellVectorAtCorner, 0.0); CHKERRQ(ierr);
1590 }
1591
1592 /* --- Group N: Statistics Output Staging (Finest Level Only) ---
1593 * The same kind of object as the corner workspace: a named buffer the
1594 * shared kernels can address, rather than simulation state. It exists only
1595 * to carry derived statistics to the nodal and VTK paths, so a run without
1596 * statistics allocates none of it. */
1597 if (level == usermg->mglevels - 1 && FieldStatisticsIsActive(simCtx)) {
1598 ierr = DMCreateGlobalVector(user->da, &user->PostScalar); CHKERRQ(ierr);
1599 ierr = VecSet(user->PostScalar, 0.0); CHKERRQ(ierr);
1600 ierr = DMCreateLocalVector(user->da, &user->lPostScalar); CHKERRQ(ierr);
1601 ierr = VecSet(user->lPostScalar, 0.0); CHKERRQ(ierr);
1602 ierr = DMCreateGlobalVector(user->da, &user->PostScalarNodal); CHKERRQ(ierr);
1603 ierr = VecSet(user->PostScalarNodal, 0.0); CHKERRQ(ierr);
1604 ierr = DMCreateGlobalVector(user->fda, &user->PostVector); CHKERRQ(ierr);
1605 ierr = VecSet(user->PostVector, 0.0); CHKERRQ(ierr);
1606 ierr = DMCreateLocalVector(user->fda, &user->lPostVector); CHKERRQ(ierr);
1607 ierr = VecSet(user->lPostVector, 0.0); CHKERRQ(ierr);
1608 ierr = DMCreateGlobalVector(user->fda, &user->PostVectorNodal); CHKERRQ(ierr);
1609 ierr = VecSet(user->PostVectorNodal, 0.0); CHKERRQ(ierr);
1610 }
1611
1612 if(level == usermg->mglevels - 1){
1613 if(simCtx->exec_mode == EXEC_MODE_POSTPROCESSOR){
1614 LOG_ALLOW(LOCAL, LOG_DEBUG, "Post-processor mode detected. Allocating derived field vectors.\n");
1615
1616 ierr = VecDuplicate(user->P, &user->P_nodal); CHKERRQ(ierr);
1617 ierr = VecSet(user->P_nodal, 0.0); CHKERRQ(ierr);
1618
1619 ierr = VecDuplicate(user->Ucat, &user->Ucat_nodal); CHKERRQ(ierr);
1620 ierr = VecSet(user->Ucat_nodal, 0.0); CHKERRQ(ierr);
1621
1622 ierr = VecDuplicate(user->P, &user->Qcrit); CHKERRQ(ierr);
1623 ierr = VecSet(user->Qcrit, 0.0); CHKERRQ(ierr);
1624
1625 LOG_ALLOW(LOCAL, LOG_DEBUG, "Derived field vectors P_nodal, Ucat_nodal, and Qcrit created.\n");
1626
1627 if(simCtx->np>0){
1628 ierr = VecDuplicate(user->Psi, &user->Psi_nodal); CHKERRQ(ierr);
1629 ierr = VecSet(user->Psi_nodal, 0.0); CHKERRQ(ierr);
1630
1631 LOG_ALLOW(LOCAL, LOG_DEBUG, "Derived field vector Psi_nodal created for particle scalar property.\n");
1632
1633 }
1634 }else{
1635 user->P_nodal = NULL;
1636 user->Ucat_nodal = NULL;
1637 user->Qcrit = NULL;
1638 user->Psi_nodal = NULL;
1639 }
1640 }
1641
1642 }
1643}
1644
1645 /* Config-counted vectors belong to this factory too. The convergence state
1646 * allocates a Vec array whose length comes from configuration, so it is
1647 * invoked here rather than by the caller, keeping every vector this run
1648 * owns created in one place. It depends only on option-parsed configuration
1649 * and on the vectors created above. */
1650 ierr = InitializeSolutionConvergenceState(simCtx); CHKERRQ(ierr);
1651
1652 LOG_ALLOW(GLOBAL, LOG_INFO, "All simulation vectors created and initialized.\n");
1653
1654 PROFILE_FUNCTION_END;
1655 PetscFunctionReturn(0);
1656}
1657
1658#undef __FUNCT__
1659#define __FUNCT__ "RepairPeriodicNormalFaceGhosts"
1660/**
1661 * @brief Repairs the adjacent normal ghost layer for periodic face-staggered data.
1662 *
1663 * PETSc wraps every component with cell-style indexing. A face family instead
1664 * needs its adjacent normal ghost shifted by one additional physical face. With
1665 * the width-three periodic DMDA, the required value is available in the deeper
1666 * PETSc ghost at -3 or n+2. Tangential ghosts retain PETSc's native wraparound.
1667 */
1668static PetscErrorCode RepairPeriodicNormalFaceGhosts(UserCtx *user, DM dm, Vec local_vec,
1669 PetscInt dof, char face_direction,
1670 PetscBool component_staggered)
1671{
1672 DMDALocalInfo info;
1673 PetscInt xs, xe, ys, ye, zs, ze;
1674 PetscInt gxs, gxe, gys, gye, gzs, gze;
1675 PetscInt mx, my, mz;
1676
1677 PetscFunctionBeginUser;
1678 if (!face_direction && !component_staggered) PetscFunctionReturn(0);
1679
1680 PetscCall(DMDAGetLocalInfo(dm, &info));
1681 xs = info.xs; xe = info.xs + info.xm;
1682 ys = info.ys; ye = info.ys + info.ym;
1683 zs = info.zs; ze = info.zs + info.zm;
1684 gxs = info.gxs; gxe = info.gxs + info.gxm;
1685 gys = info.gys; gye = info.gys + info.gym;
1686 gzs = info.gzs; gze = info.gzs + info.gzm;
1687 mx = info.mx; my = info.my; mz = info.mz;
1688
1689 if (component_staggered) {
1690 Cmpnts ***array;
1691 PetscCall(DMDAVecGetArray(dm, local_vec, &array));
1692
1693 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && xs == 0) {
1694 PetscCheck(gxs <= -3, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
1695 "Periodic Ucont.x ghost repair requires DMDA stencil width at least 3.");
1696 for (PetscInt k = gzs; k < gze; k++) for (PetscInt j = gys; j < gye; j++)
1697 array[k][j][-1].x = array[k][j][-3].x;
1698 }
1699 if (user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC && xe == mx) {
1700 PetscCheck(gxe > mx + 2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
1701 "Periodic Ucont.x ghost repair requires DMDA stencil width at least 3.");
1702 for (PetscInt k = gzs; k < gze; k++) for (PetscInt j = gys; j < gye; j++)
1703 array[k][j][mx].x = array[k][j][mx + 2].x;
1704 }
1705 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && ys == 0) {
1706 PetscCheck(gys <= -3, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
1707 "Periodic Ucont.y ghost repair requires DMDA stencil width at least 3.");
1708 for (PetscInt k = gzs; k < gze; k++) for (PetscInt i = gxs; i < gxe; i++)
1709 array[k][-1][i].y = array[k][-3][i].y;
1710 }
1711 if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC && ye == my) {
1712 PetscCheck(gye > my + 2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
1713 "Periodic Ucont.y ghost repair requires DMDA stencil width at least 3.");
1714 for (PetscInt k = gzs; k < gze; k++) for (PetscInt i = gxs; i < gxe; i++)
1715 array[k][my][i].y = array[k][my + 2][i].y;
1716 }
1717 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && zs == 0) {
1718 PetscCheck(gzs <= -3, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
1719 "Periodic Ucont.z ghost repair requires DMDA stencil width at least 3.");
1720 for (PetscInt j = gys; j < gye; j++) for (PetscInt i = gxs; i < gxe; i++)
1721 array[-1][j][i].z = array[-3][j][i].z;
1722 }
1723 if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC && ze == mz) {
1724 PetscCheck(gze > mz + 2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
1725 "Periodic Ucont.z ghost repair requires DMDA stencil width at least 3.");
1726 for (PetscInt j = gys; j < gye; j++) for (PetscInt i = gxs; i < gxe; i++)
1727 array[mz][j][i].z = array[mz + 2][j][i].z;
1728 }
1729
1730 PetscCall(DMDAVecRestoreArray(dm, local_vec, &array));
1731 PetscFunctionReturn(0);
1732 }
1733
1734 if (dof == 1) {
1735 PetscReal ***array;
1736 PetscCall(DMDAVecGetArray(dm, local_vec, &array));
1737
1738 if (face_direction == 'i') {
1739 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && xs == 0) {
1740 PetscCheck(gxs <= -3, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
1741 "Periodic I-face ghost repair requires DMDA stencil width at least 3.");
1742 for (PetscInt k = gzs; k < gze; k++) for (PetscInt j = gys; j < gye; j++)
1743 array[k][j][-1] = array[k][j][-3];
1744 }
1745 if (user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC && xe == mx) {
1746 PetscCheck(gxe > mx + 2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
1747 "Periodic I-face ghost repair requires DMDA stencil width at least 3.");
1748 for (PetscInt k = gzs; k < gze; k++) for (PetscInt j = gys; j < gye; j++)
1749 array[k][j][mx] = array[k][j][mx + 2];
1750 }
1751 } else if (face_direction == 'j') {
1752 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && ys == 0) {
1753 PetscCheck(gys <= -3, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
1754 "Periodic J-face ghost repair requires DMDA stencil width at least 3.");
1755 for (PetscInt k = gzs; k < gze; k++) for (PetscInt i = gxs; i < gxe; i++)
1756 array[k][-1][i] = array[k][-3][i];
1757 }
1758 if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC && ye == my) {
1759 PetscCheck(gye > my + 2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
1760 "Periodic J-face ghost repair requires DMDA stencil width at least 3.");
1761 for (PetscInt k = gzs; k < gze; k++) for (PetscInt i = gxs; i < gxe; i++)
1762 array[k][my][i] = array[k][my + 2][i];
1763 }
1764 } else if (face_direction == 'k') {
1765 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && zs == 0) {
1766 PetscCheck(gzs <= -3, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
1767 "Periodic K-face ghost repair requires DMDA stencil width at least 3.");
1768 for (PetscInt j = gys; j < gye; j++) for (PetscInt i = gxs; i < gxe; i++)
1769 array[-1][j][i] = array[-3][j][i];
1770 }
1771 if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC && ze == mz) {
1772 PetscCheck(gze > mz + 2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
1773 "Periodic K-face ghost repair requires DMDA stencil width at least 3.");
1774 for (PetscInt j = gys; j < gye; j++) for (PetscInt i = gxs; i < gxe; i++)
1775 array[mz][j][i] = array[mz + 2][j][i];
1776 }
1777 }
1778
1779 PetscCall(DMDAVecRestoreArray(dm, local_vec, &array));
1780 } else {
1781 Cmpnts ***array;
1782 PetscCall(DMDAVecGetArray(dm, local_vec, &array));
1783
1784 if (face_direction == 'i') {
1785 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && xs == 0) {
1786 PetscCheck(gxs <= -3, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
1787 "Periodic I-face ghost repair requires DMDA stencil width at least 3.");
1788 for (PetscInt k = gzs; k < gze; k++) for (PetscInt j = gys; j < gye; j++)
1789 array[k][j][-1] = array[k][j][-3];
1790 }
1791 if (user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC && xe == mx) {
1792 PetscCheck(gxe > mx + 2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
1793 "Periodic I-face ghost repair requires DMDA stencil width at least 3.");
1794 for (PetscInt k = gzs; k < gze; k++) for (PetscInt j = gys; j < gye; j++)
1795 array[k][j][mx] = array[k][j][mx + 2];
1796 }
1797 } else if (face_direction == 'j') {
1798 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && ys == 0) {
1799 PetscCheck(gys <= -3, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
1800 "Periodic J-face ghost repair requires DMDA stencil width at least 3.");
1801 for (PetscInt k = gzs; k < gze; k++) for (PetscInt i = gxs; i < gxe; i++)
1802 array[k][-1][i] = array[k][-3][i];
1803 }
1804 if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC && ye == my) {
1805 PetscCheck(gye > my + 2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
1806 "Periodic J-face ghost repair requires DMDA stencil width at least 3.");
1807 for (PetscInt k = gzs; k < gze; k++) for (PetscInt i = gxs; i < gxe; i++)
1808 array[k][my][i] = array[k][my + 2][i];
1809 }
1810 } else if (face_direction == 'k') {
1811 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && zs == 0) {
1812 PetscCheck(gzs <= -3, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
1813 "Periodic K-face ghost repair requires DMDA stencil width at least 3.");
1814 for (PetscInt j = gys; j < gye; j++) for (PetscInt i = gxs; i < gxe; i++)
1815 array[-1][j][i] = array[-3][j][i];
1816 }
1817 if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC && ze == mz) {
1818 PetscCheck(gze > mz + 2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
1819 "Periodic K-face ghost repair requires DMDA stencil width at least 3.");
1820 for (PetscInt j = gys; j < gye; j++) for (PetscInt i = gxs; i < gxe; i++)
1821 array[mz][j][i] = array[mz + 2][j][i];
1822 }
1823 }
1824
1825 PetscCall(DMDAVecRestoreArray(dm, local_vec, &array));
1826 }
1827
1828 PetscFunctionReturn(0);
1829}
1830
1831#undef __FUNCT__
1832#define __FUNCT__ "UpdateLocalGhosts"
1833/**
1834 * @brief Updates a catalogued field's local ghost representation.
1835 * @details Resolves the typed field view, performs the existing PETSc scatter,
1836 * and applies the established periodic normal-face repair when needed.
1837 */
1838PetscErrorCode UpdateLocalGhosts(UserCtx *user, FieldId field_id)
1839{
1840 PetscErrorCode ierr;
1841 PetscMPIInt rank;
1842 FieldView field_view;
1843 const char *field_name;
1844 Vec globalVec;
1845 Vec localVec;
1846 DM dm;
1847 PetscInt dof;
1848 char face_direction = '\0';
1849 PetscBool component_staggered = PETSC_FALSE;
1850
1851 PetscFunctionBeginUser; // Use User version for application code
1852 PROFILE_FUNCTION_BEGIN;
1853 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
1854 ierr = FieldGetView(user, field_id, &field_view); CHKERRQ(ierr);
1855 PetscCheck((field_view.descriptor->capabilities & FIELD_CAPABILITY_GHOST_UPDATE) != 0u,
1856 PETSC_COMM_SELF, PETSC_ERR_SUP,
1857 "Field '%s' does not support ghost updates.",
1858 field_view.descriptor->canonical_name);
1859 field_name = field_view.descriptor->canonical_name;
1860 globalVec = field_view.global_vec;
1861 localVec = field_view.local_vec;
1862 dm = field_view.dm;
1863 dof = field_view.descriptor->dof;
1864
1865 switch (field_view.descriptor->sync_class) {
1866 case FIELD_SYNC_I_FACE:
1867 face_direction = 'i';
1868 break;
1869 case FIELD_SYNC_J_FACE:
1870 face_direction = 'j';
1871 break;
1872 case FIELD_SYNC_K_FACE:
1873 face_direction = 'k';
1874 break;
1875 case FIELD_SYNC_COMPONENT_STAGGERED:
1876 component_staggered = PETSC_TRUE;
1877 break;
1878 case FIELD_SYNC_STANDARD:
1879 break;
1880 default:
1881 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB,
1882 "Field '%s' has an invalid ghost synchronization class.", field_name);
1883 }
1884
1885 LOG_ALLOW(GLOBAL, LOG_INFO, "Rank %d: Starting ghost update for field '%s'.\n", rank, field_name);
1886
1887 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Identified components for '%s': DM=%p, GlobalVec=%p, LocalVec=%p.\n",
1888 rank, field_name, (void*)dm, (void*)globalVec, (void*)localVec);
1889
1890 // --- 3. Optional Debugging: Norm Before Update ---
1891 // Use your logging convention check
1892 // if (get_log_level() >= LOG_LEVEL_DEBUG && is_function_allowed("UpdateLocalGhosts")) { // Example check
1893 if(get_log_level() == LOG_DEBUG && is_function_allowed(__func__)){
1894 PetscReal norm_global_before;
1895 ierr = VecNorm(globalVec, NORM_INFINITY, &norm_global_before); CHKERRQ(ierr);
1896 LOG_ALLOW(GLOBAL, LOG_INFO,"Max norm '%s' (Global) BEFORE Ghost Update: %g\n", field_name, norm_global_before);
1897 // Optional: Norm of local vector before update (might contain old ghost values)
1898 // PetscReal norm_local_before;
1899 // ierr = VecNorm(localVec, NORM_INFINITY, &norm_local_before); CHKERRQ(ierr);
1900 // LOG_ALLOW(GLOBAL, LOG_DEBUG,"Max norm '%s' (Local) BEFORE Ghost Update: %g\n", field_name, norm_local_before);
1901 }
1902
1903 // --- 4. Perform the Global-to-Local Transfer (Ghost Update) ---
1904 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Calling DMGlobalToLocalBegin/End for '%s'.\n", rank, field_name);
1905 ierr = DMGlobalToLocalBegin(dm, globalVec, INSERT_VALUES, localVec); CHKERRQ(ierr);
1906 ierr = DMGlobalToLocalEnd(dm, globalVec, INSERT_VALUES, localVec); CHKERRQ(ierr);
1907 ierr = RepairPeriodicNormalFaceGhosts(user, dm, localVec, dof, face_direction,
1908 component_staggered); CHKERRQ(ierr);
1909 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Completed DMGlobalToLocalBegin/End for '%s'.\n", rank, field_name);
1910
1911 // --- 5. Optional Debugging: Norm After Update ---
1912 // Use your logging convention check
1913 // if (get_log_level() >= LOG_LEVEL_DEBUG && is_function_allowed("UpdateLocalGhosts")) { // Example check
1914 if(get_log_level() == LOG_DEBUG && is_function_allowed(__func__)){ // Using your specific check
1915 PetscReal norm_local_after;
1916 ierr = VecNorm(localVec, NORM_INFINITY, &norm_local_after); CHKERRQ(ierr);
1917 LOG_ALLOW(GLOBAL, LOG_INFO,"Max norm '%s' (Local) AFTER Ghost Update: %g\n", field_name, norm_local_after);
1918
1919 // --- 6. Optional Debugging: Specific Point Checks (Example for Ucat on Rank 0/1) ---
1920 // (Keep this conditional if it's only for specific debug scenarios)
1921 if (field_id == FIELD_ID_UCAT) { // Only do detailed checks for Ucat for now
1922 PetscMPIInt rank_test;
1923 MPI_Comm_rank(PETSC_COMM_WORLD, &rank_test);
1924
1925 // Get Local Info needed for indexing checks
1926 DMDALocalInfo info_check;
1927 ierr = DMDAGetLocalInfo(dm, &info_check); CHKERRQ(ierr); // Use the correct dm
1928
1929 // Buffer for array pointer
1930 Cmpnts ***lUcat_arr_test = NULL;
1931 PetscErrorCode ierr_test = 0;
1932
1933 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Testing '%s' access immediately after ghost update...\n", rank_test, field_name);
1934 ierr_test = DMDAVecGetArrayDOFRead(dm, localVec, &lUcat_arr_test); // Use correct dm and localVec
1935
1936 if (ierr_test) {
1937 LOG_ALLOW(LOCAL, LOG_ERROR, "Rank %d: ERROR %d getting '%s' array after ghost update!\n", rank_test, ierr_test, field_name);
1938 } else if (!lUcat_arr_test) {
1939 LOG_ALLOW(LOCAL, LOG_ERROR, "Rank %d: ERROR NULL pointer getting '%s' array after ghost update!\n", rank_test, field_name);
1940 }
1941 else {
1942 // Check owned interior point (e.g., first interior point)
1943 PetscInt k_int = info_check.zs + (info_check.zm > 1 ? 1 : 0); // Global k index (at least zs+1 if possible)
1944 PetscInt j_int = info_check.ys + (info_check.ym > 1 ? 1 : 0); // Global j index
1945 PetscInt i_int = info_check.xs + (info_check.xm > 1 ? 1 : 0); // Global i index
1946 // Ensure indices are within global bounds if domain is very small
1947 //if (k_int >= info_check.mz-1) k_int = info_check.mz-2; if (k_int < 1) k_int = 1;
1948 //if (j_int >= info_check.my-1) j_int = info_check.my-2; if (j_int < 1) j_int = 1;
1949 // if (i_int >= info_check.mx-1) i_int = info_check.mx-2; if (i_int < 1) i_int = 1;
1950 // clamp k_int to [1 .. mz-2]
1951 if (k_int >= info_check.mz - 1) {
1952 k_int = info_check.mz - 2;
1953 }
1954 if (k_int < 1) {
1955 k_int = 1;
1956 }
1957
1958 // clamp j_int to [1 .. my-2]
1959 if (j_int >= info_check.my - 1) {
1960 j_int = info_check.my - 2;
1961 }
1962 if (j_int < 1) {
1963 j_int = 1;
1964 }
1965
1966 // clamp i_int to [1 .. mx-2]
1967 if (i_int >= info_check.mx - 1) {
1968 i_int = info_check.mx - 2;
1969 }
1970 if (i_int < 1) {
1971 i_int = 1;
1972 }
1973
1974 // Only attempt read if indices are actually owned (relevant for multi-rank)
1975 if (k_int >= info_check.zs && k_int < info_check.zs + info_check.zm &&
1976 j_int >= info_check.ys && j_int < info_check.ys + info_check.ym &&
1977 i_int >= info_check.xs && i_int < info_check.xs + info_check.xm)
1978 {
1979 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Attempting test read OWNED INTERIOR [%d][%d][%d] (Global)\n", rank_test, k_int, j_int, i_int);
1980 Cmpnts test_val_owned_interior = lUcat_arr_test[k_int][j_int][i_int];
1981 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: SUCCESS reading owned interior: x=%g\n", rank_test, test_val_owned_interior.x);
1982 } else {
1983 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Skipping interior test read for non-owned index [%d][%d][%d].\n", rank_test, k_int, j_int, i_int);
1984 }
1985
1986
1987 // Check owned boundary point (e.g., first owned point)
1988 PetscInt k_bnd = info_check.zs; // Global k index
1989 PetscInt j_bnd = info_check.ys; // Global j index
1990 PetscInt i_bnd = info_check.xs; // Global i index
1991 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Attempting test read OWNED BOUNDARY [%d][%d][%d] (Global)\n", rank_test, k_bnd, j_bnd, i_bnd);
1992 Cmpnts test_val_owned_boundary = lUcat_arr_test[k_bnd][j_bnd][i_bnd];
1993 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: SUCCESS reading owned boundary: x=%g\n", rank_test, test_val_owned_boundary.x);
1994
1995
1996 // Check ghost point (e.g., one layer below in k, if applicable)
1997 if (info_check.zs > 0) { // Only if there's a rank below
1998 PetscInt k_ghost = info_check.zs - 1;
1999 PetscInt j_ghost = info_check.ys; // Use start of owned y, simple example
2000 PetscInt i_ghost = info_check.xs; // Use start of owned x, simple example
2001 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Attempting test read GHOST [%d][%d][%d] (Global)\n", rank_test, k_ghost, j_ghost, i_ghost);
2002 Cmpnts test_val_ghost = lUcat_arr_test[k_ghost][j_ghost][i_ghost];
2003 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: SUCCESS reading ghost: x=%g\n", rank_test, test_val_ghost.x);
2004 } else {
2005 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Skipping ghost test read (zs=0).\n", rank_test);
2006 }
2007
2008 // Restore the array
2009 ierr_test = DMDAVecRestoreArrayDOFRead(dm, localVec, &lUcat_arr_test);
2010 if(ierr_test){ LOG_ALLOW(LOCAL, LOG_ERROR, "Rank %d: ERROR %d restoring '%s' array after test read!\n", rank_test, ierr_test, field_name); }
2011 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Finished testing '%s' access.\n", rank_test, field_name);
2012 }
2013 } // end if Ucat
2014 } // end debug logging check
2015
2016 LOG_ALLOW(GLOBAL, LOG_INFO, "Rank %d: Completed ghost update for field '%s'.\n", rank, field_name);
2017 PROFILE_FUNCTION_END;
2018 PetscFunctionReturn(0);
2019}
2020
2021#undef __FUNCT__
2022#define __FUNCT__ "SetupBoundaryConditions"
2023/**
2024 * @brief Internal helper implementation: `SetupBoundaryConditions()`.
2025 * @details Local to this translation unit.
2026 */
2027PetscErrorCode SetupBoundaryConditions(SimCtx *simCtx)
2028{
2029 PetscErrorCode ierr;
2030 PetscFunctionBeginUser;
2031
2032 PROFILE_FUNCTION_BEGIN;
2033
2034 LOG_ALLOW(GLOBAL,LOG_INFO, "--- Setting up Boundary Conditions ---\n");
2035 // --- Phase 1: Parse and initialize BC configuration for all blocks ---
2036 LOG_ALLOW(GLOBAL,LOG_INFO,"Parsing BC configuration files and initializing boundary condition data structures.\n");
2037 UserCtx *user_finest = simCtx->usermg.mgctx[simCtx->usermg.mglevels-1].user;
2038 for (PetscInt bi = 0; bi < simCtx->block_number; bi++) {
2039 LOG_ALLOW(GLOBAL,LOG_DEBUG, " -> Processing Block %d:\n", bi);
2040
2041 // --- Generate the filename for the current block ---
2042 const char *current_bc_filename = simCtx->bcs_files[bi];
2043 LOG_ALLOW(GLOBAL,LOG_DEBUG," -> Processing Block %d using config file '%s'\n", bi, current_bc_filename);
2044 // This will populate user_finest[bi].boundary_faces
2045
2046 //ierr = ParseAllBoundaryConditions(&user_finest[bi],current_bc_filename); CHKERRQ(ierr);
2047
2048 ierr = BoundarySystem_Initialize(&user_finest[bi], current_bc_filename); CHKERRQ(ierr);
2049 }
2050
2051 // Propogate BC Configuration to coarser levels.
2052 ierr = PropagateBoundaryConfigToCoarserLevels(simCtx); CHKERRQ(ierr);
2053
2054 // Validate the geometric contract before any metric consumes periodic geometry.
2055 for (PetscInt level = simCtx->usermg.mglevels - 1; level >= 0; level--) {
2056 UserCtx *level_users = simCtx->usermg.mgctx[level].user;
2057 for (PetscInt bi = 0; bi < simCtx->block_number; bi++) {
2058 ierr = ValidatePeriodicGeometry(&level_users[bi]); CHKERRQ(ierr);
2059 }
2060 }
2061
2062 // --- Calculate Grid Metrics (requires BC configuration) ---
2063 // NOTE: This MUST be called here (after BC initialization but before inlet/outlet calculations) because:
2064 // 1. Periodic BC corrections in metric calculations need boundary_faces data to be populated
2065 // 2. Inlet/Outlet area calculations (below) require computed metrics (Csi, Eta, Zet) to be available
2066 // Previously this was in SetupGridAndSolvers, but that caused metrics to be computed without BC info.
2067 LOG_ALLOW(GLOBAL,LOG_INFO,"Computing grid metrics with boundary condition information.\n");
2068 ierr = CalculateAllGridMetrics(simCtx); CHKERRQ(ierr);
2069
2070 // --- Phase 2: Calculate inlet/outlet properties (requires computed metrics) ---
2071 LOG_ALLOW(GLOBAL,LOG_INFO,"Calculating inlet and outlet face properties.\n");
2072 for (PetscInt bi = 0; bi < simCtx->block_number; bi++) {
2073 // Call the function to calculate the center of the inlet face & the inlet area, which may be used to calculate Boundary values.
2074 ierr = CalculateInletProperties(&user_finest[bi]); CHKERRQ(ierr);
2075
2076 // Call the function to calculate the center of the outlet face & the outlet area, which may be used to calculate Boundary values.
2077 ierr = CalculateOutletProperties(&user_finest[bi]); CHKERRQ(ierr);
2078 }
2079
2080 LOG_ALLOW(GLOBAL,LOG_INFO, "--- Boundary Conditions setup complete ---\n");
2081
2082
2083 PROFILE_FUNCTION_END;
2084 PetscFunctionReturn(0);
2085}
2086
2087/**
2088 * @brief Internal helper implementation: `Allocate3DArrayScalar()`.
2089 * @details Local to this translation unit.
2090 */
2091PetscErrorCode Allocate3DArrayScalar(PetscReal ****array, PetscInt nz, PetscInt ny, PetscInt nx)
2092{
2093 PetscErrorCode ierr;
2094 PetscReal ***data;
2095 PetscReal *dataContiguous;
2096 PetscInt k, j;
2097
2098 PetscFunctionBegin;
2099 /* Step 1: Allocate memory for an array of nz layer pointers (zero-initialized) */
2100 ierr = PetscCalloc1(nz, &data); CHKERRQ(ierr);
2101
2102 /* Step 2: Allocate memory for all row pointers (nz * ny pointers) */
2103 ierr = PetscCalloc1(nz * ny, &data[0]); CHKERRQ(ierr);
2104 for (k = 1; k < nz; k++) {
2105 data[k] = data[0] + k * ny;
2106 }
2107
2108 /* Step 3: Allocate one contiguous block for all data elements (nz*ny*nx) */
2109 ierr = PetscCalloc1(nz * ny * nx, &dataContiguous); CHKERRQ(ierr);
2110
2111 /* Build the 3D pointer structure: each row pointer gets the correct segment of data */
2112 for (k = 0; k < nz; k++) {
2113 for (j = 0; j < ny; j++) {
2114 data[k][j] = dataContiguous + (k * ny + j) * nx;
2115 /* Memory is already zeroed by PetscCalloc1, so no manual initialization is needed */
2116 }
2117 }
2118 *array = data;
2119 PetscFunctionReturn(0);
2120}
2121
2122/**
2123 * @brief Internal helper implementation: `Deallocate3DArrayScalar()`.
2124 * @details Local to this translation unit.
2125 */
2126PetscErrorCode Deallocate3DArrayScalar(PetscReal ***array, PetscInt nz, PetscInt ny)
2127{
2128 PetscErrorCode ierr;
2129 (void)nz;
2130 (void)ny;
2131
2132 PetscFunctionBegin;
2133 if (!array || !array[0] || !array[0][0] ) { // Added more robust check
2134 LOG_ALLOW(GLOBAL, LOG_WARNING, "Deallocate3DArrayScalar called with potentially unallocated or NULL array.\n");
2135 if (array) {
2136 if (array[0]) { // Check if row pointers might exist
2137 // Cannot safely access array[0][0] if array[0] might be invalid/freed
2138 // Standard deallocation below assumes valid pointers.
2139 ierr = PetscFree(array[0]); CHKERRQ(ierr); // Free row pointers if they exist
2140 }
2141 ierr = PetscFree(array); CHKERRQ(ierr); // Free layer pointers if they exist
2142 }
2143 PetscFunctionReturn(0);
2144 }
2145
2146 // --- Standard Deallocation (assuming valid allocation) ---
2147
2148 /* 1. Free the contiguous block of PetscReal values.
2149 The starting address was stored in array[0][0]. */
2150 ierr = PetscFree(array[0][0]); CHKERRQ(ierr); // Free the ACTUAL DATA
2151
2152 /* 2. Free the contiguous block of row pointers.
2153 The starting address was stored in array[0]. */
2154 ierr = PetscFree(array[0]); CHKERRQ(ierr); // Free the ROW POINTERS
2155
2156 /* 3. Free the layer pointer array.
2157 The starting address is 'array' itself. */
2158 ierr = PetscFree(array); CHKERRQ(ierr); // Free the LAYER POINTERS
2159
2160 PetscFunctionReturn(0);
2161}
2162
2163/**
2164 * @brief Implementation of \ref Allocate3DArrayVector().
2165 * @details Full API contract (arguments, ownership, side effects) is documented with
2166 * the header declaration in `include/setup.h`.
2167 * @see Allocate3DArrayVector()
2168 */
2169PetscErrorCode Allocate3DArrayVector(Cmpnts ****array, PetscInt nz, PetscInt ny, PetscInt nx)
2170{
2171 PetscErrorCode ierr;
2172 Cmpnts ***data;
2173 Cmpnts *dataContiguous;
2174 PetscInt k, j;
2175 PetscMPIInt rank;
2176
2177 PetscFunctionBegin;
2178
2179 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
2180
2181 /* Step 1: Allocate memory for nz layer pointers (zeroed) */
2182 ierr = PetscCalloc1(nz, &data); CHKERRQ(ierr);
2183
2184 LOG_ALLOW(LOCAL,LOG_DEBUG," [Rank %d] memory allocated for outermost layer (%d k-layer pointers).\n",rank,nz);
2185
2186 /* Step 2: Allocate memory for all row pointers (nz * ny pointers) */
2187 ierr = PetscCalloc1(nz * ny, &data[0]); CHKERRQ(ierr);
2188 for (k = 1; k < nz; k++) {
2189 data[k] = data[0] + k * ny;
2190 }
2191
2192 LOG_ALLOW(LOCAL,LOG_DEBUG,"[Rank %d] memory allocated for %dx%d row pointers.\n",rank,nz,ny);
2193
2194 /* Step 3: Allocate one contiguous block for nz*ny*nx Cmpnts structures (zeroed) */
2195 ierr = PetscCalloc1(nz * ny * nx, &dataContiguous); CHKERRQ(ierr);
2196
2197 LOG_ALLOW(GLOBAL,LOG_DEBUG,"[Rank %d] memory allocated for contigous block of %dx%dx%d Cmpnts structures).\n",rank,nz,ny,nx);
2198
2199 /* Build the 3D pointer structure for vector data */
2200 for (k = 0; k < nz; k++) {
2201 for (j = 0; j < ny; j++) {
2202 data[k][j] = dataContiguous + (k * ny + j) * nx;
2203 /* The PetscCalloc1 call has already initialized each Cmpnts to zero. */
2204 }
2205 }
2206
2207 LOG_ALLOW(GLOBAL,LOG_DEBUG,"[Rank %d] 3D pointer structure for vector data created. \n",rank);
2208
2209 *array = data;
2210 PetscFunctionReturn(0);
2211}
2212
2213/**
2214 * @brief Implementation of \ref Deallocate3DArrayVector().
2215 * @details Full API contract (arguments, ownership, side effects) is documented with
2216 * the header declaration in `include/setup.h`.
2217 * @see Deallocate3DArrayVector()
2218 */
2219 PetscErrorCode Deallocate3DArrayVector(Cmpnts ***array, PetscInt nz, PetscInt ny)
2220{
2221 PetscErrorCode ierr;
2222 (void)nz;
2223 (void)ny;
2224
2225 PetscFunctionBegin;
2226 // If array is NULL or hasn't been allocated properly, just return.
2227 if (!array || !array[0] || !array[0][0] ) {
2228 LOG_ALLOW(GLOBAL, LOG_WARNING, "Deallocate3DArrayVector called with potentially unallocated or NULL array.\n");
2229 // Attempt to free what might exist, but be cautious
2230 if (array) {
2231 if (array[0]) { // Check if row pointers were allocated
2232 // We don't have a direct pointer to the contiguous data block
2233 // saved separately in this allocation scheme. The allocation relies
2234 // on array[0][0] pointing to it. If array[0] was freed first,
2235 // accessing array[0][0] is unsafe.
2236 // The allocation scheme where the contiguous data block is not
2237 // stored separately makes safe deallocation tricky if freeing
2238 // happens out of order or if parts are NULL.
2239
2240 // A SAFER ALLOCATION/DEALLOCATION would store the data pointer separately.
2241 // Given the current allocation scheme, the order MUST be:
2242 // 1. Free the data block (pointed to by array[0][0])
2243 // 2. Free the row pointer block (pointed to by array[0])
2244 // 3. Free the layer pointer block (pointed to by array)
2245
2246 // Let's assume the allocation was successful and pointers are valid.
2247 // Get pointer to the contiguous data block *before* freeing row pointers
2248 Cmpnts *dataContiguous = array[0][0];
2249 ierr = PetscFree(dataContiguous); CHKERRQ(ierr); // Free data block
2250
2251 // Now free the row pointers block
2252 ierr = PetscFree(array[0]); CHKERRQ(ierr); // Free row pointers
2253
2254 }
2255 // Finally, free the array of layer pointers
2256 ierr = PetscFree(array); CHKERRQ(ierr);
2257 }
2258 PetscFunctionReturn(0); // Return gracefully if input was NULL initially
2259 }
2260
2261
2262 // --- Standard Deallocation (assuming valid allocation) ---
2263
2264 /* 1. Free the contiguous block of Cmpnts structures.
2265 The starting address was stored in array[0][0] by Allocate3DArrayVector. */
2266 ierr = PetscFree(array[0][0]); CHKERRQ(ierr); // Free the ACTUAL DATA
2267
2268 /* 2. Free the contiguous block of row pointers.
2269 The starting address was stored in array[0]. */
2270 ierr = PetscFree(array[0]); CHKERRQ(ierr); // Free the ROW POINTERS
2271
2272 /* 3. Free the layer pointer array.
2273 The starting address is 'array' itself. */
2274 ierr = PetscFree(array); CHKERRQ(ierr); // Free the LAYER POINTERS
2275
2276 PetscFunctionReturn(0);
2277}
2278
2279#undef __FUNCT__
2280#define __FUNCT__ "GetOwnedCellRange"
2281/**
2282 * @brief Internal helper implementation: `GetOwnedCellRange()`.
2283 * @details Local to this translation unit.
2284 */
2285PetscErrorCode GetOwnedCellRange(const DMDALocalInfo *info_nodes,
2286 PetscInt dim,
2287 PetscInt *xs_cell_global_out,
2288 PetscInt *xm_cell_local_out)
2289{
2290 PetscErrorCode ierr = 0; // Standard PETSc error code, not explicitly set here but good practice.
2291 PetscInt xs_node_global_rank; // Global index of the first node owned by this rank in the specified dimension.
2292 PetscInt num_nodes_owned_rank; // Number of nodes owned by this rank in this dimension (local count, excluding ghosts).
2293 PetscInt GlobalNodesInDim_from_info; // Total number of DA points in this dimension, from DMDALocalInfo.
2294
2295 PetscFunctionBeginUser;
2296
2297 // --- 1. Input Validation ---
2298 if (!info_nodes || !xs_cell_global_out || !xm_cell_local_out) {
2299 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Null pointer passed to GetOwnedCellRange.");
2300 }
2301
2302 // --- 2. Extract Node Ownership and Global Dimension Information from DMDALocalInfo ---
2303 if (dim == 0) { // I-direction
2304 xs_node_global_rank = info_nodes->xs;
2305 num_nodes_owned_rank = info_nodes->xm;
2306 GlobalNodesInDim_from_info = info_nodes->mx;
2307 } else if (dim == 1) { // J-direction
2308 xs_node_global_rank = info_nodes->ys;
2309 num_nodes_owned_rank = info_nodes->ym;
2310 GlobalNodesInDim_from_info = info_nodes->my;
2311 } else if (dim == 2) { // K-direction
2312 xs_node_global_rank = info_nodes->zs;
2313 num_nodes_owned_rank = info_nodes->zm;
2314 GlobalNodesInDim_from_info = info_nodes->mz;
2315 } else {
2316 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d in GetOwnedCellRange. Must be 0, 1, or 2.", dim);
2317 }
2318
2319 // --- 3. Correct for User-Defined Ghost Node ---
2320 // Per the function's contract (@warning), the DA size includes an extra, non-physical
2321 // node. We subtract 1 to get the true number of physical nodes for cell calculations.
2322 const PetscInt physical_nodes_in_dim = GlobalNodesInDim_from_info - 1;
2323
2324 // --- 4. Handle Edge Cases for Physical Domain Size ---
2325 // If the physical domain has 0 or 1 node, no cells can be formed.
2326 if (physical_nodes_in_dim <= 1) {
2327 *xs_cell_global_out = xs_node_global_rank; // Still report the rank's starting node
2328 *xm_cell_local_out = 0; // But 0 cells
2329 PetscFunctionReturn(0);
2330 }
2331
2332 // --- 5. Determine Cell Ownership Based on Corrected Node Ownership ---
2333 // The first cell this rank *could* define has its origin at the first node this rank owns.
2334 *xs_cell_global_out = xs_node_global_rank;
2335
2336 // If the rank owns no nodes in this dimension, it can't form any cell origins.
2337 if (num_nodes_owned_rank == 0) {
2338 *xm_cell_local_out = 0;
2339 } else {
2340 // --- BUG FIX APPLIED HERE ---
2341 // The previous logic incorrectly assumed a cell's end node (N_{k+1}) must be on the
2342 // same rank as its origin node (N_k). The correct logic is to find the intersection
2343 // between the nodes this rank owns and the nodes that are valid origins globally.
2344
2345 // The first node owned by the rank is its first potential origin.
2346 PetscInt first_owned_origin = xs_node_global_rank;
2347
2348 // The absolute last node owned by this rank. Any node up to and including this one
2349 // is a potential cell origin from this rank's perspective.
2350 PetscInt last_node_owned_by_rank = xs_node_global_rank + num_nodes_owned_rank - 1;
2351
2352 // The absolute last node in the entire PHYSICAL domain that can serve as a cell origin.
2353 // If there are `N` physical nodes (0 to N-1), this index is `N-2`.
2354 PetscInt last_possible_origin_global_idx = physical_nodes_in_dim - 2;
2355
2356 // The actual last origin this rank can provide is the *minimum* of what it owns
2357 // and what is globally possible. This correctly handles both ranks in the middle of
2358 // the domain and the very last rank.
2359 PetscInt actual_last_origin_this_rank_can_form = PetscMin(last_node_owned_by_rank, last_possible_origin_global_idx);
2360
2361 // If the first potential origin this rank owns is already beyond the actual last
2362 // origin it can form, then this rank forms no valid cell origins. This happens if
2363 // the rank only owns the very last physical node.
2364 if (first_owned_origin > actual_last_origin_this_rank_can_form) {
2365 *xm_cell_local_out = 0;
2366 } else {
2367 // The number of cells is the count of valid origins this rank owns.
2368 // (Count = Last Index - First Index + 1)
2369 *xm_cell_local_out = actual_last_origin_this_rank_can_form - first_owned_origin + 1;
2370 }
2371 }
2372
2373 PetscFunctionReturn(ierr);
2374}
2375
2376#undef __FUNCT__
2377#define __FUNCT__ "ComputeAndStoreNeighborRanks"
2378/**
2379 * @brief Internal helper implementation: `ComputeAndStoreNeighborRanks()`.
2380 * @details Local to this translation unit.
2381 */
2382PetscErrorCode ComputeAndStoreNeighborRanks(UserCtx *user)
2383{
2384 PetscErrorCode ierr;
2385 PetscMPIInt rank;
2386 PetscMPIInt size; // MPI communicator size
2387 const PetscMPIInt *neighbor_ranks_ptr; // Pointer to raw neighbor data from PETSc
2388
2389 PetscFunctionBeginUser;
2390 PROFILE_FUNCTION_BEGIN;
2391 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
2392 ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size); CHKERRQ(ierr); // Get MPI size for validation
2393
2394 LOG_ALLOW(GLOBAL, LOG_INFO, "Rank %d: Computing DMDA neighbor ranks.\n", rank);
2395
2396 if (!user || !user->da) {
2397 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "UserCtx or user->da is NULL in ComputeAndStoreNeighborRanks.");
2398 }
2399
2400 // Get the neighbor information from the DMDA
2401 // neighbor_ranks_ptr will point to an internal PETSc array of 27 ranks.
2402 ierr = DMDAGetNeighbors(user->da, &neighbor_ranks_ptr); CHKERRQ(ierr);
2403
2404 // Log the raw values from DMDAGetNeighbors for boundary-relevant directions for debugging
2405 LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG, "[Rank %d]Raw DMDAGetNeighbors: xm_raw=%d, xp_raw=%d, ym_raw=%d, yp_raw=%d, zm_raw=%d, zp_raw=%d. MPI_PROC_NULL is %d.\n",
2406 rank,
2407 neighbor_ranks_ptr[12], neighbor_ranks_ptr[14],
2408 neighbor_ranks_ptr[10], neighbor_ranks_ptr[16],
2409 neighbor_ranks_ptr[4], neighbor_ranks_ptr[22],
2410 (int)MPI_PROC_NULL);
2411
2412 // PETSc standard indices for 3D face neighbors from the 27-point stencil:
2413 // Index = k_offset*9 + j_offset*3 + i_offset (where offsets -1,0,1 map to 0,1,2)
2414 // Center: (i_off=1, j_off=1, k_off=1) => 1*9 + 1*3 + 1 = 13
2415 // X-min: (i_off=0, j_off=1, k_off=1) => 1*9 + 1*3 + 0 = 12
2416 // X-plus: (i_off=2, j_off=1, k_off=1) => 1*9 + 1*3 + 2 = 14
2417 // Y-min: (i_off=1, j_off=0, k_off=1) => 1*9 + 0*3 + 1 = 10
2418 // Y-plus: (i_off=1, j_off=2, k_off=1) => 1*9 + 2*3 + 1 = 16
2419 // Z-min: (i_off=1, j_off=1, k_off=0) => 0*9 + 1*3 + 1 = 4
2420 // Z-plus: (i_off=1, j_off=1, k_off=2) => 2*9 + 1*3 + 1 = 22
2421
2422 if (neighbor_ranks_ptr[13] != rank) {
2423 LOG_ALLOW(GLOBAL, LOG_WARNING, "Rank %d: DMDAGetNeighbors center index (13) is %d, expected current rank %d. Neighbor indexing might be non-standard or DMDA small.\n",
2424 rank, neighbor_ranks_ptr[13], rank);
2425 // This warning is important. If the center isn't the current rank, the offsets are likely wrong.
2426 // However, PETSc should ensure this unless the DM is too small for a 3x3x3 stencil.
2427 }
2428
2429 // Assign and sanitize each neighbor rank
2430 PetscMPIInt temp_neighbor;
2431
2432 temp_neighbor = neighbor_ranks_ptr[12]; // xm
2433 if (temp_neighbor < 0 || temp_neighbor >= size) {
2434 LOG_ALLOW(GLOBAL, LOG_WARNING, "[Rank %d] Correcting invalid xm neighbor %d to MPI_PROC_NULL (%d).\n", rank, temp_neighbor, (int)MPI_PROC_NULL);
2435 user->neighbors.rank_xm = MPI_PROC_NULL;
2436 } else {
2437 user->neighbors.rank_xm = temp_neighbor;
2438 }
2439
2440 temp_neighbor = neighbor_ranks_ptr[14]; // xp
2441 if (temp_neighbor < 0 || temp_neighbor >= size) {
2442 LOG_ALLOW(GLOBAL, LOG_WARNING, "[Rank %d] Correcting invalid xp neighbor %d to MPI_PROC_NULL (%d).\n", rank, temp_neighbor, (int)MPI_PROC_NULL);
2443 user->neighbors.rank_xp = MPI_PROC_NULL;
2444 } else {
2445 user->neighbors.rank_xp = temp_neighbor;
2446 }
2447
2448 temp_neighbor = neighbor_ranks_ptr[10]; // ym
2449 if (temp_neighbor < 0 || temp_neighbor >= size) {
2450 LOG_ALLOW(GLOBAL, LOG_WARNING, "[Rank %d] Correcting invalid ym neighbor %d to MPI_PROC_NULL (%d).\n", rank, temp_neighbor, (int)MPI_PROC_NULL);
2451 user->neighbors.rank_ym = MPI_PROC_NULL;
2452 } else {
2453 user->neighbors.rank_ym = temp_neighbor;
2454 }
2455
2456 temp_neighbor = neighbor_ranks_ptr[16]; // yp
2457 if (temp_neighbor < 0 || temp_neighbor >= size) {
2458 // The log for index 16 was "zm" in your output, should be yp
2459 LOG_ALLOW(GLOBAL, LOG_WARNING, "[Rank %d] Correcting invalid yp neighbor (raw index 16) %d to MPI_PROC_NULL (%d).\n", rank, temp_neighbor, (int)MPI_PROC_NULL);
2460 user->neighbors.rank_yp = MPI_PROC_NULL;
2461 } else {
2462 user->neighbors.rank_yp = temp_neighbor;
2463 }
2464
2465 temp_neighbor = neighbor_ranks_ptr[4]; // zm
2466 if (temp_neighbor < 0 || temp_neighbor >= size) {
2467 LOG_ALLOW(GLOBAL, LOG_WARNING, "[Rank %d] Correcting invalid zm neighbor %d to MPI_PROC_NULL (%d).\n", rank, temp_neighbor, (int)MPI_PROC_NULL);
2468 user->neighbors.rank_zm = MPI_PROC_NULL;
2469 } else {
2470 user->neighbors.rank_zm = temp_neighbor;
2471 }
2472
2473 temp_neighbor = neighbor_ranks_ptr[22]; // zp
2474 if (temp_neighbor < 0 || temp_neighbor >= size) {
2475 LOG_ALLOW(GLOBAL, LOG_WARNING, "[Rank %d] Correcting invalid zp neighbor %d to MPI_PROC_NULL (%d).\n", rank, temp_neighbor, (int)MPI_PROC_NULL);
2476 user->neighbors.rank_zp = MPI_PROC_NULL;
2477 } else {
2478 user->neighbors.rank_zp = temp_neighbor;
2479 }
2480
2481 LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG, "[Rank %d] Stored user->neighbors: xm=%d, xp=%d, ym=%d, yp=%d, zm=%d, zp=%d\n", rank,
2482 user->neighbors.rank_xm, user->neighbors.rank_xp,
2483 user->neighbors.rank_ym, user->neighbors.rank_yp,
2484 user->neighbors.rank_zm, user->neighbors.rank_zp);
2485 PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT); // Ensure logs are flushed
2486
2487 // Note: neighbor_ranks_ptr memory is managed by PETSc, do not free it.
2488 PROFILE_FUNCTION_END;
2489 PetscFunctionReturn(0);
2490}
2491
2492#undef __FUNCT__
2493#define __FUNCT__ "SetDMDAProcLayout"
2494/**
2495 * @brief Internal helper implementation: `SetDMDAProcLayout()`.
2496 * @details Local to this translation unit.
2497 */
2498PetscErrorCode SetDMDAProcLayout(DM dm, UserCtx *user)
2499{
2500 PetscErrorCode ierr;
2501 PetscMPIInt size, rank;
2502 PetscInt px = PETSC_DECIDE, py = PETSC_DECIDE, pz = PETSC_DECIDE;
2503 PetscBool px_set = PETSC_FALSE, py_set = PETSC_FALSE, pz_set = PETSC_FALSE;
2504 SimCtx *simCtx = user->simCtx;
2505
2506 // Set no.of processors in direction 1
2507 if(simCtx->da_procs_x) {
2508 px_set = PETSC_TRUE;
2509 px = simCtx->da_procs_x;
2510 }
2511 // Set no.of processors in direction 2
2512 if(simCtx->da_procs_y) {
2513 py_set = PETSC_TRUE;
2514 py = simCtx->da_procs_y;
2515 }
2516 // Set no.of processors in direction 1
2517 if(simCtx->da_procs_z) {
2518 pz_set = PETSC_TRUE;
2519 pz = simCtx->da_procs_z;
2520 }
2521
2522 PetscFunctionBeginUser;
2523 PROFILE_FUNCTION_BEGIN;
2524 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size); CHKERRQ(ierr);
2525 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank); CHKERRQ(ierr);
2526 LOG_ALLOW(GLOBAL, LOG_INFO, "Rank %d: Configuring DMDA processor layout for %d total processes.\n", rank, size);
2527
2528 // --- Validate User Input (Optional but Recommended) ---
2529 // Check if specified processor counts multiply to the total MPI size
2530 if (px_set && py_set && pz_set) {
2531 if (px * py * pz != size) {
2532 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_INCOMP,
2533 "Specified processor layout %d x %d x %d = %d does not match MPI size %d",
2534 px, py, pz, px * py * pz, size);
2535 }
2536 LOG_ALLOW(GLOBAL, LOG_INFO, "Using specified processor layout: %d x %d x %d\n", px, py, pz);
2537 } else if (px_set || py_set || pz_set) {
2538 // If only some are set, PETSC_DECIDE will be used for others
2539 LOG_ALLOW(GLOBAL, LOG_INFO, "Using partially specified processor layout: %d x %d x %d (PETSC_DECIDE for unspecified)\n", px, py, pz);
2540 } else {
2541 LOG_ALLOW(GLOBAL, LOG_INFO, "Using fully automatic processor layout (PETSC_DECIDE x PETSC_DECIDE x PETSC_DECIDE)\n");
2542 }
2543 // Additional checks: Ensure px, py, pz are positive if set
2544 if ((px_set && px <= 0) || (py_set && py <= 0) || (pz_set && pz <= 0)) {
2545 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Specified processor counts must be positive.");
2546 }
2547
2548
2549 // --- Apply the layout to the DMDA ---
2550 ierr = DMDASetNumProcs(dm, px, py, pz); CHKERRQ(ierr);
2551 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Rank %d: DMDASetNumProcs called with px=%d, py=%d, pz=%d.\n", rank, px, py, pz);
2552
2553 // --- Store the values in UserCtx (Optional) ---
2554 // Note: If PETSC_DECIDE was used, PETSc calculates the actual values during DMSetUp.
2555 // We store the *requested* values here. To get the *actual* values used,
2556 // you would need to call DMDAGetInfo after DMSetUp.
2557 /*
2558 if (user) {
2559 user->procs_x = px;
2560 user->procs_y = py;
2561 user->procs_z = pz;
2562 }
2563 */
2564 PROFILE_FUNCTION_END;
2565 PetscFunctionReturn(0);
2566}
2567
2568#undef __FUNCT__
2569#define __FUNCT__ "SetupDomainRankInfo"
2570/**
2571 * @brief Implementation of \ref SetupDomainRankInfo().
2572 * @details Full API contract (arguments, ownership, side effects) is documented with
2573 * the header declaration in `include/setup.h`.
2574 * @see SetupDomainRankInfo()
2575 */
2576PetscErrorCode SetupDomainRankInfo(SimCtx *simCtx)
2577{
2578 PetscErrorCode ierr;
2579 PetscInt nblk = simCtx->block_number;
2580 PetscInt size = simCtx->size;
2581 BoundingBox *final_bboxlist = NULL;
2582
2583 PetscFunctionBeginUser;
2584 PROFILE_FUNCTION_BEGIN;
2585
2586 LOG_ALLOW(GLOBAL, LOG_INFO, "Starting full rank communication setup for %d block(s).\n", nblk);
2587
2588 UserCtx *user_finest = simCtx->usermg.mgctx[simCtx->usermg.mglevels - 1].user;
2589
2590 // --- Step 1: Compute neighbor ranks (unchanged) ---
2591 for (int bi = 0; bi < nblk; bi++) {
2592 ierr = ComputeAndStoreNeighborRanks(&user_finest[bi]); CHKERRQ(ierr);
2593 }
2594 LOG_ALLOW(GLOBAL, LOG_INFO, "Neighbor ranks computed and stored for all blocks.\n");
2595
2596 // --- Step 2: Allocate the final, unified list on ALL ranks ---
2597 // Every rank will build this list in parallel.
2598 ierr = PetscMalloc1(size * nblk, &final_bboxlist); CHKERRQ(ierr);
2599
2600 // --- Step 3: Loop through each block, gather then broadcast its bbox list ---
2601 for (int bi = 0; bi < nblk; bi++) {
2602 // This is a temporary pointer for the current block's list.
2603 BoundingBox *block_bboxlist = NULL;
2604
2605 LOG_ALLOW(GLOBAL, LOG_INFO, "Processing bounding boxes for block %d...\n", bi);
2606
2607 // A) GATHER: On rank 0, block_bboxlist is allocated and filled. On others, it's NULL.
2608 ierr = GatherAllBoundingBoxes(&user_finest[bi], &block_bboxlist); CHKERRQ(ierr);
2609 LOG_ALLOW(GLOBAL, LOG_DEBUG, " -> Gather complete for block %d.\n", bi);
2610
2611 // B) BROADCAST: On non-root ranks, block_bboxlist is allocated. Then, the data
2612 // from rank 0 is broadcast to all ranks. After this call, ALL ranks have
2613 // an identical, complete copy of the bounding boxes for the current block.
2614 ierr = BroadcastAllBoundingBoxes(&user_finest[bi], &block_bboxlist); CHKERRQ(ierr);
2615 LOG_ALLOW(GLOBAL, LOG_DEBUG, " -> Broadcast complete for block %d.\n", bi);
2616
2617 // C) ASSEMBLE: Every rank now copies the data for this block into the
2618 // correct segment of its final, unified list.
2619 for (int r = 0; r < size; r++) {
2620 // The layout is [r0b0, r1b0, ..., r(size-1)b0, r0b1, r1b1, ...]
2621 final_bboxlist[bi * size + r] = block_bboxlist[r];
2622 }
2623 LOG_ALLOW(GLOBAL, LOG_DEBUG, " -> Assembly into final list complete for block %d.\n", bi);
2624
2625 // D) CLEANUP: Free the temporary list for this block on ALL ranks before the next iteration.
2626 // Your helper functions use malloc, so we must use free.
2627 free(block_bboxlist);
2628 }
2629
2630 // --- Step 4: Assign the final pointer and run the last setup step ---
2631 simCtx->bboxlist = final_bboxlist;
2632 LOG_ALLOW(GLOBAL, LOG_INFO, "Final unified bboxlist created on all ranks and stored in SimCtx.\n");
2633
2634 ierr = SetupDomainCellDecompositionMap(&user_finest[0]); CHKERRQ(ierr);
2635 LOG_ALLOW(GLOBAL, LOG_INFO, "Domain Cell Composition set and broadcasted.\n");
2636
2637 LOG_ALLOW(GLOBAL, LOG_INFO, "SetupDomainRankInfo: Completed successfully.\n");
2638
2639 PROFILE_FUNCTION_END;
2640 PetscFunctionReturn(0);
2641}
2642
2643#undef __FUNCT__
2644#define __FUNCT__ "Contra2Cart"
2645/**
2646 * @brief Internal helper implementation: `Contra2Cart()`.
2647 * @details Local to this translation unit.
2648 */
2649PetscErrorCode Contra2Cart(UserCtx *user)
2650{
2651 PetscErrorCode ierr;
2652 DMDALocalInfo info;
2653 Cmpnts ***lcsi_arr, ***leta_arr, ***lzet_arr; // Local metric arrays
2654 Cmpnts ***lucont_arr; // Local contravariant velocity array
2655 Cmpnts ***gucat_arr; // Global Cartesian velocity array
2656 PetscReal ***lnvert_arr; // Local Nvert array
2657 PetscReal ***laj_arr; // Local Jacobian Determinant inverse array
2658
2659 PetscFunctionBeginUser;
2660 PROFILE_FUNCTION_BEGIN;
2661 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Starting Contravariant-to-Cartesian velocity transformation.\n");
2662
2663 // --- 1. Get DMDA Info and Check for Valid Inputs ---
2664 // All inputs (lUcont, lCsi, etc.) and outputs (Ucat) are on DMs from the UserCtx.
2665 // We get local info from fda, which governs the layout of most arrays here.
2666 ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
2667 if (!user->lUcont || !user->lCsi || !user->lEta || !user->lZet || !user->lNvert || !user->Ucat) {
2668 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Contra2Cart requires lUcont, lCsi/Eta/Zet, lNvert, and Ucat to be non-NULL.");
2669 }
2670
2671
2672 // --- 2. Get Read-Only Array Access to Local Input Vectors (with ghosts) ---
2673 ierr = DMDAVecGetArrayRead(user->fda, user->lUcont, &lucont_arr); CHKERRQ(ierr);
2674 ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, &lcsi_arr); CHKERRQ(ierr);
2675 ierr = DMDAVecGetArrayRead(user->fda, user->lEta, &leta_arr); CHKERRQ(ierr);
2676 ierr = DMDAVecGetArrayRead(user->fda, user->lZet, &lzet_arr); CHKERRQ(ierr);
2677 ierr = DMDAVecGetArrayRead(user->da, user->lNvert, &lnvert_arr); CHKERRQ(ierr);
2678 ierr = DMDAVecGetArrayRead(user->da, user->lAj, &laj_arr); CHKERRQ(ierr);
2679
2680 // --- 3. Get Write-Only Array Access to the Global Output Vector ---
2681 // We compute for local owned cells and write into the global vector.
2682 // PETSc handles mapping the global indices to the correct local memory locations.
2683 ierr = DMDAVecGetArray(user->fda, user->Ucat, &gucat_arr); CHKERRQ(ierr);
2684
2685
2686 // --- 4. Define Loop Bounds for INTERIOR Cells ---
2687 // We use adjusted bounds to avoid calculating Ucat on the physical domain boundaries,
2688 // as these are typically set explicitly by boundary condition functions.
2689 // The stencils use indices like i-1, j-1, k-1, so we must start loops at least at index 1.
2690 PetscInt i_start = (info.xs == 0) ? info.xs + 1 : info.xs;
2691 PetscInt i_end = (info.xs + info.xm == info.mx) ? info.xs + info.xm - 1 : info.xs + info.xm;
2692
2693 PetscInt j_start = (info.ys == 0) ? info.ys + 1 : info.ys;
2694 PetscInt j_end = (info.ys + info.ym == info.my) ? info.ys + info.ym - 1 : info.ys + info.ym;
2695
2696 PetscInt k_start = (info.zs == 0) ? info.zs + 1 : info.zs;
2697 PetscInt k_end = (info.zs + info.zm == info.mz) ? info.zs + info.zm - 1 : info.zs + info.zm;
2698
2699 // --- 5. Main Computation Loop ---
2700 // Loops over the GLOBAL indices of interior cells owned by this rank.
2701 for (PetscInt k_cell = k_start; k_cell < k_end; ++k_cell) {
2702 for (PetscInt j_cell = j_start; j_cell < j_end; ++j_cell) {
2703 for (PetscInt i_cell = i_start; i_cell < i_end; ++i_cell) {
2704
2705 // Check if the cell is a fluid cell (not solid/blanked)
2706 // if (lnvert_arr[k_cell][j_cell][i_cell] > 0.1) continue; // Skip solid/blanked cells
2707
2708 // Transformation matrix [mat] is the metric tensor at the cell center,
2709 // estimated by averaging metrics from adjacent faces.
2710 PetscReal mat[3][3];
2711
2712 // PetscReal aj_center = laj_arr[k_cell+1][j_cell+1][i_cell+1];
2713
2714 mat[0][0] = 0.5 * (lcsi_arr[k_cell][j_cell][i_cell-1].x + lcsi_arr[k_cell][j_cell][i_cell].x); //* aj_center;
2715 mat[0][1] = 0.5 * (lcsi_arr[k_cell][j_cell][i_cell-1].y + lcsi_arr[k_cell][j_cell][i_cell].y); //* aj_center;
2716 mat[0][2] = 0.5 * (lcsi_arr[k_cell][j_cell][i_cell-1].z + lcsi_arr[k_cell][j_cell][i_cell].z); //* aj_center;
2717
2718 mat[1][0] = 0.5 * (leta_arr[k_cell][j_cell-1][i_cell].x + leta_arr[k_cell][j_cell][i_cell].x); //* aj_center;
2719 mat[1][1] = 0.5 * (leta_arr[k_cell][j_cell-1][i_cell].y + leta_arr[k_cell][j_cell][i_cell].y); //* aj_center;
2720 mat[1][2] = 0.5 * (leta_arr[k_cell][j_cell-1][i_cell].z + leta_arr[k_cell][j_cell][i_cell].z); //* aj_center;
2721
2722 mat[2][0] = 0.5 * (lzet_arr[k_cell-1][j_cell][i_cell].x + lzet_arr[k_cell][j_cell][i_cell].x); //* aj_center;
2723 mat[2][1] = 0.5 * (lzet_arr[k_cell-1][j_cell][i_cell].y + lzet_arr[k_cell][j_cell][i_cell].y); //* aj_center;
2724 mat[2][2] = 0.5 * (lzet_arr[k_cell-1][j_cell][i_cell].z + lzet_arr[k_cell][j_cell][i_cell].z); //* aj_center;
2725
2726 // Contravariant velocity vector `q` at the cell center,
2727 // estimated by averaging face-based contravariant velocities.
2728 PetscReal q[3];
2729 q[0] = 0.5 * (lucont_arr[k_cell][j_cell][i_cell-1].x + lucont_arr[k_cell][j_cell][i_cell].x); // U¹ at cell center
2730 q[1] = 0.5 * (lucont_arr[k_cell][j_cell-1][i_cell].y + lucont_arr[k_cell][j_cell][i_cell].y); // U² at cell center
2731 q[2] = 0.5 * (lucont_arr[k_cell-1][j_cell][i_cell].z + lucont_arr[k_cell][j_cell][i_cell].z); // U³ at cell center
2732
2733 // Solve the 3x3 system `mat * ucat = q` using Cramer's rule.
2734 PetscReal det = mat[0][0] * (mat[1][1] * mat[2][2] - mat[1][2] * mat[2][1]) -
2735 mat[0][1] * (mat[1][0] * mat[2][2] - mat[1][2] * mat[2][0]) +
2736 mat[0][2] * (mat[1][0] * mat[2][1] - mat[1][1] * mat[2][0]);
2737
2738 if (PetscAbsReal(det) < 1.0e-18) {
2739 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FLOP_COUNT, "Transformation matrix determinant is near zero at cell (%d,%d,%d) \n", i_cell, j_cell, k_cell);
2740 }
2741
2742 PetscReal det_inv = 1.0 / det;
2743
2744 PetscReal det0 = q[0] * (mat[1][1] * mat[2][2] - mat[1][2] * mat[2][1]) -
2745 q[1] * (mat[0][1] * mat[2][2] - mat[0][2] * mat[2][1]) +
2746 q[2] * (mat[0][1] * mat[1][2] - mat[0][2] * mat[1][1]);
2747
2748 PetscReal det1 = -q[0] * (mat[1][0] * mat[2][2] - mat[1][2] * mat[2][0]) +
2749 q[1] * (mat[0][0] * mat[2][2] - mat[0][2] * mat[2][0]) -
2750 q[2] * (mat[0][0] * mat[1][2] - mat[0][2] * mat[1][0]);
2751
2752 PetscReal det2 = q[0] * (mat[1][0] * mat[2][1] - mat[1][1] * mat[2][0]) -
2753 q[1] * (mat[0][0] * mat[2][1] - mat[0][1] * mat[2][0]) +
2754 q[2] * (mat[0][0] * mat[1][1] - mat[0][1] * mat[1][0]);
2755
2756 // Store computed Cartesian velocity in the GLOBAL Ucat array at the
2757 // array index corresponding to the cell's origin node.
2758 gucat_arr[k_cell][j_cell][i_cell].x = det0 * det_inv;
2759 gucat_arr[k_cell][j_cell][i_cell].y = det1 * det_inv;
2760 gucat_arr[k_cell][j_cell][i_cell].z = det2 * det_inv;
2761 }
2762 }
2763 }
2764
2765 // --- 6. Restore Array Access ---
2766 ierr = DMDAVecRestoreArrayRead(user->fda, user->lUcont, &lucont_arr); CHKERRQ(ierr);
2767 ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, &lcsi_arr); CHKERRQ(ierr);
2768 ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, &leta_arr); CHKERRQ(ierr);
2769 ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, &lzet_arr); CHKERRQ(ierr);
2770 ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, &lnvert_arr); CHKERRQ(ierr);
2771 ierr = DMDAVecRestoreArrayRead(user->da, user->lAj, &laj_arr); CHKERRQ(ierr);
2772 ierr = DMDAVecRestoreArray(user->fda, user->Ucat, &gucat_arr); CHKERRQ(ierr);
2773
2774 LOG_ALLOW(GLOBAL, LOG_INFO, "Completed Contravariant-to-Cartesian velocity transformation. \n");
2775 PROFILE_FUNCTION_END;
2776 PetscFunctionReturn(0);
2777}
2778
2779#undef __FUNCT__
2780#define __FUNCT__ "Cart2Contra"
2781/**
2782 * @brief Convert a spatially varying Cartesian velocity field to contravariant fluxes.
2783 */
2784PetscErrorCode Cart2Contra(UserCtx *user)
2785{
2786 PetscErrorCode ierr;
2787 DMDALocalInfo info;
2788 const Cmpnts ***ucat_arr, ***csi_arr, ***eta_arr, ***zet_arr;
2789 Cmpnts ***ucont_arr;
2790
2791 PetscFunctionBeginUser;
2792 PROFILE_FUNCTION_BEGIN;
2793
2794 ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
2795 ierr = DMDAVecGetArrayRead(user->fda, user->lUcat, &ucat_arr); CHKERRQ(ierr);
2796 ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, &csi_arr); CHKERRQ(ierr);
2797 ierr = DMDAVecGetArrayRead(user->fda, user->lEta, &eta_arr); CHKERRQ(ierr);
2798 ierr = DMDAVecGetArrayRead(user->fda, user->lZet, &zet_arr); CHKERRQ(ierr);
2799 ierr = DMDAVecGetArray(user->fda, user->Ucont, &ucont_arr); CHKERRQ(ierr);
2800
2801 const PetscInt i_start = PetscMax(info.xs, 1);
2802 const PetscInt j_start = PetscMax(info.ys, 1);
2803 const PetscInt k_start = PetscMax(info.zs, 1);
2804 const PetscInt i_end = PetscMin(info.xs + info.xm, info.mx - 1);
2805 const PetscInt j_end = PetscMin(info.ys + info.ym, info.my - 1);
2806 const PetscInt k_end = PetscMin(info.zs + info.zm, info.mz - 1);
2807
2808 for (PetscInt k = k_start; k < k_end; k++) {
2809 for (PetscInt j = j_start; j < j_end; j++) {
2810 for (PetscInt i = i_start; i < i_end; i++) {
2811 const Cmpnts u_xi = {
2812 0.5 * (ucat_arr[k][j][i].x + ucat_arr[k][j][i + 1].x),
2813 0.5 * (ucat_arr[k][j][i].y + ucat_arr[k][j][i + 1].y),
2814 0.5 * (ucat_arr[k][j][i].z + ucat_arr[k][j][i + 1].z)
2815 };
2816 const Cmpnts u_eta = {
2817 0.5 * (ucat_arr[k][j][i].x + ucat_arr[k][j + 1][i].x),
2818 0.5 * (ucat_arr[k][j][i].y + ucat_arr[k][j + 1][i].y),
2819 0.5 * (ucat_arr[k][j][i].z + ucat_arr[k][j + 1][i].z)
2820 };
2821 const Cmpnts u_zeta = {
2822 0.5 * (ucat_arr[k][j][i].x + ucat_arr[k + 1][j][i].x),
2823 0.5 * (ucat_arr[k][j][i].y + ucat_arr[k + 1][j][i].y),
2824 0.5 * (ucat_arr[k][j][i].z + ucat_arr[k + 1][j][i].z)
2825 };
2826 ucont_arr[k][j][i].x = csi_arr[k][j][i].x * u_xi.x + csi_arr[k][j][i].y * u_xi.y + csi_arr[k][j][i].z * u_xi.z;
2827 ucont_arr[k][j][i].y = eta_arr[k][j][i].x * u_eta.x + eta_arr[k][j][i].y * u_eta.y + eta_arr[k][j][i].z * u_eta.z;
2828 ucont_arr[k][j][i].z = zet_arr[k][j][i].x * u_zeta.x + zet_arr[k][j][i].y * u_zeta.y + zet_arr[k][j][i].z * u_zeta.z;
2829 }
2830 }
2831 }
2832
2833 ierr = DMDAVecRestoreArray(user->fda, user->Ucont, &ucont_arr); CHKERRQ(ierr);
2834 ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, &zet_arr); CHKERRQ(ierr);
2835 ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, &eta_arr); CHKERRQ(ierr);
2836 ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, &csi_arr); CHKERRQ(ierr);
2837 ierr = DMDAVecRestoreArrayRead(user->fda, user->lUcat, &ucat_arr); CHKERRQ(ierr);
2838
2839 PROFILE_FUNCTION_END;
2840 PetscFunctionReturn(0);
2841}
2842
2843#undef __FUNCT__
2844#define __FUNCT__ "UniformCart2Contra"
2845/*
2846 * Converts a uniform Cartesian velocity (u,v,w) to contravariant fluxes in Ucont.
2847 * Computes the dot product of the physical velocity with each face-area vector:
2848 * U^xi = csi · (u,v,w), U^eta = eta · (u,v,w), U^zeta = zet · (u,v,w).
2849 * Writes to all owned nodes (xs..xe, ys..ye, zs..ze); boundary ghosts are
2850 * overwritten later by ApplyBoundaryConditions.
2851 * The public header owns the rendered API contract.
2852 */
2853PetscErrorCode UniformCart2Contra(UserCtx *user, PetscReal u, PetscReal v, PetscReal w)
2854{
2855 PetscErrorCode ierr;
2856 PetscFunctionBeginUser;
2857 PROFILE_FUNCTION_BEGIN;
2858
2859 DMDALocalInfo info;
2860 Cmpnts ***ucont_arr;
2861 const Cmpnts ***csi_arr, ***eta_arr, ***zet_arr;
2862
2863 ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
2864 ierr = DMDAVecGetArray(user->fda, user->Ucont, &ucont_arr); CHKERRQ(ierr);
2865 ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, &csi_arr); CHKERRQ(ierr);
2866 ierr = DMDAVecGetArrayRead(user->fda, user->lEta, &eta_arr); CHKERRQ(ierr);
2867 ierr = DMDAVecGetArrayRead(user->fda, user->lZet, &zet_arr); CHKERRQ(ierr);
2868
2869 const PetscInt xs = info.xs, xe = info.xs + info.xm;
2870 const PetscInt ys = info.ys, ye = info.ys + info.ym;
2871 const PetscInt zs = info.zs, ze = info.zs + info.zm;
2872
2873 for (PetscInt k = zs; k < ze; k++) {
2874 for (PetscInt j = ys; j < ye; j++) {
2875 for (PetscInt i = xs; i < xe; i++) {
2876 ucont_arr[k][j][i].x = csi_arr[k][j][i].x * u + csi_arr[k][j][i].y * v + csi_arr[k][j][i].z * w;
2877 ucont_arr[k][j][i].y = eta_arr[k][j][i].x * u + eta_arr[k][j][i].y * v + eta_arr[k][j][i].z * w;
2878 ucont_arr[k][j][i].z = zet_arr[k][j][i].x * u + zet_arr[k][j][i].y * v + zet_arr[k][j][i].z * w;
2879 }
2880 }
2881 }
2882
2883 ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, &zet_arr); CHKERRQ(ierr);
2884 ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, &eta_arr); CHKERRQ(ierr);
2885 ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, &csi_arr); CHKERRQ(ierr);
2886 ierr = DMDAVecRestoreArray(user->fda, user->Ucont, &ucont_arr); CHKERRQ(ierr);
2887
2888 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Cart2Contra: set Ucont from Cartesian (%.3f, %.3f, %.3f).\n",
2889 (double)u, (double)v, (double)w);
2890 PROFILE_FUNCTION_END;
2891 PetscFunctionReturn(0);
2892}
2893
2894#undef __FUNCT__
2895#define __FUNCT__ "SetupDomainCellDecompositionMap"
2896/**
2897 * @brief Internal helper implementation: `SetupDomainCellDecompositionMap()`.
2898 * @details Local to this translation unit.
2899 */
2900PetscErrorCode SetupDomainCellDecompositionMap(UserCtx *user)
2901{
2902 PetscErrorCode ierr;
2903 DMDALocalInfo local_node_info;
2904 RankCellInfo my_cell_info;
2905 PetscMPIInt rank, size;
2906
2907 PetscFunctionBeginUser;
2908 PROFILE_FUNCTION_BEGIN;
2909
2910 // --- 1. Input Validation and MPI Info ---
2911 if (!user) {
2912 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "UserCtx pointer is NULL in SetupDomainCellDecompositionMap.");
2913 }
2914 if (!user->da) {
2915 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "user->da is not initialized in SetupDomainCellDecompositionMap.");
2916 }
2917
2918 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
2919 ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size); CHKERRQ(ierr);
2920
2921 LOG_ALLOW(GLOBAL, LOG_INFO, "Setting up domain cell decomposition map for %d ranks.\n", size);
2922
2923 // --- 2. Determine Local Cell Ownership ---
2924 // Get the local node ownership information from the primary DMDA.
2925 ierr = DMDAGetLocalInfo(user->da, &local_node_info); CHKERRQ(ierr);
2926
2927 // Use the robust helper function to convert node ownership to cell ownership.
2928 // A cell's index is defined by its origin node.
2929
2930 ierr = GetOwnedCellRange(&local_node_info, 0, &my_cell_info.xs_cell, &my_cell_info.xm_cell); CHKERRQ(ierr);
2931 ierr = GetOwnedCellRange(&local_node_info, 1, &my_cell_info.ys_cell, &my_cell_info.ym_cell); CHKERRQ(ierr);
2932 ierr = GetOwnedCellRange(&local_node_info, 2, &my_cell_info.zs_cell, &my_cell_info.zm_cell); CHKERRQ(ierr);
2933
2934 // Log the calculated local ownership for debugging purposes.
2935 LOG_ALLOW(LOCAL, LOG_DEBUG, "[Rank %d] Owns cells: i[%d, %d), j[%d, %d), k[%d, %d)\n",
2936 rank, my_cell_info.xs_cell, my_cell_info.xs_cell + my_cell_info.xm_cell,
2937 my_cell_info.ys_cell, my_cell_info.ys_cell + my_cell_info.ym_cell,
2938 my_cell_info.zs_cell, my_cell_info.zs_cell + my_cell_info.zm_cell);
2939
2940 // --- 3. Allocate and Distribute the Global Map ---
2941 // Allocate memory for the global map that will hold information from all ranks.
2942 ierr = PetscMalloc1(size, &user->RankCellInfoMap); CHKERRQ(ierr);
2943
2944 // Perform the collective communication to gather the `RankCellInfo` struct from every rank.
2945 // Each rank sends its `my_cell_info` and receives the complete array in `user->RankCellInfoMap`.
2946 // We use MPI_BYTE to ensure portability across different systems and struct padding.
2947 ierr = MPI_Allgather(&my_cell_info, sizeof(RankCellInfo), MPI_BYTE,
2948 user->RankCellInfoMap, sizeof(RankCellInfo), MPI_BYTE,
2949 PETSC_COMM_WORLD); CHKERRQ(ierr);
2950
2951 LOG_ALLOW(GLOBAL, LOG_INFO, "Domain cell decomposition map created and distributed successfully.\n");
2952
2953 PROFILE_FUNCTION_END;
2954 PetscFunctionReturn(0);
2955}
2956
2957#undef __FUNCT__
2958#define __FUNCT__ "BinarySearchInt64"
2959/**
2960 * @brief Implementation of \ref BinarySearchInt64().
2961 * @details Full API contract (arguments, ownership, side effects) is documented with
2962 * the header declaration in `include/setup.h`.
2963 * @see BinarySearchInt64()
2964 */
2965PetscErrorCode BinarySearchInt64(PetscInt n, const PetscInt64 arr[], PetscInt64 key, PetscBool *found)
2966{
2967 PetscInt low = 0, high = n - 1;
2968
2969 PetscFunctionBeginUser;
2970 PROFILE_FUNCTION_BEGIN;
2971
2972 // --- 1. Input Validation ---
2973 if (!found) {
2974 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Output pointer 'found' is NULL in PetscBinarySearchInt64.");
2975 }
2976 if (n > 0 && !arr) {
2977 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Input array 'arr' is NULL for n > 0.");
2978 }
2979
2980 // Initialize output
2981 *found = PETSC_FALSE;
2982
2983 // --- 2. Binary Search Algorithm ---
2984 while (low <= high) {
2985 // Use this form to prevent potential integer overflow on very large arrays
2986 PetscInt mid = low + (high - low) / 2;
2987
2988 if (arr[mid] == key) {
2989 *found = PETSC_TRUE; // Key found!
2990 break; // Exit the loop
2991 }
2992
2993 if (arr[mid] < key) {
2994 low = mid + 1; // Search in the right half
2995 } else {
2996 high = mid - 1; // Search in the left half
2997 }
2998 }
2999
3000 PROFILE_FUNCTION_END;
3001 PetscFunctionReturn(0);
3002}
3003
3004
3005/**
3006 * @brief Convert logical indices into the flattened global index used by setup helpers.
3007 */
3008static PetscInt Gidx(PetscInt i, PetscInt j, PetscInt k, UserCtx *user)
3009{
3010 PetscInt nidx;
3011 DMDALocalInfo info = user->info;
3012
3013 PetscInt mx = info.mx, my = info.my;
3014
3015 AO ao;
3016 DMDAGetAO(user->da, &ao);
3017 nidx=i+j*mx+k*mx*my;
3018
3019 AOApplicationToPetsc(ao,1,&nidx);
3020
3021 return (nidx);
3022}
3023
3024
3025#undef __FUNCT__
3026#define __FUNCT__ "ComputeDivergence"
3027/**
3028 * @brief Implementation of \ref ComputeDivergence().
3029 * @details Full API contract (arguments, ownership, side effects) is documented with
3030 * the header declaration in `include/setup.h`.
3031 * @see ComputeDivergence()
3032 */
3033
3034PetscErrorCode ComputeDivergence(UserCtx *user )
3035{
3036 DM da = user->da, fda = user->fda;
3037 DMDALocalInfo info = user->info;
3038
3039 PetscInt ti = user->simCtx->step;
3040
3041 PetscInt xs = info.xs, xe = info.xs + info.xm;
3042 PetscInt ys = info.ys, ye = info.ys + info.ym;
3043 PetscInt zs = info.zs, ze = info.zs + info.zm;
3044 PetscInt mx = info.mx, my = info.my, mz = info.mz;
3045
3046 PetscInt lxs, lys, lzs, lxe, lye, lze;
3047 PetscInt i, j, k;
3048
3049 Vec Div;
3050 PetscReal ***div, ***aj, ***nvert,***p;
3051 Cmpnts ***ucont;
3052 PetscReal maxdiv;
3053
3054 lxs = xs; lxe = xe;
3055 lys = ys; lye = ye;
3056 lzs = zs; lze = ze;
3057
3058 if (xs==0) lxs = xs+1;
3059 if (ys==0) lys = ys+1;
3060 if (zs==0) lzs = zs+1;
3061
3062 if (xe==mx) lxe = xe-1;
3063 if (ye==my) lye = ye-1;
3064 if (ze==mz) lze = ze-1;
3065
3066 PetscFunctionBeginUser;
3067 PROFILE_FUNCTION_BEGIN;
3068
3069 DMDAVecGetArray(fda,user->lUcont, &ucont);
3070 DMDAVecGetArray(da, user->lAj, &aj);
3071 VecDuplicate(user->P, &Div);
3072 DMDAVecGetArray(da, Div, &div);
3073 DMDAVecGetArray(da, user->lNvert, &nvert);
3074 DMDAVecGetArray(da, user->P, &p);
3075 for (k=lzs; k<lze; k++) {
3076 for (j=lys; j<lye; j++){
3077 for (i=lxs; i<lxe; i++) {
3078 if (k==10 && j==10 && i==1){
3079 LOG_ALLOW(LOCAL,LOG_INFO,"Pressure[10][10][1] = %f | Pressure[10][10][0] = %f \n ",p[k][j][i],p[k][j][i-1]);
3080 }
3081
3082 if (k==10 && j==10 && i==mx-3)
3083 LOG_ALLOW(LOCAL,LOG_INFO,"Pressure[10][10][%d] = %f | Pressure[10][10][%d] = %f \n ",mx-2,p[k][j][mx-2],mx-1,p[k][j][mx-1]);
3084 }
3085 }
3086 }
3087 DMDAVecRestoreArray(da, user->P, &p);
3088
3089
3090 for (k=lzs; k<lze; k++) {
3091 for (j=lys; j<lye; j++) {
3092 for (i=lxs; i<lxe; i++) {
3093 maxdiv = fabs((ucont[k][j][i].x - ucont[k][j][i-1].x +
3094 ucont[k][j][i].y - ucont[k][j-1][i].y +
3095 ucont[k][j][i].z - ucont[k-1][j][i].z)*aj[k][j][i]);
3096 if (nvert[k][j][i] + nvert[k+1][j][i] + nvert[k-1][j][i] +
3097 nvert[k][j+1][i] + nvert[k][j-1][i] +
3098 nvert[k][j][i+1] + nvert[k][j][i-1] > 0.1) maxdiv = 0.;
3099 div[k][j][i] = maxdiv;
3100
3101 }
3102 }
3103 }
3104
3105 if (zs==0) {
3106 k=0;
3107 for (j=ys; j<ye; j++) {
3108 for (i=xs; i<xe; i++) {
3109 div[k][j][i] = 0.;
3110 }
3111 }
3112 }
3113
3114 if (ze == mz) {
3115 k=mz-1;
3116 for (j=ys; j<ye; j++) {
3117 for (i=xs; i<xe; i++) {
3118 div[k][j][i] = 0.;
3119 }
3120 }
3121 }
3122
3123 if (xs==0) {
3124 i=0;
3125 for (k=zs; k<ze; k++) {
3126 for (j=ys; j<ye; j++) {
3127 div[k][j][i] = 0.;
3128 }
3129 }
3130 }
3131
3132 if (xe==mx) {
3133 i=mx-1;
3134 for (k=zs; k<ze; k++) {
3135 for (j=ys; j<ye; j++) {
3136 div[k][j][i] = 0;
3137 }
3138 }
3139 }
3140
3141 if (ys==0) {
3142 j=0;
3143 for (k=zs; k<ze; k++) {
3144 for (i=xs; i<xe; i++) {
3145 div[k][j][i] = 0.;
3146 }
3147 }
3148 }
3149
3150 if (ye==my) {
3151 j=my-1;
3152 for (k=zs; k<ze; k++) {
3153 for (i=xs; i<xe; i++) {
3154 div[k][j][i] = 0.;
3155 }
3156 }
3157 }
3158 DMDAVecRestoreArray(da, Div, &div);
3159 PetscInt MaxFlatIndex;
3160
3161 VecMax(Div, &MaxFlatIndex, &maxdiv);
3162
3163 LOG_ALLOW(GLOBAL,LOG_INFO,"[Step %d]] The Maximum Divergence is %e at flat index %d.\n",ti,maxdiv,MaxFlatIndex);
3164
3165 user->simCtx->MaxDivFlatArg = MaxFlatIndex;
3166 user->simCtx->MaxDiv = maxdiv;
3167
3168 for (k=zs; k<ze; k++) {
3169 for (j=ys; j<ye; j++) {
3170 for (i=xs; i<xe; i++) {
3171 if (Gidx(i,j,k,user) == MaxFlatIndex) {
3172 LOG_ALLOW(GLOBAL,LOG_INFO,"[Step %d] The Maximum Divergence(%e) is at location [%d][%d][%d]. \n", ti, maxdiv,k,j,i);
3173 user->simCtx->MaxDivz = k;
3174 user->simCtx->MaxDivy = j;
3175 user->simCtx->MaxDivx = i;
3176 }
3177 }
3178 }
3179 }
3180
3181
3182 DMDAVecRestoreArray(da, user->lNvert, &nvert);
3183 DMDAVecRestoreArray(fda, user->lUcont, &ucont);
3184 DMDAVecRestoreArray(da, user->lAj, &aj);
3185 VecDestroy(&Div);
3186
3187 PROFILE_FUNCTION_END;
3188 PetscFunctionReturn(0);
3189}
3190
3191#undef __FUNCT__
3192#define __FUNCT__ "InitializeRandomGenerators"
3193
3194/**
3195 * @brief Implementation of \ref InitializeRandomGenerators().
3196 * @details Full API contract (arguments, ownership, side effects) is documented with
3197 * the header declaration in `include/setup.h`.
3198 * @see InitializeRandomGenerators()
3199 */
3200PetscErrorCode InitializeRandomGenerators(UserCtx* user, PetscRandom *randx, PetscRandom *randy, PetscRandom *randz) {
3201 PetscErrorCode ierr; // Error code for PETSc functions
3202 PetscMPIInt rank;
3203 PetscFunctionBeginUser;
3204 PROFILE_FUNCTION_BEGIN;
3205 MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
3206
3207 // Initialize RNG for x-coordinate
3208 ierr = PetscRandomCreate(PETSC_COMM_SELF, randx); CHKERRQ(ierr);
3209 ierr = PetscRandomSetType((*randx), PETSCRAND48); CHKERRQ(ierr);
3210 ierr = PetscRandomSetInterval(*randx, user->bbox.min_coords.x, user->bbox.max_coords.x); CHKERRQ(ierr);
3211 ierr = PetscRandomSetSeed(*randx, rank + 12345); CHKERRQ(ierr); // Unique seed per rank
3212 ierr = PetscRandomSeed(*randx); CHKERRQ(ierr);
3213 LOG_ALLOW_SYNC(LOCAL,LOG_VERBOSE, "[Rank %d]Initialized RNG for X-axis.\n",rank);
3214
3215 // Initialize RNG for y-coordinate
3216 ierr = PetscRandomCreate(PETSC_COMM_SELF, randy); CHKERRQ(ierr);
3217 ierr = PetscRandomSetType((*randy), PETSCRAND48); CHKERRQ(ierr);
3218 ierr = PetscRandomSetInterval(*randy, user->bbox.min_coords.y, user->bbox.max_coords.y); CHKERRQ(ierr);
3219 ierr = PetscRandomSetSeed(*randy, rank + 67890); CHKERRQ(ierr); // Unique seed per rank
3220 ierr = PetscRandomSeed(*randy); CHKERRQ(ierr);
3221 LOG_ALLOW_SYNC(LOCAL,LOG_VERBOSE, "[Rank %d]Initialized RNG for Y-axis.\n",rank);
3222
3223 // Initialize RNG for z-coordinate
3224 ierr = PetscRandomCreate(PETSC_COMM_SELF, randz); CHKERRQ(ierr);
3225 ierr = PetscRandomSetType((*randz), PETSCRAND48); CHKERRQ(ierr);
3226 ierr = PetscRandomSetInterval(*randz, user->bbox.min_coords.z, user->bbox.max_coords.z); CHKERRQ(ierr);
3227 ierr = PetscRandomSetSeed(*randz, rank + 54321); CHKERRQ(ierr); // Unique seed per rank
3228 ierr = PetscRandomSeed(*randz); CHKERRQ(ierr);
3229 LOG_ALLOW_SYNC(LOCAL,LOG_VERBOSE, "[Rank %d]Initialized RNG for Z-axis.\n",rank);
3230
3231 PROFILE_FUNCTION_END;
3232 PetscFunctionReturn(0);
3233}
3234
3235#undef __FUNCT__
3236#define __FUNCT__ "InitializeLogicalSpaceRNGs"
3237/**
3238 * @brief Internal helper implementation: `InitializeLogicalSpaceRNGs()`.
3239 * @details Local to this translation unit.
3240 */
3241PetscErrorCode InitializeLogicalSpaceRNGs(PetscRandom *rand_logic_i, PetscRandom *rand_logic_j, PetscRandom *rand_logic_k) {
3242 PetscErrorCode ierr;
3243 PetscMPIInt rank;
3244 PetscFunctionBeginUser;
3245
3246 PROFILE_FUNCTION_BEGIN;
3247
3248 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
3249
3250 // --- RNG for i-logical dimension ---
3251 ierr = PetscRandomCreate(PETSC_COMM_SELF, rand_logic_i); CHKERRQ(ierr);
3252 ierr = PetscRandomSetType((*rand_logic_i), PETSCRAND48); CHKERRQ(ierr);
3253 ierr = PetscRandomSetInterval(*rand_logic_i, 0.0, 1.0); CHKERRQ(ierr); // Key change: [0,1)
3254 ierr = PetscRandomSetSeed(*rand_logic_i, rank + 202401); CHKERRQ(ierr); // Unique seed
3255 ierr = PetscRandomSeed(*rand_logic_i); CHKERRQ(ierr);
3256 LOG_ALLOW(LOCAL,LOG_VERBOSE, "[Rank %d] Initialized RNG for i-logical dimension [0,1).\n",rank);
3257
3258 // --- RNG for j-logical dimension ---
3259 ierr = PetscRandomCreate(PETSC_COMM_SELF, rand_logic_j); CHKERRQ(ierr);
3260 ierr = PetscRandomSetType((*rand_logic_j), PETSCRAND48); CHKERRQ(ierr);
3261 ierr = PetscRandomSetInterval(*rand_logic_j, 0.0, 1.0); CHKERRQ(ierr); // Key change: [0,1)
3262 ierr = PetscRandomSetSeed(*rand_logic_j, rank + 202402); CHKERRQ(ierr);
3263 ierr = PetscRandomSeed(*rand_logic_j); CHKERRQ(ierr);
3264 LOG_ALLOW(LOCAL,LOG_VERBOSE, "[Rank %d] Initialized RNG for j-logical dimension [0,1).\n",rank);
3265
3266 // --- RNG for k-logical dimension ---
3267 ierr = PetscRandomCreate(PETSC_COMM_SELF, rand_logic_k); CHKERRQ(ierr);
3268 ierr = PetscRandomSetType((*rand_logic_k), PETSCRAND48); CHKERRQ(ierr);
3269 ierr = PetscRandomSetInterval(*rand_logic_k, 0.0, 1.0); CHKERRQ(ierr); // Key change: [0,1)
3270 ierr = PetscRandomSetSeed(*rand_logic_k, rank + 202403); CHKERRQ(ierr);
3271 ierr = PetscRandomSeed(*rand_logic_k); CHKERRQ(ierr);
3272 LOG_ALLOW(LOCAL,LOG_VERBOSE, "[Rank %d] Initialized RNG for k-logical dimension [0,1).\n",rank);
3273
3274
3275 PROFILE_FUNCTION_END;
3276 PetscFunctionReturn(0);
3277}
3278
3279#undef __FUNCT__
3280#define __FUNCT__ "InitializeBrownianRNG"
3281/**
3282 * @brief Internal helper implementation: `InitializeBrownianRNG()`.
3283 * @details Local to this translation unit.
3284 */
3285PetscErrorCode InitializeBrownianRNG(SimCtx *simCtx) {
3286 PetscErrorCode ierr;
3287 PetscMPIInt rank;
3288
3289 PetscFunctionBeginUser;
3290 PROFILE_FUNCTION_BEGIN;
3291
3292 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
3293
3294 // 1. Create the generator (stored in SimCtx, not UserCtx, as it is global physics)
3295 ierr = PetscRandomCreate(PETSC_COMM_WORLD, &simCtx->BrownianMotionRNG); CHKERRQ(ierr);
3296 ierr = PetscRandomSetType(simCtx->BrownianMotionRNG, PETSCRAND48); CHKERRQ(ierr);
3297
3298 // 2. CRITICAL: Set interval to [0, 1).
3299 // This is required for the Gaussian math to work.
3300 ierr = PetscRandomSetInterval(simCtx->BrownianMotionRNG, 0.0, 1.0); CHKERRQ(ierr);
3301
3302 // 3. Seed based on Rank to ensure spatial randomness
3303 // Multiplying by a large prime helps separate the streams significantly
3304 unsigned long seed = (unsigned long)rank * 987654321 + (unsigned long)time(NULL);
3305 ierr = PetscRandomSetSeed(simCtx->BrownianMotionRNG, seed); CHKERRQ(ierr);
3306 ierr = PetscRandomSeed(simCtx->BrownianMotionRNG); CHKERRQ(ierr);
3307
3308 LOG_ALLOW(LOCAL, LOG_VERBOSE, "[Rank %d] Initialized Brownian Physics RNG.\n", rank);
3309
3310 PROFILE_FUNCTION_END;
3311 PetscFunctionReturn(0);
3312}
3313
3314/////////////// DERIVATIVE CALCULATION HELPERS ///////////////
3315
3316#undef __FUNCT__
3317#define __FUNCT__ "TransformScalarDerivativesToPhysical"
3318/**
3319 * @brief Implementation of \ref TransformScalarDerivativesToPhysical().
3320 * @details Full API contract (arguments, ownership, side effects) is documented with
3321 * the header declaration in `include/setup.h`.
3322 * @see TransformScalarDerivativesToPhysical()
3323 */
3324 void TransformScalarDerivativesToPhysical(PetscReal jacobian,
3325 Cmpnts csi_metrics,
3326 Cmpnts eta_metrics,
3327 Cmpnts zet_metrics,
3328 PetscReal dPhi_dcsi,
3329 PetscReal dPhi_deta,
3330 PetscReal dPhi_dzet,
3331 Cmpnts *gradPhi)
3332{
3333 // Gradient X component
3334 gradPhi->x = jacobian * (dPhi_dcsi * csi_metrics.x + dPhi_deta * eta_metrics.x + dPhi_dzet * zet_metrics.x);
3335
3336 // Gradient Y component
3337 gradPhi->y = jacobian * (dPhi_dcsi * csi_metrics.y + dPhi_deta * eta_metrics.y + dPhi_dzet * zet_metrics.y);
3338
3339 // Gradient Z component
3340 gradPhi->z = jacobian * (dPhi_dcsi * csi_metrics.z + dPhi_deta * eta_metrics.z + dPhi_dzet * zet_metrics.z);
3341}
3342
3343#undef __FUNCT__
3344#define __FUNCT__ "TransformDerivativesToPhysical"
3345/**
3346 * @brief Transform contravariant vector derivatives into physical Cartesian derivatives.
3347 */
3348static void TransformDerivativesToPhysical(PetscReal jacobian, Cmpnts csi_metrics, Cmpnts eta_metrics, Cmpnts zet_metrics,
3349 Cmpnts deriv_csi, Cmpnts deriv_eta, Cmpnts deriv_zet,
3350 Cmpnts *dudx, Cmpnts *dvdx, Cmpnts *dwdx)
3351{
3352 // Derivatives of the first component (u)
3353 dudx->x = jacobian * (deriv_csi.x * csi_metrics.x + deriv_eta.x * eta_metrics.x + deriv_zet.x * zet_metrics.x);
3354 dudx->y = jacobian * (deriv_csi.x * csi_metrics.y + deriv_eta.x * eta_metrics.y + deriv_zet.x * zet_metrics.y);
3355 dudx->z = jacobian * (deriv_csi.x * csi_metrics.z + deriv_eta.x * eta_metrics.z + deriv_zet.x * zet_metrics.z);
3356 // Derivatives of the second component (v)
3357 dvdx->x = jacobian * (deriv_csi.y * csi_metrics.x + deriv_eta.y * eta_metrics.x + deriv_zet.y * zet_metrics.x);
3358 dvdx->y = jacobian * (deriv_csi.y * csi_metrics.y + deriv_eta.y * eta_metrics.y + deriv_zet.y * zet_metrics.y);
3359 dvdx->z = jacobian * (deriv_csi.y * csi_metrics.z + deriv_eta.y * eta_metrics.z + deriv_zet.y * zet_metrics.z);
3360 // Derivatives of the third component (w)
3361 dwdx->x = jacobian * (deriv_csi.z * csi_metrics.x + deriv_eta.z * eta_metrics.x + deriv_zet.z * zet_metrics.x);
3362 dwdx->y = jacobian * (deriv_csi.z * csi_metrics.y + deriv_eta.z * eta_metrics.y + deriv_zet.z * zet_metrics.y);
3363 dwdx->z = jacobian * (deriv_csi.z * csi_metrics.z + deriv_eta.z * eta_metrics.z + deriv_zet.z * zet_metrics.z);
3364}
3365
3366#undef __FUNCT__
3367#define __FUNCT__ "ComputeScalarFieldDerivatives"
3368/**
3369 * @brief Internal helper implementation: `ComputeScalarFieldDerivatives()`.
3370 * @details Local to this translation unit.
3371 */
3372PetscErrorCode ComputeScalarFieldDerivatives(UserCtx *user, PetscInt i, PetscInt j, PetscInt k,
3373 PetscReal ***field_data, Cmpnts *grad)
3374{
3375 PetscErrorCode ierr;
3376 Cmpnts ***csi, ***eta, ***zet;
3377 PetscReal ***jac;
3378 PetscReal d_csi, d_eta, d_zet;
3379
3380 PetscFunctionBeginUser;
3381
3382 // 1. Get read-only access to metrics
3383 ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, &csi); CHKERRQ(ierr);
3384 ierr = DMDAVecGetArrayRead(user->fda, user->lEta, &eta); CHKERRQ(ierr);
3385 ierr = DMDAVecGetArrayRead(user->fda, user->lZet, &zet); CHKERRQ(ierr);
3386 ierr = DMDAVecGetArrayRead(user->da, user->lAj, &jac); CHKERRQ(ierr);
3387
3388 // 2. Compute derivatives in computational space (Central Difference)
3389 // Assumes ghosts are available at i+/-1
3390 d_csi = 0.5 * (field_data[k][j][i+1] - field_data[k][j][i-1]);
3391 d_eta = 0.5 * (field_data[k][j+1][i] - field_data[k][j-1][i]);
3392 d_zet = 0.5 * (field_data[k+1][j][i] - field_data[k-1][j][i]);
3393
3394 // 3. Transform to physical space
3395 TransformScalarDerivativesToPhysical(jac[k][j][i],
3396 csi[k][j][i], eta[k][j][i], zet[k][j][i],
3397 d_csi, d_eta, d_zet,
3398 grad);
3399
3400 // 4. Restore arrays
3401 ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, &csi); CHKERRQ(ierr);
3402 ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, &eta); CHKERRQ(ierr);
3403 ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, &zet); CHKERRQ(ierr);
3404 ierr = DMDAVecRestoreArrayRead(user->da, user->lAj, &jac); CHKERRQ(ierr);
3405
3406 PetscFunctionReturn(0);
3407}
3408
3409#undef __FUNCT__
3410#define __FUNCT__ "ComputeVectorFieldDerivatives"
3411/**
3412 * @brief Internal helper implementation: `ComputeVectorFieldDerivatives()`.
3413 * @details Local to this translation unit.
3414 */
3415PetscErrorCode ComputeVectorFieldDerivatives(UserCtx *user, PetscInt i, PetscInt j, PetscInt k, Cmpnts ***field_data,
3416 Cmpnts *dudx, Cmpnts *dvdx, Cmpnts *dwdx)
3417{
3418 PetscErrorCode ierr;
3419 Cmpnts ***csi, ***eta, ***zet;
3420 PetscReal ***jac;
3421 PetscFunctionBeginUser;
3422
3423 // 1. Get read-only access to the necessary metric data arrays
3424 ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, &csi); CHKERRQ(ierr);
3425 ierr = DMDAVecGetArrayRead(user->fda, user->lEta, &eta); CHKERRQ(ierr);
3426 ierr = DMDAVecGetArrayRead(user->fda, user->lZet, &zet); CHKERRQ(ierr);
3427 ierr = DMDAVecGetArrayRead(user->da, user->lAj, &jac); CHKERRQ(ierr);
3428
3429 // 2. Calculate derivatives in computational space using central differencing
3430 Cmpnts deriv_csi, deriv_eta, deriv_zet;
3431 deriv_csi.x = (field_data[k][j][i+1].x - field_data[k][j][i-1].x) * 0.5;
3432 deriv_csi.y = (field_data[k][j][i+1].y - field_data[k][j][i-1].y) * 0.5;
3433 deriv_csi.z = (field_data[k][j][i+1].z - field_data[k][j][i-1].z) * 0.5;
3434
3435 deriv_eta.x = (field_data[k][j+1][i].x - field_data[k][j-1][i].x) * 0.5;
3436 deriv_eta.y = (field_data[k][j+1][i].y - field_data[k][j-1][i].y) * 0.5;
3437 deriv_eta.z = (field_data[k][j+1][i].z - field_data[k][j-1][i].z) * 0.5;
3438
3439 deriv_zet.x = (field_data[k+1][j][i].x - field_data[k-1][j][i].x) * 0.5;
3440 deriv_zet.y = (field_data[k+1][j][i].y - field_data[k-1][j][i].y) * 0.5;
3441 deriv_zet.z = (field_data[k+1][j][i].z - field_data[k-1][j][i].z) * 0.5;
3442
3443 // 3. Transform derivatives to physical space
3444 TransformDerivativesToPhysical(jac[k][j][i], csi[k][j][i], eta[k][j][i], zet[k][j][i],
3445 deriv_csi, deriv_eta, deriv_zet,
3446 dudx, dvdx, dwdx);
3447
3448 // 4. Restore access to the PETSc data arrays
3449 ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, &csi); CHKERRQ(ierr);
3450 ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, &eta); CHKERRQ(ierr);
3451 ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, &zet); CHKERRQ(ierr);
3452 ierr = DMDAVecRestoreArrayRead(user->da, user->lAj, &jac); CHKERRQ(ierr);
3453
3454 PetscFunctionReturn(0);
3455}
3456
3457//================================================================================
3458//
3459// MEMORY CLEANUP FUNCTIONS
3460//
3461//================================================================================
3462
3463#undef __FUNCT__
3464#define __FUNCT__ "DestroyUserVectors"
3465/**
3466 * @brief Internal helper implementation: `DestroyUserVectors()`.
3467 * @details Local to this translation unit.
3468 */
3469PetscErrorCode DestroyUserVectors(UserCtx *user)
3470{
3471 PetscErrorCode ierr;
3472 PetscFunctionBeginUser;
3473
3474 // --- Group A: Primary Flow Fields (Always allocated at all levels) ---
3475 if (user->Ucont) { ierr = VecDestroy(&user->Ucont); CHKERRQ(ierr); }
3476 if (user->lUcont) { ierr = VecDestroy(&user->lUcont); CHKERRQ(ierr); }
3477 if (user->Ucat) { ierr = VecDestroy(&user->Ucat); CHKERRQ(ierr); }
3478 if (user->lUcat) { ierr = VecDestroy(&user->lUcat); CHKERRQ(ierr); }
3479 if (user->P) { ierr = VecDestroy(&user->P); CHKERRQ(ierr); }
3480 if (user->lP) { ierr = VecDestroy(&user->lP); CHKERRQ(ierr); }
3481 if (user->Nvert) { ierr = VecDestroy(&user->Nvert); CHKERRQ(ierr); }
3482 if (user->lNvert) { ierr = VecDestroy(&user->lNvert); CHKERRQ(ierr); }
3483
3484 // --- Group A2: Derived Flow Fields (Conditional) ---
3485 if(user->Diffusivity) {ierr = VecDestroy(&user->Diffusivity); CHKERRQ(ierr);}
3486 if(user->lDiffusivity){ierr = VecDestroy(&user->lDiffusivity); CHKERRQ(ierr);}
3487 if(user->DiffusivityGradient){ierr = VecDestroy(&user->DiffusivityGradient); CHKERRQ(ierr);}
3488 if(user->lDiffusivityGradient){ierr = VecDestroy(&user->lDiffusivityGradient); CHKERRQ(ierr);}
3489
3490 // --- Group B: Solver Work Vectors (All levels) ---
3491 if (user->Phi) { ierr = VecDestroy(&user->Phi); CHKERRQ(ierr); }
3492 if (user->lPhi) { ierr = VecDestroy(&user->lPhi); CHKERRQ(ierr); }
3493
3494 // --- Group C: Time-Stepping Vectors (Finest level only) ---
3495 if (user->Ucont_o) { ierr = VecDestroy(&user->Ucont_o); CHKERRQ(ierr); }
3496 if (user->Ucont_rm1) { ierr = VecDestroy(&user->Ucont_rm1); CHKERRQ(ierr); }
3497 if (user->Ucat_o) { ierr = VecDestroy(&user->Ucat_o); CHKERRQ(ierr); }
3498 if (user->P_o) { ierr = VecDestroy(&user->P_o); CHKERRQ(ierr); }
3499 if (user->Nvert_o) { ierr = VecDestroy(&user->Nvert_o); CHKERRQ(ierr); }
3500 if (user->lUcont_o) { ierr = VecDestroy(&user->lUcont_o); CHKERRQ(ierr); }
3501 if (user->lUcont_rm1) { ierr = VecDestroy(&user->lUcont_rm1); CHKERRQ(ierr); }
3502 if (user->lNvert_o) { ierr = VecDestroy(&user->lNvert_o); CHKERRQ(ierr); }
3503
3504 // --- Group D: Grid Metrics - Face Centered (All levels) ---
3505 if (user->Csi) { ierr = VecDestroy(&user->Csi); CHKERRQ(ierr); }
3506 if (user->Eta) { ierr = VecDestroy(&user->Eta); CHKERRQ(ierr); }
3507 if (user->Zet) { ierr = VecDestroy(&user->Zet); CHKERRQ(ierr); }
3508 if (user->Aj) { ierr = VecDestroy(&user->Aj); CHKERRQ(ierr); }
3509 if (user->lCsi) { ierr = VecDestroy(&user->lCsi); CHKERRQ(ierr); }
3510 if (user->lEta) { ierr = VecDestroy(&user->lEta); CHKERRQ(ierr); }
3511 if (user->lZet) { ierr = VecDestroy(&user->lZet); CHKERRQ(ierr); }
3512 if (user->lAj) { ierr = VecDestroy(&user->lAj); CHKERRQ(ierr); }
3513
3514 // --- Group E: Grid Metrics - Face Centered (All levels) ---
3515 if (user->ICsi) { ierr = VecDestroy(&user->ICsi); CHKERRQ(ierr); }
3516 if (user->IEta) { ierr = VecDestroy(&user->IEta); CHKERRQ(ierr); }
3517 if (user->IZet) { ierr = VecDestroy(&user->IZet); CHKERRQ(ierr); }
3518 if (user->JCsi) { ierr = VecDestroy(&user->JCsi); CHKERRQ(ierr); }
3519 if (user->JEta) { ierr = VecDestroy(&user->JEta); CHKERRQ(ierr); }
3520 if (user->JZet) { ierr = VecDestroy(&user->JZet); CHKERRQ(ierr); }
3521 if (user->KCsi) { ierr = VecDestroy(&user->KCsi); CHKERRQ(ierr); }
3522 if (user->KEta) { ierr = VecDestroy(&user->KEta); CHKERRQ(ierr); }
3523 if (user->KZet) { ierr = VecDestroy(&user->KZet); CHKERRQ(ierr); }
3524 if (user->IAj) { ierr = VecDestroy(&user->IAj); CHKERRQ(ierr); }
3525 if (user->JAj) { ierr = VecDestroy(&user->JAj); CHKERRQ(ierr); }
3526 if (user->KAj) { ierr = VecDestroy(&user->KAj); CHKERRQ(ierr); }
3527 if (user->lICsi) { ierr = VecDestroy(&user->lICsi); CHKERRQ(ierr); }
3528 if (user->lIEta) { ierr = VecDestroy(&user->lIEta); CHKERRQ(ierr); }
3529 if (user->lIZet) { ierr = VecDestroy(&user->lIZet); CHKERRQ(ierr); }
3530 if (user->lJCsi) { ierr = VecDestroy(&user->lJCsi); CHKERRQ(ierr); }
3531 if (user->lJEta) { ierr = VecDestroy(&user->lJEta); CHKERRQ(ierr); }
3532 if (user->lJZet) { ierr = VecDestroy(&user->lJZet); CHKERRQ(ierr); }
3533 if (user->lKCsi) { ierr = VecDestroy(&user->lKCsi); CHKERRQ(ierr); }
3534 if (user->lKEta) { ierr = VecDestroy(&user->lKEta); CHKERRQ(ierr); }
3535 if (user->lKZet) { ierr = VecDestroy(&user->lKZet); CHKERRQ(ierr); }
3536 if (user->lIAj) { ierr = VecDestroy(&user->lIAj); CHKERRQ(ierr); }
3537 if (user->lJAj) { ierr = VecDestroy(&user->lJAj); CHKERRQ(ierr); }
3538 if (user->lKAj) { ierr = VecDestroy(&user->lKAj); CHKERRQ(ierr); }
3539
3540 // --- Group F: Cell/Face Coordinates and Grid Spacing (All levels) ---
3541 if (user->Cent) { ierr = VecDestroy(&user->Cent); CHKERRQ(ierr); }
3542 if (user->lCent) { ierr = VecDestroy(&user->lCent); CHKERRQ(ierr); }
3543 if (user->GridSpace) { ierr = VecDestroy(&user->GridSpace); CHKERRQ(ierr); }
3544 if (user->lGridSpace) { ierr = VecDestroy(&user->lGridSpace); CHKERRQ(ierr); }
3545 if (user->Centx) { ierr = VecDestroy(&user->Centx); CHKERRQ(ierr); }
3546 if (user->Centy) { ierr = VecDestroy(&user->Centy); CHKERRQ(ierr); }
3547 if (user->Centz) { ierr = VecDestroy(&user->Centz); CHKERRQ(ierr); }
3548 if (user->lCentx) { ierr = VecDestroy(&user->lCentx); CHKERRQ(ierr); }
3549 if (user->lCenty) { ierr = VecDestroy(&user->lCenty); CHKERRQ(ierr); }
3550 if (user->lCentz) { ierr = VecDestroy(&user->lCentz); CHKERRQ(ierr); }
3551
3552 // --- Group G: Turbulence Model Vectors (Finest level, conditional on les/rans) ---
3553 if (user->Nu_t) { ierr = VecDestroy(&user->Nu_t); CHKERRQ(ierr); }
3554 if (user->lNu_t) { ierr = VecDestroy(&user->lNu_t); CHKERRQ(ierr); }
3555 if (user->CS) { ierr = VecDestroy(&user->CS); CHKERRQ(ierr); }
3556 if (user->lCs) { ierr = VecDestroy(&user->lCs); CHKERRQ(ierr); }
3557 if (user->lFriction_Velocity) { ierr = VecDestroy(&user->lFriction_Velocity); CHKERRQ(ierr); }
3558 if (user->K_Omega) { ierr = VecDestroy(&user->K_Omega); CHKERRQ(ierr); }
3559 if (user->lK_Omega) { ierr = VecDestroy(&user->lK_Omega); CHKERRQ(ierr); }
3560 if (user->K_Omega_o) { ierr = VecDestroy(&user->K_Omega_o); CHKERRQ(ierr); }
3561 if (user->lK_Omega_o) { ierr = VecDestroy(&user->lK_Omega_o); CHKERRQ(ierr); }
3562
3563 // --- Group H: Particle Vectors (Finest level, conditional on np > 0) ---
3564 if (user->ParticleCount) { ierr = VecDestroy(&user->ParticleCount); CHKERRQ(ierr); }
3565 if (user->lParticleCount) { ierr = VecDestroy(&user->lParticleCount); CHKERRQ(ierr); }
3566 if (user->Psi) { ierr = VecDestroy(&user->Psi); CHKERRQ(ierr); }
3567 if (user->lPsi) { ierr = VecDestroy(&user->lPsi); CHKERRQ(ierr); }
3568
3569 // --- Group I: Boundary Condition Vectors (All levels) ---
3570 if (user->Bcs.Ubcs) { ierr = VecDestroy(&user->Bcs.Ubcs); CHKERRQ(ierr); }
3571 if (user->Bcs.Uch) { ierr = VecDestroy(&user->Bcs.Uch); CHKERRQ(ierr); }
3572
3573 // --- Group J: Post-Processing Vectors (Finest level, postprocessor mode) ---
3574 if (user->P_nodal) { ierr = VecDestroy(&user->P_nodal); CHKERRQ(ierr); }
3575 if (user->Ucat_nodal) { ierr = VecDestroy(&user->Ucat_nodal); CHKERRQ(ierr); }
3576 if (user->Qcrit) { ierr = VecDestroy(&user->Qcrit); CHKERRQ(ierr); }
3577 if (user->Psi_nodal) { ierr = VecDestroy(&user->Psi_nodal); CHKERRQ(ierr); }
3578
3579 // --- Group K: Interpolation Vectors (Lazy allocation) ---
3580 if (user->fieldStatisticsStorage) {
3581 SimCtx *stats_ctx = user->simCtx;
3582 const PetscInt window_count = stats_ctx ? stats_ctx->fieldStatisticsWindowCount : 0;
3583 for (PetscInt w = 0; w < window_count; ++w) {
3584 ierr = PicurvWindowStorageDestroy(&user->fieldStatisticsStorage[w]); CHKERRQ(ierr);
3585 }
3586 ierr = PetscFree(user->fieldStatisticsStorage); CHKERRQ(ierr);
3587 user->fieldStatisticsStorage = NULL;
3588 }
3589 if (user->CellScalarAtCorner) { ierr = VecDestroy(&user->CellScalarAtCorner); CHKERRQ(ierr); }
3590 if (user->lCellScalarAtCorner) { ierr = VecDestroy(&user->lCellScalarAtCorner); CHKERRQ(ierr); }
3591 if (user->CellVectorAtCorner) { ierr = VecDestroy(&user->CellVectorAtCorner); CHKERRQ(ierr); }
3592 if (user->PostScalar) { ierr = VecDestroy(&user->PostScalar); CHKERRQ(ierr); }
3593 if (user->lPostScalar) { ierr = VecDestroy(&user->lPostScalar); CHKERRQ(ierr); }
3594 if (user->PostScalarNodal) { ierr = VecDestroy(&user->PostScalarNodal); CHKERRQ(ierr); }
3595 if (user->PostVector) { ierr = VecDestroy(&user->PostVector); CHKERRQ(ierr); }
3596 if (user->lPostVector) { ierr = VecDestroy(&user->lPostVector); CHKERRQ(ierr); }
3597 if (user->PostVectorNodal) { ierr = VecDestroy(&user->PostVectorNodal); CHKERRQ(ierr); }
3598 if (user->lCellVectorAtCorner) { ierr = VecDestroy(&user->lCellVectorAtCorner); CHKERRQ(ierr); }
3599
3600 // --- Group L: Implicit Solver Temporary Vectors (Destroyed after use, but check anyway) ---
3601 if (user->Rhs) { ierr = VecDestroy(&user->Rhs); CHKERRQ(ierr); }
3602 if (user->dUcont) { ierr = VecDestroy(&user->dUcont); CHKERRQ(ierr); }
3603 if (user->pUcont) { ierr = VecDestroy(&user->pUcont); CHKERRQ(ierr); }
3604
3605 // --- Group M: Poisson Solver Vectors (Destroyed after solve, but check anyway) ---
3606 if (user->B) { ierr = VecDestroy(&user->B); CHKERRQ(ierr); }
3607 if (user->R) { ierr = VecDestroy(&user->R); CHKERRQ(ierr); }
3608
3609 LOG_ALLOW(LOCAL, LOG_DEBUG, "All vectors destroyed for UserCtx.\n");
3610 PetscFunctionReturn(0);
3611}
3612#undef __FUNCT__
3613#define __FUNCT__ "DestroyUserContext"
3614/**
3615 * @brief Internal helper implementation: `DestroyUserContext()`.
3616 * @details Local to this translation unit.
3617 */
3618PetscErrorCode DestroyUserContext(UserCtx *user)
3619{
3620 PetscErrorCode ierr;
3621 PetscFunctionBeginUser;
3622
3623 if (!user) {
3624 LOG_ALLOW(LOCAL, LOG_WARNING, "DestroyUserContext called with NULL user pointer.\n");
3625 PetscFunctionReturn(0);
3626 }
3627
3628 LOG_ALLOW(LOCAL, LOG_INFO, "Destroying UserCtx at level %d...\n", user->thislevel);
3629
3630 // --- Step 1: Destroy Boundary Condition System ---
3631 // This handles all BC handlers and their private data.
3632 ierr = BoundarySystem_Destroy(user); CHKERRQ(ierr);
3633 LOG_ALLOW(LOCAL, LOG_DEBUG, " Boundary system destroyed.\n");
3634
3635 // --- Step 2: Destroy All Vectors ---
3636 // Handles ~74 Vec objects with proper NULL checking.
3637 ierr = DestroyUserVectors(user); CHKERRQ(ierr);
3638 LOG_ALLOW(LOCAL, LOG_DEBUG, " All vectors destroyed.\n");
3639
3640 // --- Step 3: Destroy Matrix and Solver Objects ---
3641 // Destroy pressure-Poisson matrices and solver.
3642 if (user->A) {
3643 ierr = MatDestroy(&user->A); CHKERRQ(ierr);
3644 LOG_ALLOW(LOCAL, LOG_DEBUG, " Matrix A destroyed.\n");
3645 }
3646 if (user->C) {
3647 ierr = MatDestroy(&user->C); CHKERRQ(ierr);
3648 LOG_ALLOW(LOCAL, LOG_DEBUG, " Matrix C destroyed.\n");
3649 }
3650 if (user->MR) {
3651 ierr = MatDestroy(&user->MR); CHKERRQ(ierr);
3652 LOG_ALLOW(LOCAL, LOG_DEBUG, " Matrix MR destroyed.\n");
3653 }
3654 if (user->MP) {
3655 ierr = MatDestroy(&user->MP); CHKERRQ(ierr);
3656 LOG_ALLOW(LOCAL, LOG_DEBUG, " Matrix MP destroyed.\n");
3657 }
3658 if (user->ksp) {
3659 ierr = KSPDestroy(&user->ksp); CHKERRQ(ierr);
3660 LOG_ALLOW(LOCAL, LOG_DEBUG, " KSP solver destroyed.\n");
3661 }
3662 if (user->nullsp) {
3663 ierr = MatNullSpaceDestroy(&user->nullsp); CHKERRQ(ierr);
3664 LOG_ALLOW(LOCAL, LOG_DEBUG, " MatNullSpace destroyed.\n");
3665 }
3666
3667 // --- Step 4: Destroy Application Ordering ---
3668 if (user->ao) {
3669 ierr = AODestroy(&user->ao); CHKERRQ(ierr);
3670 LOG_ALLOW(LOCAL, LOG_DEBUG, " AO destroyed.\n");
3671 }
3672
3673 // --- Step 5: Destroy DM Objects ---
3674 // Destroy in reverse order of dependency: post_swarm, swarm, fda6, fda2, fda, da
3675 if (user->post_swarm) {
3676 ierr = DMDestroy(&user->post_swarm); CHKERRQ(ierr);
3677 LOG_ALLOW(LOCAL, LOG_DEBUG, " post_swarm DM destroyed.\n");
3678 }
3679 if (user->swarm) {
3680 ierr = DMDestroy(&user->swarm); CHKERRQ(ierr);
3681 LOG_ALLOW(LOCAL, LOG_DEBUG, " swarm DM destroyed.\n");
3682 }
3683 if (user->fda6) {
3684 ierr = DMDestroy(&user->fda6); CHKERRQ(ierr);
3685 LOG_ALLOW(LOCAL, LOG_DEBUG, " fda6 DM destroyed.\n");
3686 }
3687 if (user->fda2) {
3688 ierr = DMDestroy(&user->fda2); CHKERRQ(ierr);
3689 LOG_ALLOW(LOCAL, LOG_DEBUG, " fda2 DM destroyed.\n");
3690 }
3691 if (user->da) {
3692 ierr = DMDestroy(&user->da); CHKERRQ(ierr);
3693 LOG_ALLOW(LOCAL, LOG_DEBUG, " da DM destroyed.\n");
3694 }
3695
3696 // --- Step 6: Free PetscMalloc'd Arrays ---
3697 // Free arrays allocated with PetscMalloc1
3698 if (user->RankCellInfoMap) {
3699 ierr = PetscFree(user->RankCellInfoMap); CHKERRQ(ierr);
3700 user->RankCellInfoMap = NULL;
3701 LOG_ALLOW(LOCAL, LOG_DEBUG, " RankCellInfoMap freed.\n");
3702 }
3703 if (user->KSKE) {
3704 ierr = PetscFree(user->KSKE); CHKERRQ(ierr);
3705 user->KSKE = NULL;
3706 LOG_ALLOW(LOCAL, LOG_DEBUG, " KSKE array freed.\n");
3707 }
3708
3709 LOG_ALLOW(LOCAL, LOG_INFO, "UserCtx at level %d fully destroyed.\n", user->thislevel);
3710 PetscFunctionReturn(0);
3711}
3712
3713#undef __FUNCT__
3714#define __FUNCT__ "FinalizeSimulation"
3715/**
3716 * @brief Implementation of \ref FinalizeSimulation().
3717 * @details Full API contract (arguments, ownership, side effects) is documented with
3718 * the header declaration in `include/setup.h`.
3719 * @see FinalizeSimulation()
3720 */
3721PetscErrorCode FinalizeSimulation(SimCtx *simCtx)
3722{
3723 PetscErrorCode ierr;
3724 PetscFunctionBeginUser;
3725
3726 if (!simCtx) {
3727 LOG_ALLOW(GLOBAL, LOG_WARNING, "FinalizeSimulation called with NULL SimCtx pointer.\n");
3728 PetscFunctionReturn(0);
3729 }
3730
3731 LOG_ALLOW(GLOBAL, LOG_INFO, "========================================\n");
3732 LOG_ALLOW(GLOBAL, LOG_INFO, "Beginning simulation memory cleanup...\n");
3733 LOG_ALLOW(GLOBAL, LOG_INFO, "========================================\n");
3734
3735 // ============================================================================
3736 // PHASE 1: DESTROY MULTIGRID HIERARCHY (All UserCtx structures)
3737 // ============================================================================
3738
3739 ierr = DestroySolutionConvergenceState(simCtx); CHKERRQ(ierr);
3740 ierr = DestroyFieldStatisticsConfig(simCtx); CHKERRQ(ierr);
3741
3742 if (simCtx->usermg.mgctx) {
3743 LOG_ALLOW(GLOBAL, LOG_INFO, "Destroying multigrid hierarchy (%d levels)...\n",
3744 simCtx->usermg.mglevels);
3745
3746 // Destroy each UserCtx from finest to coarsest (reverse order is safer)
3747 for (PetscInt level = simCtx->usermg.mglevels - 1; level >= 0; level--) {
3748 UserCtx *user = simCtx->usermg.mgctx[level].user;
3749 if (user) {
3750 LOG_ALLOW(LOCAL, LOG_INFO, " Destroying level %d of %d...\n",
3751 level, simCtx->usermg.mglevels - 1);
3752 ierr = DestroyUserContext(user); CHKERRQ(ierr);
3753
3754 // Free the UserCtx structure itself
3755 ierr = PetscFree(user); CHKERRQ(ierr);
3756 simCtx->usermg.mgctx[level].user = NULL;
3757 }
3758
3759 // Destroy the MGCtx-level packer DM
3760 if (simCtx->usermg.mgctx[level].packer) {
3761 ierr = DMDestroy(&simCtx->usermg.mgctx[level].packer); CHKERRQ(ierr);
3762 LOG_ALLOW(LOCAL, LOG_DEBUG, " MGCtx[%d].packer destroyed.\n", level);
3763 }
3764 }
3765
3766 // Free the MGCtx array itself
3767 ierr = PetscFree(simCtx->usermg.mgctx); CHKERRQ(ierr);
3768 simCtx->usermg.mgctx = NULL;
3769 LOG_ALLOW(GLOBAL, LOG_INFO, "All multigrid levels destroyed.\n");
3770 }
3771
3772 // ============================================================================
3773 // PHASE 2: DESTROY USERMG-LEVEL OBJECTS
3774 // ============================================================================
3775
3776 if (simCtx->usermg.packer) {
3777 ierr = DMDestroy(&simCtx->usermg.packer); CHKERRQ(ierr);
3778 LOG_ALLOW(LOCAL, LOG_DEBUG, "UserMG.packer DM destroyed.\n");
3779 }
3780
3781 if (simCtx->usermg.snespacker) {
3782 ierr = SNESDestroy(&simCtx->usermg.snespacker); CHKERRQ(ierr);
3783 LOG_ALLOW(LOCAL, LOG_DEBUG, "UserMG.snespacker SNES destroyed.\n");
3784 }
3785
3786 // ============================================================================
3787 // PHASE 3: DESTROY SIMCTX-LEVEL OBJECTS
3788 // ============================================================================
3789
3790 LOG_ALLOW(GLOBAL, LOG_INFO, "Destroying SimCtx-level objects...\n");
3791
3792 // --- PetscViewer for logging ---
3793 if (simCtx->logviewer) {
3794 ierr = PetscViewerDestroy(&simCtx->logviewer); CHKERRQ(ierr);
3795 LOG_ALLOW(LOCAL, LOG_DEBUG, " logviewer destroyed.\n");
3796 }
3797
3798 // --- Particle System DM ---
3799 if (simCtx->dm_swarm) {
3800 ierr = DMDestroy(&simCtx->dm_swarm); CHKERRQ(ierr);
3801 LOG_ALLOW(LOCAL, LOG_DEBUG, " dm_swarm destroyed.\n");
3802 }
3803
3804 // --- BoundingBox List (Array of BoundingBox structs) ---
3805 if (simCtx->bboxlist) {
3806 ierr = PetscFree(simCtx->bboxlist); CHKERRQ(ierr);
3807 simCtx->bboxlist = NULL;
3808 LOG_ALLOW(LOCAL, LOG_DEBUG, " bboxlist freed.\n");
3809 }
3810
3811 // --- Boundary Condition Files (Array of strings) ---
3812 if (simCtx->bcs_files) {
3813 for (PetscInt i = 0; i < simCtx->num_bcs_files; i++) {
3814 if (simCtx->bcs_files[i]) {
3815 ierr = PetscFree(simCtx->bcs_files[i]); CHKERRQ(ierr);
3816 }
3817 }
3818 ierr = PetscFree(simCtx->bcs_files); CHKERRQ(ierr);
3819 simCtx->bcs_files = NULL;
3820 LOG_ALLOW(LOCAL, LOG_DEBUG, " bcs_files array freed (%d files).\n", simCtx->num_bcs_files);
3821 }
3822
3823 // --- Brownian Motion RNG ---
3824 if (simCtx->BrownianMotionRNG) {
3825 ierr = PetscRandomDestroy(&simCtx->BrownianMotionRNG); CHKERRQ(ierr);
3826 LOG_ALLOW(LOCAL, LOG_DEBUG, " BrownianMotionRNG destroyed.\n");
3827 }
3828 // --- Post-Processing Parameters ---
3829 // pps is allocated with PetscNew and contains only static char arrays and basic types.
3830 // No internal dynamic allocations need to be freed.
3831 if (simCtx->pps) {
3832 ierr = PetscFree(simCtx->pps); CHKERRQ(ierr);
3833 simCtx->pps = NULL;
3834 LOG_ALLOW(LOCAL, LOG_DEBUG, " PostProcessParams freed.\n");
3835 }
3836
3837 // --- IBM/FSI Objects ---
3838 // Note: These are initialized to NULL and currently have no dedicated destroy functions.
3839 // If these modules are extended with cleanup routines, call them here.
3840 if (simCtx->ibm != NULL) {
3841 LOG_ALLOW(GLOBAL, LOG_WARNING, " WARNING: simCtx->ibm is non-NULL but no destroy function exists. Potential memory leak.\n");
3842 }
3843 if (simCtx->ibmv != NULL) {
3844 LOG_ALLOW(GLOBAL, LOG_WARNING, " WARNING: simCtx->ibmv is non-NULL but no destroy function exists. Potential memory leak.\n");
3845 }
3846 if (simCtx->fsi != NULL) {
3847 LOG_ALLOW(GLOBAL, LOG_WARNING, " WARNING: simCtx->fsi is non-NULL but no destroy function exists. Potential memory leak.\n");
3848 }
3849
3850 // --- Logging Allowed Functions (Array of strings) ---
3851 // Note: The logging system maintains its own copy via set_allowed_functions(),
3852 // so freeing simCtx->allowedFuncs will NOT affect LOG_ALLOW functionality.
3853 if (simCtx->allowedFuncs) {
3854 for (PetscInt i = 0; i < simCtx->nAllowed; i++) {
3855 if (simCtx->allowedFuncs[i]) {
3856 ierr = PetscFree(simCtx->allowedFuncs[i]); CHKERRQ(ierr);
3857 }
3858 }
3859 ierr = PetscFree(simCtx->allowedFuncs); CHKERRQ(ierr);
3860 simCtx->allowedFuncs = NULL;
3861 LOG_ALLOW(LOCAL, LOG_DEBUG, " allowedFuncs array freed (%d functions).\n", simCtx->nAllowed);
3862 }
3863
3864 // --- Profiling Critical Functions (Array of strings) ---
3865 if (simCtx->profilingSelectedFuncs) {
3866 for (PetscInt i = 0; i < simCtx->nProfilingSelectedFuncs; i++) {
3867 if (simCtx->profilingSelectedFuncs[i]) {
3868 ierr = PetscFree(simCtx->profilingSelectedFuncs[i]); CHKERRQ(ierr);
3869 }
3870 }
3871 ierr = PetscFree(simCtx->profilingSelectedFuncs); CHKERRQ(ierr);
3872 simCtx->profilingSelectedFuncs = NULL;
3873 LOG_ALLOW(LOCAL, LOG_DEBUG, " profilingSelectedFuncs array freed (%d functions).\n", simCtx->nProfilingSelectedFuncs);
3874 }
3875
3876 // ============================================================================
3877 // PHASE 4: FINAL SUMMARY
3878 // ============================================================================
3879
3880 LOG_ALLOW(GLOBAL, LOG_INFO, "========================================\n");
3881 LOG_ALLOW(GLOBAL, LOG_INFO, "Simulation cleanup completed successfully.\n");
3882 LOG_ALLOW(GLOBAL, LOG_INFO, "All PETSc objects have been destroyed.\n");
3883 LOG_ALLOW(GLOBAL, LOG_INFO, "========================================\n");
3884
3885 ierr = PetscFree(simCtx); CHKERRQ(ierr);
3886 PetscFunctionReturn(0);
3887}
BoundarySystem_Initialize
PetscErrorCode BoundarySystem_Initialize(UserCtx *user, const char *bcs_filename)
Initializes the entire boundary system.
Definition Boundaries.c:850
PropagateBoundaryConfigToCoarserLevels
PetscErrorCode PropagateBoundaryConfigToCoarserLevels(SimCtx *simCtx)
Propagates boundary condition configuration from finest to all coarser multigrid levels.
Definition Boundaries.c:947
BoundarySystem_Destroy
PetscErrorCode BoundarySystem_Destroy(UserCtx *user)
Cleans up and destroys all boundary system resources.
Definition Boundaries.c:1490
CalculateAllGridMetrics
PetscErrorCode CalculateAllGridMetrics(SimCtx *simCtx)
Orchestrates the calculation of all grid metrics.
Definition Metric.c:1940
FIELD_CAPABILITY_GHOST_UPDATE
@ FIELD_CAPABILITY_GHOST_UPDATE
Definition field_catalog.h:110
FIELD_SYNC_STANDARD
@ FIELD_SYNC_STANDARD
Definition field_catalog.h:90
FIELD_SYNC_K_FACE
@ FIELD_SYNC_K_FACE
Definition field_catalog.h:93
FIELD_SYNC_J_FACE
@ FIELD_SYNC_J_FACE
Definition field_catalog.h:92
FIELD_SYNC_COMPONENT_STAGGERED
@ FIELD_SYNC_COMPONENT_STAGGERED
Definition field_catalog.h:94
FIELD_SYNC_I_FACE
@ FIELD_SYNC_I_FACE
Definition field_catalog.h:91
FieldDescriptor::capabilities
unsigned int capabilities
Definition field_catalog.h:129
FieldView::descriptor
const FieldDescriptor * descriptor
Definition field_catalog.h:136
FieldGetView
PetscErrorCode FieldGetView(UserCtx *user, FieldId field_id, FieldView *view)
Resolve the existing DM and global/local vectors for one field.
Definition field_catalog.c:259
FieldView::local_vec
Vec local_vec
Definition field_catalog.h:139
FieldView::dm
DM dm
Definition field_catalog.h:137
FieldDescriptor::dof
PetscInt dof
Definition field_catalog.h:124
FieldView::global_vec
Vec global_vec
Definition field_catalog.h:138
FieldDescriptor::canonical_name
const char * canonical_name
Definition field_catalog.h:121
FieldDescriptor::sync_class
FieldSyncClass sync_class
Definition field_catalog.h:127
FieldId
FieldId
Compile-time identity for a catalogued Eulerian field.
Definition field_catalog.h:23
FIELD_ID_UCAT
@ FIELD_ID_UCAT
Definition field_catalog.h:26
FieldView
Non-owning runtime objects resolved for one field and UserCtx.
Definition field_catalog.h:135
DefineAllGridDimensions
PetscErrorCode DefineAllGridDimensions(SimCtx *simCtx)
Orchestrates the parsing and setting of grid dimensions for all blocks.
Definition grid.c:57
CalculateOutletProperties
PetscErrorCode CalculateOutletProperties(UserCtx *user)
Calculates the center and area of the primary OUTLET face.
Definition grid.c:1164
BroadcastAllBoundingBoxes
PetscErrorCode BroadcastAllBoundingBoxes(UserCtx *user, BoundingBox **bboxlist)
Broadcasts the bounding box information collected on rank 0 to all other ranks.
Definition grid.c:1061
ValidatePeriodicGeometry
PetscErrorCode ValidatePeriodicGeometry(UserCtx *user)
Validates that configured geometric periodic seams match by translation.
Definition grid.c:421
InitializeAllGridDMs
PetscErrorCode InitializeAllGridDMs(SimCtx *simCtx)
Orchestrates the creation of DMDA objects for every block and multigrid level.
Definition grid.c:276
AssignAllGridCoordinates
PetscErrorCode AssignAllGridCoordinates(SimCtx *simCtx)
Orchestrates the assignment of physical coordinates to all DMDA objects.
Definition grid.c:358
CalculateInletProperties
PetscErrorCode CalculateInletProperties(UserCtx *user)
Calculates the center and area of the primary INLET face.
Definition grid.c:1111
GatherAllBoundingBoxes
PetscErrorCode GatherAllBoundingBoxes(UserCtx *user, BoundingBox **allBBoxes)
Gathers local bounding boxes from all MPI processes to rank 0.
Definition grid.c:999
ParsePostProcessingSettings
PetscErrorCode ParsePostProcessingSettings(SimCtx *simCtx)
Initializes post-processing settings from a config file and command-line overrides.
Definition io.c:3076
ParseScalingInformation
PetscErrorCode ParseScalingInformation(SimCtx *simCtx)
Parses physical scaling parameters from command-line options.
Definition io.c:3242
VerifyPathExistence
PetscErrorCode VerifyPathExistence(const char *path, PetscBool is_dir, PetscBool is_optional, const char *description, PetscBool *exists)
A parallel-safe helper to verify the existence of a generic file or directory path.
Definition io.c:1128
set_allowed_functions
void set_allowed_functions(const char **functionList, int count)
Sets the global list of function names that are allowed to log.
Definition logging.c:155
is_function_allowed
PetscBool is_function_allowed(const char *functionName)
Checks if a given function is in the allow-list.
Definition logging.c:186
LOG_ALLOW_SYNC
#define LOG_ALLOW_SYNC(scope, level, fmt,...)
Synchronized logging macro that checks both the log level and whether the calling function is in the ...
Definition logging.h:253
LOCAL
#define LOCAL
Logging scope definitions for controlling message output.
Definition logging.h:45
GLOBAL
#define GLOBAL
Scope for global logging across all processes.
Definition logging.h:46
LOG_ALLOW
#define LOG_ALLOW(scope, level, fmt,...)
Logging macro that checks both the log level and whether the calling function is in the allowed-funct...
Definition logging.h:200
print_log_level
PetscErrorCode print_log_level(void)
Prints the current logging level to the console.
Definition logging.c:119
PROFILE_FUNCTION_END
#define PROFILE_FUNCTION_END
Marks the end of a profiled code block.
Definition logging.h:859
LOG
#define LOG(scope, level, fmt,...)
Logging macro for PETSc-based applications with scope control.
Definition logging.h:84
LoadAllowedFunctionsFromFile
PetscErrorCode LoadAllowedFunctionsFromFile(const char filename[], char ***funcsOut, PetscInt *nOut)
Load function names from a text file.
Definition logging.c:598
get_log_level
LogLevel get_log_level()
Retrieves the current logging level from the environment variable LOG_LEVEL.
Definition logging.c:87
ProfilingInitialize
PetscErrorCode ProfilingInitialize(SimCtx *simCtx)
Initializes the custom profiling system using configuration from SimCtx.
Definition logging.c:1928
LOG_ERROR
@ LOG_ERROR
Critical errors that may halt the program.
Definition logging.h:29
LOG_INFO
@ LOG_INFO
Informational messages about program execution.
Definition logging.h:31
LOG_WARNING
@ LOG_WARNING
Non-critical issues that warrant attention.
Definition logging.h:30
LOG_DEBUG
@ LOG_DEBUG
Detailed debugging information.
Definition logging.h:32
LOG_VERBOSE
@ LOG_VERBOSE
Extremely detailed logs, typically for development use only.
Definition logging.h:34
PROFILE_FUNCTION_BEGIN
#define PROFILE_FUNCTION_BEGIN
Marks the beginning of a profiled code block (typically a function).
Definition logging.h:850
ParticleInitializationToString
const char * ParticleInitializationToString(ParticleInitializationType ParticleInitialization)
Returns the canonical log token for a particle-initialization mode.
Definition logging.c:724
ComputeVectorFieldDerivatives
PetscErrorCode ComputeVectorFieldDerivatives(UserCtx *user, PetscInt i, PetscInt j, PetscInt k, Cmpnts ***field_data, Cmpnts *dudx, Cmpnts *dvdx, Cmpnts *dwdx)
Internal helper implementation: ComputeVectorFieldDerivatives().
Definition setup.c:3415
DestroyUserContext
PetscErrorCode DestroyUserContext(UserCtx *user)
Internal helper implementation: DestroyUserContext().
Definition setup.c:3618
GetOwnedCellRange
PetscErrorCode GetOwnedCellRange(const DMDALocalInfo *info_nodes, PetscInt dim, PetscInt *xs_cell_global_out, PetscInt *xm_cell_local_out)
Internal helper implementation: GetOwnedCellRange().
Definition setup.c:2285
SetupDomainRankInfo
PetscErrorCode SetupDomainRankInfo(SimCtx *simCtx)
Implementation of SetupDomainRankInfo().
Definition setup.c:2576
UniformCart2Contra
PetscErrorCode UniformCart2Contra(UserCtx *user, PetscReal u, PetscReal v, PetscReal w)
Populate contravariant fluxes from one uniform Cartesian velocity.
Definition setup.c:2853
InitializeRandomGenerators
PetscErrorCode InitializeRandomGenerators(UserCtx *user, PetscRandom *randx, PetscRandom *randy, PetscRandom *randz)
Implementation of InitializeRandomGenerators().
Definition setup.c:3200
Deallocate3DArrayVector
PetscErrorCode Deallocate3DArrayVector(Cmpnts ***array, PetscInt nz, PetscInt ny)
Implementation of Deallocate3DArrayVector().
Definition setup.c:2219
SetupGridAndSolvers
PetscErrorCode SetupGridAndSolvers(SimCtx *simCtx)
Implementation of SetupGridAndSolvers().
Definition setup.c:1364
SetupSolverParameters
static PetscErrorCode SetupSolverParameters(SimCtx *simCtx)
Configure solver tolerances, operators, and time-integration controls from the context.
Definition setup.c:1332
InitializeBrownianRNG
PetscErrorCode InitializeBrownianRNG(SimCtx *simCtx)
Internal helper implementation: InitializeBrownianRNG().
Definition setup.c:3285
Gidx
static PetscInt Gidx(PetscInt i, PetscInt j, PetscInt k, UserCtx *user)
Convert logical indices into the flattened global index used by setup helpers.
Definition setup.c:3008
SetupSimulationEnvironment
PetscErrorCode SetupSimulationEnvironment(SimCtx *simCtx)
Internal helper implementation: SetupSimulationEnvironment().
Definition setup.c:1063
CreateAndInitializeAllVectors
PetscErrorCode CreateAndInitializeAllVectors(SimCtx *simCtx)
Internal helper implementation: CreateAndInitializeAllVectors().
Definition setup.c:1400
ComputeAndStoreNeighborRanks
PetscErrorCode ComputeAndStoreNeighborRanks(UserCtx *user)
Internal helper implementation: ComputeAndStoreNeighborRanks().
Definition setup.c:2382
Contra2Cart
PetscErrorCode Contra2Cart(UserCtx *user)
Internal helper implementation: Contra2Cart().
Definition setup.c:2649
TransformScalarDerivativesToPhysical
void TransformScalarDerivativesToPhysical(PetscReal jacobian, Cmpnts csi_metrics, Cmpnts eta_metrics, Cmpnts zet_metrics, PetscReal dPhi_dcsi, PetscReal dPhi_deta, PetscReal dPhi_dzet, Cmpnts *gradPhi)
Implementation of TransformScalarDerivativesToPhysical().
Definition setup.c:3324
PetscMkdirRecursive
static PetscErrorCode PetscMkdirRecursive(const char *path)
Create a directory path recursively using PETSc-compatible error handling.
Definition setup.c:1011
DestroySolutionConvergenceState
PetscErrorCode DestroySolutionConvergenceState(SimCtx *simCtx)
Implementation of DestroySolutionConvergenceState().
Definition setup.c:101
Allocate3DArrayScalar
PetscErrorCode Allocate3DArrayScalar(PetscReal ****array, PetscInt nz, PetscInt ny, PetscInt nx)
Internal helper implementation: Allocate3DArrayScalar().
Definition setup.c:2091
CreateSimulationContext
PetscErrorCode CreateSimulationContext(int argc, char **argv, SimCtx **p_simCtx)
Implementation of CreateSimulationContext().
Definition setup.c:160
InitializeSolutionConvergenceState
PetscErrorCode InitializeSolutionConvergenceState(SimCtx *simCtx)
Implementation of InitializeSolutionConvergenceState().
Definition setup.c:49
SetDMDAProcLayout
PetscErrorCode SetDMDAProcLayout(DM dm, UserCtx *user)
Internal helper implementation: SetDMDAProcLayout().
Definition setup.c:2498
RepairPeriodicNormalFaceGhosts
static PetscErrorCode RepairPeriodicNormalFaceGhosts(UserCtx *user, DM dm, Vec local_vec, PetscInt dof, char face_direction, PetscBool component_staggered)
Repairs the adjacent normal ghost layer for periodic face-staggered data.
Definition setup.c:1668
InitializeLogicalSpaceRNGs
PetscErrorCode InitializeLogicalSpaceRNGs(PetscRandom *rand_logic_i, PetscRandom *rand_logic_j, PetscRandom *rand_logic_k)
Internal helper implementation: InitializeLogicalSpaceRNGs().
Definition setup.c:3241
ComputeScalarFieldDerivatives
PetscErrorCode ComputeScalarFieldDerivatives(UserCtx *user, PetscInt i, PetscInt j, PetscInt k, PetscReal ***field_data, Cmpnts *grad)
Internal helper implementation: ComputeScalarFieldDerivatives().
Definition setup.c:3372
ComputeDivergence
PetscErrorCode ComputeDivergence(UserCtx *user)
Implementation of ComputeDivergence().
Definition setup.c:3034
UpdateLocalGhosts
PetscErrorCode UpdateLocalGhosts(UserCtx *user, FieldId field_id)
Updates a catalogued field's local ghost representation.
Definition setup.c:1838
BinarySearchInt64
PetscErrorCode BinarySearchInt64(PetscInt n, const PetscInt64 arr[], PetscInt64 key, PetscBool *found)
Implementation of BinarySearchInt64().
Definition setup.c:2965
AllocateContextHierarchy
static PetscErrorCode AllocateContextHierarchy(SimCtx *simCtx)
Allocate the user-context objects required by every multigrid level.
Definition setup.c:1230
Cart2Contra
PetscErrorCode Cart2Contra(UserCtx *user)
Convert a spatially varying Cartesian velocity field to contravariant fluxes.
Definition setup.c:2784
DestroyUserVectors
PetscErrorCode DestroyUserVectors(UserCtx *user)
Internal helper implementation: DestroyUserVectors().
Definition setup.c:3469
Allocate3DArrayVector
PetscErrorCode Allocate3DArrayVector(Cmpnts ****array, PetscInt nz, PetscInt ny, PetscInt nx)
Implementation of Allocate3DArrayVector().
Definition setup.c:2169
SetupBoundaryConditions
PetscErrorCode SetupBoundaryConditions(SimCtx *simCtx)
Internal helper implementation: SetupBoundaryConditions().
Definition setup.c:2027
TransformDerivativesToPhysical
static void TransformDerivativesToPhysical(PetscReal jacobian, Cmpnts csi_metrics, Cmpnts eta_metrics, Cmpnts zet_metrics, Cmpnts deriv_csi, Cmpnts deriv_eta, Cmpnts deriv_zet, Cmpnts *dudx, Cmpnts *dvdx, Cmpnts *dwdx)
Transform contravariant vector derivatives into physical Cartesian derivatives.
Definition setup.c:3348
__FUNCT__
#define __FUNCT__
Definition setup.c:152
SetupDomainCellDecompositionMap
PetscErrorCode SetupDomainCellDecompositionMap(UserCtx *user)
Internal helper implementation: SetupDomainCellDecompositionMap().
Definition setup.c:2900
FinalizeSimulation
PetscErrorCode FinalizeSimulation(SimCtx *simCtx)
Implementation of FinalizeSimulation().
Definition setup.c:3721
Deallocate3DArrayScalar
PetscErrorCode Deallocate3DArrayScalar(PetscReal ***array, PetscInt nz, PetscInt ny)
Internal helper implementation: Deallocate3DArrayScalar().
Definition setup.c:2126
RuntimeWalltimeGuardParsePositiveSeconds
PetscBool RuntimeWalltimeGuardParsePositiveSeconds(const char *text, PetscReal *seconds_out)
Implementation of RuntimeWalltimeGuardParsePositiveSeconds().
Definition setup.c:20
statistics_accumulator.h
Per-window PETSc accumulator storage and pointwise application.
PicurvWindowStorageCreate
PetscErrorCode PicurvWindowStorageCreate(UserCtx *user, const PicurvWindowDefinition *definition, PicurvWindowStorage *storage)
Allocates the accumulator state one window owns on one block.
Definition statistics_accumulator.c:124
PicurvWindowStorageDestroy
PetscErrorCode PicurvWindowStorageDestroy(PicurvWindowStorage *storage)
Releases accumulator state previously created for one window.
Definition statistics_accumulator.c:204
statistics_config.h
Control ingress for the field-statistics pipeline.
DestroyFieldStatisticsConfig
PetscErrorCode DestroyFieldStatisticsConfig(SimCtx *simCtx)
Releases the window definitions resolved by ParseFieldStatisticsConfig().
Definition statistics_config.c:358
ParseFieldStatisticsConfig
PetscErrorCode ParseFieldStatisticsConfig(SimCtx *simCtx)
Resolves field-statistics configuration from the control file.
Definition statistics_config.c:291
PicurvWindow::definition
PicurvWindowDefinition definition
Definition statistics_window.h:76
FieldStatisticsIsActive
PetscBool FieldStatisticsIsActive(const struct SimCtx *simCtx)
Reports whether this run has live field-statistics state.
RankNeighbors::rank_zm
PetscMPIInt rank_zm
Definition variables.h:199
LESModelType
LESModelType
Identifies the six logical faces of a structured computational block.
Definition variables.h:520
NO_LES_MODEL
@ NO_LES_MODEL
Definition variables.h:521
SimCtx::icVelocityPhysical
PetscReal icVelocityPhysical
Definition variables.h:759
SimCtx::MHV
PetscInt MHV
Definition variables.h:732
UserCtx::lFriction_Velocity
Vec lFriction_Velocity
Definition variables.h:935
UserCtx::lDiffusivityGradient
Vec lDiffusivityGradient
Definition variables.h:943
UserCtx::isc
PetscInt isc
Definition variables.h:924
UserMG::packer
DM packer
Definition variables.h:582
SimCtx::turbine
PetscInt turbine
Definition variables.h:732
SimCtx::mom_nk_monitor_history
PetscBool mom_nk_monitor_history
Definition variables.h:752
SimCtx::fishcyl
PetscInt fishcyl
Definition variables.h:732
SimCtx::clark
PetscInt clark
Definition variables.h:822
UserCtx::lPostScalar
Vec lPostScalar
Definition variables.h:959
SimCtx::fieldStatisticsWindowCount
PetscInt fieldStatisticsWindowCount
Definition variables.h:770
PostProcessParams::statistics_output_prefix
char statistics_output_prefix[256]
basename for CSV output, e.g.
Definition variables.h:617
SimCtx::movefsi
PetscInt movefsi
Definition variables.h:726
UserCtx::lCent
Vec lCent
Definition variables.h:974
PERIODIC
@ PERIODIC
Definition variables.h:292
UserCtx::GridSpace
Vec GridSpace
Definition variables.h:974
SimCtx::continueMode
PetscBool continueMode
Definition variables.h:712
SimCtx::moveframe
PetscInt moveframe
Definition variables.h:727
UserCtx::P_nodal
Vec P_nodal
Definition variables.h:1001
UserCtx::JCsi
Vec JCsi
Definition variables.h:978
UserCtx::KAj
Vec KAj
Definition variables.h:979
SimCtx::TwoD
PetscInt TwoD
Definition variables.h:727
SimCtx::pseudo_periodic
PetscInt pseudo_periodic
Definition variables.h:791
MGCtx::user
UserCtx * user
Definition variables.h:571
SimCtx::fish_c
PetscInt fish_c
Definition variables.h:732
RankCellInfo::ys_cell
PetscInt ys_cell
Definition variables.h:204
SimCtx::dgf_z
PetscInt dgf_z
Definition variables.h:728
UserCtx::JEta
Vec JEta
Definition variables.h:978
SimCtx::poisson_tol
PetscReal poisson_tol
Definition variables.h:741
UserCtx::Zet
Vec Zet
Definition variables.h:974
UserCtx::Rhs
Vec Rhs
Definition variables.h:947
SimCtx::profilingFinalSummary
PetscBool profilingFinalSummary
Definition variables.h:868
PostProcessParams::particle_output_prefix
char particle_output_prefix[256]
Definition variables.h:612
RankCellInfo::xs_cell
PetscInt xs_cell
Definition variables.h:204
SimCtx::schmidt_number
PetscReal schmidt_number
Definition variables.h:787
SimCtx::rank
PetscMPIInt rank
Definition variables.h:698
UserCtx::mglevels
PetscInt mglevels
Definition variables.h:988
SimCtx::profilingTimestepFile
char profilingTimestepFile[PETSC_MAX_PATH_LEN]
Definition variables.h:867
SimCtx::fish
PetscInt fish
Definition variables.h:732
SimCtx::LV
PetscInt LV
Definition variables.h:732
SimCtx::angle
PetscReal angle
Definition variables.h:782
SimCtx::Turbulent_schmidt_number
PetscReal Turbulent_schmidt_number
Definition variables.h:787
UserCtx::boundary_faces
BoundaryFaceConfig boundary_faces[6]
Definition variables.h:931
RankNeighbors::rank_yp
PetscMPIInt rank_yp
Definition variables.h:198
SearchMetricsState::searchLocatedCount
PetscInt64 searchLocatedCount
Definition variables.h:241
SimCtx::thin
PetscInt thin
Definition variables.h:727
UserCtx::nullsp
MatNullSpace nullsp
Definition variables.h:965
SimCtx::statisticsConsoleOutputFreq
PetscInt statisticsConsoleOutputFreq
Definition variables.h:772
SimCtx::grid1d
PetscInt grid1d
Definition variables.h:790
SimCtx::block_number
PetscInt block_number
Definition variables.h:790
UserCtx::lIEta
Vec lIEta
Definition variables.h:977
UserCtx::KSKE
PetscInt * KSKE
Definition variables.h:966
SimCtx::mom_rtol
PetscReal mom_rtol
Definition variables.h:738
SearchMetricsState::searchLostCount
PetscInt64 searchLostCount
Definition variables.h:242
SimCtx::da_procs_z
PetscInt da_procs_z
Definition variables.h:796
SimCtx::blkpbc
PetscInt blkpbc
Definition variables.h:791
SimCtx::sediment
PetscInt sediment
Definition variables.h:726
SimCtx::targetVolumetricFlux
PetscReal targetVolumetricFlux
Definition variables.h:807
UserCtx::solutionConvergencePeriodicPRef
Vec * solutionConvergencePeriodicPRef
Definition variables.h:949
SimCtx::walltimeGuardActive
PetscBool walltimeGuardActive
Definition variables.h:870
UserMG::snespacker
SNES snespacker
Definition variables.h:583
UserCtx::lIZet
Vec lIZet
Definition variables.h:977
UserCtx::user_f
UserCtx * user_f
Definition variables.h:989
SimCtx::mom_last_lambda_max
PetscReal mom_last_lambda_max
Definition variables.h:751
SimCtx::channelz
PetscInt channelz
Definition variables.h:733
UserCtx::lNvert
Vec lNvert
Definition variables.h:939
UserCtx::Phi
Vec Phi
Definition variables.h:939
SimCtx::walltimeGuardWarmupTotalSeconds
PetscReal walltimeGuardWarmupTotalSeconds
Definition variables.h:878
SimCtx::forceScalingFactor
PetscReal forceScalingFactor
Definition variables.h:801
SimCtx::pseudo_cfl_reduction_factor
PetscReal pseudo_cfl_reduction_factor
Definition variables.h:745
SimCtx::initialConditionMode
InitialConditionMode initialConditionMode
Definition variables.h:754
UserCtx::simCtx
SimCtx * simCtx
Back-pointer to the master simulation context.
Definition variables.h:909
SimCtx::rans
PetscInt rans
Definition variables.h:821
ParticleInitializationType
ParticleInitializationType
Enumerator to identify the particle initialization strategy.
Definition variables.h:551
PARTICLE_INIT_SURFACE_RANDOM
@ PARTICLE_INIT_SURFACE_RANDOM
Random placement on the inlet face.
Definition variables.h:552
SimCtx::StartTime
PetscReal StartTime
Definition variables.h:709
SimCtx::dgf_az
PetscInt dgf_az
Definition variables.h:728
SimCtx::solutionConvergenceMeanSpeedHistory
PetscReal * solutionConvergenceMeanSpeedHistory
Definition variables.h:766
SimCtx::FluxOutSum
PetscReal FluxOutSum
Definition variables.h:799
RankNeighbors::rank_ym
PetscMPIInt rank_ym
Definition variables.h:198
SimCtx::walltimeGuardHasEWMA
PetscBool walltimeGuardHasEWMA
Definition variables.h:880
SimCtx::CMy_c
PetscReal CMy_c
Definition variables.h:783
UserCtx::K_Omega_o
Vec K_Omega_o
Definition variables.h:982
UserCtx::IZet
Vec IZet
Definition variables.h:977
SimCtx::flowDirection
FlowDirection flowDirection
Definition variables.h:758
RankNeighbors::rank_xp
PetscMPIInt rank_xp
Definition variables.h:197
UserCtx::Centz
Vec Centz
Definition variables.h:975
SimCtx::runtimeMemoryLogEnabled
PetscBool runtimeMemoryLogEnabled
Enable the rank-reduced runtime memory log.
Definition variables.h:883
PostProcessParams::output_prefix
char output_prefix[256]
Definition variables.h:609
SimCtx::bcs_files
char ** bcs_files
Definition variables.h:798
SimCtx::boundaryVelocityCorrection
PetscReal boundaryVelocityCorrection
Definition variables.h:814
SimCtx::max_angle
PetscReal max_angle
Definition variables.h:782
UserCtx::IEta
Vec IEta
Definition variables.h:977
SimCtx::min_pseudo_cfl
PetscReal min_pseudo_cfl
Definition variables.h:746
SearchMetricsState::boundaryClampCount
PetscInt64 boundaryClampCount
Definition variables.h:248
UserCtx::ksc
PetscInt ksc
Definition variables.h:924
SimCtx::particlesLostLastStep
PetscInt particlesLostLastStep
Definition variables.h:834
SimCtx::duplicate
PetscInt duplicate
Definition variables.h:849
SimCtx::tiout
PetscInt tiout
Definition variables.h:707
UserCtx::lZet
Vec lZet
Definition variables.h:974
UserCtx::assignedA
PetscBool assignedA
Definition variables.h:970
SimCtx::usermg
UserMG usermg
Definition variables.h:852
SimCtx::walltimeGuardMinSeconds
PetscReal walltimeGuardMinSeconds
Definition variables.h:873
SimCtx::allowedFile
char allowedFile[PETSC_MAX_PATH_LEN]
Definition variables.h:853
UserCtx::Csi
Vec Csi
Definition variables.h:974
UserCtx::solutionConvergencePeriodicUcatRef
Vec * solutionConvergencePeriodicUcatRef
Definition variables.h:948
SimCtx::da_procs_y
PetscInt da_procs_y
Definition variables.h:796
SearchMetricsState::traversalStepsSum
PetscInt64 traversalStepsSum
Definition variables.h:243
UserCtx::K_Omega
Vec K_Omega
Definition variables.h:982
SimCtx::mom_last_converged
PetscBool mom_last_converged
Definition variables.h:750
SimCtx::testfilter_1d
PetscInt testfilter_1d
Definition variables.h:824
SimCtx::psrc_x
PetscReal psrc_x
Definition variables.h:784
SimCtx::ren
PetscReal ren
Definition variables.h:744
SimCtx::Crotz
PetscReal Crotz
Definition variables.h:794
UserCtx::post_swarm
DM post_swarm
Definition variables.h:1000
UserCtx::lUcont_rm1
Vec lUcont_rm1
Definition variables.h:947
SimCtx::mixed
PetscInt mixed
Definition variables.h:822
RankCellInfo::zm_cell
PetscInt zm_cell
Definition variables.h:205
SimCtx::solutionConvergenceSamplesRecorded
PetscInt solutionConvergenceSamplesRecorded
Definition variables.h:765
UserCtx::lIAj
Vec lIAj
Definition variables.h:977
BoundingBox::max_coords
Cmpnts max_coords
Maximum x, y, z coordinates of the bounding box.
Definition variables.h:173
RankCellInfo::zs_cell
PetscInt zs_cell
Definition variables.h:204
SimCtx::drivenFluxTargetLatched
PetscBool drivenFluxTargetLatched
Definition variables.h:812
SimCtx::ibmv
IBMVNodes * ibmv
Definition variables.h:846
UserCtx::_this
PetscInt _this
Definition variables.h:924
UserCtx::lKEta
Vec lKEta
Definition variables.h:979
SearchMetricsState::searchPopulation
PetscInt64 searchPopulation
Definition variables.h:240
SimCtx::output_dir
char output_dir[PETSC_MAX_PATH_LEN]
Definition variables.h:717
SimCtx::solutionConvergenceEnabled
PetscBool solutionConvergenceEnabled
Definition variables.h:761
SimCtx::solutionConvergenceMeanKEHistory
PetscReal * solutionConvergenceMeanKEHistory
Definition variables.h:767
SimCtx::walltimeGuardLatestStepSeconds
PetscReal walltimeGuardLatestStepSeconds
Definition variables.h:882
SimCtx::dt
PetscReal dt
Definition variables.h:710
SimCtx::runtimeMemoryLogFile
char runtimeMemoryLogFile[PETSC_MAX_PATH_LEN]
File name written under log_dir.
Definition variables.h:884
SimCtx::runtimeMemoryLogStarted
PetscBool runtimeMemoryLogStarted
True after rank 0 writes the log header.
Definition variables.h:885
SimCtx::occupiedCellCount
PetscInt occupiedCellCount
Definition variables.h:838
SimCtx::StepsToRun
PetscInt StepsToRun
Definition variables.h:706
SimCtx::profilingTimestepMode
char profilingTimestepMode[32]
Definition variables.h:866
SimCtx::k_periodic
PetscInt k_periodic
Definition variables.h:791
SimCtx::inletprofile
PetscInt inletprofile
Definition variables.h:790
UserCtx::Ucat_nodal
Vec Ucat_nodal
Definition variables.h:1002
UserCtx::neighbors
RankNeighbors neighbors
Definition variables.h:923
SimCtx::bulkVelocityCorrection
PetscReal bulkVelocityCorrection
Definition variables.h:813
SimCtx::cdisy
PetscReal cdisy
Definition variables.h:744
SimCtx::mom_atol
PetscReal mom_atol
Definition variables.h:738
UserCtx::lPsi
Vec lPsi
Definition variables.h:997
SimCtx::rstart_fsi
PetscBool rstart_fsi
Definition variables.h:848
SearchMetricsState::currentSettlementPass
PetscInt currentSettlementPass
Definition variables.h:252
SimCtx::np
PetscInt np
Definition variables.h:827
UserCtx::PostScalarNodal
Vec PostScalarNodal
Definition variables.h:959
SimCtx::fieldStatisticsEnabled
PetscBool fieldStatisticsEnabled
Definition variables.h:769
UserCtx::jsc
PetscInt jsc
Definition variables.h:924
MGCtx::thislevel
PetscInt thislevel
Definition variables.h:572
SimCtx::no_pseudo_cfl_backtrack
PetscBool no_pseudo_cfl_backtrack
Definition variables.h:748
SimCtx::C_IEM
PetscReal C_IEM
Definition variables.h:842
UserCtx::DiffusivityGradient
Vec DiffusivityGradient
Definition variables.h:943
UserCtx::lJCsi
Vec lJCsi
Definition variables.h:978
SimCtx::ccc
PetscInt ccc
Definition variables.h:817
UserCtx::lCs
Vec lCs
Definition variables.h:982
SimCtx::ratio
PetscReal ratio
Definition variables.h:818
SimCtx::mg_idx
PetscInt mg_idx
Definition variables.h:739
UserCtx::Ucont
Vec Ucont
Definition variables.h:939
SimCtx::StartStep
PetscInt StartStep
Definition variables.h:705
SimCtx::mg_MAX_IT
PetscInt mg_MAX_IT
Definition variables.h:739
BoundingBox::min_coords
Cmpnts min_coords
Minimum x, y, z coordinates of the bounding box.
Definition variables.h:172
SimCtx::OnlySetup
PetscBool OnlySetup
Definition variables.h:711
SimCtx::rotatefsi
PetscInt rotatefsi
Definition variables.h:726
UserCtx::PostScalar
Vec PostScalar
Definition variables.h:959
BCS::Ubcs
Vec Ubcs
Physical Cartesian velocity at boundary faces. Full 3D array but only boundary-face entries are meani...
Definition variables.h:123
MOMENTUM_SOLVER_DUALTIME_PICARD_JAMESON_RK
@ MOMENTUM_SOLVER_DUALTIME_PICARD_JAMESON_RK
Definition variables.h:536
MOMENTUM_SOLVER_EXPLICIT_RK
@ MOMENTUM_SOLVER_EXPLICIT_RK
Definition variables.h:535
MOMENTUM_SOLVER_NEWTON_KRYLOV
@ MOMENTUM_SOLVER_NEWTON_KRYLOV
Definition variables.h:537
SimCtx::solutionConvergencePeriodSteps
PetscInt solutionConvergencePeriodSteps
Definition variables.h:763
SimCtx::cdisz
PetscReal cdisz
Definition variables.h:744
UserCtx::Qcrit
Vec Qcrit
Definition variables.h:1003
Cmpnts::x
PetscScalar x
Definition variables.h:103
UserCtx::JZet
Vec JZet
Definition variables.h:978
UserCtx::CellScalarAtCorner
Vec CellScalarAtCorner
Definition variables.h:952
SearchMetricsState::reSearchCount
PetscInt64 reSearchCount
Definition variables.h:244
SimCtx::dgf_x
PetscInt dgf_x
Definition variables.h:728
SimCtx::current_io_directory
char * current_io_directory
Definition variables.h:720
SimCtx::pizza
PetscInt pizza
Definition variables.h:732
SimCtx::MaxDiv
PetscReal MaxDiv
Definition variables.h:859
UserCtx::Centx
Vec Centx
Definition variables.h:975
UserCtx::Bcs
BCS Bcs
Definition variables.h:934
SimCtx::grid_file
char grid_file[PETSC_MAX_PATH_LEN]
Definition variables.h:795
UserCtx::lPhi
Vec lPhi
Definition variables.h:939
SimCtx::max_cs
PetscReal max_cs
Definition variables.h:823
UserCtx::lParticleCount
Vec lParticleCount
Definition variables.h:996
SimCtx::invicid
PetscInt invicid
Definition variables.h:727
SimCtx::allowedFuncs
char ** allowedFuncs
Definition variables.h:855
RankCellInfo::xm_cell
PetscInt xm_cell
Definition variables.h:205
PostProcessParams::statistics_pipeline
char statistics_pipeline[1024]
e.g.
Definition variables.h:616
UserCtx::lUcont_o
Vec lUcont_o
Definition variables.h:946
SearchMetricsState::bboxGuessFallbackCount
PetscInt64 bboxGuessFallbackCount
Definition variables.h:250
SimCtx::interpolationMethod
InterpolationMethod interpolationMethod
Definition variables.h:832
UserCtx::RankCellInfoMap
RankCellInfo * RankCellInfoMap
Definition variables.h:995
SimCtx::psrc_z
PetscReal psrc_z
Point source location for PARTICLE_INIT_POINT_SOURCE.
Definition variables.h:784
UserCtx::CellVectorAtCorner
Vec CellVectorAtCorner
Definition variables.h:953
SimCtx::mg_poItr
PetscInt mg_poItr
Definition variables.h:739
SimCtx::STRONG_COUPLING
PetscInt STRONG_COUPLING
Definition variables.h:742
RankCellInfo::ym_cell
PetscInt ym_cell
Definition variables.h:205
SimCtx::verificationScalar
VerificationScalarConfig verificationScalar
Definition variables.h:778
SimCtx::max_pseudo_cfl
PetscReal max_pseudo_cfl
Definition variables.h:746
UserCtx::Ucat_o
Vec Ucat_o
Definition variables.h:946
SimCtx::MaxDivx
PetscInt MaxDivx
Definition variables.h:860
UserCtx::user_c
UserCtx * user_c
Definition variables.h:989
SimCtx::poisson
PetscInt poisson
Definition variables.h:740
SimCtx::k_homo_filter
PetscInt k_homo_filter
Definition variables.h:824
SimCtx::profilingSelectedFuncsFile
char profilingSelectedFuncsFile[PETSC_MAX_PATH_LEN]
Definition variables.h:862
SimCtx::MaxDivy
PetscInt MaxDivy
Definition variables.h:860
SimCtx::NumberOfBodies
PetscInt NumberOfBodies
Definition variables.h:781
SimCtx::particleRestartMode
char particleRestartMode[16]
Definition variables.h:833
SimCtx::Ogrid
PetscInt Ogrid
Definition variables.h:790
SearchMetricsState::bboxGuessSuccessCount
PetscInt64 bboxGuessSuccessCount
Definition variables.h:249
SimCtx::MaxDivz
PetscInt MaxDivz
Definition variables.h:860
SimCtx::bboxlist
BoundingBox * bboxlist
Definition variables.h:830
SimCtx::j_homo_filter
PetscInt j_homo_filter
Definition variables.h:824
UserCtx::lKZet
Vec lKZet
Definition variables.h:979
SimCtx::fieldStatisticsWindows
struct PicurvWindow * fieldStatisticsWindows
Definition variables.h:771
UserCtx::Eta
Vec Eta
Definition variables.h:974
SimCtx::eel
PetscInt eel
Definition variables.h:732
SimCtx::log_dir
char log_dir[PETSC_MAX_PATH_LEN]
Definition variables.h:718
SimCtx::MaxDivFlatArg
PetscInt MaxDivFlatArg
Definition variables.h:860
UserCtx::lNu_t
Vec lNu_t
Definition variables.h:982
SimCtx::FluxInSum
PetscReal FluxInSum
Definition variables.h:799
RankNeighbors::rank_xm
PetscMPIInt rank_xm
Definition variables.h:197
UserCtx::Nu_t
Vec Nu_t
Definition variables.h:982
SimCtx::walltimeGuardCompletedSteps
PetscInt walltimeGuardCompletedSteps
Definition variables.h:877
SearchMetricsState::maxParticlePassDepth
PetscInt64 maxParticlePassDepth
Definition variables.h:251
PostProcessParams::source_dir
char source_dir[PETSC_MAX_PATH_LEN]
Definition variables.h:598
UserCtx::lCellScalarAtCorner
Vec lCellScalarAtCorner
Definition variables.h:952
UserCtx::lJEta
Vec lJEta
Definition variables.h:978
UserCtx::lCsi
Vec lCsi
Definition variables.h:974
SimCtx::CMz_c
PetscReal CMz_c
Definition variables.h:783
UserCtx::lGridSpace
Vec lGridSpace
Definition variables.h:974
SearchMetricsState::maxTraversalSteps
PetscInt64 maxTraversalSteps
Definition variables.h:245
SimCtx::generate_grid
PetscBool generate_grid
Definition variables.h:792
UserCtx::thislevel
PetscInt thislevel
Definition variables.h:988
SimCtx::AnalyticalUniformVelocity
Cmpnts AnalyticalUniformVelocity
Definition variables.h:760
SimCtx::eulerianSource
char eulerianSource[PETSC_MAX_PATH_LEN]
Definition variables.h:715
SimCtx::imp_stol
PetscReal imp_stol
Definition variables.h:738
SimCtx::nAllowed
PetscInt nAllowed
Definition variables.h:856
SimCtx::walltimeGuardEnabled
PetscBool walltimeGuardEnabled
Definition variables.h:869
SimCtx::checkpointGeometryHashReady
PetscBool checkpointGeometryHashReady
Definition variables.h:722
UserCtx::ICsi
Vec ICsi
Definition variables.h:977
SimCtx::wall_roughness_height
PetscReal wall_roughness_height
Definition variables.h:786
SimCtx::useProfilingSelectedFuncsCfg
PetscBool useProfilingSelectedFuncsCfg
Definition variables.h:863
SimCtx::walltimeGuardWarmupSteps
PetscInt walltimeGuardWarmupSteps
Definition variables.h:871
SimCtx::ParticleInitialization
ParticleInitializationType ParticleInitialization
Definition variables.h:831
SimCtx::mom_dt_jameson_residual_norm_noise_allowance_factor
PetscReal mom_dt_jameson_residual_norm_noise_allowance_factor
Definition variables.h:747
Cmpnts::z
PetscScalar z
Definition variables.h:103
UserCtx::pUcont
Vec pUcont
Definition variables.h:947
InterpolationMethod
InterpolationMethod
Selects the grid-to-particle interpolation method.
Definition variables.h:564
INTERP_TRILINEAR
@ INTERP_TRILINEAR
Definition variables.h:565
UserCtx::lK_Omega_o
Vec lK_Omega_o
Definition variables.h:982
UserCtx::lKCsi
Vec lKCsi
Definition variables.h:979
UserCtx::Ucat
Vec Ucat
Definition variables.h:939
UserCtx::ParticleCount
Vec ParticleCount
Definition variables.h:996
SimCtx::Const_CS
PetscReal Const_CS
Definition variables.h:823
UserCtx::lK_Omega
Vec lK_Omega
Definition variables.h:982
UserCtx::Ucont_o
Vec Ucont_o
Definition variables.h:946
SimCtx::i_homo_filter
PetscInt i_homo_filter
Definition variables.h:824
UserCtx::lCenty
Vec lCenty
Definition variables.h:976
SimCtx::wallfunction
PetscInt wallfunction
Definition variables.h:822
SimCtx::drivingForceStep
PetscInt drivingForceStep
Definition variables.h:806
SimCtx::rheology
PetscInt rheology
Definition variables.h:726
SimCtx::Flux_in
PetscReal Flux_in
Definition variables.h:782
UserCtx::PostVectorNodal
Vec PostVectorNodal
Definition variables.h:960
SimCtx::runtimeMemoryLogHasPrevious
PetscBool runtimeMemoryLogHasPrevious
True after the first process-memory sample.
Definition variables.h:886
UserMG::mglevels
PetscInt mglevels
Definition variables.h:578
SimCtx::profilingSelectedFuncs
char ** profilingSelectedFuncs
Definition variables.h:864
SimCtx::cdisx
PetscReal cdisx
Definition variables.h:744
SimCtx::dgf_ax
PetscInt dgf_ax
Definition variables.h:728
SimCtx::mglevels
PetscInt mglevels
Definition variables.h:739
MGCtx::packer
DM packer
Definition variables.h:573
SimCtx::num_bcs_files
PetscInt num_bcs_files
Definition variables.h:797
SimCtx::dm_swarm
DM dm_swarm
Definition variables.h:829
SimCtx::useCfg
PetscBool useCfg
Definition variables.h:854
SimCtx::psrc_y
PetscReal psrc_y
Definition variables.h:784
SimCtx::readFields
PetscBool readFields
Definition variables.h:828
SimCtx::solutionConvergenceWindowSteps
PetscInt solutionConvergenceWindowSteps
Definition variables.h:764
SimCtx::central
PetscInt central
Definition variables.h:742
SimCtx::Fluxsum
PetscReal Fluxsum
Definition variables.h:799
UserCtx::lJZet
Vec lJZet
Definition variables.h:978
UserCtx::Nvert_o
Vec Nvert_o
Definition variables.h:946
FlowDirection
FlowDirection
Primary flow direction for streamwise IC and Poiseuille modes.
Definition variables.h:272
FLOW_DIR_UNSET
@ FLOW_DIR_UNSET
Definition variables.h:279
SimCtx::pseudo_cfl_growth_factor
PetscReal pseudo_cfl_growth_factor
Definition variables.h:745
PostProcessParams::outputParticles
PetscBool outputParticles
Definition variables.h:604
SimCtx::Croty
PetscReal Croty
Definition variables.h:794
SimCtx::particlesLostCumulative
PetscInt particlesLostCumulative
Definition variables.h:835
SimCtx::nProfilingSelectedFuncs
PetscInt nProfilingSelectedFuncs
Definition variables.h:865
UserCtx::IAj
Vec IAj
Definition variables.h:977
UserCtx::lCellVectorAtCorner
Vec lCellVectorAtCorner
Definition variables.h:953
SimCtx::particlesMigratedLastStep
PetscInt particlesMigratedLastStep
Definition variables.h:837
SimCtx::initialConditionDirectory
char initialConditionDirectory[PETSC_MAX_PATH_LEN]
Definition variables.h:756
SimCtx::grid_rotation_angle
PetscReal grid_rotation_angle
Definition variables.h:793
SimCtx::dynamic_freq
PetscInt dynamic_freq
Definition variables.h:822
UserCtx::fieldStatisticsStorage
struct PicurvWindowStorage * fieldStatisticsStorage
Definition variables.h:962
UserCtx::Psi_nodal
Vec Psi_nodal
Definition variables.h:1004
SimCtx::AnalyticalSolutionType
char AnalyticalSolutionType[PETSC_MAX_PATH_LEN]
Definition variables.h:729
SimCtx::da_procs_x
PetscInt da_procs_x
Definition variables.h:796
UserCtx::JAj
Vec JAj
Definition variables.h:978
SimCtx::U_bc
PetscReal U_bc
Definition variables.h:816
UserCtx::lPostVector
Vec lPostVector
Definition variables.h:960
SimCtx::walltimeGuardWarmupAverageSeconds
PetscReal walltimeGuardWarmupAverageSeconds
Definition variables.h:879
UserCtx::KEta
Vec KEta
Definition variables.h:979
InitialConditionMode
InitialConditionMode
Selects the algorithm used to populate a fresh Eulerian velocity field.
Definition variables.h:151
IC_MODE_FILE
@ IC_MODE_FILE
Definition variables.h:156
IC_MODE_ZERO
@ IC_MODE_ZERO
Definition variables.h:152
SimCtx::particleConsoleOutputFreq
PetscInt particleConsoleOutputFreq
Definition variables.h:708
SimCtx::InitialConstantContra
Cmpnts InitialConstantContra
Definition variables.h:757
UserCtx::lCentx
Vec lCentx
Definition variables.h:976
RankNeighbors::rank_zp
PetscMPIInt rank_zp
Definition variables.h:199
UserCtx::Ucont_rm1
Vec Ucont_rm1
Definition variables.h:947
SimCtx::searchMetrics
SearchMetricsState searchMetrics
Definition variables.h:840
SimCtx::i_periodic
PetscInt i_periodic
Definition variables.h:791
SimCtx::mom_resid_rtol
PetscReal mom_resid_rtol
Definition variables.h:738
SimCtx::checkpointGeometrySHA256
char checkpointGeometrySHA256[65]
Definition variables.h:721
UserCtx::lUcont
Vec lUcont
Definition variables.h:939
SimCtx::runtimeMemoryLogPreviousProcessMB
PetscReal runtimeMemoryLogPreviousProcessMB
Previous local process memory sample in MB.
Definition variables.h:887
SimCtx::step
PetscInt step
Definition variables.h:703
UserCtx::Diffusivity
Vec Diffusivity
Definition variables.h:942
SimCtx::walltimeGuardEWMASeconds
PetscReal walltimeGuardEWMASeconds
Definition variables.h:881
SimCtx::AreaOutSum
PetscReal AreaOutSum
Definition variables.h:815
SimCtx::dgf_ay
PetscInt dgf_ay
Definition variables.h:728
SimCtx::mom_max_pseudo_steps
PetscInt mom_max_pseudo_steps
Definition variables.h:737
UserCtx::lAj
Vec lAj
Definition variables.h:974
SimCtx::BrownianMotionRNG
PetscRandom BrownianMotionRNG
Definition variables.h:841
UserCtx::lICsi
Vec lICsi
Definition variables.h:977
SimCtx::testfilter_ik
PetscInt testfilter_ik
Definition variables.h:824
UserCtx::info
DMDALocalInfo info
Definition variables.h:918
UserCtx::dUcont
Vec dUcont
Definition variables.h:947
SimCtx::hydro
PetscInt hydro
Definition variables.h:732
UserCtx::lUcat
Vec lUcat
Definition variables.h:939
SimCtx::pps
PostProcessParams * pps
Definition variables.h:890
SimCtx::migrationPassesLastStep
PetscInt migrationPassesLastStep
Definition variables.h:836
Cmpnts::y
PetscScalar y
Definition variables.h:103
InitialConditionField
InitialConditionField
Selects the authoritative velocity representation in a staged file IC.
Definition variables.h:160
IC_FIELD_UCONT
@ IC_FIELD_UCONT
Definition variables.h:162
IC_FIELD_UCAT
@ IC_FIELD_UCAT
Definition variables.h:161
SimCtx::size
PetscMPIInt size
Definition variables.h:699
EXEC_MODE_SOLVER
@ EXEC_MODE_SOLVER
Definition variables.h:668
EXEC_MODE_POSTPROCESSOR
@ EXEC_MODE_POSTPROCESSOR
Definition variables.h:669
SimCtx::_io_context_buffer
char _io_context_buffer[PETSC_MAX_PATH_LEN]
Definition variables.h:719
SimCtx::walltimeGuardLimitSeconds
PetscReal walltimeGuardLimitSeconds
Definition variables.h:876
UserCtx::lEta
Vec lEta
Definition variables.h:974
UserCtx::ksp
KSP ksp
Definition variables.h:965
UserCtx::PostVector
Vec PostVector
Definition variables.h:960
SimCtx::ps_ksp_pic_monitor_true_residual
PetscBool ps_ksp_pic_monitor_true_residual
Definition variables.h:753
UserCtx::KZet
Vec KZet
Definition variables.h:979
UserCtx::Cent
Vec Cent
Definition variables.h:974
SimCtx::walltimeGuardEstimatorAlpha
PetscReal walltimeGuardEstimatorAlpha
Definition variables.h:874
SimCtx::les
PetscInt les
Definition variables.h:821
UserCtx::Nvert
Vec Nvert
Definition variables.h:939
UserCtx::KCsi
Vec KCsi
Definition variables.h:979
UserMG::mgctx
MGCtx * mgctx
Definition variables.h:581
SOLUTION_CONVERGENCE_TRANSIENT
@ SOLUTION_CONVERGENCE_TRANSIENT
Definition variables.h:547
SOLUTION_CONVERGENCE_PERIODIC_DETERMINISTIC
@ SOLUTION_CONVERGENCE_PERIODIC_DETERMINISTIC
Definition variables.h:545
SOLUTION_CONVERGENCE_STATISTICAL_STEADY
@ SOLUTION_CONVERGENCE_STATISTICAL_STEADY
Definition variables.h:546
SOLUTION_CONVERGENCE_STEADY_DETERMINISTIC
@ SOLUTION_CONVERGENCE_STEADY_DETERMINISTIC
Definition variables.h:544
SimCtx::mg_preItr
PetscInt mg_preItr
Definition variables.h:739
UserCtx::lDiffusivity
Vec lDiffusivity
Definition variables.h:942
UserCtx::lNvert_o
Vec lNvert_o
Definition variables.h:946
SimCtx::mom_ratio_ema_alpha
PetscReal mom_ratio_ema_alpha
Definition variables.h:749
BoundaryFaceConfig::mathematical_type
BCType mathematical_type
Definition variables.h:368
UserCtx::Centy
Vec Centy
Definition variables.h:975
SimCtx::solutionConvergenceMode
SolutionConvergenceMode solutionConvergenceMode
Definition variables.h:762
SimCtx::logviewer
PetscViewer logviewer
Definition variables.h:713
UserCtx::multinullspace
PetscBool multinullspace
Definition variables.h:967
UserCtx::lCentz
Vec lCentz
Definition variables.h:976
SearchMetricsState::searchAttempts
PetscInt64 searchAttempts
Definition variables.h:239
SimCtx::initialConditionField
InitialConditionField initialConditionField
Definition variables.h:755
SimCtx::exec_mode
ExecutionMode exec_mode
Definition variables.h:714
SearchMetricsState::tieBreakCount
PetscInt64 tieBreakCount
Definition variables.h:247
SimCtx::mom_resid_atol
PetscReal mom_resid_atol
Definition variables.h:738
UserCtx::lJAj
Vec lJAj
Definition variables.h:978
UserCtx::bbox
BoundingBox bbox
Definition variables.h:922
SimCtx::restartHistoryAvailable
PetscBool restartHistoryAvailable
Definition variables.h:723
SimCtx::cop
PetscInt cop
Definition variables.h:732
SimCtx::ti
PetscReal ti
Definition variables.h:704
SimCtx::walltimeGuardMultiplier
PetscReal walltimeGuardMultiplier
Definition variables.h:872
SimCtx::Pipe
PetscInt Pipe
Definition variables.h:732
SimCtx::rotateframe
PetscInt rotateframe
Definition variables.h:727
SimCtx::ibm
IBMNodes * ibm
Definition variables.h:845
SimCtx::AreaInSum
PetscReal AreaInSum
Definition variables.h:815
SimCtx::mom_solver_type
MomentumSolverType mom_solver_type
Definition variables.h:736
SimCtx::summationRHS
PetscReal summationRHS
Definition variables.h:858
UserCtx::lKAj
Vec lKAj
Definition variables.h:979
SimCtx::immersed
PetscInt immersed
Definition variables.h:726
SearchMetricsState::maxTraversalFailCount
PetscInt64 maxTraversalFailCount
Definition variables.h:246
SimCtx::PostprocessingControlFile
char PostprocessingControlFile[PETSC_MAX_PATH_LEN]
Definition variables.h:889
SimCtx::restart_dir
char restart_dir[PETSC_MAX_PATH_LEN]
Definition variables.h:716
SimCtx::verificationDiffusivity
VerificationDiffusivityConfig verificationDiffusivity
Definition variables.h:777
SimCtx::blank
PetscInt blank
Definition variables.h:727
SimCtx::dgf_y
PetscInt dgf_y
Definition variables.h:728
SimCtx::walltimeGuardJobStartEpochSeconds
PetscReal walltimeGuardJobStartEpochSeconds
Definition variables.h:875
SimCtx::pseudo_cfl
PetscReal pseudo_cfl
Definition variables.h:744
SimCtx::LoggingFrequency
PetscInt LoggingFrequency
Definition variables.h:857
SimCtx::CMx_c
PetscReal CMx_c
Definition variables.h:783
SimCtx::drivingForceMagnitude
PetscReal drivingForceMagnitude
Definition variables.h:801
UserCtx::Psi
Vec Psi
Definition variables.h:997
SimCtx::particleLoadImbalance
PetscReal particleLoadImbalance
Definition variables.h:839
UserCtx::P_o
Vec P_o
Definition variables.h:946
SimCtx::fieldStatisticsContinue
PetscBool fieldStatisticsContinue
Definition variables.h:776
BCS::Uch
Vec Uch
Characteristic velocity for boundary conditions.
Definition variables.h:124
BC_FACE_NEG_X
@ BC_FACE_NEG_X
Definition variables.h:262
BC_FACE_POS_Z
@ BC_FACE_POS_Z
Definition variables.h:264
BC_FACE_POS_Y
@ BC_FACE_POS_Y
Definition variables.h:263
BC_FACE_NEG_Z
@ BC_FACE_NEG_Z
Definition variables.h:264
BC_FACE_POS_X
@ BC_FACE_POS_X
Definition variables.h:262
BC_FACE_NEG_Y
@ BC_FACE_NEG_Y
Definition variables.h:263
SimCtx::j_periodic
PetscInt j_periodic
Definition variables.h:791
SimCtx::wing
PetscInt wing
Definition variables.h:732
SimCtx::fsi
FSInfo * fsi
Definition variables.h:847
BoundingBox
Defines a 3D axis-aligned bounding box.
Definition variables.h:171
Cmpnts
A 3D point or vector with PetscScalar components.
Definition variables.h:102
MGCtx
Context for Multigrid operations.
Definition variables.h:570
PostProcessParams
Holds all configuration parameters for a post-processing run.
Definition variables.h:596
RankCellInfo
A lean struct to hold the global cell ownership range for a single MPI rank.
Definition variables.h:203
SimCtx
The master context for the entire simulation.
Definition variables.h:695
UserCtx
User-defined context containing data specific to a single computational grid level.
Definition variables.h:906
UserMG
User-level context for managing the entire multigrid hierarchy.
Definition variables.h:577
Generated by doxygen 1.9.8