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 "setup.h"
8
9#undef __FUNCT__
10#define __FUNCT__ "CreateSimulationContext"
11
12/**
13 * @brief Allocates and populates the master SimulationContext object.
14 *
15 * This function serves as the single, authoritative entry point for all
16 * simulation configuration. It merges the setup logic from both the legacy
17 * FSI/IBM solver and the modern particle solver into a unified, robust
18 * process.
19 *
20 * The function follows a strict sequence:
21 * 1. **Allocate Context & Set Defaults:** It first allocates the `SimulationContext`
22 * and populates every field with a sane, hardcoded default value. This
23 * ensures the simulation starts from a known, predictable state.
24 * 2. **Configure Logging System:** It configures the custom logging framework. It
25 * parses the `-func_config_file` option to load a list of function names
26 * allowed to produce log output. This configuration (the file path and the
27 * list of function names) is stored within the `SimulationContext` for
28 * later reference and cleanup.
29 * 3. **Parse All Options:** It performs a comprehensive pass of `PetscOptionsGet...`
30 * calls for every possible command-line flag, overriding the default
31 * values set in step 1.
32 * 4. **Log Summary:** After all options are parsed, it uses the now-active
33 * logging system to print a summary of the key simulation parameters.
34 *
35 * @param[in] argc Argument count passed from `main`.
36 * @param[in] argv Argument vector passed from `main`.
37 * @param[out] p_simCtx On success, this will point to the newly created and
38 * fully configured `SimulationContext` pointer. The caller
39 * is responsible for eventually destroying this object by
40 * calling `FinalizeSimulation()`.
41 *
42 * @return PetscErrorCode Returns 0 on success, or a non-zero PETSc error code on failure.
43 */
44PetscErrorCode CreateSimulationContext(int argc, char **argv, SimCtx **p_simCtx)
45{
46 PetscErrorCode ierr;
47 SimCtx *simCtx;
48 char control_filename[PETSC_MAX_PATH_LEN] = ""; // Temporary placeholder for control file name.
49 PetscBool control_flg; // Temporary placeholder for control file tag existence check flag.
50
51 PetscFunctionBeginUser;
52
53 PROFILE_FUNCTION_BEGIN;
54
55 // === 1. Allocate the Context Struct and Set ALL Defaults ==================
56 ierr = PetscNew(p_simCtx); CHKERRQ(ierr);
57 simCtx = *p_simCtx;
58
59 // --- Group 1: Parallelism & MPI ---
60 simCtx->rank = 0; simCtx->size = 1;
61
62 // --- Group 2: Simulation Control, Time, and I/O ---
63 simCtx->step = 0; simCtx->ti = 0.0; simCtx->StartStep = 0; simCtx->StepsToRun = 10;
64 simCtx->tiout = 10; simCtx->StartTime = 0.0; simCtx->dt = 0.001;
65 simCtx->OnlySetup = PETSC_FALSE;
66 simCtx->logviewer = NULL; simCtx->OutputFreq = simCtx->tiout;
67 strcpy(simCtx->eulerianSource,"solve");
68 strcpy(simCtx->restart_dir,"results");
69 strcpy(simCtx->output_dir,"results");
70 strcpy(simCtx->log_dir,"logs");
71 strcpy(simCtx->euler_subdir,"eulerian");
72 strcpy(simCtx->particle_subdir,"particles");
73 simCtx->_io_context_buffer[0] = '\0';
74 simCtx->current_io_directory = NULL;
75
76 // --- Group 3: High-Level Physics & Model Selection Flags ---
77 simCtx->immersed = 0; simCtx->movefsi = 0; simCtx->rotatefsi = 0;
78 simCtx->sediment = 0; simCtx->rheology = 0; simCtx->invicid = 0;
79 simCtx->TwoD = 0; simCtx->thin = 0; simCtx->moveframe = 0;
80 simCtx->rotateframe = 0; simCtx->blank = 0;
81 simCtx->dgf_x = 0; simCtx->dgf_y = 1; simCtx->dgf_z = 0;
82 simCtx->dgf_ax = 1; simCtx->dgf_ay = 0; simCtx->dgf_az = 0;
83 simCtx->st = 1.0;
84
85 // --- Group 4: Specific Simulation Case Flags --- (DEPRICATED)
86 simCtx->cop=0; simCtx->fish=0; simCtx->fish_c=0; simCtx->fishcyl=0;
87 simCtx->eel=0; simCtx->pizza=0; simCtx->turbine=0; simCtx->Pipe=0;
88 simCtx->wing=0; simCtx->hydro=0; simCtx->MHV=0; simCtx->LV=0;
89
90 // --- Group 5: Solver & Numerics Parameters ---
91 simCtx->implicit = 0; simCtx->implicit_type = 0; simCtx->imp_MAX_IT = 50;
92 simCtx->imp_atol = 1e-7; simCtx->imp_rtol = 1e-4; simCtx->imp_stol = 1.e-8;
93 simCtx->mglevels = 3; simCtx->mg_MAX_IT = 30; simCtx->mg_idx = 1;
94 simCtx->mg_preItr = 1; simCtx->mg_poItr = 1;
95 simCtx->poisson = 0; simCtx->poisson_tol = 5.e-9;
96 simCtx->STRONG_COUPLING = 0;simCtx->central=0;
97 simCtx->ren = 100.0; simCtx->cfl = 0.1; simCtx->vnn = 0.1;
98 simCtx->cdisx = 0.0; simCtx->cdisy = 0.0; simCtx->cdisz = 0.0;
99 simCtx->FieldInitialization = 0;
100 simCtx->InitialConstantContra.x = 0.0;
101 simCtx->InitialConstantContra.y = 0.0;
102 simCtx->InitialConstantContra.z = 0.0;
103
104 // --- Group 6: Physical & Geometric Parameters ---
105 simCtx->NumberOfBodies = 1; simCtx->Flux_in = 1.0; simCtx->angle = 0.0;
106 simCtx->max_angle = -54. * 3.1415926 / 180.;
107 simCtx->CMx_c=0.0; simCtx->CMy_c=0.0; simCtx->CMz_c=0.0;
108
109 // --- Group 7: Grid, Domain, and Boundary Condition Settings ---
110 simCtx->block_number = 1; simCtx->inletprofile = 1;
111 simCtx->grid1d = 0; simCtx->Ogrid = 0; simCtx->channelz = 0;
112 simCtx->i_periodic = 0; simCtx->j_periodic = 0; simCtx->k_periodic = 0;
113 simCtx->blkpbc = 10; simCtx->pseudo_periodic = 0;
114 strcpy(simCtx->grid_file, "config/grid.run");
115 simCtx->generate_grid = PETSC_FALSE;
116 simCtx->da_procs_x = PETSC_DECIDE;
117 simCtx->da_procs_y = PETSC_DECIDE;
118 simCtx->da_procs_z = PETSC_DECIDE;
119 simCtx->grid_rotation_angle = 0.0;
120 simCtx->Croty = 0.0; simCtx->Crotz = 0.0;
121 simCtx->num_bcs_files = 1;
122 ierr = PetscMalloc1(1, &simCtx->bcs_files); CHKERRQ(ierr);
123 ierr = PetscStrallocpy("config/bcs.run", &simCtx->bcs_files[0]); CHKERRQ(ierr);
124 simCtx->FluxInSum = 0.0; simCtx->FluxOutSum = 0.0; simCtx->Fluxsum = 0.0;
125 simCtx->AreaInSum = 0.0; simCtx->AreaOutSum = 0.0;
126 simCtx->U_bc = 0.0; simCtx->ccc = 0;
127 simCtx->ratio = 0.0;
128
129
130 // --- Group 8: Turbulence Modeling (LES/RANS) ---
131 simCtx->averaging = PETSC_FALSE; simCtx->les = 0; simCtx->rans = 0;
132 simCtx->wallfunction = 0; simCtx->mixed = 0; simCtx->clark = 0;
133 simCtx->dynamic_freq = 1; simCtx->max_cs = 0.5;
134 simCtx->Const_CS = 0.03;
135 simCtx->testfilter_ik = 0; simCtx->testfilter_1d = 0;
136 simCtx->i_homo_filter = 0; simCtx->j_homo_filter = 0; simCtx->k_homo_filter = 0;
137
138 // --- Group 9: Particle / DMSwarm Data & Settings ---
139 simCtx->np = 0; simCtx->readFields = PETSC_FALSE;
140 simCtx->dm_swarm = NULL; simCtx->bboxlist = NULL;
141 simCtx->ParticleInitialization = 0;
142 strcpy(simCtx->particleRestartMode,"load");
143
144 // --- Group 10: Immersed Boundary & FSI Data Object Pointers ---
145 simCtx->ibm = NULL; simCtx->ibmv = NULL; simCtx->fsi = NULL;
146 simCtx->rstart_fsi = PETSC_FALSE; simCtx->duplicate = 0;
147
148 // --- Group 11: Logging and Custom Configuration ---
149 strcpy(simCtx->allowedFile, "config/whitelist.run");
150 simCtx->useCfg = PETSC_FALSE;
151 simCtx->allowedFuncs = NULL;
152 simCtx->nAllowed = 0;
153 simCtx->LoggingFrequency = 10;
154 simCtx->summationRHS = 0.0;
155 simCtx->MaxDiv = 0.0;
156 simCtx->MaxDivFlatArg = 0; simCtx->MaxDivx = 0; simCtx->MaxDivy = 0; simCtx->MaxDivz = 0;
157 strcpy(simCtx->criticalFuncsFile, "config/profile.run");
158 simCtx->useCriticalFuncsCfg = PETSC_FALSE;
159 simCtx->criticalFuncs = NULL;
160 simCtx->nCriticalFuncs = 0;
161 // --- Group 11: Post-Processing Information ---
162 strcpy(simCtx->PostprocessingControlFile, "config/post.run");
163 ierr = PetscNew(&simCtx->pps); CHKERRQ(ierr);
164
165 // === 2. Get MPI Info and Handle Config File =============================
166 // -- Group 1: Parallelism & MPI Information
167 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &simCtx->rank); CHKERRQ(ierr);
168 ierr = MPI_Comm_size(PETSC_COMM_WORLD, &simCtx->size); CHKERRQ(ierr);
169
170 // First, check if the -control_file argument was provided by the user/script.
171 ierr = PetscOptionsGetString(NULL, NULL, "-control_file", control_filename, sizeof(control_filename), &control_flg); CHKERRQ(ierr);
172
173 // If the flag is NOT present or the filename is empty, abort with a helpful error.
174 if (!control_flg || strlen(control_filename) == 0) {
175 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
176 "\n\n*** MANDATORY ARGUMENT MISSING ***\n"
177 "The -control_file argument was not provided.\n"
178 "This program must be launched with a configuration file.\n"
179 "Example: mpiexec -n 4 ./picsolver -control_file /path/to/your/config.control\n"
180 "This is typically handled automatically by the 'pic-flow' script.\n");
181 }
182
183 // At this point, we have a valid filename. Attempt to load it.
184 LOG(GLOBAL, LOG_INFO, "Loading mandatory configuration from: %s\n", control_filename);
185 ierr = PetscOptionsInsertFile(PETSC_COMM_WORLD, NULL, control_filename, PETSC_FALSE);
186 if (ierr == PETSC_ERR_FILE_OPEN) {
187 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_FILE_OPEN, "The specified control file was not found or could not be opened: %s", control_filename);
188 }
189 CHKERRQ(ierr);
190
191 // === 3. A Configure Logging System ========================================
192 // This logic determines the logging configuration and STORES it in simCtx for
193 // later reference and cleanup.
194 ierr = PetscOptionsGetString(NULL, NULL, "-whitelist_config_file", simCtx->allowedFile, PETSC_MAX_PATH_LEN, &simCtx->useCfg); CHKERRQ(ierr);
195
196 if (simCtx->useCfg) {
197 ierr = LoadAllowedFunctionsFromFile(simCtx->allowedFile, &simCtx->allowedFuncs, &simCtx->nAllowed);
198 if (ierr) {
199 // Use direct PetscPrintf as logging system isn't fully active yet.
200 PetscPrintf(PETSC_COMM_SELF, "[%s] WARNING: Failed to load allowed functions from '%s'. Falling back to default list.\n", __func__, simCtx->allowedFile);
201 simCtx->useCfg = PETSC_FALSE; // Mark as failed.
202 ierr = 0; // Clear the error to allow fallback.
203 }
204 }
205 if (!simCtx->useCfg) {
206 // Fallback to default logging functions if no file was used or if loading failed.
207 simCtx->nAllowed = 2;
208 ierr = PetscMalloc1(simCtx->nAllowed, &simCtx->allowedFuncs); CHKERRQ(ierr);
209 ierr = PetscStrallocpy("main", &simCtx->allowedFuncs[0]); CHKERRQ(ierr);
210 ierr = PetscStrallocpy("CreateSimulationContext", &simCtx->allowedFuncs[1]); CHKERRQ(ierr);
211 }
212
213 // Activate the configuration by passing it to the logging module's setup function.
214 set_allowed_functions((const char**)simCtx->allowedFuncs, (size_t)simCtx->nAllowed);
215
216 // Now that the logger is configured, we can use it.
217 LOG_ALLOW_SYNC(LOCAL, LOG_INFO, "Context created. Initializing on rank %d of %d.\n", simCtx->rank, simCtx->size);
218 print_log_level(); // This will now correctly reflect the LOG_LEVEL environment variable.
219
220 // === 3.B Configure Profiling System ========================================
221 ierr = PetscOptionsGetString(NULL, NULL, "-profile_config_file", simCtx->criticalFuncsFile, PETSC_MAX_PATH_LEN, &simCtx->useCriticalFuncsCfg); CHKERRQ(ierr);
222 if (simCtx->useCriticalFuncsCfg) {
223 ierr = LoadAllowedFunctionsFromFile(simCtx->criticalFuncsFile, &simCtx->criticalFuncs, &simCtx->nCriticalFuncs);
224 if (ierr) {
225 PetscPrintf(PETSC_COMM_SELF, "[%s] WARNING: Failed to load critical profiling functions from '%s'. Falling back to default list.\n", __func__, simCtx->criticalFuncsFile);
226 simCtx->useCriticalFuncsCfg = PETSC_FALSE;
227 ierr = 0;
228 }
229 }
230 if (!simCtx->useCriticalFuncsCfg) {
231 // Fallback to a hardcoded default list if no file or loading failed
232 simCtx->nCriticalFuncs = 4;
233 ierr = PetscMalloc1(simCtx->nCriticalFuncs, &simCtx->criticalFuncs); CHKERRQ(ierr);
234 ierr = PetscStrallocpy("FlowSolver", &simCtx->criticalFuncs[0]); CHKERRQ(ierr);
235 ierr = PetscStrallocpy("AdvanceSimulation", &simCtx->criticalFuncs[1]); CHKERRQ(ierr);
236 ierr = PetscStrallocpy("LocateAllParticlesInGrid", &simCtx->criticalFuncs[2]); CHKERRQ(ierr);
237 ierr = PetscStrallocpy("InterpolateAllFieldsToSwarm", &simCtx->criticalFuncs[3]); CHKERRQ(ierr);
238 }
239
240 // Initialize the profiling system with the current updated simulation context.
241 ierr = ProfilingInitialize(simCtx); CHKERRQ(ierr);
242
243 // === 4. Parse All Command Line Options ==================================
244 LOG_ALLOW(GLOBAL, LOG_INFO, "Parsing command-line options...\n");
245
246 // --- Group 2
247 LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 2: Simulation Control,Time and I/O.\n");
248 // Read the physical time to start from.
249 // The default is already 0.0, so this will only be non-zero if the user provides it.
250 ierr = PetscOptionsGetInt(NULL, NULL, "-start_step", &simCtx->StartStep, NULL); CHKERRQ(ierr);
251 ierr = PetscOptionsGetInt(NULL,NULL, "-totalsteps", &simCtx->StepsToRun, NULL); CHKERRQ(ierr);
252 ierr = PetscOptionsGetBool(NULL, NULL, "-only_setup", &simCtx->OnlySetup, NULL); CHKERRQ(ierr);
253 ierr = PetscOptionsGetReal(NULL, NULL, "-dt", &simCtx->dt, NULL); CHKERRQ(ierr);
254 ierr = PetscOptionsGetInt(NULL, NULL, "-tio", &simCtx->tiout, NULL); CHKERRQ(ierr);
255 ierr = PetscOptionsGetString(NULL,NULL,"-euler_field_source",simCtx->eulerianSource,64,NULL);CHKERRQ(ierr);
256 ierr = PetscOptionsGetString(NULL,NULL,"-output_dir",&simCtx->output_dir,sizeof(simCtx->output_dir),NULL);CHKERRQ(ierr);
257 ierr = PetscOptionsGetString(NULL,NULL,"-restart_dir",&simCtx->restart_dir,sizeof(simCtx->restart_dir),NULL);CHKERRQ(ierr);
258 ierr = PetscOptionsGetString(NULL,NULL,"-log_dir",&simCtx->log_dir,sizeof(simCtx->log_dir),NULL);CHKERRQ(ierr);
259 ierr = PetscOptionsGetString(NULL,NULL,"-euler_subdir",&simCtx->euler_subdir,sizeof(simCtx->euler_subdir),NULL);CHKERRQ(ierr);
260 ierr = PetscOptionsGetString(NULL,NULL,"-particle_subdir",&simCtx->particle_subdir,sizeof(simCtx->particle_subdir),NULL);CHKERRQ(ierr);
261
262 simCtx->OutputFreq = simCtx->tiout; // backward compatibility related redundancy.
263 if(strcmp(simCtx->eulerianSource,"solve")!= 0 && strcmp(simCtx->eulerianSource,"load") != 0){
264 SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Invalid value for -euler_field_source. Must be 'load' or 'solve'. You provided '%s'.",simCtx->eulerianSource);
265 }
266
267 // --- Group 3
268 LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 3: High-Level Physics & Model Selection Flags\n");
269 ierr = PetscOptionsGetInt(NULL, NULL, "-imm", &simCtx->immersed, NULL); CHKERRQ(ierr);
270 ierr = PetscOptionsGetInt(NULL, NULL, "-fsi", &simCtx->movefsi, NULL); CHKERRQ(ierr);
271 ierr = PetscOptionsGetInt(NULL, NULL, "-rfsi", &simCtx->rotatefsi, NULL); CHKERRQ(ierr);
272 ierr = PetscOptionsGetInt(NULL, NULL, "-sediment", &simCtx->sediment, NULL); CHKERRQ(ierr);
273 ierr = PetscOptionsGetInt(NULL, NULL, "-rheology", &simCtx->rheology, NULL); CHKERRQ(ierr);
274 ierr = PetscOptionsGetInt(NULL, NULL, "-inv", &simCtx->invicid, NULL); CHKERRQ(ierr);
275 ierr = PetscOptionsGetInt(NULL, NULL, "-TwoD", &simCtx->TwoD, NULL); CHKERRQ(ierr);
276 ierr = PetscOptionsGetInt(NULL, NULL, "-thin", &simCtx->thin, NULL); CHKERRQ(ierr);
277 ierr = PetscOptionsGetInt(NULL, NULL, "-mframe", &simCtx->moveframe, NULL); CHKERRQ(ierr);
278 ierr = PetscOptionsGetInt(NULL, NULL, "-rframe", &simCtx->rotateframe, NULL); CHKERRQ(ierr);
279 ierr = PetscOptionsGetInt(NULL, NULL, "-blk", &simCtx->blank, NULL); CHKERRQ(ierr);
280 ierr = PetscOptionsGetInt(NULL, NULL, "-dgf_z", &simCtx->dgf_z, NULL); CHKERRQ(ierr);
281 ierr = PetscOptionsGetInt(NULL, NULL, "-dgf_y", &simCtx->dgf_y, NULL); CHKERRQ(ierr);
282 ierr = PetscOptionsGetInt(NULL, NULL, "-dgf_x", &simCtx->dgf_x, NULL); CHKERRQ(ierr);
283 ierr = PetscOptionsGetInt(NULL, NULL, "-dgf_az", &simCtx->dgf_az, NULL); CHKERRQ(ierr);
284 ierr = PetscOptionsGetInt(NULL, NULL, "-dgf_ay", &simCtx->dgf_ay, NULL); CHKERRQ(ierr);
285 ierr = PetscOptionsGetInt(NULL, NULL, "-dgf_ax", &simCtx->dgf_ax, NULL); CHKERRQ(ierr);
286
287 // --- Group 4
288 LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 4: Specific Simulation Case Flags \n");
289 ierr = PetscOptionsGetInt(NULL, NULL, "-cop", &simCtx->cop, NULL); CHKERRQ(ierr);
290 ierr = PetscOptionsGetInt(NULL, NULL, "-fish", &simCtx->fish, NULL); CHKERRQ(ierr);
291 ierr = PetscOptionsGetInt(NULL, NULL, "-pizza", &simCtx->pizza, NULL); CHKERRQ(ierr);
292 ierr = PetscOptionsGetInt(NULL, NULL, "-turbine", &simCtx->turbine, NULL); CHKERRQ(ierr);
293 ierr = PetscOptionsGetInt(NULL, NULL, "-fishcyl", &simCtx->fishcyl, NULL); CHKERRQ(ierr);
294 ierr = PetscOptionsGetInt(NULL, NULL, "-eel", &simCtx->eel, NULL); CHKERRQ(ierr);
295 ierr = PetscOptionsGetInt(NULL, NULL, "-cstart", &simCtx->fish_c, NULL); CHKERRQ(ierr);
296 ierr = PetscOptionsGetInt(NULL, NULL, "-wing", &simCtx->wing, NULL); CHKERRQ(ierr);
297 ierr = PetscOptionsGetInt(NULL, NULL, "-mhv", &simCtx->MHV, NULL); CHKERRQ(ierr);
298 ierr = PetscOptionsGetInt(NULL, NULL, "-hydro", &simCtx->hydro, NULL); CHKERRQ(ierr);
299 ierr = PetscOptionsGetInt(NULL, NULL, "-lv", &simCtx->LV, NULL); CHKERRQ(ierr);
300 ierr = PetscOptionsGetInt(NULL, NULL, "-Pipe", &simCtx->Pipe, NULL); CHKERRQ(ierr);
301
302 // --- Group 5
303 LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 5: Solver & Numerics Parameters \n");
304 ierr = PetscOptionsGetInt(NULL, NULL, "-imp", &simCtx->implicit, NULL); CHKERRQ(ierr);
305 ierr = PetscOptionsGetInt(NULL, NULL, "-imp_type", &simCtx->implicit_type, NULL); CHKERRQ(ierr);
306 ierr = PetscOptionsGetInt(NULL, NULL, "-imp_MAX_IT", &simCtx->imp_MAX_IT, NULL); CHKERRQ(ierr);
307 ierr = PetscOptionsGetReal(NULL, NULL, "-imp_atol", &simCtx->imp_atol, NULL); CHKERRQ(ierr);
308 ierr = PetscOptionsGetReal(NULL, NULL, "-imp_rtol", &simCtx->imp_rtol, NULL); CHKERRQ(ierr);
309 ierr = PetscOptionsGetReal(NULL, NULL, "-imp_stol", &simCtx->imp_stol, NULL); CHKERRQ(ierr);
310 ierr = PetscOptionsGetInt(NULL, NULL, "-central", &simCtx->central, NULL); CHKERRQ(ierr);
311
312 // --- Multigrid Options ---
313 ierr = PetscOptionsGetInt(NULL, NULL, "-mg_level", &simCtx->mglevels, NULL); CHKERRQ(ierr);
314 ierr = PetscOptionsGetInt(NULL, NULL, "-mg_max_it", &simCtx->mg_MAX_IT, NULL); CHKERRQ(ierr);
315 ierr = PetscOptionsGetInt(NULL, NULL, "-mg_idx", &simCtx->mg_idx, NULL); CHKERRQ(ierr);
316 ierr = PetscOptionsGetInt(NULL, NULL, "-mg_pre_it", &simCtx->mg_preItr, NULL); CHKERRQ(ierr);
317 ierr = PetscOptionsGetInt(NULL, NULL, "-mg_post_it", &simCtx->mg_poItr, NULL); CHKERRQ(ierr);
318
319 // --- Other Solver Options ---
320 ierr = PetscOptionsGetInt(NULL, NULL, "-poisson", &simCtx->poisson, NULL); CHKERRQ(ierr);
321 ierr = PetscOptionsGetReal(NULL, NULL, "-poisson_tol", &simCtx->poisson_tol, NULL); CHKERRQ(ierr);
322 ierr = PetscOptionsGetInt(NULL, NULL, "-str", &simCtx->STRONG_COUPLING, NULL); CHKERRQ(ierr);
323 ierr = PetscOptionsGetReal(NULL, NULL, "-ren", &simCtx->ren, NULL); CHKERRQ(ierr);
324 ierr = PetscOptionsGetReal(NULL, NULL, "-cfl", &simCtx->cfl, NULL); CHKERRQ(ierr);
325 ierr = PetscOptionsGetReal(NULL, NULL, "-vnn", &simCtx->vnn, NULL); CHKERRQ(ierr);
326 ierr = PetscOptionsGetInt(NULL, NULL, "-finit", &simCtx->FieldInitialization, NULL); CHKERRQ(ierr);
327 ierr = PetscOptionsGetReal(NULL, NULL, "-ucont_x", &simCtx->InitialConstantContra.x, NULL); CHKERRQ(ierr);
328 ierr = PetscOptionsGetReal(NULL, NULL, "-ucont_y", &simCtx->InitialConstantContra.y, NULL); CHKERRQ(ierr);
329 ierr = PetscOptionsGetReal(NULL, NULL, "-ucont_z", &simCtx->InitialConstantContra.z, NULL); CHKERRQ(ierr);
330 // NOTE: cdisx,cdisy,cdisz haven't been parsed, add if necessary.
331
332 // --- Group 6
333 LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 6: Physical & Geometric Parameters \n");
334 ierr = PetscOptionsGetInt(NULL, NULL, "-no_of_bodies", &simCtx->NumberOfBodies, NULL); CHKERRQ(ierr);
335 // NOTE: angle is not parsed in the original code, it set programmatically. We will follow that.
336 // NOTE: max_angle is calculated based on other flags (like MHV) in the legacy code.
337 // We will defer that logic to a later setup stage and not parse them directly.
338 // The Scaling Information is calculated here
339 ierr = ParseScalingInformation(simCtx); CHKERRQ(ierr);
340
341 // --- Group 7
342 LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 7: Grid, Domain, and Boundary Condition Settings \n");
343 ierr = PetscOptionsGetInt(NULL, NULL, "-nblk", &simCtx->block_number, NULL); CHKERRQ(ierr); // This is also a modern option
344 ierr = PetscOptionsGetInt(NULL, NULL, "-inlet", &simCtx->inletprofile, NULL); CHKERRQ(ierr);
345 ierr = PetscOptionsGetInt(NULL, NULL, "-Ogrid", &simCtx->Ogrid, NULL); CHKERRQ(ierr);
346 // NOTE: channelz was not parsed, likely set programmatically. We will omit its parsing call.
347 ierr = PetscOptionsGetInt(NULL, NULL, "-grid1d", &simCtx->grid1d, NULL); CHKERRQ(ierr);
348 ierr = PetscOptionsGetBool(NULL, NULL, "-grid", &simCtx->generate_grid, NULL); CHKERRQ(ierr);
349 ierr = PetscOptionsGetString(NULL, NULL, "-grid_file", simCtx->grid_file, PETSC_MAX_PATH_LEN, NULL); CHKERRQ(ierr);
350 ierr = PetscOptionsGetInt(NULL, NULL, "-da_processors_x", &simCtx->da_procs_x, NULL); CHKERRQ(ierr);
351 ierr = PetscOptionsGetInt(NULL, NULL, "-da_processors_y", &simCtx->da_procs_y, NULL); CHKERRQ(ierr);
352 ierr = PetscOptionsGetInt(NULL, NULL, "-da_processors_z", &simCtx->da_procs_z, NULL); CHKERRQ(ierr);
353 ierr = PetscOptionsGetInt(NULL, NULL, "-i_periodic", &simCtx->i_periodic, NULL); CHKERRQ(ierr);
354 ierr = PetscOptionsGetInt(NULL, NULL, "-j_periodic", &simCtx->j_periodic, NULL); CHKERRQ(ierr);
355 ierr = PetscOptionsGetInt(NULL, NULL, "-k_periodic", &simCtx->k_periodic, NULL); CHKERRQ(ierr);
356 ierr = PetscOptionsGetInt(NULL, NULL, "-pbc_domain", &simCtx->blkpbc, NULL); CHKERRQ(ierr);
357 // NOTE: pseudo_periodic was not parsed. We will omit its parsing call.
358 ierr = PetscOptionsGetReal(NULL, NULL, "-grid_rotation_angle", &simCtx->grid_rotation_angle, NULL); CHKERRQ(ierr);
359 ierr = PetscOptionsGetReal(NULL, NULL, "-Croty", &simCtx->Croty, NULL); CHKERRQ(ierr);
360 ierr = PetscOptionsGetReal(NULL, NULL, "-Crotz", &simCtx->Crotz, NULL); CHKERRQ(ierr);
361 PetscBool bcs_flg;
362 char file_list_str[PETSC_MAX_PATH_LEN * 10]; // Buffer for comma-separated list
363
364 ierr = PetscOptionsGetString(NULL, NULL, "-bcs_files", file_list_str, sizeof(file_list_str), &bcs_flg); CHKERRQ(ierr);
365 ierr = PetscOptionsGetReal(NULL, NULL, "-U_bc", &simCtx->U_bc, NULL); CHKERRQ(ierr);
366
367 if (bcs_flg) {
368 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Found -bcs_files option, overriding default.\n");
369
370 // A. Clean up the default memory we allocated in Phase 1.
371 ierr = PetscFree(simCtx->bcs_files[0]); CHKERRQ(ierr);
372 ierr = PetscFree(simCtx->bcs_files); CHKERRQ(ierr);
373 simCtx->num_bcs_files = 0;
374 simCtx->bcs_files = NULL;
375
376 // B. Parse the user-provided comma-separated list.
377 char *token;
378 char *str_copy;
379 ierr = PetscStrallocpy(file_list_str, &str_copy); CHKERRQ(ierr);
380
381 // First pass: count the number of files.
382 token = strtok(str_copy, ",");
383 while (token) {
384 simCtx->num_bcs_files++;
385 token = strtok(NULL, ",");
386 }
387 ierr = PetscFree(str_copy); CHKERRQ(ierr);
388
389 // Second pass: allocate memory and store the filenames.
390 ierr = PetscMalloc1(simCtx->num_bcs_files, &simCtx->bcs_files); CHKERRQ(ierr);
391 ierr = PetscStrallocpy(file_list_str, &str_copy); CHKERRQ(ierr);
392 token = strtok(str_copy, ",");
393 for (PetscInt i = 0; i < simCtx->num_bcs_files; i++) {
394 ierr = PetscStrallocpy(token, &simCtx->bcs_files[i]); CHKERRQ(ierr);
395 token = strtok(NULL, ",");
396 }
397 ierr = PetscFree(str_copy); CHKERRQ(ierr);
398 }
399
400
401 // --- Group 8
402 LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 8: Turbulence Modeling (LES/RANS) \n");
403 ierr = PetscOptionsGetInt(NULL, NULL, "-les", &simCtx->les, NULL); CHKERRQ(ierr);
404 ierr = PetscOptionsGetInt(NULL, NULL, "-rans", &simCtx->rans, NULL); CHKERRQ(ierr);
405 ierr = PetscOptionsGetInt(NULL, NULL, "-wallfunction", &simCtx->wallfunction, NULL); CHKERRQ(ierr);
406 ierr = PetscOptionsGetInt(NULL, NULL, "-mixed", &simCtx->mixed, NULL); CHKERRQ(ierr);
407 ierr = PetscOptionsGetInt(NULL, NULL, "-clark", &simCtx->clark, NULL); CHKERRQ(ierr);
408 ierr = PetscOptionsGetInt(NULL, NULL, "-dynamic_freq", &simCtx->dynamic_freq, NULL); CHKERRQ(ierr);
409 ierr = PetscOptionsGetReal(NULL, NULL, "-max_cs", &simCtx->max_cs, NULL); CHKERRQ(ierr);
410 ierr = PetscOptionsGetReal(NULL, NULL, "-const_cs", &simCtx->Const_CS, NULL); CHKERRQ(ierr);
411 ierr = PetscOptionsGetInt(NULL, NULL, "-testfilter_ik", &simCtx->testfilter_ik, NULL); CHKERRQ(ierr);
412 ierr = PetscOptionsGetInt(NULL, NULL, "-testfilter_1d", &simCtx->testfilter_1d, NULL); CHKERRQ(ierr);
413 ierr = PetscOptionsGetInt(NULL, NULL, "-i_homo_filter", &simCtx->i_homo_filter, NULL); CHKERRQ(ierr);
414 ierr = PetscOptionsGetInt(NULL, NULL, "-j_homo_filter", &simCtx->j_homo_filter, NULL); CHKERRQ(ierr);
415 ierr = PetscOptionsGetInt(NULL, NULL, "-k_homo_filter", &simCtx->k_homo_filter, NULL); CHKERRQ(ierr);
416 ierr = PetscOptionsGetBool(NULL, NULL, "-averaging", &simCtx->averaging, NULL); CHKERRQ(ierr);
417
418 // --- Group 9
419 LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 9: Particle / DMSwarm Data & Settings \n");
420 ierr = PetscOptionsGetInt(NULL, NULL, "-numParticles", &simCtx->np, NULL); CHKERRQ(ierr);
421 ierr = PetscOptionsGetBool(NULL, NULL, "-read_fields", &simCtx->readFields, NULL); CHKERRQ(ierr);
422 ierr = PetscOptionsGetInt(NULL, NULL, "-pinit", &simCtx->ParticleInitialization, NULL); CHKERRQ(ierr);
423 ierr = PetscOptionsGetString(NULL,NULL,"-particle_restart_mode",simCtx->particleRestartMode,sizeof(simCtx->particleRestartMode),NULL); CHKERRQ(ierr);
424 // Validation for Particle Restart Mode
425 if (strcmp(simCtx->particleRestartMode, "load") != 0 && strcmp(simCtx->particleRestartMode, "init") != 0) {
426 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Invalid value for -particle_restart_mode. Must be 'load' or 'init'. You provided '%s'.", simCtx->particleRestartMode);
427 }
428 // --- Group 10
429 LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 10: Immersed Boundary & FSI Data Object Pointers \n");
430 ierr = PetscOptionsGetBool(NULL, NULL, "-rs_fsi", &simCtx->rstart_fsi, NULL); CHKERRQ(ierr);
431 ierr = PetscOptionsGetInt(NULL, NULL, "-duplicate", &simCtx->duplicate, NULL); CHKERRQ(ierr);
432
433 // --- Group 11
434 LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 11: Top-Level Managers & Custom Configuration \n");
435 ierr = PetscOptionsGetInt(NULL, NULL, "-logfreq", &simCtx->LoggingFrequency, NULL); CHKERRQ(ierr);
436
437 if (simCtx->num_bcs_files != simCtx->block_number) {
438 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);
439 }
440
441 // --- Group 12
442 LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 12: Post-Processing Information.\n");
443 // This logic determines the Post Processing configuration and STORES it in simCtx for later reference and cleanup.
444 ierr = PetscOptionsGetString(NULL,NULL,"-postprocessing_config_file",simCtx->PostprocessingControlFile,PETSC_MAX_PATH_LEN,NULL); CHKERRQ(ierr);
445 ierr = PetscNew(&simCtx->pps); CHKERRQ(ierr);
446 ierr = ParsePostProcessingSettings(simCtx);
447
448 // === 5. Dependent Parameter Calculations ================================
449 // Some parameters depend on others, so we calculate them here.
450 simCtx->StartTime = (PetscReal)simCtx->StartStep*simCtx->dt;
451 simCtx->ti = simCtx->StartTime;
452 simCtx->step = simCtx->StartStep;
453
454 // === 5. Log Summary and Finalize Setup ==================================
455 LOG_ALLOW(GLOBAL, LOG_DEBUG, "-- Console Output Functions [Total : %d] : --\n", simCtx->nAllowed);
456 for (PetscInt i = 0; i < simCtx->nAllowed; ++i) {
457 LOG_ALLOW(GLOBAL, LOG_DEBUG, " [%2d] «%s»\n", i, simCtx->allowedFuncs[i]);
458 }
459
460 LOG_ALLOW(GLOBAL, LOG_INFO, "Configuration complete. Key parameters:\n");
461 LOG_ALLOW(GLOBAL, LOG_INFO, " - Run mode: %s\n", simCtx->OnlySetup ? "SETUP ONLY" : "Full Simulation");
462 LOG_ALLOW(GLOBAL, LOG_INFO, " - Time steps: %d (from %d to %d)\n", simCtx->StepsToRun, simCtx->StartStep, simCtx->StartStep + simCtx->StepsToRun);
463 LOG_ALLOW(GLOBAL, LOG_INFO, " - Time step size (dt): %g\n", simCtx->dt);
464 LOG_ALLOW(GLOBAL, LOG_INFO, " - Immersed Boundary: %s\n", simCtx->immersed ? "ENABLED" : "DISABLED");
465 LOG_ALLOW(GLOBAL, LOG_INFO, " - Particles: %d\n", simCtx->np);
466 if (simCtx->StartStep > 0 && simCtx->np > 0) {
467 LOG_ALLOW(GLOBAL, LOG_INFO, " - Particle Restart Mode: %s\n", simCtx->particleRestartMode);
468 }
469
470 // --- Initialize PETSc's internal performance logging stage ---
471 ierr = PetscLogDefaultBegin(); CHKERRQ(ierr); // REDUNDANT but safe.
472
473 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Finished CreateSimulationContext successfully on rank %d.\n", simCtx->rank);
474
475 PROFILE_FUNCTION_END;
476 PetscFunctionReturn(0);
477}
478
479#undef __FUNCT__
480#define __FUNCT__ "PetscMkdirRecursive"
481/**
482 * @brief Creates a directory path recursively, similar to `mkdir -p`.
483 *
484 * This function is designed to be called by Rank 0 only. It takes a full path,
485 * parses it, and creates each directory component in sequence.
486 *
487 * @param[in] path The full directory path to create.
488 *
489 * @return PetscErrorCode
490 */
491static PetscErrorCode PetscMkdirRecursive(const char *path)
492{
493 PetscErrorCode ierr;
494 char tmp_path[PETSC_MAX_PATH_LEN];
495 char *p = NULL;
496 size_t len;
497 PetscBool exists;
498
499 PetscFunctionBeginUser;
500
501 // Create a mutable copy of the path
502 len = strlen(path);
503 if (len >= sizeof(tmp_path)) {
504 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Path is too long to process: %s", path);
505 }
506 strcpy(tmp_path, path);
507
508 // If the path ends with a separator, remove it
509 if (tmp_path[len - 1] == '/') {
510 tmp_path[len - 1] = 0;
511 }
512
513 // Iterate through the path, creating each directory level
514 for (p = tmp_path + 1; *p; p++) {
515 if (*p == '/') {
516 *p = 0; // Temporarily terminate the string
517
518 // Check if this directory level exists
519 ierr = PetscTestDirectory(tmp_path, 'r', &exists); CHKERRQ(ierr);
520 if (!exists) {
521 ierr = PetscMkdir(tmp_path); CHKERRQ(ierr);
522 }
523
524 *p = '/'; // Restore the separator
525 }
526 }
527
528 // Create the final, full directory path
529 ierr = PetscTestDirectory(tmp_path, 'r', &exists); CHKERRQ(ierr);
530 if (!exists) {
531 ierr = PetscMkdir(tmp_path); CHKERRQ(ierr);
532 }
533
534 PetscFunctionReturn(0);
535}
536
537#undef __FUNCT__
538#define __FUNCT__ "SetupSimulationEnvironment"
539/**
540 * @brief Verifies and prepares the complete I/O environment for a simulation run.
541 *
542 * This function performs a comprehensive series of checks and setup actions to
543 * ensure a valid and clean environment. It is parallel-safe; all filesystem
544 * operations and checks are performed by Rank 0, with collective error handling.
545 *
546 * The function's responsibilities include:
547 * 1. **Checking Mandatory Inputs:** Verifies existence of grid and BCs files.
548 * 2. **Checking Optional Inputs:** Warns if optional config files (whitelist, profile) are missing.
549 * 3. **Validating Run Mode Paths:** Ensures `restart_dir` or post-processing source directories exist when needed.
550 * 4. **Preparing Log Directory:** Creates the log directory and cleans it of previous logs.
551 * 5. **Preparing Output Directories:** Creates the main output directory and its required subdirectories.
552 *
553 * @param[in] simCtx The fully configured SimulationContext object.
554 *
555 * @return PetscErrorCode Returns 0 on success, or a non-zero error code if a
556 * mandatory file/directory is missing or a critical operation fails.
557 */
558PetscErrorCode SetupSimulationEnvironment(SimCtx *simCtx)
559{
560 PetscErrorCode ierr;
561 PetscMPIInt rank;
562 PetscBool exists;
563
564 PetscFunctionBeginUser;
565 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
566
567 LOG_ALLOW(GLOBAL, LOG_INFO, "--- Setting up simulation environment ---\n");
568
569 /* =====================================================================
570 * Phase 1: Check for all required and optional INPUT files.
571 * ===================================================================== */
572 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Phase 1: Verifying input files...\n");
573
574 // --- Mandatory Inputs ---
575 if (!simCtx->generate_grid) {
576 ierr = VerifyPathExistence(simCtx->grid_file, PETSC_FALSE, PETSC_FALSE, "Grid file", &exists); CHKERRQ(ierr);
577 }
578 for (PetscInt i = 0; i < simCtx->num_bcs_files; i++) {
579 char desc[128];
580 ierr = PetscSNPrintf(desc, sizeof(desc), "BCS file #%d", i + 1); CHKERRQ(ierr);
581 ierr = VerifyPathExistence(simCtx->bcs_files[i], PETSC_FALSE, PETSC_FALSE, desc, &exists); CHKERRQ(ierr);
582 }
583
584 // --- Optional Inputs (these produce warnings if missing) ---
585 if (simCtx->useCfg) {
586 ierr = VerifyPathExistence(simCtx->allowedFile, PETSC_FALSE, PETSC_TRUE, "Whitelist config file", &exists); CHKERRQ(ierr);
587 }
588 if (simCtx->useCriticalFuncsCfg) {
589 ierr = VerifyPathExistence(simCtx->criticalFuncsFile, PETSC_FALSE, PETSC_TRUE, "Profiling config file", &exists); CHKERRQ(ierr);
590 }
591 if (simCtx->exec_mode == EXEC_MODE_POSTPROCESSOR) {
592 ierr = VerifyPathExistence(simCtx->PostprocessingControlFile, PETSC_FALSE, PETSC_TRUE, "Post-processing control file", &exists); CHKERRQ(ierr);
593 }
594
595
596 /* =====================================================================
597 * Phase 2: Validate directories specific to the execution mode.
598 * ===================================================================== */
599 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Phase 2: Verifying execution mode directories...\n");
600 // The data source directory must exist if we intend to load any data from it.
601 // This is true if:
602 // 1. We are restarting from a previous time step (StartStep > 0), which implies
603 // loading Eulerian fields and/or particle fields.
604 // 2. We are starting from t=0 but are explicitly told to load the initial
605 // Eulerian fields from a file (eulerianSource == "load").
606 if (simCtx->StartStep > 0 || strcmp(simCtx->eulerianSource,"load")== 0){ // If this is a restart run
607 ierr = VerifyPathExistence(simCtx->restart_dir, PETSC_TRUE, PETSC_FALSE, "Restart source directory", &exists); CHKERRQ(ierr);
608 }
609 if (simCtx->exec_mode == EXEC_MODE_POSTPROCESSOR) {
610 ierr = VerifyPathExistence(simCtx->pps->source_dir, PETSC_TRUE, PETSC_FALSE, "Post-processing source directory", &exists); CHKERRQ(ierr);
611 }
612
613 /* =====================================================================
614 * Phase 3: Create and prepare all OUTPUT directories.
615 * ===================================================================== */
616 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Phase 3: Preparing output directories...\n");
617
618 if (rank == 0){
619 if(simCtx->exec_mode == EXEC_MODE_SOLVER){
620 // --- Prepare Log Directory ---
621 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Creating/cleaning log directory: %s\n", simCtx->log_dir);
622 ierr = PetscRMTree(simCtx->log_dir); // Wipes the directory and its contents
623 if (ierr) { /* Ignore file-not-found error, but fail on others */
624 PetscError(PETSC_COMM_SELF, __LINE__, __FUNCT__, __FILE__, ierr, PETSC_ERROR_INITIAL, "Could not remove existing log directory '%s'. Check permissions.", simCtx->log_dir);
625 }
626 ierr = PetscMkdir(simCtx->log_dir); CHKERRQ(ierr);
627
628 // --- Prepare Output Directories ---
629 char path_buffer[PETSC_MAX_PATH_LEN];
630
631 // 1. Check/Create the main output directory
632 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Verifying main output directory: %s\n", simCtx->output_dir);
633 ierr = PetscTestDirectory(simCtx->output_dir, 'r', &exists); CHKERRQ(ierr);
634 if (!exists) {
635 LOG_ALLOW(GLOBAL, LOG_INFO, "Output directory not found. Creating: %s\n", simCtx->output_dir);
636 ierr = PetscMkdir(simCtx->output_dir); CHKERRQ(ierr);
637 }
638
639 // 2. Check/Create the Eulerian subdirectory
640 ierr = PetscSNPrintf(path_buffer, sizeof(path_buffer), "%s/%s", simCtx->output_dir, simCtx->euler_subdir); CHKERRQ(ierr);
641 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Verifying Eulerian subdirectory: %s\n", path_buffer);
642 ierr = PetscTestDirectory(path_buffer, 'r', &exists); CHKERRQ(ierr);
643 if (!exists) {
644 LOG_ALLOW(GLOBAL, LOG_INFO, "Eulerian subdirectory not found. Creating: %s\n", path_buffer);
645 ierr = PetscMkdir(path_buffer); CHKERRQ(ierr);
646 }
647
648 // 3. Check/Create the Particle subdirectory if needed
649 if (simCtx->np > 0) {
650 ierr = PetscSNPrintf(path_buffer, sizeof(path_buffer), "%s/%s", simCtx->output_dir, simCtx->particle_subdir); CHKERRQ(ierr);
651 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Verifying Particle subdirectory: %s\n", path_buffer);
652 ierr = PetscTestDirectory(path_buffer, 'r', &exists); CHKERRQ(ierr);
653 if (!exists) {
654 LOG_ALLOW(GLOBAL, LOG_INFO, "Particle subdirectory not found. Creating: %s\n", path_buffer);
655 ierr = PetscMkdir(path_buffer); CHKERRQ(ierr);
656 }
657 }
658 } else if(simCtx->exec_mode == EXEC_MODE_POSTPROCESSOR){
659 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Preparing post-processing output directories ...\n");
660
661 PostProcessParams *pps = simCtx->pps;
662 char path_buffer[PETSC_MAX_PATH_LEN];
663
664 const char *last_slash_euler = strrchr(pps->output_prefix, '/');
665 if(last_slash_euler){
666 size_t dir_len = last_slash_euler - pps->output_prefix;
667 if(dir_len > 0){
668 if(dir_len >= sizeof(path_buffer)) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Post-processing output prefix path is too long.");
669 strncpy(path_buffer, pps->output_prefix, dir_len);
670 path_buffer[dir_len] = '\0';
671
672 ierr = PetscTestDirectory(path_buffer, 'r', &exists); CHKERRQ(ierr);
673 if (!exists){
674 LOG_ALLOW(GLOBAL, LOG_INFO, "Creating post-processing Eulerian output directory: %s\n", path_buffer);
675 ierr = PetscMkdirRecursive(path_buffer); CHKERRQ(ierr);
676 }
677 }
678 }
679
680 // Particle output directory
681 if(pps->outputParticles){
682 const char *last_slash_particle = strrchr(pps->particle_output_prefix, '/');
683 if(last_slash_particle){
684 size_t dir_len = last_slash_particle - pps->particle_output_prefix;
685 if(dir_len > 0){
686 if(dir_len > sizeof(path_buffer)) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Post-processing particle output prefix path is too long.");
687 strncpy(path_buffer, pps->particle_output_prefix, dir_len);
688 path_buffer[dir_len] = '\0';
689
690 ierr = PetscTestDirectory(path_buffer, 'r', &exists); CHKERRQ(ierr);
691
692 if (!exists){
693 LOG_ALLOW(GLOBAL, LOG_INFO, "Creating post-processing Particle output directory: %s\n", path_buffer);
694 ierr = PetscMkdirRecursive(path_buffer); CHKERRQ(ierr);
695 }
696 }
697 }
698 }
699 }
700 }
701
702 // Synchronize all processes before proceeding
703 ierr = MPI_Barrier(PETSC_COMM_WORLD); CHKERRMPI(ierr);
704
705 LOG_ALLOW(GLOBAL, LOG_INFO, "--- Environment setup complete ---\n");
706
707 PetscFunctionReturn(0);
708}
709
710#undef __FUNCT__
711#define __FUNCT__ "AllocateContextHeirarchy"
712/**
713 * @brief Allocates the memory for the UserMG and UserCtx hierarchy.
714 *
715 * This function performs the foundational memory allocation for the solver's
716 * data structures. It reads the number of multigrid levels and grid blocks
717 * from the master SimulationContext, then allocates the arrays for MGCtx and
718 * UserCtx.
719 *
720 * It performs three critical tasks:
721 * 1. Allocates the MGCtx array within the UserMG struct.
722 * 2. For each level, allocates the UserCtx array for all blocks.
723 * 3. Sets the `simCtx` back-pointer in every single UserCtx, linking it to the
724 * master configuration.
725 *
726 * @param simCtx The master SimulationContext, which contains configuration and
727 * will store the resulting allocated hierarchy in its `usermg` member.
728 * @return PetscErrorCode
729 */
730static PetscErrorCode AllocateContextHierarchy(SimCtx *simCtx)
731{
732 PetscErrorCode ierr;
733 UserMG *usermg = &simCtx->usermg;
734 MGCtx *mgctx;
735 PetscInt nblk = simCtx->block_number;
736 PetscBool found;
737 PetscFunctionBeginUser;
738 PROFILE_FUNCTION_BEGIN;
739
740 LOG_ALLOW(GLOBAL, LOG_INFO, "Allocating context hierarchy for %d levels and %d blocks...\n", simCtx->mglevels, nblk);
741
742 // Store the number of levels in the UserMG struct itself
743 usermg->mglevels = simCtx->mglevels;
744
745 // --- 1. Allocate the array of MGCtx structs ---
746 ierr = PetscMalloc(usermg->mglevels * sizeof(MGCtx), &usermg->mgctx); CHKERRQ(ierr);
747 mgctx = usermg->mgctx;
748 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Allocated MGCtx array of size %d.\n", simCtx->rank, usermg->mglevels);
749
750 // --- 2. Parse semi-coarsening options (logic from MG_Initial) ---
751 // These flags determine if a dimension is coarsened in the multigrid hierarchy.
752 PetscInt *isc, *jsc, *ksc;
753 ierr = PetscMalloc3(nblk, &isc, nblk, &jsc, nblk, &ksc); CHKERRQ(ierr);
754 // Set defaults to FALSE (full coarsening)
755 for (PetscInt i = 0; i < nblk; ++i) {
756 isc[i] = 0; jsc[i] = 0; ksc[i] = 0;
757 }
758
759// Use a temporary variable for the 'count' argument to the parsing function.
760 // This protects the original 'nblk' which is needed for the loop bounds.
761 PetscInt n_opts_found = nblk;
762 ierr = PetscOptionsGetIntArray(NULL, NULL, "-mg_i_semi", isc, &n_opts_found, &found); CHKERRQ(ierr);
763
764 n_opts_found = nblk; // Reset the temp variable before the next call
765 ierr = PetscOptionsGetIntArray(NULL, NULL, "-mg_j_semi", jsc, &n_opts_found, &found); CHKERRQ(ierr);
766
767 n_opts_found = nblk; // Reset the temp variable before the next call
768 ierr = PetscOptionsGetIntArray(NULL, NULL, "-mg_k_semi", ksc, &n_opts_found, &found); CHKERRQ(ierr);
769
770 // --- 3. Loop over levels and blocks to allocate UserCtx arrays ---
771 for (PetscInt level = 0; level < simCtx->mglevels; level++) {
772
773 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Setting up MG Level %d...\n", simCtx->rank, level);
774 // Allocate the array of UserCtx structs for this level
775 ierr = PetscMalloc(nblk * sizeof(UserCtx), &mgctx[level].user); CHKERRQ(ierr);
776 // It's good practice to zero out the memory to avoid uninitialized values
777 ierr = PetscMemzero(mgctx[level].user, nblk * sizeof(UserCtx)); CHKERRQ(ierr);
778 mgctx[level].thislevel = level;
779
780 for (PetscInt bi = 0; bi < nblk; bi++) {
781 UserCtx *currentUser = &mgctx[level].user[bi];
782 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Initializing UserCtx for Level %d, Block %d.\n", simCtx->rank, level, bi);
783
784 // --- CRITICAL STEP: Set the back-pointer to the master context ---
785 currentUser->simCtx = simCtx;
786
787 // Initialize other per-context values
788 currentUser->thislevel = level;
789 currentUser->_this = bi; //
790 currentUser->mglevels = usermg->mglevels;
791
792 // Assign semi-coarsening flags
793 currentUser->isc = isc[bi];
794 currentUser->jsc = jsc[bi];
795 currentUser->ksc = ksc[bi];
796
797 // Link to finer/coarser contexts for multigrid operations
798 if (level > 0) {
799 currentUser->user_c = &mgctx[level-1].user[bi];
800 LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG, "Rank %d: -> Linked to coarser context (user_c).\n", simCtx->rank);
801 }
802 if (level < usermg->mglevels - 1) {
803 currentUser->user_f = &mgctx[level+1].user[bi];
804 LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG, "Rank %d: -> Linked to finer context (user_f).\n", simCtx->rank);
805 }
806 }
807 }
808
809 // Log a summary of the parsed flags on each rank.
810 if (get_log_level() >= LOG_DEBUG && nblk > 0) {
811 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Final semi-coarsening configuration view:\n", simCtx->rank);
812 for (PetscInt bi = 0; bi < nblk; ++bi) {
813 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]);
814 }
815 }
816
817 // Clean up temporary arrays
818 ierr = PetscFree3(isc, jsc, ksc); CHKERRQ(ierr);
819
820 LOG_ALLOW(GLOBAL, LOG_INFO, "Context hierarchy allocation complete.\n");
821 PROFILE_FUNCTION_END;
822 PetscFunctionReturn(0);
823}
824
825#undef __FUNCT__
826#define __FUNCT__ "SetupSolverParameters"
827static PetscErrorCode SetupSolverParameters(SimCtx *simCtx){
828
829 PetscFunctionBeginUser;
830 PROFILE_FUNCTION_BEGIN;
831
832 LOG_ALLOW(GLOBAL,LOG_INFO, " -- Setting up solver parameters -- .\n");
833
834 UserMG *usermg = &simCtx->usermg;
835 MGCtx *mgctx = usermg->mgctx;
836 PetscInt nblk = simCtx->block_number;
837
838 for (PetscInt level = usermg->mglevels-1; level >=0; level--) {
839 for (PetscInt bi = 0; bi < nblk; bi++) {
840 UserCtx *user = &mgctx[level].user[bi];
841 LOG_ALLOW_SYNC(LOCAL, LOG_DEBUG, "Rank %d: Setting up parameters for level %d, block %d\n", simCtx->rank, level, bi);
842
843 user->assignedA = PETSC_FALSE;
844 user->multinullspace = PETSC_FALSE;
845 }
846 }
847 PROFILE_FUNCTION_END;
848 PetscFunctionReturn(0);
849}
850
851#undef __FUNCT__
852#define __FUNCT__ "SetupGridAndSolvers"
853/**
854 * @brief The main orchestrator for setting up all grid-related components.
855 * (Implementation of the orchestrator function itself)
856 */
857PetscErrorCode SetupGridAndSolvers(SimCtx *simCtx)
858{
859 PetscErrorCode ierr;
860 PetscFunctionBeginUser;
861
862 PROFILE_FUNCTION_BEGIN;
863
864 LOG_ALLOW(GLOBAL, LOG_INFO, "--- Starting Grid and Solvers Setup ---\n");
865
866 // Phase 1: Allocate the UserMG and UserCtx hierarchy
867 ierr = AllocateContextHierarchy(simCtx); CHKERRQ(ierr);
868
869 // [ The next phases will be added here ]
870 ierr = DefineAllGridDimensions(simCtx); CHKERRQ(ierr);
871 ierr = InitializeAllGridDMs(simCtx); CHKERRQ(ierr);
872 ierr = AssignAllGridCoordinates(simCtx);
873 ierr = CreateAndInitializeAllVectors(simCtx); CHKERRQ(ierr);
874 ierr = SetupSolverParameters(simCtx); CHKERRQ(ierr);
875 ierr = CalculateAllGridMetrics(simCtx); CHKERRQ(ierr);
876
877 LOG_ALLOW(GLOBAL, LOG_INFO, "--- Grid and Solvers Setup Complete ---\n");
878
879 PROFILE_FUNCTION_END;
880 PetscFunctionReturn(0);
881}
882
883
884#undef __FUNCT__
885#define __FUNCT__ "CreateAndInitializeAllVectors"
886/**
887 * @brief Creates and initializes all PETSc Vec objects for all fields.
888 *
889 * This function iterates through every UserCtx in the multigrid and multi-block
890 * hierarchy. For each context, it creates the comprehensive set of global and
891 * local PETSc Vecs required by the flow solver (e.g., Ucont, P, Nvert, metrics,
892 * turbulence fields, etc.). Each vector is initialized to zero.
893 *
894 * @param simCtx The master SimCtx, containing the configured UserCtx hierarchy.
895 * @return PetscErrorCode
896 */
897PetscErrorCode CreateAndInitializeAllVectors(SimCtx *simCtx)
898{
899 PetscErrorCode ierr;
900 UserMG *usermg = &simCtx->usermg;
901 MGCtx *mgctx = usermg->mgctx;
902 PetscInt nblk = simCtx->block_number;
903
904 PetscFunctionBeginUser;
905
906 PROFILE_FUNCTION_BEGIN;
907
908 LOG_ALLOW(GLOBAL, LOG_INFO, "Creating and initializing all simulation vectors...\n");
909
910 for (PetscInt level = usermg->mglevels-1; level >=0; level--) {
911 for (PetscInt bi = 0; bi < nblk; bi++) {
912 UserCtx *user = &mgctx[level].user[bi];
913
914 if(!user->da || !user->fda) {
915 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "DMs not properly initialized in UserCtx before vector creation.");
916 }
917
918 LOG_ALLOW_SYNC(LOCAL, LOG_DEBUG, "Rank %d: Creating vectors for level %d, block %d\n", simCtx->rank, level, bi);
919
920 // --- Group A: Primary Flow Fields (Global and Local) ---
921 // These are the core solution variables.
922 ierr = DMCreateGlobalVector(user->fda, &user->Ucont); CHKERRQ(ierr); ierr = VecSet(user->Ucont, 0.0); CHKERRQ(ierr);
923 ierr = DMCreateGlobalVector(user->fda, &user->Ucat); CHKERRQ(ierr); ierr = VecSet(user->Ucat, 0.0); CHKERRQ(ierr);
924 ierr = DMCreateGlobalVector(user->da, &user->P); CHKERRQ(ierr); ierr = VecSet(user->P, 0.0); CHKERRQ(ierr);
925 ierr = DMCreateGlobalVector(user->da, &user->Nvert); CHKERRQ(ierr); ierr = VecSet(user->Nvert, 0.0); CHKERRQ(ierr);
926
927 ierr = DMCreateLocalVector(user->fda, &user->lUcont); CHKERRQ(ierr); ierr = VecSet(user->lUcont, 0.0); CHKERRQ(ierr);
928 ierr = DMCreateLocalVector(user->fda, &user->lUcat); CHKERRQ(ierr); ierr = VecSet(user->lUcat, 0.0); CHKERRQ(ierr);
929 ierr = DMCreateLocalVector(user->da, &user->lP); CHKERRQ(ierr); ierr = VecSet(user->lP, 0.0); CHKERRQ(ierr);
930 ierr = DMCreateLocalVector(user->da, &user->lNvert); CHKERRQ(ierr); ierr = VecSet(user->lNvert, 0.0); CHKERRQ(ierr);
931
932 // -- Group B: Solver Work Vectors (Global and Local) ---
933 ierr = VecDuplicate(user->P, &user->Phi); CHKERRQ(ierr); ierr = VecSet(user->Phi, 0.0); CHKERRQ(ierr);
934 ierr = VecDuplicate(user->lP, &user->lPhi); CHKERRQ(ierr); ierr = VecSet(user->lPhi, 0.0); CHKERRQ(ierr);
935
936 // --- Group C: Time-Stepping & Workspace Fields (Finest Level Only) ---
937 if (level == usermg->mglevels - 1) {
938 ierr = VecDuplicate(user->Ucont, &user->Ucont_o); CHKERRQ(ierr); ierr = VecSet(user->Ucont_o, 0.0); CHKERRQ(ierr);
939 ierr = VecDuplicate(user->Ucont, &user->Ucont_rm1); CHKERRQ(ierr); ierr = VecSet(user->Ucont_rm1, 0.0); CHKERRQ(ierr);
940 ierr = VecDuplicate(user->Ucat, &user->Ucat_o); CHKERRQ(ierr); ierr = VecSet(user->Ucat_o, 0.0); CHKERRQ(ierr);
941 ierr = VecDuplicate(user->P, &user->P_o); CHKERRQ(ierr); ierr = VecSet(user->P_o, 0.0); CHKERRQ(ierr);
942 ierr = VecDuplicate(user->lUcont, &user->lUcont_o); CHKERRQ(ierr); ierr = VecSet(user->lUcont_o, 0.0); CHKERRQ(ierr);
943 ierr = VecDuplicate(user->lUcont, &user->lUcont_rm1); CHKERRQ(ierr); ierr = VecSet(user->lUcont_rm1, 0.0); CHKERRQ(ierr);
944 ierr = DMCreateLocalVector(user->da, &user->lNvert_o); CHKERRQ(ierr); ierr = VecSet(user->lNvert_o, 0.0); CHKERRQ(ierr);
945 ierr = VecDuplicate(user->Nvert, &user->Nvert_o); CHKERRQ(ierr); ierr = VecSet(user->Nvert_o, 0.0); CHKERRQ(ierr);
946 }
947
948 // --- Group D: Grid Metrics (Cell-Centered) ---
949 ierr = DMCreateGlobalVector(user->fda, &user->Csi); CHKERRQ(ierr); ierr = VecSet(user->Csi, 0.0); CHKERRQ(ierr);
950 ierr = VecDuplicate(user->Csi, &user->Eta); CHKERRQ(ierr); ierr = VecSet(user->Eta, 0.0); CHKERRQ(ierr);
951 ierr = VecDuplicate(user->Csi, &user->Zet); CHKERRQ(ierr); ierr = VecSet(user->Zet, 0.0); CHKERRQ(ierr);
952 ierr = DMCreateGlobalVector(user->da, &user->Aj); CHKERRQ(ierr); ierr = VecSet(user->Aj, 0.0); CHKERRQ(ierr);
953
954 ierr = DMCreateLocalVector(user->fda, &user->lCsi); CHKERRQ(ierr); ierr = VecSet(user->lCsi, 0.0); CHKERRQ(ierr);
955 ierr = VecDuplicate(user->lCsi, &user->lEta); CHKERRQ(ierr); ierr = VecSet(user->lEta, 0.0); CHKERRQ(ierr);
956 ierr = VecDuplicate(user->lCsi, &user->lZet); CHKERRQ(ierr); ierr = VecSet(user->lZet, 0.0); CHKERRQ(ierr);
957 ierr = DMCreateLocalVector(user->da, &user->lAj); CHKERRQ(ierr); ierr = VecSet(user->lAj, 0.0); CHKERRQ(ierr);
958
959
960 // --- Group E: Grid Metrics (Face-Centered) ---
961 // Vector metrics are duplicated from Csi (DOF=3, fda-based)
962 ierr = VecDuplicate(user->Csi, &user->ICsi); CHKERRQ(ierr); ierr = VecSet(user->ICsi, 0.0); CHKERRQ(ierr);
963 ierr = VecDuplicate(user->Csi, &user->IEta); CHKERRQ(ierr); ierr = VecSet(user->IEta, 0.0); CHKERRQ(ierr);
964 ierr = VecDuplicate(user->Csi, &user->IZet); CHKERRQ(ierr); ierr = VecSet(user->IZet, 0.0); CHKERRQ(ierr);
965 ierr = VecDuplicate(user->Csi, &user->JCsi); CHKERRQ(ierr); ierr = VecSet(user->JCsi, 0.0); CHKERRQ(ierr);
966 ierr = VecDuplicate(user->Csi, &user->JEta); CHKERRQ(ierr); ierr = VecSet(user->JEta, 0.0); CHKERRQ(ierr);
967 ierr = VecDuplicate(user->Csi, &user->JZet); CHKERRQ(ierr); ierr = VecSet(user->JZet, 0.0); CHKERRQ(ierr);
968 ierr = VecDuplicate(user->Csi, &user->KCsi); CHKERRQ(ierr); ierr = VecSet(user->KCsi, 0.0); CHKERRQ(ierr);
969 ierr = VecDuplicate(user->Csi, &user->KEta); CHKERRQ(ierr); ierr = VecSet(user->KEta, 0.0); CHKERRQ(ierr);
970 ierr = VecDuplicate(user->Csi, &user->KZet); CHKERRQ(ierr); ierr = VecSet(user->KZet, 0.0); CHKERRQ(ierr);
971 // Scalar metrics are duplicated from Aj (DOF=1, da-based)
972 ierr = VecDuplicate(user->Aj, &user->IAj); CHKERRQ(ierr); ierr = VecSet(user->IAj, 0.0); CHKERRQ(ierr);
973 ierr = VecDuplicate(user->Aj, &user->JAj); CHKERRQ(ierr); ierr = VecSet(user->JAj, 0.0); CHKERRQ(ierr);
974 ierr = VecDuplicate(user->Aj, &user->KAj); CHKERRQ(ierr); ierr = VecSet(user->KAj, 0.0); CHKERRQ(ierr);
975
976 ierr = VecDuplicate(user->lCsi, &user->lICsi); CHKERRQ(ierr); ierr = VecSet(user->lICsi, 0.0); CHKERRQ(ierr);
977 ierr = VecDuplicate(user->lCsi, &user->lIEta); CHKERRQ(ierr); ierr = VecSet(user->lIEta, 0.0); CHKERRQ(ierr);
978 ierr = VecDuplicate(user->lCsi, &user->lIZet); CHKERRQ(ierr); ierr = VecSet(user->lIZet, 0.0); CHKERRQ(ierr);
979 ierr = VecDuplicate(user->lCsi, &user->lJCsi); CHKERRQ(ierr); ierr = VecSet(user->lJCsi, 0.0); CHKERRQ(ierr);
980 ierr = VecDuplicate(user->lCsi, &user->lJEta); CHKERRQ(ierr); ierr = VecSet(user->lJEta, 0.0); CHKERRQ(ierr);
981 ierr = VecDuplicate(user->lCsi, &user->lJZet); CHKERRQ(ierr); ierr = VecSet(user->lJZet, 0.0); CHKERRQ(ierr);
982 ierr = VecDuplicate(user->lCsi, &user->lKCsi); CHKERRQ(ierr); ierr = VecSet(user->lKCsi, 0.0); CHKERRQ(ierr);
983 ierr = VecDuplicate(user->lCsi, &user->lKEta); CHKERRQ(ierr); ierr = VecSet(user->lKEta, 0.0); CHKERRQ(ierr);
984 ierr = VecDuplicate(user->lCsi, &user->lKZet); CHKERRQ(ierr); ierr = VecSet(user->lKZet, 0.0); CHKERRQ(ierr);
985
986 ierr = VecDuplicate(user->lAj, &user->lIAj); CHKERRQ(ierr); ierr = VecSet(user->lIAj, 0.0); CHKERRQ(ierr);
987 ierr = VecDuplicate(user->lAj, &user->lJAj); CHKERRQ(ierr); ierr = VecSet(user->lJAj, 0.0); CHKERRQ(ierr);
988 ierr = VecDuplicate(user->lAj, &user->lKAj); CHKERRQ(ierr); ierr = VecSet(user->lKAj, 0.0); CHKERRQ(ierr);
989
990 // --- Group F: Cell/Face Center Coordinates and Grid Spacing ---
991 ierr = DMCreateGlobalVector(user->fda, &user->Cent); CHKERRQ(ierr); ierr = VecSet(user->Cent, 0.0); CHKERRQ(ierr);
992 ierr = DMCreateLocalVector(user->fda, &user->lCent); CHKERRQ(ierr); ierr = VecSet(user->lCent, 0.0); CHKERRQ(ierr);
993
994 ierr = VecDuplicate(user->Cent, &user->GridSpace); CHKERRQ(ierr); ierr = VecSet(user->GridSpace, 0.0); CHKERRQ(ierr);
995 ierr = VecDuplicate(user->lCent, &user->lGridSpace); CHKERRQ(ierr); ierr = VecSet(user->lGridSpace, 0.0); CHKERRQ(ierr);
996
997 // Face-center coordinate vectors are LOCAL to hold calculated values before scattering
998 ierr = VecDuplicate(user->lCent, &user->Centx); CHKERRQ(ierr); ierr = VecSet(user->Centx, 0.0); CHKERRQ(ierr);
999 ierr = VecDuplicate(user->lCent, &user->Centy); CHKERRQ(ierr); ierr = VecSet(user->Centy, 0.0); CHKERRQ(ierr);
1000 ierr = VecDuplicate(user->lCent, &user->Centz); CHKERRQ(ierr); ierr = VecSet(user->Centz, 0.0); CHKERRQ(ierr);
1001
1002 if(level == usermg->mglevels -1){
1003 // --- Group G: Turbulence Models (Finest Level Only) ---
1004 if (simCtx->les || simCtx->rans) {
1005 ierr = DMCreateGlobalVector(user->da, &user->Nu_t); CHKERRQ(ierr); ierr = VecSet(user->Nu_t, 0.0); CHKERRQ(ierr);
1006 ierr = DMCreateLocalVector(user->da, &user->lNu_t); CHKERRQ(ierr); ierr = VecSet(user->lNu_t, 0.0); CHKERRQ(ierr);
1007 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Turbulence viscosity (Nu_t) vectors created for LES/RANS model.\n");
1008 if(simCtx->les){
1009 ierr = DMCreateGlobalVector(user->da,&user->CS); CHKERRQ(ierr); ierr = VecSet(user->CS,0.0); CHKERRQ(ierr);
1010 ierr = DMCreateLocalVector(user->da,&user->lCs); CHKERRQ(ierr); ierr = VecSet(user->lCs,0.0); CHKERRQ(ierr);
1011 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Smagorinsky constant (CS) vectors created for LES model.\n");
1012 }
1013 // Add K_Omega etc. here as needed
1014
1015 // Note: Add any other vectors from the legacy MG_Initial here as needed.
1016 // For example: Rhs, Forcing, turbulence Vecs (K_Omega, Nu_t)...
1017
1018 }
1019 // --- Group H: Particle Methods
1020 if(simCtx->np>0){
1021 ierr = DMCreateGlobalVector(user->da,&user->ParticleCount); CHKERRQ(ierr); ierr = VecSet(user->ParticleCount,0.0); CHKERRQ(ierr);
1022 ierr = DMCreateLocalVector(user->da,&user->lParticleCount); CHKERRQ(ierr); ierr = VecSet(user->lParticleCount,0.0); CHKERRQ(ierr);
1023 // Scalar field to hold particle scalar property (e.g., temperature, concentration)
1024 ierr = DMCreateGlobalVector(user->da,&user->Psi); CHKERRQ(ierr); ierr = VecSet(user->Psi,0.0); CHKERRQ(ierr);
1025 ierr = DMCreateLocalVector(user->da,&user->lPsi); CHKERRQ(ierr); ierr = VecSet(user->lPsi,0.0); CHKERRQ(ierr);
1026 LOG_ALLOW(GLOBAL,LOG_DEBUG,"ParticleCount & Scalar(Psi) created for %d particles.\n",simCtx->np);
1027 }
1028 }
1029 // --- Group I: Boundary Condition vectors needed by the legacy FormBCS ---
1030 ierr = DMCreateGlobalVector(user->fda, &user->Bcs.Ubcs); CHKERRQ(ierr);
1031 ierr = VecSet(user->Bcs.Ubcs, 0.0); CHKERRQ(ierr);
1032 ierr = DMCreateGlobalVector(user->fda, &user->Bcs.Uch); CHKERRQ(ierr);
1033 ierr = VecSet(user->Bcs.Uch, 0.0); CHKERRQ(ierr);
1034
1035 if(level == usermg->mglevels - 1){
1036 if(simCtx->exec_mode == EXEC_MODE_POSTPROCESSOR){
1037 LOG_ALLOW(LOCAL, LOG_DEBUG, "Post-processor mode detected. Allocating derived field vectors.\n");
1038
1039 ierr = VecDuplicate(user->P, &user->P_nodal); CHKERRQ(ierr);
1040 ierr = VecSet(user->P_nodal, 0.0); CHKERRQ(ierr);
1041
1042 ierr = VecDuplicate(user->Ucat, &user->Ucat_nodal); CHKERRQ(ierr);
1043 ierr = VecSet(user->Ucat_nodal, 0.0); CHKERRQ(ierr);
1044
1045 ierr = VecDuplicate(user->P, &user->Qcrit); CHKERRQ(ierr);
1046 ierr = VecSet(user->Qcrit, 0.0); CHKERRQ(ierr);
1047
1048 LOG_ALLOW(LOCAL, LOG_DEBUG, "Derived field vectors P_nodal, Ucat_nodal, and Qcrit created.\n");
1049
1050 if(simCtx->np>0){
1051 ierr = VecDuplicate(user->Psi, &user->Psi_nodal); CHKERRQ(ierr);
1052 ierr = VecSet(user->Psi_nodal, 0.0); CHKERRQ(ierr);
1053
1054 LOG_ALLOW(LOCAL, LOG_DEBUG, "Derived field vector Psi_nodal created for particle scalar property.\n");
1055
1056 }
1057 }else{
1058 user->P_nodal = NULL;
1059 user->Ucat_nodal = NULL;
1060 user->Qcrit = NULL;
1061 user->Psi_nodal = NULL;
1062 }
1063 }
1064
1065 }
1066}
1067
1068 LOG_ALLOW(GLOBAL, LOG_INFO, "All simulation vectors created and initialized.\n");
1069
1070 PROFILE_FUNCTION_END;
1071 PetscFunctionReturn(0);
1072}
1073
1074#undef __FUNCT__
1075#define __FUNCT__ "UpdateLocalGhosts"
1076/**
1077 * @brief Updates the local vector (including ghost points) from its corresponding global vector.
1078 *
1079 * This function identifies the correct global vector, local vector, and DM based on the
1080 * provided fieldName and performs the standard PETSc DMGlobalToLocalBegin/End sequence.
1081 * Includes optional debugging output (max norms before/after).
1082 *
1083 * @param user The UserCtx structure containing the vectors and DMs.
1084 * @param fieldName The name of the field to update ("Ucat", "Ucont", "P", "Nvert", etc.).
1085 *
1086 * @return PetscErrorCode 0 on success, non-zero on failure.
1087 *
1088 * @note This function assumes the global vector associated with fieldName has already
1089 * been populated with the desired data (including any boundary conditions).
1090 */
1091PetscErrorCode UpdateLocalGhosts(UserCtx* user, const char *fieldName)
1092{
1093 PetscErrorCode ierr;
1094 PetscMPIInt rank;
1095 Vec globalVec = NULL;
1096 Vec localVec = NULL;
1097 DM dm = NULL; // The DM associated with this field pair
1098
1099 PetscFunctionBeginUser; // Use User version for application code
1100 PROFILE_FUNCTION_BEGIN;
1101 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
1102 LOG_ALLOW(GLOBAL, LOG_INFO, "Rank %d: Starting ghost update for field '%s'.\n", rank, fieldName);
1103
1104 // --- 1. Identify the correct Vectors and DM ---
1105 if (strcmp(fieldName, "Ucat") == 0) {
1106 globalVec = user->Ucat;
1107 localVec = user->lUcat;
1108 dm = user->fda;
1109 } else if (strcmp(fieldName, "Ucont") == 0) {
1110 globalVec = user->Ucont;
1111 localVec = user->lUcont;
1112 dm = user->fda;
1113 } else if (strcmp(fieldName, "P") == 0) {
1114 globalVec = user->P;
1115 localVec = user->lP;
1116 dm = user->da;
1117 } else if (strcmp(fieldName, "Csi") == 0) {
1118 globalVec = user->Csi;
1119 localVec = user->lCsi;
1120 dm = user->fda;
1121 } else if (strcmp(fieldName, "Eta") == 0) {
1122 globalVec = user->Eta;
1123 localVec = user->lEta;
1124 dm = user->fda;
1125 } else if (strcmp(fieldName, "Zet") == 0) {
1126 globalVec = user->Zet;
1127 localVec = user->lZet;
1128 dm = user->fda;
1129 }else if (strcmp(fieldName, "Nvert") == 0) {
1130 globalVec = user->Nvert;
1131 localVec = user->lNvert;
1132 dm = user->da;
1133 // Add other fields as needed
1134 } else if (strcmp(fieldName, "Aj") == 0) {
1135 globalVec = user->Aj;
1136 localVec = user->lAj;
1137 dm = user->da;
1138 } else if (strcmp(fieldName, "Cent") == 0) {
1139 globalVec = user->Cent;
1140 localVec = user->lCent;
1141 dm = user->fda;
1142 }else if (strcmp(fieldName, "GridSpace") == 0) {
1143 globalVec = user->GridSpace;
1144 localVec = user->lGridSpace;
1145 dm = user->fda;
1146 }else if (strcmp(fieldName,"ICsi") == 0){
1147 globalVec = user->ICsi;
1148 localVec = user->lICsi;
1149 dm = user->fda;
1150 }else if (strcmp(fieldName,"IEta") == 0){
1151 globalVec = user->IEta;
1152 localVec = user->lIEta;
1153 dm = user->fda;
1154 }else if (strcmp(fieldName,"IZet") == 0){
1155 globalVec = user->IZet;
1156 localVec = user->lIZet;
1157 dm = user->fda;
1158 }else if (strcmp(fieldName,"JCsi") == 0){
1159 globalVec = user->JCsi;
1160 localVec = user->lJCsi;
1161 dm = user->fda;
1162 }else if (strcmp(fieldName,"JEta") == 0){
1163 globalVec = user->JEta;
1164 localVec = user->lJEta;
1165 dm = user->fda;
1166 }else if (strcmp(fieldName,"JZet") == 0){
1167 globalVec = user->JZet;
1168 localVec = user->lJZet;
1169 dm = user->fda;
1170 }else if (strcmp(fieldName,"KCsi") == 0){
1171 globalVec = user->KCsi;
1172 localVec = user->lKCsi;
1173 dm = user->fda;
1174 }else if (strcmp(fieldName,"KEta") == 0){
1175 globalVec = user->KEta;
1176 localVec = user->lKEta;
1177 dm = user->fda;
1178 }else if (strcmp(fieldName,"KZet") == 0){
1179 globalVec = user->KZet;
1180 localVec = user->lKZet;
1181 dm = user->fda;
1182 }else if (strcmp(fieldName,"IAj") == 0){
1183 globalVec = user->IAj;
1184 localVec = user->lIAj;
1185 dm = user->da;
1186 }else if (strcmp(fieldName,"JAj") == 0){
1187 globalVec = user->JAj;
1188 localVec = user->lJAj;
1189 dm = user->da;
1190 }else if (strcmp(fieldName,"KAj") == 0){
1191 globalVec = user->KAj;
1192 localVec = user->lKAj;
1193 dm = user->da;
1194 }else if (strcmp(fieldName,"Phi") == 0){ // Pressure correction term.
1195 globalVec = user->Phi;
1196 localVec = user->lPhi;
1197 dm = user->da;
1198 }else if (strcmp(fieldName,"Psi") == 0){ // Particle scalar property.
1199 globalVec = user->Psi;
1200 localVec = user->lPsi;
1201 dm = user->da;
1202 }else {
1203 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Field '%s' not recognized for ghost update.", fieldName);
1204 }
1205
1206 // --- 2. Check if components were found ---
1207 if (!globalVec) {
1208 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Global vector for field '%s' is NULL.", fieldName);
1209 }
1210 if (!localVec) {
1211 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Local vector for field '%s' is NULL.", fieldName);
1212 }
1213 if (!dm) {
1214 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "DM for field '%s' is NULL.", fieldName);
1215 }
1216
1217 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Identified components for '%s': DM=%p, GlobalVec=%p, LocalVec=%p.\n",
1218 rank, fieldName, (void*)dm, (void*)globalVec, (void*)localVec);
1219
1220 // --- 3. Optional Debugging: Norm Before Update ---
1221 // Use your logging convention check
1222 // if (get_log_level() >= LOG_LEVEL_DEBUG && is_function_allowed("UpdateLocalGhosts")) { // Example check
1223 if(get_log_level() == LOG_DEBUG && is_function_allowed(__func__)){
1224 PetscReal norm_global_before;
1225 ierr = VecNorm(globalVec, NORM_INFINITY, &norm_global_before); CHKERRQ(ierr);
1226 LOG_ALLOW(GLOBAL, LOG_INFO,"Max norm '%s' (Global) BEFORE Ghost Update: %g\n", fieldName, norm_global_before);
1227 // Optional: Norm of local vector before update (might contain old ghost values)
1228 // PetscReal norm_local_before;
1229 // ierr = VecNorm(localVec, NORM_INFINITY, &norm_local_before); CHKERRQ(ierr);
1230 // LOG_ALLOW(GLOBAL, LOG_DEBUG,"Max norm '%s' (Local) BEFORE Ghost Update: %g\n", fieldName, norm_local_before);
1231 }
1232
1233 // --- 4. Perform the Global-to-Local Transfer (Ghost Update) ---
1234 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Calling DMGlobalToLocalBegin/End for '%s'.\n", rank, fieldName);
1235 ierr = DMGlobalToLocalBegin(dm, globalVec, INSERT_VALUES, localVec); CHKERRQ(ierr);
1236 ierr = DMGlobalToLocalEnd(dm, globalVec, INSERT_VALUES, localVec); CHKERRQ(ierr);
1237 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Completed DMGlobalToLocalBegin/End for '%s'.\n", rank, fieldName);
1238
1239 // --- 5. Optional Debugging: Norm After Update ---
1240 // Use your logging convention check
1241 // if (get_log_level() >= LOG_LEVEL_DEBUG && is_function_allowed("UpdateLocalGhosts")) { // Example check
1242 if(get_log_level() == LOG_DEBUG && is_function_allowed(__func__)){ // Using your specific check
1243 PetscReal norm_local_after;
1244 ierr = VecNorm(localVec, NORM_INFINITY, &norm_local_after); CHKERRQ(ierr);
1245 LOG_ALLOW(GLOBAL, LOG_INFO,"Max norm '%s' (Local) AFTER Ghost Update: %g\n", fieldName, norm_local_after);
1246
1247 // --- 6. Optional Debugging: Specific Point Checks (Example for Ucat on Rank 0/1) ---
1248 // (Keep this conditional if it's only for specific debug scenarios)
1249 if (strcmp(fieldName, "Ucat") == 0) { // Only do detailed checks for Ucat for now
1250 PetscMPIInt rank_test;
1251 MPI_Comm_rank(PETSC_COMM_WORLD, &rank_test);
1252
1253 // Get Local Info needed for indexing checks
1254 DMDALocalInfo info_check;
1255 ierr = DMDAGetLocalInfo(dm, &info_check); CHKERRQ(ierr); // Use the correct dm
1256
1257 // Buffer for array pointer
1258 Cmpnts ***lUcat_arr_test = NULL;
1259 PetscErrorCode ierr_test = 0;
1260
1261 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Testing '%s' access immediately after ghost update...\n", rank_test, fieldName);
1262 ierr_test = DMDAVecGetArrayDOFRead(dm, localVec, &lUcat_arr_test); // Use correct dm and localVec
1263
1264 if (ierr_test) {
1265 LOG_ALLOW(LOCAL, LOG_ERROR, "Rank %d: ERROR %d getting '%s' array after ghost update!\n", rank_test, ierr_test, fieldName);
1266 } else if (!lUcat_arr_test) {
1267 LOG_ALLOW(LOCAL, LOG_ERROR, "Rank %d: ERROR NULL pointer getting '%s' array after ghost update!\n", rank_test, fieldName);
1268 }
1269 else {
1270 // Check owned interior point (e.g., first interior point)
1271 PetscInt k_int = info_check.zs + (info_check.zm > 1 ? 1 : 0); // Global k index (at least zs+1 if possible)
1272 PetscInt j_int = info_check.ys + (info_check.ym > 1 ? 1 : 0); // Global j index
1273 PetscInt i_int = info_check.xs + (info_check.xm > 1 ? 1 : 0); // Global i index
1274 // Ensure indices are within global bounds if domain is very small
1275 //if (k_int >= info_check.mz-1) k_int = info_check.mz-2; if (k_int < 1) k_int = 1;
1276 //if (j_int >= info_check.my-1) j_int = info_check.my-2; if (j_int < 1) j_int = 1;
1277 // if (i_int >= info_check.mx-1) i_int = info_check.mx-2; if (i_int < 1) i_int = 1;
1278 // clamp k_int to [1 .. mz-2]
1279 if (k_int >= info_check.mz - 1) {
1280 k_int = info_check.mz - 2;
1281 }
1282 if (k_int < 1) {
1283 k_int = 1;
1284 }
1285
1286 // clamp j_int to [1 .. my-2]
1287 if (j_int >= info_check.my - 1) {
1288 j_int = info_check.my - 2;
1289 }
1290 if (j_int < 1) {
1291 j_int = 1;
1292 }
1293
1294 // clamp i_int to [1 .. mx-2]
1295 if (i_int >= info_check.mx - 1) {
1296 i_int = info_check.mx - 2;
1297 }
1298 if (i_int < 1) {
1299 i_int = 1;
1300 }
1301
1302 // Only attempt read if indices are actually owned (relevant for multi-rank)
1303 if (k_int >= info_check.zs && k_int < info_check.zs + info_check.zm &&
1304 j_int >= info_check.ys && j_int < info_check.ys + info_check.ym &&
1305 i_int >= info_check.xs && i_int < info_check.xs + info_check.xm)
1306 {
1307 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);
1308 Cmpnts test_val_owned_interior = lUcat_arr_test[k_int][j_int][i_int];
1309 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: SUCCESS reading owned interior: x=%g\n", rank_test, test_val_owned_interior.x);
1310 } else {
1311 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);
1312 }
1313
1314
1315 // Check owned boundary point (e.g., first owned point)
1316 PetscInt k_bnd = info_check.zs; // Global k index
1317 PetscInt j_bnd = info_check.ys; // Global j index
1318 PetscInt i_bnd = info_check.xs; // Global i index
1319 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);
1320 Cmpnts test_val_owned_boundary = lUcat_arr_test[k_bnd][j_bnd][i_bnd];
1321 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: SUCCESS reading owned boundary: x=%g\n", rank_test, test_val_owned_boundary.x);
1322
1323
1324 // Check ghost point (e.g., one layer below in k, if applicable)
1325 if (info_check.zs > 0) { // Only if there's a rank below
1326 PetscInt k_ghost = info_check.zs - 1;
1327 PetscInt j_ghost = info_check.ys; // Use start of owned y, simple example
1328 PetscInt i_ghost = info_check.xs; // Use start of owned x, simple example
1329 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Attempting test read GHOST [%d][%d][%d] (Global)\n", rank_test, k_ghost, j_ghost, i_ghost);
1330 Cmpnts test_val_ghost = lUcat_arr_test[k_ghost][j_ghost][i_ghost];
1331 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: SUCCESS reading ghost: x=%g\n", rank_test, test_val_ghost.x);
1332 } else {
1333 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Skipping ghost test read (zs=0).\n", rank_test);
1334 }
1335
1336 // Restore the array
1337 ierr_test = DMDAVecRestoreArrayDOFRead(dm, localVec, &lUcat_arr_test);
1338 if(ierr_test){ LOG_ALLOW(LOCAL, LOG_ERROR, "Rank %d: ERROR %d restoring '%s' array after test read!\n", rank_test, ierr_test, fieldName); }
1339 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Finished testing '%s' access.\n", rank_test, fieldName);
1340 }
1341 } // end if Ucat
1342 } // end debug logging check
1343
1344 LOG_ALLOW(GLOBAL, LOG_INFO, "Rank %d: Completed ghost update for field '%s'.\n", rank, fieldName);
1345 PROFILE_FUNCTION_END;
1346 PetscFunctionReturn(0);
1347}
1348
1349#undef __FUNCT__
1350#define __FUNCT__ "SetupBoundaryConditions"
1351/**
1352 * @brief (Orchestrator) Sets up all boundary conditions for the simulation.
1353 */
1354PetscErrorCode SetupBoundaryConditions(SimCtx *simCtx)
1355{
1356 PetscErrorCode ierr;
1357 PetscFunctionBeginUser;
1358
1359 PROFILE_FUNCTION_BEGIN;
1360
1361 LOG_ALLOW(GLOBAL,LOG_INFO, "--- Setting up Boundary Conditions ---\n");
1362
1363 // --- Parse and Adapt for each block on the finest level ---
1364 LOG_ALLOW(GLOBAL,LOG_INFO,"Parsing BC configuration file and adapting to legacy system for finest grid.\n");
1365 UserCtx *user_finest = simCtx->usermg.mgctx[simCtx->usermg.mglevels-1].user;
1366 for (PetscInt bi = 0; bi < simCtx->block_number; bi++) {
1367 LOG_ALLOW(GLOBAL,LOG_DEBUG, " -> Processing Block %d:\n", bi);
1368
1369 // --- Generate the filename for the current block ---
1370 const char *current_bc_filename = simCtx->bcs_files[bi];
1371 LOG_ALLOW(GLOBAL,LOG_DEBUG," -> Processing Block %d using config file '%s'\n", bi, current_bc_filename);
1372 // This will populate user_finest[bi].boundary_faces
1373 ierr = ParseAllBoundaryConditions(&user_finest[bi],current_bc_filename); CHKERRQ(ierr);
1374
1375 // Call the adapter to translate into the legacy format
1376 ierr = TranslateModernBCsToLegacy(&user_finest[bi]); CHKERRQ(ierr);
1377
1378 // Call the function to calculate the center of the inlet face, which may be used to calculate Boundary values.
1379 ierr = CalculateInletCenter(&user_finest[bi]); CHKERRQ(ierr);
1380 }
1381
1382 LOG_ALLOW(GLOBAL,LOG_INFO, "--- Boundary Conditions setup complete ---\n");
1383
1384
1385 PROFILE_FUNCTION_END;
1386 PetscFunctionReturn(0);
1387}
1388
1389/**
1390 * @brief Allocates a 3D array of PetscReal values using PetscCalloc.
1391 *
1392 * This function dynamically allocates memory for a 3D array of PetscReal values
1393 * with dimensions nz (layers) x ny (rows) x nx (columns). It uses PetscCalloc1
1394 * to ensure the memory is zero-initialized.
1395 *
1396 * The allocation is done in three steps:
1397 * 1. Allocate an array of nz pointers (one for each layer).
1398 * 2. Allocate a contiguous block for nz*ny row pointers and assign each layer’s row pointers.
1399 * 3. Allocate a contiguous block for all nz*ny*nx PetscReal values.
1400 *
1401 * This setup allows the array to be accessed as array[k][j][i], and the memory
1402 * for the data is contiguous, which improves cache efficiency.
1403 *
1404 * @param[out] array Pointer to the 3D array to be allocated.
1405 * @param[in] nz Number of layers (z-direction).
1406 * @param[in] ny Number of rows (y-direction).
1407 * @param[in] nx Number of columns (x-direction).
1408 *
1409 * @return PetscErrorCode 0 on success, nonzero on failure.
1410 */
1411PetscErrorCode Allocate3DArrayScalar(PetscReal ****array, PetscInt nz, PetscInt ny, PetscInt nx)
1412{
1413 PetscErrorCode ierr;
1414 PetscReal ***data;
1415 PetscReal *dataContiguous;
1416 PetscInt k, j;
1417
1418 PetscFunctionBegin;
1419 /* Step 1: Allocate memory for an array of nz layer pointers (zero-initialized) */
1420 ierr = PetscCalloc1(nz, &data); CHKERRQ(ierr);
1421
1422 /* Step 2: Allocate memory for all row pointers (nz * ny pointers) */
1423 ierr = PetscCalloc1(nz * ny, &data[0]); CHKERRQ(ierr);
1424 for (k = 1; k < nz; k++) {
1425 data[k] = data[0] + k * ny;
1426 }
1427
1428 /* Step 3: Allocate one contiguous block for all data elements (nz*ny*nx) */
1429 ierr = PetscCalloc1(nz * ny * nx, &dataContiguous); CHKERRQ(ierr);
1430
1431 /* Build the 3D pointer structure: each row pointer gets the correct segment of data */
1432 for (k = 0; k < nz; k++) {
1433 for (j = 0; j < ny; j++) {
1434 data[k][j] = dataContiguous + (k * ny + j) * nx;
1435 /* Memory is already zeroed by PetscCalloc1, so no manual initialization is needed */
1436 }
1437 }
1438 *array = data;
1439 PetscFunctionReturn(0);
1440}
1441
1442/**
1443 * @brief Deallocates a 3D array of PetscReal values allocated by Allocate3DArrayScalar.
1444 *
1445 * This function frees the memory allocated for a 3D array of PetscReal values.
1446 * It assumes the memory was allocated using Allocate3DArrayScalar, which allocated
1447 * three separate memory blocks: one for the contiguous data, one for the row pointers,
1448 * and one for the layer pointers.
1449 *
1450 * @param[in] array Pointer to the 3D array to be deallocated.
1451 * @param[in] nz Number of layers (z-direction).
1452 * @param[in] ny Number of rows (y-direction).
1453 *
1454 * @return PetscErrorCode 0 on success, nonzero on failure.
1455 */
1456PetscErrorCode Deallocate3DArrayScalar(PetscReal ***array, PetscInt nz, PetscInt ny)
1457{
1458 PetscErrorCode ierr;
1459
1460 PetscFunctionBegin;
1461 if (!array || !array[0] || !array[0][0] ) { // Added more robust check
1462 LOG_ALLOW(GLOBAL, LOG_WARNING, "Deallocate3DArrayScalar called with potentially unallocated or NULL array.\n");
1463 if (array) {
1464 if (array[0]) { // Check if row pointers might exist
1465 // Cannot safely access array[0][0] if array[0] might be invalid/freed
1466 // Standard deallocation below assumes valid pointers.
1467 ierr = PetscFree(array[0]); CHKERRQ(ierr); // Free row pointers if they exist
1468 }
1469 ierr = PetscFree(array); CHKERRQ(ierr); // Free layer pointers if they exist
1470 }
1471 PetscFunctionReturn(0);
1472 }
1473
1474 // --- Standard Deallocation (assuming valid allocation) ---
1475
1476 /* 1. Free the contiguous block of PetscReal values.
1477 The starting address was stored in array[0][0]. */
1478 ierr = PetscFree(array[0][0]); CHKERRQ(ierr); // Free the ACTUAL DATA
1479
1480 /* 2. Free the contiguous block of row pointers.
1481 The starting address was stored in array[0]. */
1482 ierr = PetscFree(array[0]); CHKERRQ(ierr); // Free the ROW POINTERS
1483
1484 /* 3. Free the layer pointer array.
1485 The starting address is 'array' itself. */
1486 ierr = PetscFree(array); CHKERRQ(ierr); // Free the LAYER POINTERS
1487
1488 PetscFunctionReturn(0);
1489}
1490
1491/**
1492 * @brief Allocates a 3D array of Cmpnts structures using PetscCalloc.
1493 *
1494 * This function dynamically allocates memory for a 3D array of Cmpnts (vector) structures
1495 * with dimensions nz (layers) x ny (rows) x nx (columns). It uses PetscCalloc1 to ensure
1496 * that all allocated memory is zero-initialized.
1497 *
1498 * The allocation procedure is similar to Allocate3DArrayScalar:
1499 * 1. Allocate an array of nz pointers (one for each layer).
1500 * 2. Allocate a contiguous block for nz*ny row pointers.
1501 * 3. Allocate one contiguous block for nz*ny*nx Cmpnts structures.
1502 *
1503 * After allocation, the array can be accessed as array[k][j][i] and each element
1504 * (a Cmpnts structure) will have its x, y, and z fields initialized to 0.0.
1505 *
1506 * @param[out] array Pointer to the 3D array to be allocated.
1507 * @param[in] nz Number of layers in the z-direction.
1508 * @param[in] ny Number of rows in the y-direction.
1509 * @param[in] nx Number of columns in the x-direction.
1510 *
1511 * @return PetscErrorCode 0 on success, nonzero on failure.
1512 */
1513PetscErrorCode Allocate3DArrayVector(Cmpnts ****array, PetscInt nz, PetscInt ny, PetscInt nx)
1514{
1515 PetscErrorCode ierr;
1516 Cmpnts ***data;
1517 Cmpnts *dataContiguous;
1518 PetscInt k, j;
1519 PetscMPIInt rank;
1520
1521 PetscFunctionBegin;
1522
1523 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
1524
1525 /* Step 1: Allocate memory for nz layer pointers (zeroed) */
1526 ierr = PetscCalloc1(nz, &data); CHKERRQ(ierr);
1527
1528 LOG_ALLOW(LOCAL,LOG_DEBUG," [Rank %d] memory allocated for outermost layer (%d k-layer pointers).\n",rank,nz);
1529
1530 /* Step 2: Allocate memory for all row pointers (nz * ny pointers) */
1531 ierr = PetscCalloc1(nz * ny, &data[0]); CHKERRQ(ierr);
1532 for (k = 1; k < nz; k++) {
1533 data[k] = data[0] + k * ny;
1534 }
1535
1536 LOG_ALLOW(LOCAL,LOG_DEBUG,"[Rank %d] memory allocated for %dx%d row pointers.\n",rank,nz,ny);
1537
1538 /* Step 3: Allocate one contiguous block for nz*ny*nx Cmpnts structures (zeroed) */
1539 ierr = PetscCalloc1(nz * ny * nx, &dataContiguous); CHKERRQ(ierr);
1540
1541 LOG_ALLOW(GLOBAL,LOG_DEBUG,"[Rank %d] memory allocated for contigous block of %dx%dx%d Cmpnts structures).\n",rank,nz,ny,nx);
1542
1543 /* Build the 3D pointer structure for vector data */
1544 for (k = 0; k < nz; k++) {
1545 for (j = 0; j < ny; j++) {
1546 data[k][j] = dataContiguous + (k * ny + j) * nx;
1547 /* The PetscCalloc1 call has already initialized each Cmpnts to zero. */
1548 }
1549 }
1550
1551 LOG_ALLOW(GLOBAL,LOG_DEBUG,"[Rank %d] 3D pointer structure for vector data created. \n",rank);
1552
1553 *array = data;
1554 PetscFunctionReturn(0);
1555}
1556
1557/**
1558 * @brief Deallocates a 3D array of Cmpnts structures allocated by Allocate3DArrayVector.
1559 *
1560 * This function frees the memory allocated for a 3D array of Cmpnts structures.
1561 * It assumes the memory was allocated using Allocate3DArrayVector, which created three
1562 * separate memory blocks: one for the contiguous vector data, one for the row pointers,
1563 * and one for the layer pointers.
1564 *
1565 * @param[in] array Pointer to the 3D array to be deallocated.
1566 * @param[in] nz Number of layers in the z-direction.
1567 * @param[in] ny Number of rows in the y-direction.
1568 *
1569 * @return PetscErrorCode 0 on success, nonzero on failure.
1570 */
1571 PetscErrorCode Deallocate3DArrayVector(Cmpnts ***array, PetscInt nz, PetscInt ny)
1572{
1573 PetscErrorCode ierr;
1574
1575 PetscFunctionBegin;
1576 // If array is NULL or hasn't been allocated properly, just return.
1577 if (!array || !array[0] || !array[0][0] ) {
1578 LOG_ALLOW(GLOBAL, LOG_WARNING, "Deallocate3DArrayVector called with potentially unallocated or NULL array.\n");
1579 // Attempt to free what might exist, but be cautious
1580 if (array) {
1581 if (array[0]) { // Check if row pointers were allocated
1582 // We don't have a direct pointer to the contiguous data block
1583 // saved separately in this allocation scheme. The allocation relies
1584 // on array[0][0] pointing to it. If array[0] was freed first,
1585 // accessing array[0][0] is unsafe.
1586 // The allocation scheme where the contiguous data block is not
1587 // stored separately makes safe deallocation tricky if freeing
1588 // happens out of order or if parts are NULL.
1589
1590 // A SAFER ALLOCATION/DEALLOCATION would store the data pointer separately.
1591 // Given the current allocation scheme, the order MUST be:
1592 // 1. Free the data block (pointed to by array[0][0])
1593 // 2. Free the row pointer block (pointed to by array[0])
1594 // 3. Free the layer pointer block (pointed to by array)
1595
1596 // Let's assume the allocation was successful and pointers are valid.
1597 // Get pointer to the contiguous data block *before* freeing row pointers
1598 Cmpnts *dataContiguous = array[0][0];
1599 ierr = PetscFree(dataContiguous); CHKERRQ(ierr); // Free data block
1600
1601 // Now free the row pointers block
1602 ierr = PetscFree(array[0]); CHKERRQ(ierr); // Free row pointers
1603
1604 }
1605 // Finally, free the array of layer pointers
1606 ierr = PetscFree(array); CHKERRQ(ierr);
1607 }
1608 PetscFunctionReturn(0); // Return gracefully if input was NULL initially
1609 }
1610
1611
1612 // --- Standard Deallocation (assuming valid allocation) ---
1613
1614 /* 1. Free the contiguous block of Cmpnts structures.
1615 The starting address was stored in array[0][0] by Allocate3DArrayVector. */
1616 ierr = PetscFree(array[0][0]); CHKERRQ(ierr); // Free the ACTUAL DATA
1617
1618 /* 2. Free the contiguous block of row pointers.
1619 The starting address was stored in array[0]. */
1620 ierr = PetscFree(array[0]); CHKERRQ(ierr); // Free the ROW POINTERS
1621
1622 /* 3. Free the layer pointer array.
1623 The starting address is 'array' itself. */
1624 ierr = PetscFree(array); CHKERRQ(ierr); // Free the LAYER POINTERS
1625
1626 PetscFunctionReturn(0);
1627}
1628
1629#undef __FUNCT__
1630#define __FUNCT__ "GetOwnedCellRange"
1631/**
1632 * @brief Gets the global starting index of cells owned by this rank and the number of such cells.
1633 *
1634 * A cell's global index is considered the same as its origin node's global index.
1635 * This function assumes a node-centered DMDA where `info_nodes` provides all necessary
1636 * information:
1637 * - `info_nodes->xs, ys, zs`: Global starting index of the first node owned by this rank (excluding ghosts).
1638 * - `info_nodes->xm, ym, zm`: Number of nodes owned by this rank in each dimension (excluding ghosts).
1639 * - `info_nodes->mx, my, mz`: Total number of global nodes in each dimension for the entire domain.
1640 *
1641 * A cell `C_k` (0-indexed) is defined by its origin node `N_k` and extends to node `N_{k+1}`.
1642 * Thus, the last node in the global domain cannot be an origin for a cell. The last possible
1643 * cell origin node index is `PhysicalGlobalNodesInDim - 2`.
1644 *
1645 * @warning THIS FUNCTION MAKES A CRITICAL ASSUMPTION: It assumes the DMDA it
1646 * receives information from was created with dimensions of
1647 * (physical_nodes_x + 1), (physical_nodes_y + 1), etc., where the `+1`
1648 * node is a "user-defined ghost" not meant for physical calculations.
1649 * It will internally subtract 1 from the global dimensions (mx, my, mz)
1650 * before calculating the number of cells.
1651 *
1652 * @param[in] info_nodes Pointer to the DMDALocalInfo struct for the current rank.
1653 * This struct contains local ownership information (xs, xm, etc.)
1654 * and global domain dimensions (mx, my, mz for nodes).
1655 * @param[in] dim The dimension for which to get the cell range (0 for i/x, 1 for j/y, 2 for k/z).
1656 * @param[out] xs_cell_global_out Pointer to store the global index of the first cell whose origin node
1657 * is owned by this rank.
1658 * @param[out] xm_cell_local_out Pointer to store the number of cells for which this rank owns the
1659 * origin node.
1660 *
1661 * @return PetscErrorCode 0 on success, or an error code on failure.
1662 */
1663PetscErrorCode GetOwnedCellRange(const DMDALocalInfo *info_nodes,
1664 PetscInt dim,
1665 PetscInt *xs_cell_global_out,
1666 PetscInt *xm_cell_local_out)
1667{
1668 PetscErrorCode ierr = 0; // Standard PETSc error code, not explicitly set here but good practice.
1669 PetscInt xs_node_global_rank; // Global index of the first node owned by this rank in the specified dimension.
1670 PetscInt num_nodes_owned_rank; // Number of nodes owned by this rank in this dimension (local count, excluding ghosts).
1671 PetscInt GlobalNodesInDim_from_info; // Total number of DA points in this dimension, from DMDALocalInfo.
1672
1673 PetscFunctionBeginUser;
1674
1675 // --- 1. Input Validation ---
1676 if (!info_nodes || !xs_cell_global_out || !xm_cell_local_out) {
1677 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Null pointer passed to GetOwnedCellRange.");
1678 }
1679
1680 // --- 2. Extract Node Ownership and Global Dimension Information from DMDALocalInfo ---
1681 if (dim == 0) { // I-direction
1682 xs_node_global_rank = info_nodes->xs;
1683 num_nodes_owned_rank = info_nodes->xm;
1684 GlobalNodesInDim_from_info = info_nodes->mx;
1685 } else if (dim == 1) { // J-direction
1686 xs_node_global_rank = info_nodes->ys;
1687 num_nodes_owned_rank = info_nodes->ym;
1688 GlobalNodesInDim_from_info = info_nodes->my;
1689 } else if (dim == 2) { // K-direction
1690 xs_node_global_rank = info_nodes->zs;
1691 num_nodes_owned_rank = info_nodes->zm;
1692 GlobalNodesInDim_from_info = info_nodes->mz;
1693 } else {
1694 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d in GetOwnedCellRange. Must be 0, 1, or 2.", dim);
1695 }
1696
1697 // --- 3. Correct for User-Defined Ghost Node ---
1698 // Per the function's contract (@warning), the DA size includes an extra, non-physical
1699 // node. We subtract 1 to get the true number of physical nodes for cell calculations.
1700 const PetscInt physical_nodes_in_dim = GlobalNodesInDim_from_info - 1;
1701
1702 // --- 4. Handle Edge Cases for Physical Domain Size ---
1703 // If the physical domain has 0 or 1 node, no cells can be formed.
1704 if (physical_nodes_in_dim <= 1) {
1705 *xs_cell_global_out = xs_node_global_rank; // Still report the rank's starting node
1706 *xm_cell_local_out = 0; // But 0 cells
1707 PetscFunctionReturn(0);
1708 }
1709
1710 // --- 5. Determine Cell Ownership Based on Corrected Node Ownership ---
1711 // The first cell this rank *could* define has its origin at the first node this rank owns.
1712 *xs_cell_global_out = xs_node_global_rank;
1713
1714 // If the rank owns no nodes in this dimension, it can't form any cell origins.
1715 if (num_nodes_owned_rank == 0) {
1716 *xm_cell_local_out = 0;
1717 } else {
1718 // --- BUG FIX APPLIED HERE ---
1719 // The previous logic incorrectly assumed a cell's end node (N_{k+1}) must be on the
1720 // same rank as its origin node (N_k). The correct logic is to find the intersection
1721 // between the nodes this rank owns and the nodes that are valid origins globally.
1722
1723 // The first node owned by the rank is its first potential origin.
1724 PetscInt first_owned_origin = xs_node_global_rank;
1725
1726 // The absolute last node owned by this rank. Any node up to and including this one
1727 // is a potential cell origin from this rank's perspective.
1728 PetscInt last_node_owned_by_rank = xs_node_global_rank + num_nodes_owned_rank - 1;
1729
1730 // The absolute last node in the entire PHYSICAL domain that can serve as a cell origin.
1731 // If there are `N` physical nodes (0 to N-1), this index is `N-2`.
1732 PetscInt last_possible_origin_global_idx = physical_nodes_in_dim - 2;
1733
1734 // The actual last origin this rank can provide is the *minimum* of what it owns
1735 // and what is globally possible. This correctly handles both ranks in the middle of
1736 // the domain and the very last rank.
1737 PetscInt actual_last_origin_this_rank_can_form = PetscMin(last_node_owned_by_rank, last_possible_origin_global_idx);
1738
1739 // If the first potential origin this rank owns is already beyond the actual last
1740 // origin it can form, then this rank forms no valid cell origins. This happens if
1741 // the rank only owns the very last physical node.
1742 if (first_owned_origin > actual_last_origin_this_rank_can_form) {
1743 *xm_cell_local_out = 0;
1744 } else {
1745 // The number of cells is the count of valid origins this rank owns.
1746 // (Count = Last Index - First Index + 1)
1747 *xm_cell_local_out = actual_last_origin_this_rank_can_form - first_owned_origin + 1;
1748 }
1749 }
1750
1751 PetscFunctionReturn(ierr);
1752}
1753
1754#undef __FUNCT__
1755#define __FUNCT__ "ComputeAndStoreNeighborRanks"
1756/**
1757 * @brief Computes and stores the Cartesian neighbor ranks for the DMDA decomposition.
1758 *
1759 * This function retrieves the neighbor information from the primary DMDA (user->da)
1760 * and stores the face neighbors (xm, xp, ym, yp, zm, zp) in the user->neighbors structure.
1761 * It assumes a standard PETSc ordering for the neighbors array returned by DMDAGetNeighbors.
1762 * If DMDAGetNeighbors returns a negative rank that is not MPI_PROC_NULL (which can happen
1763 * in some PETSc/MPI configurations for non-periodic boundaries if not fully standard),
1764 * this function will sanitize it to MPI_PROC_NULL to prevent issues.
1765 *
1766 * @param[in,out] user Pointer to the UserCtx structure where neighbor info will be stored.
1767 *
1768 * @return PetscErrorCode 0 on success, non-zero on failure.
1769 */
1770PetscErrorCode ComputeAndStoreNeighborRanks(UserCtx *user)
1771{
1772 PetscErrorCode ierr;
1773 PetscMPIInt rank;
1774 PetscMPIInt size; // MPI communicator size
1775 const PetscMPIInt *neighbor_ranks_ptr; // Pointer to raw neighbor data from PETSc
1776
1777 PetscFunctionBeginUser;
1778 PROFILE_FUNCTION_BEGIN;
1779 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
1780 ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size); CHKERRQ(ierr); // Get MPI size for validation
1781
1782 LOG_ALLOW(GLOBAL, LOG_INFO, "Rank %d: Computing DMDA neighbor ranks.\n", rank);
1783
1784 if (!user || !user->da) {
1785 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "UserCtx or user->da is NULL in ComputeAndStoreNeighborRanks.");
1786 }
1787
1788 // Get the neighbor information from the DMDA
1789 // neighbor_ranks_ptr will point to an internal PETSc array of 27 ranks.
1790 ierr = DMDAGetNeighbors(user->da, &neighbor_ranks_ptr); CHKERRQ(ierr);
1791
1792 // Log the raw values from DMDAGetNeighbors for boundary-relevant directions for debugging
1793 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",
1794 rank,
1795 neighbor_ranks_ptr[12], neighbor_ranks_ptr[14],
1796 neighbor_ranks_ptr[10], neighbor_ranks_ptr[16],
1797 neighbor_ranks_ptr[4], neighbor_ranks_ptr[22],
1798 (int)MPI_PROC_NULL);
1799
1800 // PETSc standard indices for 3D face neighbors from the 27-point stencil:
1801 // Index = k_offset*9 + j_offset*3 + i_offset (where offsets -1,0,1 map to 0,1,2)
1802 // Center: (i_off=1, j_off=1, k_off=1) => 1*9 + 1*3 + 1 = 13
1803 // X-min: (i_off=0, j_off=1, k_off=1) => 1*9 + 1*3 + 0 = 12
1804 // X-plus: (i_off=2, j_off=1, k_off=1) => 1*9 + 1*3 + 2 = 14
1805 // Y-min: (i_off=1, j_off=0, k_off=1) => 1*9 + 0*3 + 1 = 10
1806 // Y-plus: (i_off=1, j_off=2, k_off=1) => 1*9 + 2*3 + 1 = 16
1807 // Z-min: (i_off=1, j_off=1, k_off=0) => 0*9 + 1*3 + 1 = 4
1808 // Z-plus: (i_off=1, j_off=1, k_off=2) => 2*9 + 1*3 + 1 = 22
1809
1810 if (neighbor_ranks_ptr[13] != rank) {
1811 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",
1812 rank, neighbor_ranks_ptr[13], rank);
1813 // This warning is important. If the center isn't the current rank, the offsets are likely wrong.
1814 // However, PETSc should ensure this unless the DM is too small for a 3x3x3 stencil.
1815 }
1816
1817 // Assign and sanitize each neighbor rank
1818 PetscMPIInt temp_neighbor;
1819
1820 temp_neighbor = neighbor_ranks_ptr[12]; // xm
1821 if (temp_neighbor < 0 || temp_neighbor >= size) {
1822 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);
1823 user->neighbors.rank_xm = MPI_PROC_NULL;
1824 } else {
1825 user->neighbors.rank_xm = temp_neighbor;
1826 }
1827
1828 temp_neighbor = neighbor_ranks_ptr[14]; // xp
1829 if (temp_neighbor < 0 || temp_neighbor >= size) {
1830 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);
1831 user->neighbors.rank_xp = MPI_PROC_NULL;
1832 } else {
1833 user->neighbors.rank_xp = temp_neighbor;
1834 }
1835
1836 temp_neighbor = neighbor_ranks_ptr[10]; // ym
1837 if (temp_neighbor < 0 || temp_neighbor >= size) {
1838 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);
1839 user->neighbors.rank_ym = MPI_PROC_NULL;
1840 } else {
1841 user->neighbors.rank_ym = temp_neighbor;
1842 }
1843
1844 temp_neighbor = neighbor_ranks_ptr[16]; // yp
1845 if (temp_neighbor < 0 || temp_neighbor >= size) {
1846 // The log for index 16 was "zm" in your output, should be yp
1847 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);
1848 user->neighbors.rank_yp = MPI_PROC_NULL;
1849 } else {
1850 user->neighbors.rank_yp = temp_neighbor;
1851 }
1852
1853 temp_neighbor = neighbor_ranks_ptr[4]; // zm
1854 if (temp_neighbor < 0 || temp_neighbor >= size) {
1855 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);
1856 user->neighbors.rank_zm = MPI_PROC_NULL;
1857 } else {
1858 user->neighbors.rank_zm = temp_neighbor;
1859 }
1860
1861 temp_neighbor = neighbor_ranks_ptr[22]; // zp
1862 if (temp_neighbor < 0 || temp_neighbor >= size) {
1863 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);
1864 user->neighbors.rank_zp = MPI_PROC_NULL;
1865 } else {
1866 user->neighbors.rank_zp = temp_neighbor;
1867 }
1868
1869 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,
1870 user->neighbors.rank_xm, user->neighbors.rank_xp,
1871 user->neighbors.rank_ym, user->neighbors.rank_yp,
1872 user->neighbors.rank_zm, user->neighbors.rank_zp);
1873 PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT); // Ensure logs are flushed
1874
1875 // Note: neighbor_ranks_ptr memory is managed by PETSc, do not free it.
1876 PROFILE_FUNCTION_END;
1877 PetscFunctionReturn(0);
1878}
1879
1880#undef __FUNCT__
1881#define __FUNCT__ "SetDMDAProcLayout"
1882/**
1883 * @brief Sets the processor layout for a given DMDA based on PETSc options.
1884 *
1885 * Reads the desired number of processors in x, y, and z directions using
1886 * PETSc options (e.g., -dm_processors_x, -dm_processors_y, -dm_processors_z).
1887 * If an option is not provided for a direction, PETSC_DECIDE is used for that direction.
1888 * Applies the layout using DMDASetNumProcs.
1889 *
1890 * Also stores the retrieved/decided values in user->procs_x/y/z if user context is provided.
1891 *
1892 * @param dm The DMDA object to configure the layout for.
1893 * @param user Pointer to the UserCtx structure (optional, used to store layout values).
1894 *
1895 * @return PetscErrorCode 0 on success, non-zero on failure.
1896 */
1897PetscErrorCode SetDMDAProcLayout(DM dm, UserCtx *user)
1898{
1899 PetscErrorCode ierr;
1900 PetscMPIInt size, rank;
1901 PetscInt px = PETSC_DECIDE, py = PETSC_DECIDE, pz = PETSC_DECIDE;
1902 PetscBool px_set = PETSC_FALSE, py_set = PETSC_FALSE, pz_set = PETSC_FALSE;
1903 SimCtx *simCtx = user->simCtx;
1904
1905 // Set no.of processors in direction 1
1906 if(simCtx->da_procs_x) {
1907 px_set = PETSC_TRUE;
1908 px = simCtx->da_procs_x;
1909 }
1910 // Set no.of processors in direction 2
1911 if(simCtx->da_procs_y) {
1912 py_set = PETSC_TRUE;
1913 py = simCtx->da_procs_y;
1914 }
1915 // Set no.of processors in direction 1
1916 if(simCtx->da_procs_z) {
1917 pz_set = PETSC_TRUE;
1918 pz = simCtx->da_procs_z;
1919 }
1920
1921 PetscFunctionBeginUser;
1922 PROFILE_FUNCTION_BEGIN;
1923 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size); CHKERRQ(ierr);
1924 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank); CHKERRQ(ierr);
1925 LOG_ALLOW(GLOBAL, LOG_INFO, "Rank %d: Configuring DMDA processor layout for %d total processes.\n", rank, size);
1926
1927 // --- Validate User Input (Optional but Recommended) ---
1928 // Check if specified processor counts multiply to the total MPI size
1929 if (px_set && py_set && pz_set) {
1930 if (px * py * pz != size) {
1931 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_INCOMP,
1932 "Specified processor layout %d x %d x %d = %d does not match MPI size %d",
1933 px, py, pz, px * py * pz, size);
1934 }
1935 LOG_ALLOW(GLOBAL, LOG_INFO, "Using specified processor layout: %d x %d x %d\n", px, py, pz);
1936 } else if (px_set || py_set || pz_set) {
1937 // If only some are set, PETSC_DECIDE will be used for others
1938 LOG_ALLOW(GLOBAL, LOG_INFO, "Using partially specified processor layout: %d x %d x %d (PETSC_DECIDE for unspecified)\n", px, py, pz);
1939 } else {
1940 LOG_ALLOW(GLOBAL, LOG_INFO, "Using fully automatic processor layout (PETSC_DECIDE x PETSC_DECIDE x PETSC_DECIDE)\n");
1941 }
1942 // Additional checks: Ensure px, py, pz are positive if set
1943 if ((px_set && px <= 0) || (py_set && py <= 0) || (pz_set && pz <= 0)) {
1944 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Specified processor counts must be positive.");
1945 }
1946
1947
1948 // --- Apply the layout to the DMDA ---
1949 ierr = DMDASetNumProcs(dm, px, py, pz); CHKERRQ(ierr);
1950 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Rank %d: DMDASetNumProcs called with px=%d, py=%d, pz=%d.\n", rank, px, py, pz);
1951
1952 // --- Store the values in UserCtx (Optional) ---
1953 // Note: If PETSC_DECIDE was used, PETSc calculates the actual values during DMSetUp.
1954 // We store the *requested* values here. To get the *actual* values used,
1955 // you would need to call DMDAGetInfo after DMSetUp.
1956 /*
1957 if (user) {
1958 user->procs_x = px;
1959 user->procs_y = py;
1960 user->procs_z = pz;
1961 }
1962 */
1963 PROFILE_FUNCTION_END;
1964 PetscFunctionReturn(0);
1965}
1966
1967#undef __FUNCT__
1968#define __FUNCT__ "SetupDomainRankInfo"
1969/**
1970 * @brief Sets up the full rank communication infrastructure for all blocks.
1971 *
1972 * This function orchestrates the setup of the parallel communication map. It uses the
1973 * existing helper functions to sequentially gather and broadcast the bounding box
1974 * information for each computational block. All ranks participate in each step,
1975 * assembling the final, unified multi-block list in parallel.
1976 *
1977 * @param[in,out] simCtx The master simulation context. After this call,
1978 * simCtx->bboxlist and the relevant user->neighbors fields
1979 * will be populated on all ranks.
1980 */
1981PetscErrorCode SetupDomainRankInfo(SimCtx *simCtx)
1982{
1983 PetscErrorCode ierr;
1984 PetscInt nblk = simCtx->block_number;
1985 PetscInt size = simCtx->size;
1986 BoundingBox *final_bboxlist = NULL;
1987
1988 PetscFunctionBeginUser;
1989 PROFILE_FUNCTION_BEGIN;
1990
1991 LOG_ALLOW(GLOBAL, LOG_INFO, "Starting full rank communication setup for %d block(s).\n", nblk);
1992
1993 UserCtx *user_finest = simCtx->usermg.mgctx[simCtx->usermg.mglevels - 1].user;
1994
1995 // --- Step 1: Compute neighbor ranks (unchanged) ---
1996 for (int bi = 0; bi < nblk; bi++) {
1997 ierr = ComputeAndStoreNeighborRanks(&user_finest[bi]); CHKERRQ(ierr);
1998 }
1999 LOG_ALLOW(GLOBAL, LOG_INFO, "Neighbor ranks computed and stored for all blocks.\n");
2000
2001 // --- Step 2: Allocate the final, unified list on ALL ranks ---
2002 // Every rank will build this list in parallel.
2003 ierr = PetscMalloc1(size * nblk, &final_bboxlist); CHKERRQ(ierr);
2004
2005 // --- Step 3: Loop through each block, gather then broadcast its bbox list ---
2006 for (int bi = 0; bi < nblk; bi++) {
2007 // This is a temporary pointer for the current block's list.
2008 BoundingBox *block_bboxlist = NULL;
2009
2010 LOG_ALLOW(GLOBAL, LOG_INFO, "Processing bounding boxes for block %d...\n", bi);
2011
2012 // A) GATHER: On rank 0, block_bboxlist is allocated and filled. On others, it's NULL.
2013 ierr = GatherAllBoundingBoxes(&user_finest[bi], &block_bboxlist); CHKERRQ(ierr);
2014 LOG_ALLOW(GLOBAL, LOG_DEBUG, " -> Gather complete for block %d.\n", bi);
2015
2016 // B) BROADCAST: On non-root ranks, block_bboxlist is allocated. Then, the data
2017 // from rank 0 is broadcast to all ranks. After this call, ALL ranks have
2018 // an identical, complete copy of the bounding boxes for the current block.
2019 ierr = BroadcastAllBoundingBoxes(&user_finest[bi], &block_bboxlist); CHKERRQ(ierr);
2020 LOG_ALLOW(GLOBAL, LOG_DEBUG, " -> Broadcast complete for block %d.\n", bi);
2021
2022 // C) ASSEMBLE: Every rank now copies the data for this block into the
2023 // correct segment of its final, unified list.
2024 for (int r = 0; r < size; r++) {
2025 // The layout is [r0b0, r1b0, ..., r(size-1)b0, r0b1, r1b1, ...]
2026 final_bboxlist[bi * size + r] = block_bboxlist[r];
2027 }
2028 LOG_ALLOW(GLOBAL, LOG_DEBUG, " -> Assembly into final list complete for block %d.\n", bi);
2029
2030 // D) CLEANUP: Free the temporary list for this block on ALL ranks before the next iteration.
2031 // Your helper functions use malloc, so we must use free.
2032 free(block_bboxlist);
2033 }
2034
2035 // --- Step 4: Assign the final pointer and run the last setup step ---
2036 simCtx->bboxlist = final_bboxlist;
2037 LOG_ALLOW(GLOBAL, LOG_INFO, "Final unified bboxlist created on all ranks and stored in SimCtx.\n");
2038
2039 ierr = SetupDomainCellDecompositionMap(&user_finest[0]); CHKERRQ(ierr);
2040 LOG_ALLOW(GLOBAL, LOG_INFO, "Domain Cell Composition set and broadcasted.\n");
2041
2042 LOG_ALLOW(GLOBAL, LOG_INFO, "SetupDomainRankInfo: Completed successfully.\n");
2043
2044 PROFILE_FUNCTION_END;
2045 PetscFunctionReturn(0);
2046}
2047
2048#undef __FUNCT__
2049#define __FUNCT__ "Contra2Cart"
2050/**
2051 * @brief Reconstructs Cartesian velocity (Ucat) at cell centers from contravariant
2052 * velocity (Ucont) defined on cell faces.
2053 *
2054 * This function performs the transformation from a contravariant velocity representation
2055 * (which is natural on a curvilinear grid) to a Cartesian (x,y,z) representation.
2056 * For each interior computational cell owned by the rank, it performs the following:
2057 *
2058 * 1. It averages the contravariant velocity components (U¹, U², U³) from the
2059 * surrounding faces to get an estimate of the contravariant velocity at the cell center.
2060 * 2. It averages the metric vectors (Csi, Eta, Zet) from the surrounding faces
2061 * to get an estimate of the metric tensor at the cell center. This tensor forms
2062 * the transformation matrix.
2063 * 3. It solves the linear system `[MetricTensor] * [ucat] = [ucont]` for the
2064 * Cartesian velocity vector `ucat = (u,v,w)` using Cramer's rule.
2065 * 4. The computed Cartesian velocity is stored in the global `user->Ucat` vector.
2066 *
2067 * The function operates on local, ghosted versions of the input vectors (`user->lUcont`,
2068 * `user->lCsi`, etc.) to ensure stencils are valid across processor boundaries.
2069 *
2070 * @param[in,out] user Pointer to the UserCtx structure. The function reads from
2071 * `user->lUcont`, `user->lCsi`, `user->lEta`, `user->lZet`, `user->lNvert`
2072 * and writes to the global `user->Ucat` vector.
2073 *
2074 * @return PetscErrorCode 0 on success.
2075 *
2076 * @note
2077 * - This function should be called AFTER `user->lUcont` and all local metric vectors
2078 * (`user->lCsi`, etc.) have been populated with up-to-date ghost values via `UpdateLocalGhosts`.
2079 * - It only computes `Ucat` for interior cells (not on physical boundaries) and for
2080 * cells not marked as solid/blanked by `user->lNvert`.
2081 * - The caller is responsible for subsequently applying boundary conditions to `user->Ucat`
2082 * and calling `UpdateLocalGhosts(user, "Ucat")` to populate `user->lUcat`.
2083 */
2084PetscErrorCode Contra2Cart(UserCtx *user)
2085{
2086 PetscErrorCode ierr;
2087 DMDALocalInfo info;
2088 Cmpnts ***lcsi_arr, ***leta_arr, ***lzet_arr; // Local metric arrays
2089 Cmpnts ***lucont_arr; // Local contravariant velocity array
2090 Cmpnts ***gucat_arr; // Global Cartesian velocity array
2091 PetscReal ***lnvert_arr; // Local Nvert array
2092 PetscReal ***laj_arr; // Local Jacobian Determinant inverse array
2093
2094 PetscFunctionBeginUser;
2095 PROFILE_FUNCTION_BEGIN;
2096 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Starting Contravariant-to-Cartesian velocity transformation.\n");
2097
2098 // --- 1. Get DMDA Info and Check for Valid Inputs ---
2099 // All inputs (lUcont, lCsi, etc.) and outputs (Ucat) are on DMs from the UserCtx.
2100 // We get local info from fda, which governs the layout of most arrays here.
2101 ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
2102 if (!user->lUcont || !user->lCsi || !user->lEta || !user->lZet || !user->lNvert || !user->Ucat) {
2103 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Contra2Cart requires lUcont, lCsi/Eta/Zet, lNvert, and Ucat to be non-NULL.");
2104 }
2105
2106
2107 // --- 2. Get Read-Only Array Access to Local Input Vectors (with ghosts) ---
2108 ierr = DMDAVecGetArrayRead(user->fda, user->lUcont, &lucont_arr); CHKERRQ(ierr);
2109 ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, &lcsi_arr); CHKERRQ(ierr);
2110 ierr = DMDAVecGetArrayRead(user->fda, user->lEta, &leta_arr); CHKERRQ(ierr);
2111 ierr = DMDAVecGetArrayRead(user->fda, user->lZet, &lzet_arr); CHKERRQ(ierr);
2112 ierr = DMDAVecGetArrayRead(user->da, user->lNvert, &lnvert_arr); CHKERRQ(ierr);
2113 ierr = DMDAVecGetArrayRead(user->da, user->lAj, &laj_arr); CHKERRQ(ierr);
2114
2115 // --- 3. Get Write-Only Array Access to the Global Output Vector ---
2116 // We compute for local owned cells and write into the global vector.
2117 // PETSc handles mapping the global indices to the correct local memory locations.
2118 ierr = DMDAVecGetArray(user->fda, user->Ucat, &gucat_arr); CHKERRQ(ierr);
2119
2120
2121 // --- 4. Define Loop Bounds for INTERIOR Cells ---
2122 // We use adjusted bounds to avoid calculating Ucat on the physical domain boundaries,
2123 // as these are typically set explicitly by boundary condition functions.
2124 // The stencils use indices like i-1, j-1, k-1, so we must start loops at least at index 1.
2125 PetscInt i_start = (info.xs == 0) ? info.xs + 1 : info.xs;
2126 PetscInt i_end = (info.xs + info.xm == info.mx) ? info.xs + info.xm - 1 : info.xs + info.xm;
2127
2128 PetscInt j_start = (info.ys == 0) ? info.ys + 1 : info.ys;
2129 PetscInt j_end = (info.ys + info.ym == info.my) ? info.ys + info.ym - 1 : info.ys + info.ym;
2130
2131 PetscInt k_start = (info.zs == 0) ? info.zs + 1 : info.zs;
2132 PetscInt k_end = (info.zs + info.zm == info.mz) ? info.zs + info.zm - 1 : info.zs + info.zm;
2133
2134 // --- 5. Main Computation Loop ---
2135 // Loops over the GLOBAL indices of interior cells owned by this rank.
2136 for (PetscInt k_cell = k_start; k_cell < k_end; ++k_cell) {
2137 for (PetscInt j_cell = j_start; j_cell < j_end; ++j_cell) {
2138 for (PetscInt i_cell = i_start; i_cell < i_end; ++i_cell) {
2139
2140 // Check if the cell is a fluid cell (not solid/blanked)
2141 // if (lnvert_arr[k_cell][j_cell][i_cell] > 0.1) continue; // Skip solid/blanked cells
2142
2143 // Transformation matrix [mat] is the metric tensor at the cell center,
2144 // estimated by averaging metrics from adjacent faces.
2145 PetscReal mat[3][3];
2146
2147 // PetscReal aj_center = laj_arr[k_cell+1][j_cell+1][i_cell+1];
2148
2149 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;
2150 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;
2151 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;
2152
2153 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;
2154 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;
2155 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;
2156
2157 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;
2158 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;
2159 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;
2160
2161 // Contravariant velocity vector `q` at the cell center,
2162 // estimated by averaging face-based contravariant velocities.
2163 PetscReal q[3];
2164 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
2165 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
2166 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
2167
2168 // Solve the 3x3 system `mat * ucat = q` using Cramer's rule.
2169 PetscReal det = mat[0][0] * (mat[1][1] * mat[2][2] - mat[1][2] * mat[2][1]) -
2170 mat[0][1] * (mat[1][0] * mat[2][2] - mat[1][2] * mat[2][0]) +
2171 mat[0][2] * (mat[1][0] * mat[2][1] - mat[1][1] * mat[2][0]);
2172
2173 if (PetscAbsReal(det) < 1.0e-18) {
2174 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);
2175 }
2176
2177 PetscReal det_inv = 1.0 / det;
2178
2179 PetscReal det0 = q[0] * (mat[1][1] * mat[2][2] - mat[1][2] * mat[2][1]) -
2180 q[1] * (mat[0][1] * mat[2][2] - mat[0][2] * mat[2][1]) +
2181 q[2] * (mat[0][1] * mat[1][2] - mat[0][2] * mat[1][1]);
2182
2183 PetscReal det1 = -q[0] * (mat[1][0] * mat[2][2] - mat[1][2] * mat[2][0]) +
2184 q[1] * (mat[0][0] * mat[2][2] - mat[0][2] * mat[2][0]) -
2185 q[2] * (mat[0][0] * mat[1][2] - mat[0][2] * mat[1][0]);
2186
2187 PetscReal det2 = q[0] * (mat[1][0] * mat[2][1] - mat[1][1] * mat[2][0]) -
2188 q[1] * (mat[0][0] * mat[2][1] - mat[0][1] * mat[2][0]) +
2189 q[2] * (mat[0][0] * mat[1][1] - mat[0][1] * mat[1][0]);
2190
2191 // Store computed Cartesian velocity in the GLOBAL Ucat array at the
2192 // array index corresponding to the cell's origin node.
2193 gucat_arr[k_cell][j_cell][i_cell].x = det0 * det_inv;
2194 gucat_arr[k_cell][j_cell][i_cell].y = det1 * det_inv;
2195 gucat_arr[k_cell][j_cell][i_cell].z = det2 * det_inv;
2196 }
2197 }
2198 }
2199
2200 // --- 6. Restore Array Access ---
2201 ierr = DMDAVecRestoreArrayRead(user->fda, user->lUcont, &lucont_arr); CHKERRQ(ierr);
2202 ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, &lcsi_arr); CHKERRQ(ierr);
2203 ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, &leta_arr); CHKERRQ(ierr);
2204 ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, &lzet_arr); CHKERRQ(ierr);
2205 ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, &lnvert_arr); CHKERRQ(ierr);
2206 ierr = DMDAVecRestoreArrayRead(user->da, user->lAj, &laj_arr); CHKERRQ(ierr);
2207 ierr = DMDAVecRestoreArray(user->fda, user->Ucat, &gucat_arr); CHKERRQ(ierr);
2208
2209 LOG_ALLOW(GLOBAL, LOG_INFO, "Completed Contravariant-to-Cartesian velocity transformation. \n");
2210 PROFILE_FUNCTION_END;
2211 PetscFunctionReturn(0);
2212}
2213
2214#undef __FUNCT__
2215#define __FUNCT__ "SetupDomainCellDecompositionMap"
2216/**
2217 * @brief Creates and distributes a map of the domain's cell decomposition to all ranks.
2218 * @ingroup DomainInfo
2219 *
2220 * This function is a critical part of the simulation setup. It determines the global
2221 * cell ownership for each MPI rank and makes this information available to all
2222 * other ranks. This "decomposition map" is essential for the robust "Walk and Handoff"
2223 * particle migration strategy, allowing any rank to quickly identify the owner of a
2224 * target cell.
2225 *
2226 * The process involves:
2227 * 1. Each rank gets its own node ownership information from the DMDA.
2228 * 2. It converts this node information into cell ownership ranges using the
2229 * `GetOwnedCellRange` helper function.
2230 * 3. It participates in an `MPI_Allgather` collective operation to build a complete
2231 * array (`user->RankCellInfoMap`) containing the ownership information for every rank.
2232 *
2233 * This function should be called once during initialization after the primary DMDA
2234 * (user->da) has been set up.
2235 *
2236 * @param[in,out] user Pointer to the UserCtx structure. The function will allocate and
2237 * populate `user->RankCellInfoMap` and set `user->num_ranks`.
2238 *
2239 * @return PetscErrorCode 0 on success, or a non-zero PETSc error code on failure.
2240 * Errors can occur if input pointers are NULL or if MPI communication fails.
2241 */
2242PetscErrorCode SetupDomainCellDecompositionMap(UserCtx *user)
2243{
2244 PetscErrorCode ierr;
2245 DMDALocalInfo local_node_info;
2246 RankCellInfo my_cell_info;
2247 PetscMPIInt rank, size;
2248
2249 PetscFunctionBeginUser;
2250 PROFILE_FUNCTION_BEGIN;
2251
2252 // --- 1. Input Validation and MPI Info ---
2253 if (!user) {
2254 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "UserCtx pointer is NULL in SetupDomainCellDecompositionMap.");
2255 }
2256 if (!user->da) {
2257 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "user->da is not initialized in SetupDomainCellDecompositionMap.");
2258 }
2259
2260 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
2261 ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size); CHKERRQ(ierr);
2262
2263 LOG_ALLOW(GLOBAL, LOG_INFO, "Setting up domain cell decomposition map for %d ranks.\n", size);
2264
2265 // --- 2. Determine Local Cell Ownership ---
2266 // Get the local node ownership information from the primary DMDA.
2267 ierr = DMDAGetLocalInfo(user->da, &local_node_info); CHKERRQ(ierr);
2268
2269 // Use the robust helper function to convert node ownership to cell ownership.
2270 // A cell's index is defined by its origin node.
2271
2272 ierr = GetOwnedCellRange(&local_node_info, 0, &my_cell_info.xs_cell, &my_cell_info.xm_cell); CHKERRQ(ierr);
2273 ierr = GetOwnedCellRange(&local_node_info, 1, &my_cell_info.ys_cell, &my_cell_info.ym_cell); CHKERRQ(ierr);
2274 ierr = GetOwnedCellRange(&local_node_info, 2, &my_cell_info.zs_cell, &my_cell_info.zm_cell); CHKERRQ(ierr);
2275
2276 // Log the calculated local ownership for debugging purposes.
2277 LOG_ALLOW(LOCAL, LOG_DEBUG, "[Rank %d] Owns cells: i[%d, %d), j[%d, %d), k[%d, %d)\n",
2278 rank, my_cell_info.xs_cell, my_cell_info.xs_cell + my_cell_info.xm_cell,
2279 my_cell_info.ys_cell, my_cell_info.ys_cell + my_cell_info.ym_cell,
2280 my_cell_info.zs_cell, my_cell_info.zs_cell + my_cell_info.zm_cell);
2281
2282 // --- 3. Allocate and Distribute the Global Map ---
2283 // Allocate memory for the global map that will hold information from all ranks.
2284 ierr = PetscMalloc1(size, &user->RankCellInfoMap); CHKERRQ(ierr);
2285
2286 // Perform the collective communication to gather the `RankCellInfo` struct from every rank.
2287 // Each rank sends its `my_cell_info` and receives the complete array in `user->RankCellInfoMap`.
2288 // We use MPI_BYTE to ensure portability across different systems and struct padding.
2289 ierr = MPI_Allgather(&my_cell_info, sizeof(RankCellInfo), MPI_BYTE,
2290 user->RankCellInfoMap, sizeof(RankCellInfo), MPI_BYTE,
2291 PETSC_COMM_WORLD); CHKERRQ(ierr);
2292
2293 LOG_ALLOW(GLOBAL, LOG_INFO, "Domain cell decomposition map created and distributed successfully.\n");
2294
2295 PROFILE_FUNCTION_END;
2296 PetscFunctionReturn(0);
2297}
2298
2299#undef __FUNCT__
2300#define __FUNCT__ "BinarySearchInt64"
2301/**
2302 * @brief Performs a binary search for a key in a sorted array of PetscInt64.
2303 *
2304 * This is a standard binary search algorithm implemented as a PETSc-style helper function.
2305 * It efficiently determines if a given `key` exists within a `sorted` array.
2306 *
2307 * @param[in] n The number of elements in the array.
2308 * @param[in] arr A pointer to the sorted array of PetscInt64 values to be searched.
2309 * @param[in] key The PetscInt64 value to search for.
2310 * @param[out] found A pointer to a PetscBool that will be set to PETSC_TRUE if the key
2311 * is found, and PETSC_FALSE otherwise.
2312 *
2313 * @return PetscErrorCode 0 on success, or a non-zero PETSc error code on failure.
2314 *
2315 * @note The input array `arr` **must** be sorted in ascending order for the algorithm
2316 * to work correctly.
2317 */
2318PetscErrorCode BinarySearchInt64(PetscInt n, const PetscInt64 arr[], PetscInt64 key, PetscBool *found)
2319{
2320 PetscInt low = 0, high = n - 1;
2321
2322 PetscFunctionBeginUser;
2323 PROFILE_FUNCTION_BEGIN;
2324
2325 // --- 1. Input Validation ---
2326 if (!found) {
2327 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Output pointer 'found' is NULL in PetscBinarySearchInt64.");
2328 }
2329 if (n > 0 && !arr) {
2330 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Input array 'arr' is NULL for n > 0.");
2331 }
2332
2333 // Initialize output
2334 *found = PETSC_FALSE;
2335
2336 // --- 2. Binary Search Algorithm ---
2337 while (low <= high) {
2338 // Use this form to prevent potential integer overflow on very large arrays
2339 PetscInt mid = low + (high - low) / 2;
2340
2341 if (arr[mid] == key) {
2342 *found = PETSC_TRUE; // Key found!
2343 break; // Exit the loop
2344 }
2345
2346 if (arr[mid] < key) {
2347 low = mid + 1; // Search in the right half
2348 } else {
2349 high = mid - 1; // Search in the left half
2350 }
2351 }
2352
2353 PROFILE_FUNCTION_END;
2354 PetscFunctionReturn(0);
2355}
2356
2357
2358static PetscInt Gidx(PetscInt i, PetscInt j, PetscInt k, UserCtx *user)
2359{
2360 PetscInt nidx;
2361 DMDALocalInfo info = user->info;
2362
2363 PetscInt mx = info.mx, my = info.my;
2364
2365 AO ao;
2366 DMDAGetAO(user->da, &ao);
2367 nidx=i+j*mx+k*mx*my;
2368
2369 AOApplicationToPetsc(ao,1,&nidx);
2370
2371 return (nidx);
2372}
2373
2374
2375#undef __FUNCT__
2376#define __FUNCT__ "ComputeDivergence"
2377PetscErrorCode ComputeDivergence(UserCtx *user )
2378{
2379 DM da = user->da, fda = user->fda;
2380 DMDALocalInfo info = user->info;
2381
2382 PetscInt ti = user->simCtx->step;
2383
2384 PetscInt xs = info.xs, xe = info.xs + info.xm;
2385 PetscInt ys = info.ys, ye = info.ys + info.ym;
2386 PetscInt zs = info.zs, ze = info.zs + info.zm;
2387 PetscInt mx = info.mx, my = info.my, mz = info.mz;
2388
2389 PetscInt lxs, lys, lzs, lxe, lye, lze;
2390 PetscInt i, j, k;
2391
2392 Vec Div;
2393 PetscReal ***div, ***aj, ***nvert,***p;
2394 Cmpnts ***ucont;
2395 PetscReal maxdiv;
2396
2397 lxs = xs; lxe = xe;
2398 lys = ys; lye = ye;
2399 lzs = zs; lze = ze;
2400
2401 if (xs==0) lxs = xs+1;
2402 if (ys==0) lys = ys+1;
2403 if (zs==0) lzs = zs+1;
2404
2405 if (xe==mx) lxe = xe-1;
2406 if (ye==my) lye = ye-1;
2407 if (ze==mz) lze = ze-1;
2408
2409 PetscFunctionBeginUser;
2410 PROFILE_FUNCTION_BEGIN;
2411
2412 DMDAVecGetArray(fda,user->lUcont, &ucont);
2413 DMDAVecGetArray(da, user->lAj, &aj);
2414 VecDuplicate(user->P, &Div);
2415 DMDAVecGetArray(da, Div, &div);
2416 DMDAVecGetArray(da, user->lNvert, &nvert);
2417 DMDAVecGetArray(da, user->P, &p);
2418 for (k=lzs; k<lze; k++) {
2419 for (j=lys; j<lye; j++){
2420 for (i=lxs; i<lxe; i++) {
2421 if (k==10 && j==10 && i==1){
2422 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]);
2423 }
2424
2425 if (k==10 && j==10 && i==mx-3)
2426 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]);
2427 }
2428 }
2429 }
2430 DMDAVecRestoreArray(da, user->P, &p);
2431
2432
2433 for (k=lzs; k<lze; k++) {
2434 for (j=lys; j<lye; j++) {
2435 for (i=lxs; i<lxe; i++) {
2436 maxdiv = fabs((ucont[k][j][i].x - ucont[k][j][i-1].x +
2437 ucont[k][j][i].y - ucont[k][j-1][i].y +
2438 ucont[k][j][i].z - ucont[k-1][j][i].z)*aj[k][j][i]);
2439 if (nvert[k][j][i] + nvert[k+1][j][i] + nvert[k-1][j][i] +
2440 nvert[k][j+1][i] + nvert[k][j-1][i] +
2441 nvert[k][j][i+1] + nvert[k][j][i-1] > 0.1) maxdiv = 0.;
2442 div[k][j][i] = maxdiv;
2443
2444 }
2445 }
2446 }
2447
2448 if (zs==0) {
2449 k=0;
2450 for (j=ys; j<ye; j++) {
2451 for (i=xs; i<xe; i++) {
2452 div[k][j][i] = 0.;
2453 }
2454 }
2455 }
2456
2457 if (ze == mz) {
2458 k=mz-1;
2459 for (j=ys; j<ye; j++) {
2460 for (i=xs; i<xe; i++) {
2461 div[k][j][i] = 0.;
2462 }
2463 }
2464 }
2465
2466 if (xs==0) {
2467 i=0;
2468 for (k=zs; k<ze; k++) {
2469 for (j=ys; j<ye; j++) {
2470 div[k][j][i] = 0.;
2471 }
2472 }
2473 }
2474
2475 if (xe==mx) {
2476 i=mx-1;
2477 for (k=zs; k<ze; k++) {
2478 for (j=ys; j<ye; j++) {
2479 div[k][j][i] = 0;
2480 }
2481 }
2482 }
2483
2484 if (ys==0) {
2485 j=0;
2486 for (k=zs; k<ze; k++) {
2487 for (i=xs; i<xe; i++) {
2488 div[k][j][i] = 0.;
2489 }
2490 }
2491 }
2492
2493 if (ye==my) {
2494 j=my-1;
2495 for (k=zs; k<ze; k++) {
2496 for (i=xs; i<xe; i++) {
2497 div[k][j][i] = 0.;
2498 }
2499 }
2500 }
2501 DMDAVecRestoreArray(da, Div, &div);
2502 PetscInt MaxFlatIndex;
2503
2504 VecMax(Div, &MaxFlatIndex, &maxdiv);
2505
2506 LOG_ALLOW(GLOBAL,LOG_INFO,"[Step %d]] The Maximum Divergence is %e at flat index %d.\n",ti,maxdiv,MaxFlatIndex);
2507
2508 user->simCtx->MaxDivFlatArg = MaxFlatIndex;
2509 user->simCtx->MaxDiv = maxdiv;
2510
2511 for (k=zs; k<ze; k++) {
2512 for (j=ys; j<ye; j++) {
2513 for (i=xs; i<xe; i++) {
2514 if (Gidx(i,j,k,user) == MaxFlatIndex) {
2515 LOG_ALLOW(GLOBAL,LOG_INFO,"[Step %d] The Maximum Divergence(%e) is at location [%d][%d][%d]. \n", ti, maxdiv,k,j,i);
2516 user->simCtx->MaxDivz = k;
2517 user->simCtx->MaxDivy = j;
2518 user->simCtx->MaxDivx = i;
2519 }
2520 }
2521 }
2522 }
2523
2524
2525 DMDAVecRestoreArray(da, user->lNvert, &nvert);
2526 DMDAVecRestoreArray(fda, user->lUcont, &ucont);
2527 DMDAVecRestoreArray(da, user->lAj, &aj);
2528 VecDestroy(&Div);
2529
2530 PROFILE_FUNCTION_END;
2531 PetscFunctionReturn(0);
2532}
2533
2534#undef __FUNCT__
2535#define __FUNCT__ "InitializeRandomGenerators"
2536
2537/**
2538 * @brief Initializes random number generators for assigning particle properties.
2539 *
2540 * This function creates and configures separate PETSc random number generators for the x, y, and z coordinates.
2541 *
2542 * @param[in,out] user Pointer to the UserCtx structure containing simulation context.
2543 * @param[out] randx Pointer to store the RNG for the x-coordinate.
2544 * @param[out] randy Pointer to store the RNG for the y-coordinate.
2545 * @param[out] randz Pointer to store the RNG for the z-coordinate.
2546 *
2547 * @return PetscErrorCode Returns 0 on success, non-zero on failure.
2548 */
2549PetscErrorCode InitializeRandomGenerators(UserCtx* user, PetscRandom *randx, PetscRandom *randy, PetscRandom *randz) {
2550 PetscErrorCode ierr; // Error code for PETSc functions
2551 PetscMPIInt rank;
2552 PetscFunctionBeginUser;
2553 PROFILE_FUNCTION_BEGIN;
2554 MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
2555
2556 // Initialize RNG for x-coordinate
2557 ierr = PetscRandomCreate(PETSC_COMM_SELF, randx); CHKERRQ(ierr);
2558 ierr = PetscRandomSetType((*randx), PETSCRAND48); CHKERRQ(ierr);
2559 ierr = PetscRandomSetInterval(*randx, user->bbox.min_coords.x, user->bbox.max_coords.x); CHKERRQ(ierr);
2560 ierr = PetscRandomSetSeed(*randx, rank + 12345); CHKERRQ(ierr); // Unique seed per rank
2561 ierr = PetscRandomSeed(*randx); CHKERRQ(ierr);
2562 LOG_ALLOW_SYNC(LOCAL,LOG_VERBOSE, "[Rank %d]Initialized RNG for X-axis.\n",rank);
2563
2564 // Initialize RNG for y-coordinate
2565 ierr = PetscRandomCreate(PETSC_COMM_SELF, randy); CHKERRQ(ierr);
2566 ierr = PetscRandomSetType((*randy), PETSCRAND48); CHKERRQ(ierr);
2567 ierr = PetscRandomSetInterval(*randy, user->bbox.min_coords.y, user->bbox.max_coords.y); CHKERRQ(ierr);
2568 ierr = PetscRandomSetSeed(*randy, rank + 67890); CHKERRQ(ierr); // Unique seed per rank
2569 ierr = PetscRandomSeed(*randy); CHKERRQ(ierr);
2570 LOG_ALLOW_SYNC(LOCAL,LOG_VERBOSE, "[Rank %d]Initialized RNG for Y-axis.\n",rank);
2571
2572 // Initialize RNG for z-coordinate
2573 ierr = PetscRandomCreate(PETSC_COMM_SELF, randz); CHKERRQ(ierr);
2574 ierr = PetscRandomSetType((*randz), PETSCRAND48); CHKERRQ(ierr);
2575 ierr = PetscRandomSetInterval(*randz, user->bbox.min_coords.z, user->bbox.max_coords.z); CHKERRQ(ierr);
2576 ierr = PetscRandomSetSeed(*randz, rank + 54321); CHKERRQ(ierr); // Unique seed per rank
2577 ierr = PetscRandomSeed(*randz); CHKERRQ(ierr);
2578 LOG_ALLOW_SYNC(LOCAL,LOG_VERBOSE, "[Rank %d]Initialized RNG for Z-axis.\n",rank);
2579
2580 PROFILE_FUNCTION_END;
2581 PetscFunctionReturn(0);
2582}
2583
2584#undef __FUNCT__
2585#define __FUNCT__ "InitializeLogicalSpaceRNGs"
2586/**
2587 * @brief Initializes random number generators for logical space operations [0.0, 1.0).
2588 *
2589 * This function creates and configures three separate PETSc random number generators,
2590 * one for each logical dimension (i, j, k or xi, eta, zeta equivalent).
2591 * Each RNG is configured to produce uniformly distributed real numbers in the interval [0.0, 1.0).
2592 * These are typically used for selecting owned cells or generating intra-cell logical coordinates.
2593 *
2594 * @param[out] rand_logic_i Pointer to store the RNG for the i-logical dimension.
2595 * @param[out] rand_logic_j Pointer to store the RNG for the j-logical dimension.
2596 * @param[out] rand_logic_k Pointer to store the RNG for the k-logical dimension.
2597 *
2598 * @return PetscErrorCode Returns 0 on success, non-zero on failure.
2599 */
2600PetscErrorCode InitializeLogicalSpaceRNGs(PetscRandom *rand_logic_i, PetscRandom *rand_logic_j, PetscRandom *rand_logic_k) {
2601 PetscErrorCode ierr;
2602 PetscMPIInt rank;
2603 PetscFunctionBeginUser;
2604
2605 PROFILE_FUNCTION_BEGIN;
2606
2607 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
2608
2609 // --- RNG for i-logical dimension ---
2610 ierr = PetscRandomCreate(PETSC_COMM_SELF, rand_logic_i); CHKERRQ(ierr);
2611 ierr = PetscRandomSetType((*rand_logic_i), PETSCRAND48); CHKERRQ(ierr);
2612 ierr = PetscRandomSetInterval(*rand_logic_i, 0.0, 1.0); CHKERRQ(ierr); // Key change: [0,1)
2613 ierr = PetscRandomSetSeed(*rand_logic_i, rank + 202401); CHKERRQ(ierr); // Unique seed
2614 ierr = PetscRandomSeed(*rand_logic_i); CHKERRQ(ierr);
2615 LOG_ALLOW(LOCAL,LOG_VERBOSE, "[Rank %d] Initialized RNG for i-logical dimension [0,1).\n",rank);
2616
2617 // --- RNG for j-logical dimension ---
2618 ierr = PetscRandomCreate(PETSC_COMM_SELF, rand_logic_j); CHKERRQ(ierr);
2619 ierr = PetscRandomSetType((*rand_logic_j), PETSCRAND48); CHKERRQ(ierr);
2620 ierr = PetscRandomSetInterval(*rand_logic_j, 0.0, 1.0); CHKERRQ(ierr); // Key change: [0,1)
2621 ierr = PetscRandomSetSeed(*rand_logic_j, rank + 202402); CHKERRQ(ierr);
2622 ierr = PetscRandomSeed(*rand_logic_j); CHKERRQ(ierr);
2623 LOG_ALLOW(LOCAL,LOG_VERBOSE, "[Rank %d] Initialized RNG for j-logical dimension [0,1).\n",rank);
2624
2625 // --- RNG for k-logical dimension ---
2626 ierr = PetscRandomCreate(PETSC_COMM_SELF, rand_logic_k); CHKERRQ(ierr);
2627 ierr = PetscRandomSetType((*rand_logic_k), PETSCRAND48); CHKERRQ(ierr);
2628 ierr = PetscRandomSetInterval(*rand_logic_k, 0.0, 1.0); CHKERRQ(ierr); // Key change: [0,1)
2629 ierr = PetscRandomSetSeed(*rand_logic_k, rank + 202403); CHKERRQ(ierr);
2630 ierr = PetscRandomSeed(*rand_logic_k); CHKERRQ(ierr);
2631 LOG_ALLOW(LOCAL,LOG_VERBOSE, "[Rank %d] Initialized RNG for k-logical dimension [0,1).\n",rank);
2632
2633
2634 PROFILE_FUNCTION_END;
2635 PetscFunctionReturn(0);
2636}
2637
2638/////////////// DERIVATIVE CALCULATION HELPERS ///////////////
2639
2640/**
2641 * @brief Transforms derivatives from computational space to physical space using the chain rule.
2642 */
2643static void TransformDerivativesToPhysical(PetscReal jacobian, Cmpnts csi_metrics, Cmpnts eta_metrics, Cmpnts zet_metrics,
2644 Cmpnts deriv_csi, Cmpnts deriv_eta, Cmpnts deriv_zet,
2645 Cmpnts *dudx, Cmpnts *dvdx, Cmpnts *dwdx)
2646{
2647 // Derivatives of the first component (u)
2648 dudx->x = jacobian * (deriv_csi.x * csi_metrics.x + deriv_eta.x * eta_metrics.x + deriv_zet.x * zet_metrics.x);
2649 dudx->y = jacobian * (deriv_csi.x * csi_metrics.y + deriv_eta.x * eta_metrics.y + deriv_zet.x * zet_metrics.y);
2650 dudx->z = jacobian * (deriv_csi.x * csi_metrics.z + deriv_eta.x * eta_metrics.z + deriv_zet.x * zet_metrics.z);
2651 // Derivatives of the second component (v)
2652 dvdx->x = jacobian * (deriv_csi.y * csi_metrics.x + deriv_eta.y * eta_metrics.x + deriv_zet.y * zet_metrics.x);
2653 dvdx->y = jacobian * (deriv_csi.y * csi_metrics.y + deriv_eta.y * eta_metrics.y + deriv_zet.y * zet_metrics.y);
2654 dvdx->z = jacobian * (deriv_csi.y * csi_metrics.z + deriv_eta.y * eta_metrics.z + deriv_zet.y * zet_metrics.z);
2655 // Derivatives of the third component (w)
2656 dwdx->x = jacobian * (deriv_csi.z * csi_metrics.x + deriv_eta.z * eta_metrics.x + deriv_zet.z * zet_metrics.x);
2657 dwdx->y = jacobian * (deriv_csi.z * csi_metrics.y + deriv_eta.z * eta_metrics.y + deriv_zet.z * zet_metrics.y);
2658 dwdx->z = jacobian * (deriv_csi.z * csi_metrics.z + deriv_eta.z * eta_metrics.z + deriv_zet.z * zet_metrics.z);
2659}
2660
2661#undef __FUNCT__
2662#define __FUNCT__ "ComputeVectorFieldDerivatives"
2663/**
2664 * @brief Computes the derivatives of a cell-centered vector field at a specific grid point.
2665 *
2666 * This function orchestrates the calculation of spatial derivatives. It first computes
2667 * the derivatives in computational space (d/dcsi, d/deta, d/dzet) using a central
2668 * difference scheme and then transforms them into physical space (d/dx, d/dy, d/dz).
2669 *
2670 * @param user The user context for the current computational block.
2671 * @param i, j, k The grid indices of the cell center where derivatives are required.
2672 * @param field_data A 3D array pointer to the raw local data of the vector field (e.g., from lUcat).
2673 * @param dudx Output: A Cmpnts struct storing [du/dx, du/dy, du/dz].
2674 * @param dvdx Output: A Cmpnts struct storing [dv/dx, dv/dy, dv/dz].
2675 * @param dwdx Output: A Cmpnts struct storing [dw/dx, dw/dy, dw/dz].
2676 * @return PetscErrorCode 0 on success.
2677 */
2678PetscErrorCode ComputeVectorFieldDerivatives(UserCtx *user, PetscInt i, PetscInt j, PetscInt k, Cmpnts ***field_data,
2679 Cmpnts *dudx, Cmpnts *dvdx, Cmpnts *dwdx)
2680{
2681 PetscErrorCode ierr;
2682 Cmpnts ***csi, ***eta, ***zet;
2683 PetscReal ***jac;
2684 PetscFunctionBeginUser;
2685
2686 // 1. Get read-only access to the necessary metric data arrays
2687 ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, &csi); CHKERRQ(ierr);
2688 ierr = DMDAVecGetArrayRead(user->fda, user->lEta, &eta); CHKERRQ(ierr);
2689 ierr = DMDAVecGetArrayRead(user->fda, user->lZet, &zet); CHKERRQ(ierr);
2690 ierr = DMDAVecGetArrayRead(user->da, user->lAj, &jac); CHKERRQ(ierr);
2691
2692 // 2. Calculate derivatives in computational space using central differencing
2693 Cmpnts deriv_csi, deriv_eta, deriv_zet;
2694 deriv_csi.x = (field_data[k][j][i+1].x - field_data[k][j][i-1].x) * 0.5;
2695 deriv_csi.y = (field_data[k][j][i+1].y - field_data[k][j][i-1].y) * 0.5;
2696 deriv_csi.z = (field_data[k][j][i+1].z - field_data[k][j][i-1].z) * 0.5;
2697
2698 deriv_eta.x = (field_data[k][j+1][i].x - field_data[k][j-1][i].x) * 0.5;
2699 deriv_eta.y = (field_data[k][j+1][i].y - field_data[k][j-1][i].y) * 0.5;
2700 deriv_eta.z = (field_data[k][j+1][i].z - field_data[k][j-1][i].z) * 0.5;
2701
2702 deriv_zet.x = (field_data[k+1][j][i].x - field_data[k-1][j][i].x) * 0.5;
2703 deriv_zet.y = (field_data[k+1][j][i].y - field_data[k-1][j][i].y) * 0.5;
2704 deriv_zet.z = (field_data[k+1][j][i].z - field_data[k-1][j][i].z) * 0.5;
2705
2706 // 3. Transform derivatives to physical space
2707 TransformDerivativesToPhysical(jac[k][j][i], csi[k][j][i], eta[k][j][i], zet[k][j][i],
2708 deriv_csi, deriv_eta, deriv_zet,
2709 dudx, dvdx, dwdx);
2710
2711 // 4. Restore access to the PETSc data arrays
2712 ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, &csi); CHKERRQ(ierr);
2713 ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, &eta); CHKERRQ(ierr);
2714 ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, &zet); CHKERRQ(ierr);
2715 ierr = DMDAVecRestoreArrayRead(user->da, user->lAj, &jac); CHKERRQ(ierr);
2716
2717 PetscFunctionReturn(0);
2718}
TranslateModernBCsToLegacy
PetscErrorCode TranslateModernBCsToLegacy(UserCtx *user)
Definition Boundaries.c:635
CalculateAllGridMetrics
PetscErrorCode CalculateAllGridMetrics(SimCtx *simCtx)
Orchestrates the calculation of all grid metrics.
Definition Metric.c:1809
DefineAllGridDimensions
PetscErrorCode DefineAllGridDimensions(SimCtx *simCtx)
Orchestrates the parsing and setting of grid dimensions for all blocks.
Definition grid.c:65
BroadcastAllBoundingBoxes
PetscErrorCode BroadcastAllBoundingBoxes(UserCtx *user, BoundingBox **bboxlist)
Broadcasts the bounding box information collected on rank 0 to all other ranks.
Definition grid.c:996
CalculateInletCenter
PetscErrorCode CalculateInletCenter(UserCtx *user)
Calculates the geometric center of the primary inlet face.
Definition grid.c:1054
InitializeAllGridDMs
PetscErrorCode InitializeAllGridDMs(SimCtx *simCtx)
Orchestrates the creation of DMDA objects for every block and multigrid level.
Definition grid.c:247
AssignAllGridCoordinates
PetscErrorCode AssignAllGridCoordinates(SimCtx *simCtx)
Orchestrates the assignment of physical coordinates to all DMDA objects.
Definition grid.c:340
GatherAllBoundingBoxes
PetscErrorCode GatherAllBoundingBoxes(UserCtx *user, BoundingBox **allBBoxes)
Gathers local bounding boxes from all MPI processes to rank 0.
Definition grid.c:928
ParsePostProcessingSettings
PetscErrorCode ParsePostProcessingSettings(SimCtx *simCtx)
Initializes post-processing settings from a config file and command-line overrides.
Definition io.c:2199
ParseScalingInformation
PetscErrorCode ParseScalingInformation(SimCtx *simCtx)
Parses physical scaling parameters from command-line options.
Definition io.c:2334
ParseAllBoundaryConditions
PetscErrorCode ParseAllBoundaryConditions(UserCtx *user, const char *bcs_input_filename)
Parses the boundary conditions file to configure the type, handler, and any associated parameters for...
Definition io.c:399
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:590
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:122
is_function_allowed
PetscBool is_function_allowed(const char *functionName)
Checks if a given function is in the allow-list.
Definition logging.c:157
LOG_ALLOW_SYNC
#define LOG_ALLOW_SYNC(scope, level, fmt,...)
----— DEBUG ---------------------------------------— #define LOG_ALLOW(scope, level,...
Definition logging.h:268
LOCAL
#define LOCAL
Logging scope definitions for controlling message output.
Definition logging.h:46
GLOBAL
#define GLOBAL
Scope for global logging across all processes.
Definition logging.h:47
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:201
print_log_level
PetscErrorCode print_log_level(void)
Prints the current logging level to the console.
Definition logging.c:82
PROFILE_FUNCTION_END
#define PROFILE_FUNCTION_END
Marks the end of a profiled code block.
Definition logging.h:740
LOG
#define LOG(scope, level, fmt,...)
Logging macro for PETSc-based applications with scope control.
Definition logging.h:85
LoadAllowedFunctionsFromFile
PetscErrorCode LoadAllowedFunctionsFromFile(const char filename[], char ***funcsOut, PetscInt *nOut)
Load function names from a text file.
Definition logging.c:566
get_log_level
LogLevel get_log_level()
Retrieves the current logging level from the environment variable LOG_LEVEL.
Definition logging.c:39
ProfilingInitialize
PetscErrorCode ProfilingInitialize(SimCtx *simCtx)
Initializes the custom profiling system using configuration from SimCtx.
Definition logging.c:998
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:32
LOG_WARNING
@ LOG_WARNING
Non-critical issues that warrant attention.
Definition logging.h:30
LOG_DEBUG
@ LOG_DEBUG
Detailed debugging information.
Definition logging.h:33
LOG_VERBOSE
@ LOG_VERBOSE
Extremely detailed logs, typically for development use only.
Definition logging.h:35
PROFILE_FUNCTION_BEGIN
#define PROFILE_FUNCTION_BEGIN
Marks the beginning of a profiled code block (typically a function).
Definition logging.h:731
ComputeVectorFieldDerivatives
PetscErrorCode ComputeVectorFieldDerivatives(UserCtx *user, PetscInt i, PetscInt j, PetscInt k, Cmpnts ***field_data, Cmpnts *dudx, Cmpnts *dvdx, Cmpnts *dwdx)
Computes the derivatives of a cell-centered vector field at a specific grid point.
Definition setup.c:2678
GetOwnedCellRange
PetscErrorCode GetOwnedCellRange(const DMDALocalInfo *info_nodes, PetscInt dim, PetscInt *xs_cell_global_out, PetscInt *xm_cell_local_out)
Gets the global starting index of cells owned by this rank and the number of such cells.
Definition setup.c:1663
SetupDomainRankInfo
PetscErrorCode SetupDomainRankInfo(SimCtx *simCtx)
Sets up the full rank communication infrastructure for all blocks.
Definition setup.c:1981
InitializeRandomGenerators
PetscErrorCode InitializeRandomGenerators(UserCtx *user, PetscRandom *randx, PetscRandom *randy, PetscRandom *randz)
Initializes random number generators for assigning particle properties.
Definition setup.c:2549
Deallocate3DArrayVector
PetscErrorCode Deallocate3DArrayVector(Cmpnts ***array, PetscInt nz, PetscInt ny)
Deallocates a 3D array of Cmpnts structures allocated by Allocate3DArrayVector.
Definition setup.c:1571
SetupGridAndSolvers
PetscErrorCode SetupGridAndSolvers(SimCtx *simCtx)
The main orchestrator for setting up all grid-related components.
Definition setup.c:857
SetupSolverParameters
static PetscErrorCode SetupSolverParameters(SimCtx *simCtx)
Definition setup.c:827
Gidx
static PetscInt Gidx(PetscInt i, PetscInt j, PetscInt k, UserCtx *user)
Definition setup.c:2358
SetupSimulationEnvironment
PetscErrorCode SetupSimulationEnvironment(SimCtx *simCtx)
Verifies and prepares the complete I/O environment for a simulation run.
Definition setup.c:558
CreateAndInitializeAllVectors
PetscErrorCode CreateAndInitializeAllVectors(SimCtx *simCtx)
Creates and initializes all PETSc Vec objects for all fields.
Definition setup.c:897
ComputeAndStoreNeighborRanks
PetscErrorCode ComputeAndStoreNeighborRanks(UserCtx *user)
Computes and stores the Cartesian neighbor ranks for the DMDA decomposition.
Definition setup.c:1770
Contra2Cart
PetscErrorCode Contra2Cart(UserCtx *user)
Reconstructs Cartesian velocity (Ucat) at cell centers from contravariant velocity (Ucont) defined on...
Definition setup.c:2084
PetscMkdirRecursive
static PetscErrorCode PetscMkdirRecursive(const char *path)
Creates a directory path recursively, similar to mkdir -p.
Definition setup.c:491
Allocate3DArrayScalar
PetscErrorCode Allocate3DArrayScalar(PetscReal ****array, PetscInt nz, PetscInt ny, PetscInt nx)
Allocates a 3D array of PetscReal values using PetscCalloc.
Definition setup.c:1411
CreateSimulationContext
PetscErrorCode CreateSimulationContext(int argc, char **argv, SimCtx **p_simCtx)
Allocates and populates the master SimulationContext object.
Definition setup.c:44
SetDMDAProcLayout
PetscErrorCode SetDMDAProcLayout(DM dm, UserCtx *user)
Sets the processor layout for a given DMDA based on PETSc options.
Definition setup.c:1897
InitializeLogicalSpaceRNGs
PetscErrorCode InitializeLogicalSpaceRNGs(PetscRandom *rand_logic_i, PetscRandom *rand_logic_j, PetscRandom *rand_logic_k)
Initializes random number generators for logical space operations [0.0, 1.0).
Definition setup.c:2600
ComputeDivergence
PetscErrorCode ComputeDivergence(UserCtx *user)
Definition setup.c:2377
BinarySearchInt64
PetscErrorCode BinarySearchInt64(PetscInt n, const PetscInt64 arr[], PetscInt64 key, PetscBool *found)
Performs a binary search for a key in a sorted array of PetscInt64.
Definition setup.c:2318
AllocateContextHierarchy
static PetscErrorCode AllocateContextHierarchy(SimCtx *simCtx)
Allocates the memory for the UserMG and UserCtx hierarchy.
Definition setup.c:730
Allocate3DArrayVector
PetscErrorCode Allocate3DArrayVector(Cmpnts ****array, PetscInt nz, PetscInt ny, PetscInt nx)
Allocates a 3D array of Cmpnts structures using PetscCalloc.
Definition setup.c:1513
UpdateLocalGhosts
PetscErrorCode UpdateLocalGhosts(UserCtx *user, const char *fieldName)
Updates the local vector (including ghost points) from its corresponding global vector.
Definition setup.c:1091
SetupBoundaryConditions
PetscErrorCode SetupBoundaryConditions(SimCtx *simCtx)
(Orchestrator) Sets up all boundary conditions for the simulation.
Definition setup.c:1354
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)
Transforms derivatives from computational space to physical space using the chain rule.
Definition setup.c:2643
__FUNCT__
#define __FUNCT__
Definition setup.c:10
SetupDomainCellDecompositionMap
PetscErrorCode SetupDomainCellDecompositionMap(UserCtx *user)
Creates and distributes a map of the domain's cell decomposition to all ranks.
Definition setup.c:2242
Deallocate3DArrayScalar
PetscErrorCode Deallocate3DArrayScalar(PetscReal ***array, PetscInt nz, PetscInt ny)
Deallocates a 3D array of PetscReal values allocated by Allocate3DArrayScalar.
Definition setup.c:1456
RankNeighbors::rank_zm
PetscMPIInt rank_zm
Definition variables.h:179
SimCtx::MHV
PetscInt MHV
Definition variables.h:572
UserCtx::isc
PetscInt isc
Definition variables.h:674
SimCtx::turbine
PetscInt turbine
Definition variables.h:572
SimCtx::fishcyl
PetscInt fishcyl
Definition variables.h:572
SimCtx::clark
PetscInt clark
Definition variables.h:610
SimCtx::implicit_type
PetscInt implicit_type
Definition variables.h:575
SimCtx::movefsi
PetscInt movefsi
Definition variables.h:567
UserCtx::lCent
Vec lCent
Definition variables.h:705
UserCtx::GridSpace
Vec GridSpace
Definition variables.h:705
SimCtx::moveframe
PetscInt moveframe
Definition variables.h:568
UserCtx::P_nodal
Vec P_nodal
Definition variables.h:734
UserCtx::JCsi
Vec JCsi
Definition variables.h:708
UserCtx::KAj
Vec KAj
Definition variables.h:709
SimCtx::TwoD
PetscInt TwoD
Definition variables.h:568
SimCtx::pseudo_periodic
PetscInt pseudo_periodic
Definition variables.h:594
MGCtx::user
UserCtx * user
Definition variables.h:427
SimCtx::fish_c
PetscInt fish_c
Definition variables.h:572
RankCellInfo::ys_cell
PetscInt ys_cell
Definition variables.h:184
SimCtx::dgf_z
PetscInt dgf_z
Definition variables.h:569
UserCtx::JEta
Vec JEta
Definition variables.h:708
SimCtx::poisson_tol
PetscReal poisson_tol
Definition variables.h:579
UserCtx::Zet
Vec Zet
Definition variables.h:705
PostProcessParams::particle_output_prefix
char particle_output_prefix[256]
Definition variables.h:469
RankCellInfo::xs_cell
PetscInt xs_cell
Definition variables.h:184
SimCtx::rank
PetscMPIInt rank
Definition variables.h:541
UserCtx::mglevels
PetscInt mglevels
Definition variables.h:721
SimCtx::fish
PetscInt fish
Definition variables.h:572
SimCtx::LV
PetscInt LV
Definition variables.h:572
SimCtx::angle
PetscReal angle
Definition variables.h:588
RankNeighbors::rank_yp
PetscMPIInt rank_yp
Definition variables.h:178
SimCtx::thin
PetscInt thin
Definition variables.h:568
SimCtx::grid1d
PetscInt grid1d
Definition variables.h:593
SimCtx::eulerianSource
char eulerianSource[64]
Definition variables.h:557
SimCtx::block_number
PetscInt block_number
Definition variables.h:593
UserCtx::lIEta
Vec lIEta
Definition variables.h:707
SimCtx::da_procs_z
PetscInt da_procs_z
Definition variables.h:599
SimCtx::blkpbc
PetscInt blkpbc
Definition variables.h:594
SimCtx::sediment
PetscInt sediment
Definition variables.h:567
UserCtx::lIZet
Vec lIZet
Definition variables.h:707
UserCtx::user_f
UserCtx * user_f
Definition variables.h:722
SimCtx::channelz
PetscInt channelz
Definition variables.h:593
UserCtx::lNvert
Vec lNvert
Definition variables.h:688
UserCtx::Phi
Vec Phi
Definition variables.h:688
SimCtx::euler_subdir
char euler_subdir[PETSC_MAX_PATH_LEN]
Definition variables.h:560
UserCtx::simCtx
SimCtx * simCtx
Back-pointer to the master simulation context.
Definition variables.h:664
SimCtx::rans
PetscInt rans
Definition variables.h:609
SimCtx::StartTime
PetscReal StartTime
Definition variables.h:551
SimCtx::dgf_az
PetscInt dgf_az
Definition variables.h:569
SimCtx::FluxOutSum
PetscReal FluxOutSum
Definition variables.h:602
RankNeighbors::rank_ym
PetscMPIInt rank_ym
Definition variables.h:178
SimCtx::CMy_c
PetscReal CMy_c
Definition variables.h:589
UserCtx::IZet
Vec IZet
Definition variables.h:707
RankNeighbors::rank_xp
PetscMPIInt rank_xp
Definition variables.h:177
UserCtx::Centz
Vec Centz
Definition variables.h:706
PostProcessParams::output_prefix
char output_prefix[256]
Definition variables.h:466
SimCtx::bcs_files
char ** bcs_files
Definition variables.h:601
SimCtx::max_angle
PetscReal max_angle
Definition variables.h:588
UserCtx::IEta
Vec IEta
Definition variables.h:707
SimCtx::imp_rtol
PetscReal imp_rtol
Definition variables.h:576
UserCtx::ksc
PetscInt ksc
Definition variables.h:674
SimCtx::st
PetscReal st
Definition variables.h:581
SimCtx::duplicate
PetscInt duplicate
Definition variables.h:628
SimCtx::tiout
PetscInt tiout
Definition variables.h:550
SimCtx::implicit
PetscInt implicit
Definition variables.h:575
UserCtx::lZet
Vec lZet
Definition variables.h:705
UserCtx::assignedA
PetscBool assignedA
Definition variables.h:701
SimCtx::usermg
UserMG usermg
Definition variables.h:631
SimCtx::allowedFile
char allowedFile[PETSC_MAX_PATH_LEN]
Definition variables.h:632
UserCtx::Csi
Vec Csi
Definition variables.h:705
SimCtx::da_procs_y
PetscInt da_procs_y
Definition variables.h:599
SimCtx::imp_atol
PetscReal imp_atol
Definition variables.h:576
SimCtx::testfilter_1d
PetscInt testfilter_1d
Definition variables.h:612
SimCtx::FieldInitialization
PetscInt FieldInitialization
Definition variables.h:582
SimCtx::ren
PetscReal ren
Definition variables.h:581
SimCtx::Crotz
PetscReal Crotz
Definition variables.h:597
UserCtx::lUcont_rm1
Vec lUcont_rm1
Definition variables.h:692
SimCtx::mixed
PetscInt mixed
Definition variables.h:610
RankCellInfo::zm_cell
PetscInt zm_cell
Definition variables.h:185
UserCtx::lIAj
Vec lIAj
Definition variables.h:707
BoundingBox::max_coords
Cmpnts max_coords
Maximum x, y, z coordinates of the bounding box.
Definition variables.h:156
RankCellInfo::zs_cell
PetscInt zs_cell
Definition variables.h:184
SimCtx::ibmv
IBMVNodes * ibmv
Definition variables.h:625
UserCtx::_this
PetscInt _this
Definition variables.h:674
UserCtx::lKEta
Vec lKEta
Definition variables.h:709
SimCtx::criticalFuncsFile
char criticalFuncsFile[PETSC_MAX_PATH_LEN]
Definition variables.h:641
SimCtx::output_dir
char output_dir[PETSC_MAX_PATH_LEN]
Definition variables.h:559
SimCtx::dt
PetscReal dt
Definition variables.h:552
SimCtx::StepsToRun
PetscInt StepsToRun
Definition variables.h:549
SimCtx::k_periodic
PetscInt k_periodic
Definition variables.h:594
SimCtx::inletprofile
PetscInt inletprofile
Definition variables.h:593
UserCtx::Ucat_nodal
Vec Ucat_nodal
Definition variables.h:735
UserCtx::neighbors
RankNeighbors neighbors
Definition variables.h:673
SimCtx::cdisy
PetscReal cdisy
Definition variables.h:581
UserCtx::lPsi
Vec lPsi
Definition variables.h:730
SimCtx::criticalFuncs
char ** criticalFuncs
Definition variables.h:643
SimCtx::rstart_fsi
PetscBool rstart_fsi
Definition variables.h:627
SimCtx::np
PetscInt np
Definition variables.h:616
UserCtx::jsc
PetscInt jsc
Definition variables.h:674
MGCtx::thislevel
PetscInt thislevel
Definition variables.h:428
SimCtx::averaging
PetscBool averaging
Definition variables.h:613
UserCtx::lJCsi
Vec lJCsi
Definition variables.h:708
SimCtx::ccc
PetscInt ccc
Definition variables.h:605
UserCtx::lCs
Vec lCs
Definition variables.h:712
SimCtx::ratio
PetscReal ratio
Definition variables.h:606
SimCtx::mg_idx
PetscInt mg_idx
Definition variables.h:577
UserCtx::Ucont
Vec Ucont
Definition variables.h:688
SimCtx::StartStep
PetscInt StartStep
Definition variables.h:548
SimCtx::mg_MAX_IT
PetscInt mg_MAX_IT
Definition variables.h:577
BoundingBox::min_coords
Cmpnts min_coords
Minimum x, y, z coordinates of the bounding box.
Definition variables.h:155
SimCtx::OnlySetup
PetscBool OnlySetup
Definition variables.h:553
SimCtx::rotatefsi
PetscInt rotatefsi
Definition variables.h:567
BCS::Ubcs
Vec Ubcs
Boundary condition velocity values. (Comment: "An ugly hack, waste of memory")
Definition variables.h:121
SimCtx::cdisz
PetscReal cdisz
Definition variables.h:581
UserCtx::Qcrit
Vec Qcrit
Definition variables.h:736
Cmpnts::x
PetscScalar x
Definition variables.h:101
UserCtx::JZet
Vec JZet
Definition variables.h:708
SimCtx::dgf_x
PetscInt dgf_x
Definition variables.h:569
SimCtx::current_io_directory
char * current_io_directory
Definition variables.h:564
SimCtx::pizza
PetscInt pizza
Definition variables.h:572
SimCtx::MaxDiv
PetscReal MaxDiv
Definition variables.h:638
UserCtx::Centx
Vec Centx
Definition variables.h:706
UserCtx::Bcs
BCS Bcs
Definition variables.h:682
SimCtx::grid_file
char grid_file[PETSC_MAX_PATH_LEN]
Definition variables.h:598
UserCtx::lPhi
Vec lPhi
Definition variables.h:688
SimCtx::max_cs
PetscReal max_cs
Definition variables.h:611
UserCtx::lParticleCount
Vec lParticleCount
Definition variables.h:729
SimCtx::invicid
PetscInt invicid
Definition variables.h:568
SimCtx::allowedFuncs
char ** allowedFuncs
Definition variables.h:634
RankCellInfo::xm_cell
PetscInt xm_cell
Definition variables.h:185
UserCtx::lUcont_o
Vec lUcont_o
Definition variables.h:691
UserCtx::RankCellInfoMap
RankCellInfo * RankCellInfoMap
Definition variables.h:728
SimCtx::mg_poItr
PetscInt mg_poItr
Definition variables.h:577
SimCtx::STRONG_COUPLING
PetscInt STRONG_COUPLING
Definition variables.h:580
RankCellInfo::ym_cell
PetscInt ym_cell
Definition variables.h:185
UserCtx::Ucat_o
Vec Ucat_o
Definition variables.h:691
SimCtx::MaxDivx
PetscInt MaxDivx
Definition variables.h:639
UserCtx::user_c
UserCtx * user_c
Definition variables.h:722
SimCtx::poisson
PetscInt poisson
Definition variables.h:578
SimCtx::k_homo_filter
PetscInt k_homo_filter
Definition variables.h:612
SimCtx::MaxDivy
PetscInt MaxDivy
Definition variables.h:639
SimCtx::NumberOfBodies
PetscInt NumberOfBodies
Definition variables.h:587
SimCtx::particleRestartMode
char particleRestartMode[16]
Definition variables.h:621
SimCtx::Ogrid
PetscInt Ogrid
Definition variables.h:593
SimCtx::particle_subdir
char particle_subdir[PETSC_MAX_PATH_LEN]
Definition variables.h:561
SimCtx::MaxDivz
PetscInt MaxDivz
Definition variables.h:639
SimCtx::bboxlist
BoundingBox * bboxlist
Definition variables.h:619
SimCtx::j_homo_filter
PetscInt j_homo_filter
Definition variables.h:612
UserCtx::lKZet
Vec lKZet
Definition variables.h:709
UserCtx::Eta
Vec Eta
Definition variables.h:705
SimCtx::eel
PetscInt eel
Definition variables.h:572
SimCtx::log_dir
char log_dir[PETSC_MAX_PATH_LEN]
Definition variables.h:562
SimCtx::MaxDivFlatArg
PetscInt MaxDivFlatArg
Definition variables.h:639
UserCtx::lNu_t
Vec lNu_t
Definition variables.h:712
SimCtx::FluxInSum
PetscReal FluxInSum
Definition variables.h:602
RankNeighbors::rank_xm
PetscMPIInt rank_xm
Definition variables.h:177
UserCtx::Nu_t
Vec Nu_t
Definition variables.h:712
PostProcessParams::source_dir
char source_dir[PETSC_MAX_PATH_LEN]
Definition variables.h:454
UserCtx::lJEta
Vec lJEta
Definition variables.h:708
UserCtx::lCsi
Vec lCsi
Definition variables.h:705
SimCtx::CMz_c
PetscReal CMz_c
Definition variables.h:589
UserCtx::lGridSpace
Vec lGridSpace
Definition variables.h:705
SimCtx::generate_grid
PetscBool generate_grid
Definition variables.h:595
UserCtx::thislevel
PetscInt thislevel
Definition variables.h:721
SimCtx::imp_stol
PetscReal imp_stol
Definition variables.h:576
SimCtx::nAllowed
PetscInt nAllowed
Definition variables.h:635
UserCtx::ICsi
Vec ICsi
Definition variables.h:707
Cmpnts::z
PetscScalar z
Definition variables.h:101
SimCtx::OutputFreq
PetscInt OutputFreq
Definition variables.h:555
UserCtx::lKCsi
Vec lKCsi
Definition variables.h:709
UserCtx::Ucat
Vec Ucat
Definition variables.h:688
UserCtx::ParticleCount
Vec ParticleCount
Definition variables.h:729
SimCtx::Const_CS
PetscReal Const_CS
Definition variables.h:611
UserCtx::Ucont_o
Vec Ucont_o
Definition variables.h:691
SimCtx::i_homo_filter
PetscInt i_homo_filter
Definition variables.h:612
SimCtx::wallfunction
PetscInt wallfunction
Definition variables.h:610
SimCtx::rheology
PetscInt rheology
Definition variables.h:567
SimCtx::Flux_in
PetscReal Flux_in
Definition variables.h:588
UserMG::mglevels
PetscInt mglevels
Definition variables.h:434
SimCtx::cdisx
PetscReal cdisx
Definition variables.h:581
SimCtx::dgf_ax
PetscInt dgf_ax
Definition variables.h:569
SimCtx::cfl
PetscReal cfl
Definition variables.h:581
SimCtx::mglevels
PetscInt mglevels
Definition variables.h:577
SimCtx::num_bcs_files
PetscInt num_bcs_files
Definition variables.h:600
SimCtx::dm_swarm
DM dm_swarm
Definition variables.h:618
SimCtx::useCfg
PetscBool useCfg
Definition variables.h:633
SimCtx::readFields
PetscBool readFields
Definition variables.h:617
SimCtx::central
PetscInt central
Definition variables.h:580
SimCtx::useCriticalFuncsCfg
PetscBool useCriticalFuncsCfg
Definition variables.h:642
SimCtx::Fluxsum
PetscReal Fluxsum
Definition variables.h:602
UserCtx::lJZet
Vec lJZet
Definition variables.h:708
UserCtx::Nvert_o
Vec Nvert_o
Definition variables.h:691
PostProcessParams::outputParticles
PetscBool outputParticles
Definition variables.h:460
SimCtx::Croty
PetscReal Croty
Definition variables.h:597
UserCtx::IAj
Vec IAj
Definition variables.h:707
SimCtx::grid_rotation_angle
PetscReal grid_rotation_angle
Definition variables.h:596
SimCtx::dynamic_freq
PetscInt dynamic_freq
Definition variables.h:610
UserCtx::Psi_nodal
Vec Psi_nodal
Definition variables.h:737
SimCtx::da_procs_x
PetscInt da_procs_x
Definition variables.h:599
UserCtx::JAj
Vec JAj
Definition variables.h:708
SimCtx::U_bc
PetscReal U_bc
Definition variables.h:604
UserCtx::KEta
Vec KEta
Definition variables.h:709
SimCtx::InitialConstantContra
Cmpnts InitialConstantContra
Definition variables.h:583
RankNeighbors::rank_zp
PetscMPIInt rank_zp
Definition variables.h:179
UserCtx::Ucont_rm1
Vec Ucont_rm1
Definition variables.h:692
SimCtx::i_periodic
PetscInt i_periodic
Definition variables.h:594
UserCtx::lUcont
Vec lUcont
Definition variables.h:688
SimCtx::step
PetscInt step
Definition variables.h:546
SimCtx::AreaOutSum
PetscReal AreaOutSum
Definition variables.h:603
SimCtx::dgf_ay
PetscInt dgf_ay
Definition variables.h:569
UserCtx::lAj
Vec lAj
Definition variables.h:705
UserCtx::lICsi
Vec lICsi
Definition variables.h:707
SimCtx::testfilter_ik
PetscInt testfilter_ik
Definition variables.h:612
UserCtx::info
DMDALocalInfo info
Definition variables.h:668
SimCtx::hydro
PetscInt hydro
Definition variables.h:572
UserCtx::lUcat
Vec lUcat
Definition variables.h:688
SimCtx::pps
PostProcessParams * pps
Definition variables.h:648
Cmpnts::y
PetscScalar y
Definition variables.h:101
SimCtx::size
PetscMPIInt size
Definition variables.h:542
EXEC_MODE_SOLVER
@ EXEC_MODE_SOLVER
Definition variables.h:511
EXEC_MODE_POSTPROCESSOR
@ EXEC_MODE_POSTPROCESSOR
Definition variables.h:512
SimCtx::_io_context_buffer
char _io_context_buffer[PETSC_MAX_PATH_LEN]
Definition variables.h:563
UserCtx::lEta
Vec lEta
Definition variables.h:705
UserCtx::KZet
Vec KZet
Definition variables.h:709
UserCtx::Cent
Vec Cent
Definition variables.h:705
SimCtx::les
PetscInt les
Definition variables.h:609
UserCtx::Nvert
Vec Nvert
Definition variables.h:688
UserCtx::KCsi
Vec KCsi
Definition variables.h:709
UserMG::mgctx
MGCtx * mgctx
Definition variables.h:437
SimCtx::mg_preItr
PetscInt mg_preItr
Definition variables.h:577
UserCtx::lNvert_o
Vec lNvert_o
Definition variables.h:691
UserCtx::Centy
Vec Centy
Definition variables.h:706
SimCtx::logviewer
PetscViewer logviewer
Definition variables.h:554
UserCtx::multinullspace
PetscBool multinullspace
Definition variables.h:698
SimCtx::exec_mode
ExecutionMode exec_mode
Definition variables.h:556
SimCtx::ParticleInitialization
PetscInt ParticleInitialization
Definition variables.h:620
SimCtx::imp_MAX_IT
PetscInt imp_MAX_IT
Definition variables.h:575
UserCtx::lJAj
Vec lJAj
Definition variables.h:708
UserCtx::bbox
BoundingBox bbox
Definition variables.h:672
SimCtx::cop
PetscInt cop
Definition variables.h:572
SimCtx::ti
PetscReal ti
Definition variables.h:547
SimCtx::Pipe
PetscInt Pipe
Definition variables.h:572
SimCtx::vnn
PetscReal vnn
Definition variables.h:581
SimCtx::rotateframe
PetscInt rotateframe
Definition variables.h:568
SimCtx::ibm
IBMNodes * ibm
Definition variables.h:624
SimCtx::AreaInSum
PetscReal AreaInSum
Definition variables.h:603
SimCtx::nCriticalFuncs
PetscInt nCriticalFuncs
Definition variables.h:644
SimCtx::summationRHS
PetscReal summationRHS
Definition variables.h:637
UserCtx::lKAj
Vec lKAj
Definition variables.h:709
SimCtx::immersed
PetscInt immersed
Definition variables.h:567
SimCtx::PostprocessingControlFile
char PostprocessingControlFile[PETSC_MAX_PATH_LEN]
Definition variables.h:647
SimCtx::restart_dir
char restart_dir[PETSC_MAX_PATH_LEN]
Definition variables.h:558
SimCtx::blank
PetscInt blank
Definition variables.h:568
SimCtx::dgf_y
PetscInt dgf_y
Definition variables.h:569
SimCtx::LoggingFrequency
PetscInt LoggingFrequency
Definition variables.h:636
SimCtx::CMx_c
PetscReal CMx_c
Definition variables.h:589
UserCtx::Psi
Vec Psi
Definition variables.h:730
UserCtx::P_o
Vec P_o
Definition variables.h:691
BCS::Uch
Vec Uch
Characteristic velocity for boundary conditions.
Definition variables.h:122
SimCtx::j_periodic
PetscInt j_periodic
Definition variables.h:594
SimCtx::wing
PetscInt wing
Definition variables.h:572
SimCtx::fsi
FSInfo * fsi
Definition variables.h:626
BoundingBox
Defines a 3D axis-aligned bounding box.
Definition variables.h:154
Cmpnts
A 3D point or vector with PetscScalar components.
Definition variables.h:100
MGCtx
Context for Multigrid operations.
Definition variables.h:426
PostProcessParams
Holds all configuration parameters for a post-processing run.
Definition variables.h:452
RankCellInfo
A lean struct to hold the global cell ownership range for a single MPI rank.
Definition variables.h:183
SimCtx
The master context for the entire simulation.
Definition variables.h:538
UserCtx
User-defined context containing data specific to a single computational grid level.
Definition variables.h:661
UserMG
User-level context for managing the entire multigrid hierarchy.
Definition variables.h:433