PICurv 0.1.0
A Parallel Particle-In-Cell Solver for Curvilinear LES
Loading...
Searching...
No Matches
Macros | Functions
setup.c File Reference

// Setup code for running any simulation More...

#include <ctype.h>
#include <errno.h>
#include "setup.h"
Include dependency graph for setup.c:

Go to the source code of this file.

Macros

#define __FUNCT__   "CreateSimulationContext"
 
#define __FUNCT__   "PetscMkdirRecursive"
 
#define __FUNCT__   "SetupSimulationEnvironment"
 
#define __FUNCT__   "AllocateContextHeirarchy"
 
#define __FUNCT__   "SetupSolverParameters"
 
#define __FUNCT__   "SetupGridAndSolvers"
 
#define __FUNCT__   "CreateAndInitializeAllVectors"
 
#define __FUNCT__   "UpdateLocalGhosts"
 
#define __FUNCT__   "SetupBoundaryConditions"
 
#define __FUNCT__   "GetOwnedCellRange"
 
#define __FUNCT__   "ComputeAndStoreNeighborRanks"
 
#define __FUNCT__   "SetDMDAProcLayout"
 
#define __FUNCT__   "SetupDomainRankInfo"
 
#define __FUNCT__   "Contra2Cart"
 
#define __FUNCT__   "SetupDomainCellDecompositionMap"
 
#define __FUNCT__   "BinarySearchInt64"
 
#define __FUNCT__   "ComputeDivergence"
 
#define __FUNCT__   "InitializeRandomGenerators"
 
#define __FUNCT__   "InitializeLogicalSpaceRNGs"
 
#define __FUNCT__   "InitializeBrownianRNG"
 
#define __FUNCT__   "TransformScalarDerivativesToPhysical"
 
#define __FUNCT__   "TransformDerivativesToPhysical"
 
#define __FUNCT__   "ComputeScalarFieldDerivatives"
 
#define __FUNCT__   "ComputeVectorFieldDerivatives"
 
#define __FUNCT__   "DestroyUserVectors"
 
#define __FUNCT__   "DestroyUserContext"
 
#define __FUNCT__   "FinalizeSimulation"
 

Functions

PetscBool RuntimeWalltimeGuardParsePositiveSeconds (const char *text, PetscReal *seconds_out)
 Implementation of RuntimeWalltimeGuardParsePositiveSeconds().
 
PetscErrorCode CreateSimulationContext (int argc, char **argv, SimCtx **p_simCtx)
 Implementation of CreateSimulationContext().
 
static PetscErrorCode PetscMkdirRecursive (const char *path)
 Internal helper implementation: PetscMkdirRecursive().
 
PetscErrorCode SetupSimulationEnvironment (SimCtx *simCtx)
 Internal helper implementation: SetupSimulationEnvironment().
 
static PetscErrorCode AllocateContextHierarchy (SimCtx *simCtx)
 Internal helper implementation: AllocateContextHierarchy().
 
static PetscErrorCode SetupSolverParameters (SimCtx *simCtx)
 Internal helper implementation: SetupSolverParameters().
 
PetscErrorCode SetupGridAndSolvers (SimCtx *simCtx)
 Implementation of SetupGridAndSolvers().
 
PetscErrorCode CreateAndInitializeAllVectors (SimCtx *simCtx)
 Internal helper implementation: CreateAndInitializeAllVectors().
 
PetscErrorCode UpdateLocalGhosts (UserCtx *user, const char *fieldName)
 Internal helper implementation: UpdateLocalGhosts().
 
PetscErrorCode SetupBoundaryConditions (SimCtx *simCtx)
 Internal helper implementation: SetupBoundaryConditions().
 
PetscErrorCode Allocate3DArrayScalar (PetscReal ****array, PetscInt nz, PetscInt ny, PetscInt nx)
 Internal helper implementation: Allocate3DArrayScalar().
 
PetscErrorCode Deallocate3DArrayScalar (PetscReal ***array, PetscInt nz, PetscInt ny)
 Internal helper implementation: Deallocate3DArrayScalar().
 
PetscErrorCode Allocate3DArrayVector (Cmpnts ****array, PetscInt nz, PetscInt ny, PetscInt nx)
 Implementation of Allocate3DArrayVector().
 
PetscErrorCode Deallocate3DArrayVector (Cmpnts ***array, PetscInt nz, PetscInt ny)
 Implementation of Deallocate3DArrayVector().
 
PetscErrorCode GetOwnedCellRange (const DMDALocalInfo *info_nodes, PetscInt dim, PetscInt *xs_cell_global_out, PetscInt *xm_cell_local_out)
 Internal helper implementation: GetOwnedCellRange().
 
PetscErrorCode ComputeAndStoreNeighborRanks (UserCtx *user)
 Internal helper implementation: ComputeAndStoreNeighborRanks().
 
PetscErrorCode SetDMDAProcLayout (DM dm, UserCtx *user)
 Internal helper implementation: SetDMDAProcLayout().
 
PetscErrorCode SetupDomainRankInfo (SimCtx *simCtx)
 Implementation of SetupDomainRankInfo().
 
PetscErrorCode Contra2Cart (UserCtx *user)
 Internal helper implementation: Contra2Cart().
 
PetscErrorCode SetupDomainCellDecompositionMap (UserCtx *user)
 Internal helper implementation: SetupDomainCellDecompositionMap().
 
PetscErrorCode BinarySearchInt64 (PetscInt n, const PetscInt64 arr[], PetscInt64 key, PetscBool *found)
 Implementation of BinarySearchInt64().
 
static PetscInt Gidx (PetscInt i, PetscInt j, PetscInt k, UserCtx *user)
 Internal helper implementation: Gidx().
 
PetscErrorCode ComputeDivergence (UserCtx *user)
 Implementation of ComputeDivergence().
 
PetscErrorCode InitializeRandomGenerators (UserCtx *user, PetscRandom *randx, PetscRandom *randy, PetscRandom *randz)
 Implementation of InitializeRandomGenerators().
 
PetscErrorCode InitializeLogicalSpaceRNGs (PetscRandom *rand_logic_i, PetscRandom *rand_logic_j, PetscRandom *rand_logic_k)
 Internal helper implementation: InitializeLogicalSpaceRNGs().
 
PetscErrorCode InitializeBrownianRNG (SimCtx *simCtx)
 Internal helper implementation: InitializeBrownianRNG().
 
void TransformScalarDerivativesToPhysical (PetscReal jacobian, Cmpnts csi_metrics, Cmpnts eta_metrics, Cmpnts zet_metrics, PetscReal dPhi_dcsi, PetscReal dPhi_deta, PetscReal dPhi_dzet, Cmpnts *gradPhi)
 Implementation of TransformScalarDerivativesToPhysical().
 
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)
 Internal helper implementation: TransformDerivativesToPhysical().
 
PetscErrorCode ComputeScalarFieldDerivatives (UserCtx *user, PetscInt i, PetscInt j, PetscInt k, PetscReal ***field_data, Cmpnts *grad)
 Internal helper implementation: ComputeScalarFieldDerivatives().
 
PetscErrorCode ComputeVectorFieldDerivatives (UserCtx *user, PetscInt i, PetscInt j, PetscInt k, Cmpnts ***field_data, Cmpnts *dudx, Cmpnts *dvdx, Cmpnts *dwdx)
 Internal helper implementation: ComputeVectorFieldDerivatives().
 
PetscErrorCode DestroyUserVectors (UserCtx *user)
 Internal helper implementation: DestroyUserVectors().
 
PetscErrorCode DestroyUserContext (UserCtx *user)
 Internal helper implementation: DestroyUserContext().
 
PetscErrorCode FinalizeSimulation (SimCtx *simCtx)
 Implementation of FinalizeSimulation().
 

Detailed Description

// Setup code for running any simulation

Test program for DMSwarm interpolation using the fdf-curvIB method. Provides the setup to start any simulation with DMSwarm and DMDAs.

Definition in file setup.c.

Macro Definition Documentation

◆ __FUNCT__ [1/27]

#define __FUNCT__   "CreateSimulationContext"

Definition at line 42 of file setup.c.

◆ __FUNCT__ [2/27]

#define __FUNCT__   "PetscMkdirRecursive"

Definition at line 42 of file setup.c.

◆ __FUNCT__ [3/27]

#define __FUNCT__   "SetupSimulationEnvironment"

Definition at line 42 of file setup.c.

◆ __FUNCT__ [4/27]

#define __FUNCT__   "AllocateContextHeirarchy"

Definition at line 42 of file setup.c.

◆ __FUNCT__ [5/27]

#define __FUNCT__   "SetupSolverParameters"

Definition at line 42 of file setup.c.

◆ __FUNCT__ [6/27]

#define __FUNCT__   "SetupGridAndSolvers"

Definition at line 42 of file setup.c.

◆ __FUNCT__ [7/27]

#define __FUNCT__   "CreateAndInitializeAllVectors"

Definition at line 42 of file setup.c.

◆ __FUNCT__ [8/27]

#define __FUNCT__   "UpdateLocalGhosts"

Definition at line 42 of file setup.c.

◆ __FUNCT__ [9/27]

#define __FUNCT__   "SetupBoundaryConditions"

Definition at line 42 of file setup.c.

◆ __FUNCT__ [10/27]

#define __FUNCT__   "GetOwnedCellRange"

Definition at line 42 of file setup.c.

◆ __FUNCT__ [11/27]

#define __FUNCT__   "ComputeAndStoreNeighborRanks"

Definition at line 42 of file setup.c.

◆ __FUNCT__ [12/27]

#define __FUNCT__   "SetDMDAProcLayout"

Definition at line 42 of file setup.c.

◆ __FUNCT__ [13/27]

#define __FUNCT__   "SetupDomainRankInfo"

Definition at line 42 of file setup.c.

◆ __FUNCT__ [14/27]

#define __FUNCT__   "Contra2Cart"

Definition at line 42 of file setup.c.

◆ __FUNCT__ [15/27]

#define __FUNCT__   "SetupDomainCellDecompositionMap"

Definition at line 42 of file setup.c.

◆ __FUNCT__ [16/27]

#define __FUNCT__   "BinarySearchInt64"

Definition at line 42 of file setup.c.

◆ __FUNCT__ [17/27]

#define __FUNCT__   "ComputeDivergence"

Definition at line 42 of file setup.c.

◆ __FUNCT__ [18/27]

#define __FUNCT__   "InitializeRandomGenerators"

Definition at line 42 of file setup.c.

◆ __FUNCT__ [19/27]

#define __FUNCT__   "InitializeLogicalSpaceRNGs"

Definition at line 42 of file setup.c.

◆ __FUNCT__ [20/27]

#define __FUNCT__   "InitializeBrownianRNG"

Definition at line 42 of file setup.c.

◆ __FUNCT__ [21/27]

#define __FUNCT__   "TransformScalarDerivativesToPhysical"

Definition at line 42 of file setup.c.

◆ __FUNCT__ [22/27]

#define __FUNCT__   "TransformDerivativesToPhysical"

Definition at line 42 of file setup.c.

◆ __FUNCT__ [23/27]

#define __FUNCT__   "ComputeScalarFieldDerivatives"

Definition at line 42 of file setup.c.

◆ __FUNCT__ [24/27]

#define __FUNCT__   "ComputeVectorFieldDerivatives"

Definition at line 42 of file setup.c.

◆ __FUNCT__ [25/27]

#define __FUNCT__   "DestroyUserVectors"

Definition at line 42 of file setup.c.

◆ __FUNCT__ [26/27]

#define __FUNCT__   "DestroyUserContext"

Definition at line 42 of file setup.c.

◆ __FUNCT__ [27/27]

#define __FUNCT__   "FinalizeSimulation"

Definition at line 42 of file setup.c.

Function Documentation

◆ RuntimeWalltimeGuardParsePositiveSeconds()

PetscBool RuntimeWalltimeGuardParsePositiveSeconds ( const char *  text,
PetscReal *  seconds_out 
)

Implementation of RuntimeWalltimeGuardParsePositiveSeconds().

Parse a positive floating-point seconds value from runtime metadata.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/setup.h.

See also
RuntimeWalltimeGuardParsePositiveSeconds()

Definition at line 18 of file setup.c.

19{
20 char *endptr = NULL;
21 double parsed_value;
22
23 if (seconds_out) *seconds_out = 0.0;
24 if (!text || text[0] == '\0') return PETSC_FALSE;
25
26 errno = 0;
27 parsed_value = strtod(text, &endptr);
28 if (endptr == text || errno == ERANGE || !isfinite(parsed_value) || parsed_value <= 0.0) {
29 return PETSC_FALSE;
30 }
31
32 while (*endptr != '\0' && isspace((unsigned char)*endptr)) {
33 endptr++;
34 }
35 if (*endptr != '\0') return PETSC_FALSE;
36
37 if (seconds_out) *seconds_out = (PetscReal)parsed_value;
38 return PETSC_TRUE;
39}
Here is the caller graph for this function:

◆ CreateSimulationContext()

PetscErrorCode CreateSimulationContext ( int  argc,
char **  argv,
SimCtx **  p_simCtx 
)

Implementation of CreateSimulationContext().

Allocates and populates the master SimulationContext object.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/setup.h.

See also
CreateSimulationContext()

Definition at line 50 of file setup.c.

51{
52 PetscErrorCode ierr;
53 (void)argc;
54 (void)argv;
55 SimCtx *simCtx;
56 char control_filename[PETSC_MAX_PATH_LEN] = ""; // Temporary placeholder for control file name.
57 PetscBool control_flg; // Temporary placeholder for control file tag existence check flag.
58 PetscBool particle_console_output_freq_flg = PETSC_FALSE;
59
60 PetscFunctionBeginUser;
61
63
64 // === 1. Allocate the Context Struct and Set ALL Defaults ==================
65 ierr = PetscNew(p_simCtx); CHKERRQ(ierr);
66 simCtx = *p_simCtx;
67
68 // --- Group 1: Parallelism & MPI ---
69 simCtx->rank = 0; simCtx->size = 1;
70
71 // --- Group 2: Simulation Control, Time, and I/O ---
72 simCtx->step = 0; simCtx->ti = 0.0; simCtx->StartStep = 0; simCtx->StepsToRun = 10;
73 simCtx->tiout = 10; simCtx->particleConsoleOutputFreq = simCtx->tiout;
74 simCtx->StartTime = 0.0; simCtx->dt = 0.001;
75 simCtx->OnlySetup = PETSC_FALSE;
76 simCtx->continueMode = PETSC_FALSE;
77 simCtx->logviewer = NULL;
78 strcpy(simCtx->eulerianSource,"solve");
79 strcpy(simCtx->restart_dir,"restart");
80 strcpy(simCtx->output_dir,"output");
81 strcpy(simCtx->log_dir,"logs");
82 strcpy(simCtx->euler_subdir,"eulerian");
83 strcpy(simCtx->particle_subdir,"particles");
84 simCtx->_io_context_buffer[0] = '\0';
85 simCtx->current_io_directory = NULL;
86
87 // --- Group 3: High-Level Physics & Model Selection Flags ---
88 simCtx->immersed = 0; simCtx->movefsi = 0; simCtx->rotatefsi = 0;
89 simCtx->sediment = 0; simCtx->rheology = 0; simCtx->invicid = 0;
90 simCtx->TwoD = 0; simCtx->thin = 0; simCtx->moveframe = 0;
91 simCtx->rotateframe = 0; simCtx->blank = 0;
92 simCtx->dgf_x = 0; simCtx->dgf_y = 1; simCtx->dgf_z = 0;
93 simCtx->dgf_ax = 1; simCtx->dgf_ay = 0; simCtx->dgf_az = 0;
94 strcpy(simCtx->AnalyticalSolutionType,"TGV3D");
95
96 // --- Group 4: Specific Simulation Case Flags --- (DEPRICATED)
97 simCtx->cop=0; simCtx->fish=0; simCtx->fish_c=0; simCtx->fishcyl=0;
98 simCtx->eel=0; simCtx->pizza=0; simCtx->turbine=0; simCtx->Pipe=0;
99 simCtx->wing=0; simCtx->hydro=0; simCtx->MHV=0; simCtx->LV=0;
100 simCtx->channelz = 0;
101
102 // --- Group 5: Solver & Numerics Parameters ---
104 simCtx->mom_dt_rk4_residual_norm_noise_allowance_factor = 1.05; // New addition for divergence detection
105 simCtx->mom_atol = 1e-7; simCtx->mom_rtol = 1e-4; simCtx->imp_stol = 1.e-8;
106 simCtx->mglevels = 3; simCtx->mg_MAX_IT = 30; simCtx->mg_idx = 1;
107 simCtx->mg_preItr = 1; simCtx->mg_poItr = 1;
108 simCtx->poisson = 0; simCtx->poisson_tol = 5.e-9;
109 simCtx->STRONG_COUPLING = 0;simCtx->central=0;
110 simCtx->ren = 100.0; simCtx->pseudo_cfl = 0.1;
111 simCtx->max_pseudo_cfl = 1.0; simCtx->min_pseudo_cfl = 0.001;
112 simCtx->pseudo_cfl_reduction_factor = 1.0;
113 simCtx->pseudo_cfl_growth_factor = 1.0; //simCtx->vnn = 0.1;
114 simCtx->ps_ksp_pic_monitor_true_residual = PETSC_FALSE;
115 simCtx->cdisx = 0.0; simCtx->cdisy = 0.0; simCtx->cdisz = 0.0;
116 simCtx->FieldInitialization = 0;
117 simCtx->InitialConstantContra.x = 0.0;
118 simCtx->InitialConstantContra.y = 0.0;
119 simCtx->InitialConstantContra.z = 0.0;
120 simCtx->AnalyticalUniformVelocity.x = 0.0;
121 simCtx->AnalyticalUniformVelocity.y = 0.0;
122 simCtx->AnalyticalUniformVelocity.z = 0.0;
123 simCtx->verificationDiffusivity.enabled = PETSC_FALSE;
124 strcpy(simCtx->verificationDiffusivity.mode, "");
125 strcpy(simCtx->verificationDiffusivity.profile, "");
126 simCtx->verificationDiffusivity.gamma0 = 0.0;
127 simCtx->verificationDiffusivity.slope_x = 0.0;
128
129 // --- Group 6: Physical & Geometric Parameters ---
130 simCtx->NumberOfBodies = 1; simCtx->Flux_in = 1.0; simCtx->angle = 0.0;
131 simCtx->max_angle = -54. * 3.1415926 / 180.;
132 simCtx->CMx_c=0.0; simCtx->CMy_c=0.0; simCtx->CMz_c=0.0;
133 simCtx->wall_roughness_height = 1e-16;
134 simCtx->schmidt_number = 1.0; simCtx->Turbulent_schmidt_number = 0.7;
135
136 // --- Group 7: Grid, Domain, and Boundary Condition Settings ---
137 simCtx->block_number = 1; simCtx->inletprofile = 1;
138 simCtx->grid1d = 0; simCtx->Ogrid = 0;
139 simCtx->i_periodic = 0; simCtx->j_periodic = 0; simCtx->k_periodic = 0;
140 simCtx->blkpbc = 10; simCtx->pseudo_periodic = 0;
141 strcpy(simCtx->grid_file, "config/grid.run");
142 simCtx->generate_grid = PETSC_FALSE;
143 simCtx->da_procs_x = PETSC_DECIDE;
144 simCtx->da_procs_y = PETSC_DECIDE;
145 simCtx->da_procs_z = PETSC_DECIDE;
146 simCtx->grid_rotation_angle = 0.0;
147 simCtx->Croty = 0.0; simCtx->Crotz = 0.0;
148 simCtx->num_bcs_files = 1;
149 ierr = PetscMalloc1(1, &simCtx->bcs_files); CHKERRQ(ierr);
150 ierr = PetscStrallocpy("config/bcs.run", &simCtx->bcs_files[0]); CHKERRQ(ierr);
151 simCtx->FluxInSum = 0.0; simCtx->FluxOutSum = 0.0; simCtx->Fluxsum = 0.0;
152 simCtx->drivingForceMagnitude = 0.0, simCtx->forceScalingFactor = 1.8;
153 simCtx->targetVolumetricFlux = 0.0;
154 simCtx->AreaInSum = 0.0; simCtx->AreaOutSum = 0.0;
155 simCtx->U_bc = 0.0; simCtx->ccc = 0;
156 simCtx->ratio = 0.0;
157
158
159 // --- Group 8: Turbulence Modeling (LES/RANS) ---
160 simCtx->averaging = PETSC_FALSE; simCtx->les = NO_LES_MODEL; simCtx->rans = 0;
161 simCtx->wallfunction = 0; simCtx->mixed = 0; simCtx->clark = 0;
162 simCtx->dynamic_freq = 1; simCtx->max_cs = 0.5;
163 simCtx->Const_CS = 0.03;
164 simCtx->testfilter_ik = 0; simCtx->testfilter_1d = 0;
165 simCtx->i_homo_filter = 0; simCtx->j_homo_filter = 0; simCtx->k_homo_filter = 0;
166
167 // --- Group 9: Particle / DMSwarm Data & Settings ---
168 simCtx->np = 0; simCtx->readFields = PETSC_FALSE;
169 simCtx->dm_swarm = NULL; simCtx->bboxlist = NULL;
172 strcpy(simCtx->particleRestartMode,"load");
173 simCtx->particlesLostLastStep = 0;
174 simCtx->particlesLostCumulative = 0;
175 simCtx->particlesMigratedLastStep = 0;
176 simCtx->occupiedCellCount = 0;
177 simCtx->particleLoadImbalance = 0.0;
178 simCtx->migrationPassesLastStep = 0;
179 simCtx->searchMetrics.searchAttempts = 0;
182 simCtx->searchMetrics.searchLostCount = 0;
184 simCtx->searchMetrics.reSearchCount = 0;
187 simCtx->searchMetrics.tieBreakCount = 0;
193 simCtx->BrownianMotionRNG = NULL;
194 simCtx->C_IEM = 2.0;
195
196 // --- Group 10: Immersed Boundary & FSI Data Object Pointers ---
197 simCtx->ibm = NULL; simCtx->ibmv = NULL; simCtx->fsi = NULL;
198 simCtx->rstart_fsi = PETSC_FALSE; simCtx->duplicate = 0;
199
200 // --- Group 11: Logging and Custom Configuration ---
201 strcpy(simCtx->allowedFile, "config/whitelist.run");
202 simCtx->useCfg = PETSC_FALSE;
203 simCtx->allowedFuncs = NULL;
204 simCtx->nAllowed = 0;
205 simCtx->LoggingFrequency = 10;
206 simCtx->summationRHS = 0.0;
207 simCtx->MaxDiv = 0.0;
208 simCtx->MaxDivFlatArg = 0; simCtx->MaxDivx = 0; simCtx->MaxDivy = 0; simCtx->MaxDivz = 0;
209 strcpy(simCtx->profilingSelectedFuncsFile, "config/profile.run");
210 simCtx->useProfilingSelectedFuncsCfg = PETSC_FALSE;
211 simCtx->profilingSelectedFuncs = NULL;
212 simCtx->nProfilingSelectedFuncs = 0;
213 strcpy(simCtx->profilingTimestepMode, "selected");
214 strcpy(simCtx->profilingTimestepFile, "Profiling_Timestep_Summary.csv");
215 simCtx->profilingFinalSummary = PETSC_TRUE;
216 simCtx->walltimeGuardEnabled = PETSC_FALSE;
217 simCtx->walltimeGuardActive = PETSC_FALSE;
218 simCtx->walltimeGuardWarmupSteps = 10;
219 simCtx->walltimeGuardMultiplier = 2.0;
220 simCtx->walltimeGuardMinSeconds = 60.0;
221 simCtx->walltimeGuardEstimatorAlpha = 0.35;
223 simCtx->walltimeGuardLimitSeconds = 0.0;
224 simCtx->walltimeGuardCompletedSteps = 0;
227 simCtx->walltimeGuardHasEWMA = PETSC_FALSE;
228 simCtx->walltimeGuardEWMASeconds = 0.0;
229 simCtx->walltimeGuardLatestStepSeconds = 0.0;
230 // --- Group 11: Post-Processing Information ---
231 strcpy(simCtx->PostprocessingControlFile, "config/post.run");
232 ierr = PetscNew(&simCtx->pps); CHKERRQ(ierr);
233
234 // === 2. Get MPI Info and Handle Config File =============================
235 // -- Group 1: Parallelism & MPI Information
236 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &simCtx->rank); CHKERRQ(ierr);
237 ierr = MPI_Comm_size(PETSC_COMM_WORLD, &simCtx->size); CHKERRQ(ierr);
238
239 // First, check if the -control_file argument was provided by the user/script.
240 ierr = PetscOptionsGetString(NULL, NULL, "-control_file", control_filename, sizeof(control_filename), &control_flg); CHKERRQ(ierr);
241
242 // If the flag is NOT present or the filename is empty, abort with a helpful error.
243 if (!control_flg || strlen(control_filename) == 0) {
244 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
245 "\n\n*** MANDATORY ARGUMENT MISSING ***\n"
246 "The -control_file argument was not provided.\n"
247 "This program must be launched with a configuration file.\n"
248 "Example: mpiexec -n 4 ./simulator -control_file /path/to/your/config.control\n"
249 "This is typically handled automatically by the 'picurv' script.\n");
250 }
251
252 // At this point, we have a valid filename. Attempt to load it.
253 LOG(GLOBAL, LOG_INFO, "Loading mandatory configuration from: %s\n", control_filename);
254 ierr = PetscOptionsInsertFile(PETSC_COMM_WORLD, NULL, control_filename, PETSC_FALSE);
255 if (ierr == PETSC_ERR_FILE_OPEN) {
256 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_FILE_OPEN, "The specified control file was not found or could not be opened: %s", control_filename);
257 }
258 CHKERRQ(ierr);
259
260 // === 3. A Configure Logging System ========================================
261 // This logic determines the logging configuration and STORES it in simCtx for
262 // later reference and cleanup.
263 ierr = PetscOptionsGetString(NULL, NULL, "-whitelist_config_file", simCtx->allowedFile, PETSC_MAX_PATH_LEN, &simCtx->useCfg); CHKERRQ(ierr);
264
265 if (simCtx->useCfg) {
266 ierr = LoadAllowedFunctionsFromFile(simCtx->allowedFile, &simCtx->allowedFuncs, &simCtx->nAllowed);
267 if (ierr) {
268 // Use direct PetscPrintf as logging system isn't fully active yet.
269 PetscPrintf(PETSC_COMM_SELF, "[%s] WARNING: Failed to load allowed functions from '%s'. Falling back to default list.\n", __func__, simCtx->allowedFile);
270 simCtx->useCfg = PETSC_FALSE; // Mark as failed.
271 ierr = 0; // Clear the error to allow fallback.
272 } else if (simCtx->nAllowed == 0) {
273 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
274 "Whitelist config file '%s' is empty. Omit -whitelist_config_file to use the default allow-list, or list at least one function.",
275 simCtx->allowedFile);
276 }
277 }
278 if (!simCtx->useCfg) {
279 // Fallback to default logging functions if no file was used or if loading failed.
280 simCtx->nAllowed = 2;
281 ierr = PetscMalloc1(simCtx->nAllowed, &simCtx->allowedFuncs); CHKERRQ(ierr);
282 ierr = PetscStrallocpy("main", &simCtx->allowedFuncs[0]); CHKERRQ(ierr);
283 ierr = PetscStrallocpy("CreateSimulationContext", &simCtx->allowedFuncs[1]); CHKERRQ(ierr);
284 }
285
286 // Activate the configuration by passing it to the logging module's setup function.
287 set_allowed_functions((const char**)simCtx->allowedFuncs, (size_t)simCtx->nAllowed);
288
289 // Now that the logger is configured, we can use it.
290 LOG_ALLOW_SYNC(LOCAL, LOG_INFO, "Context created. Initializing on rank %d of %d.\n", simCtx->rank, simCtx->size);
291 print_log_level(); // This will now correctly reflect the LOG_LEVEL environment variable.
292
293 // === 3.B Configure Profiling System ========================================
294 ierr = PetscOptionsGetString(NULL, NULL, "-profiling_timestep_mode", simCtx->profilingTimestepMode, sizeof(simCtx->profilingTimestepMode), NULL); CHKERRQ(ierr);
295 ierr = PetscOptionsGetString(NULL, NULL, "-profiling_timestep_file", simCtx->profilingTimestepFile, PETSC_MAX_PATH_LEN, NULL); CHKERRQ(ierr);
296 ierr = PetscOptionsGetBool(NULL, NULL, "-profiling_final_summary", &simCtx->profilingFinalSummary, NULL); CHKERRQ(ierr);
297 if (strcmp(simCtx->profilingTimestepMode, "off") != 0 &&
298 strcmp(simCtx->profilingTimestepMode, "selected") != 0 &&
299 strcmp(simCtx->profilingTimestepMode, "all") != 0) {
300 PetscPrintf(PETSC_COMM_SELF, "[%s] WARNING: Unknown profiling timestep mode '%s'. Falling back to 'selected'.\n", __func__, simCtx->profilingTimestepMode);
301 strcpy(simCtx->profilingTimestepMode, "selected");
302 }
303
304 if (strcmp(simCtx->profilingTimestepMode, "selected") == 0) {
305 ierr = PetscOptionsGetString(NULL, NULL, "-profile_config_file", simCtx->profilingSelectedFuncsFile, PETSC_MAX_PATH_LEN, &simCtx->useProfilingSelectedFuncsCfg); CHKERRQ(ierr);
306 if (simCtx->useProfilingSelectedFuncsCfg) {
308 if (ierr) {
309 PetscPrintf(PETSC_COMM_SELF, "[%s] WARNING: Failed to load selected profiling functions from '%s'. Falling back to default list.\n", __func__, simCtx->profilingSelectedFuncsFile);
310 simCtx->useProfilingSelectedFuncsCfg = PETSC_FALSE;
311 ierr = 0;
312 }
313 }
314 if (!simCtx->useProfilingSelectedFuncsCfg) {
315 // Fallback to a hardcoded default list if no file was provided or loading failed.
316 simCtx->nProfilingSelectedFuncs = 4;
317 ierr = PetscMalloc1(simCtx->nProfilingSelectedFuncs, &simCtx->profilingSelectedFuncs); CHKERRQ(ierr);
318 ierr = PetscStrallocpy("FlowSolver", &simCtx->profilingSelectedFuncs[0]); CHKERRQ(ierr);
319 ierr = PetscStrallocpy("AdvanceSimulation", &simCtx->profilingSelectedFuncs[1]); CHKERRQ(ierr);
320 ierr = PetscStrallocpy("LocateAllParticlesInGrid", &simCtx->profilingSelectedFuncs[2]); CHKERRQ(ierr);
321 ierr = PetscStrallocpy("InterpolateAllFieldsToSwarm", &simCtx->profilingSelectedFuncs[3]); CHKERRQ(ierr);
322 }
323 }
324
325 // Initialize the profiling system with the current updated simulation context.
326 ierr = ProfilingInitialize(simCtx); CHKERRQ(ierr);
327
328 // === 4. Parse All Command Line Options ==================================
329 LOG_ALLOW(GLOBAL, LOG_INFO, "Parsing command-line options...\n");
330
331 // --- Group 2
332 LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 2: Simulation Control,Time and I/O.\n");
333 // Read the physical time to start from.
334 // The default is already 0.0, so this will only be non-zero if the user provides it.
335 ierr = PetscOptionsGetInt(NULL, NULL, "-start_step", &simCtx->StartStep, NULL); CHKERRQ(ierr);
336 ierr = PetscOptionsGetInt(NULL,NULL, "-totalsteps", &simCtx->StepsToRun, NULL); CHKERRQ(ierr);
337 ierr = PetscOptionsGetBool(NULL, NULL, "-only_setup", &simCtx->OnlySetup, NULL); CHKERRQ(ierr);
338 ierr = PetscOptionsGetBool(NULL, NULL, "-continue_mode", &simCtx->continueMode, NULL); CHKERRQ(ierr);
339 ierr = PetscOptionsGetReal(NULL, NULL, "-dt", &simCtx->dt, NULL); CHKERRQ(ierr);
340 ierr = PetscOptionsGetInt(NULL, NULL, "-tio", &simCtx->tiout, NULL); CHKERRQ(ierr);
341 ierr = PetscOptionsGetInt(NULL, NULL, "-particle_console_output_freq", &simCtx->particleConsoleOutputFreq, &particle_console_output_freq_flg); CHKERRQ(ierr);
342 if (!particle_console_output_freq_flg) {
343 simCtx->particleConsoleOutputFreq = simCtx->tiout;
344 }
345 ierr = PetscOptionsGetString(NULL,NULL,"-euler_field_source",simCtx->eulerianSource,sizeof(simCtx->eulerianSource),NULL);CHKERRQ(ierr);
346 ierr = PetscOptionsGetString(NULL,NULL,"-output_dir",simCtx->output_dir,sizeof(simCtx->output_dir),NULL);CHKERRQ(ierr);
347 ierr = PetscOptionsGetString(NULL,NULL,"-restart_dir",simCtx->restart_dir,sizeof(simCtx->restart_dir),NULL);CHKERRQ(ierr);
348 ierr = PetscOptionsGetString(NULL,NULL,"-log_dir",simCtx->log_dir,sizeof(simCtx->log_dir),NULL);CHKERRQ(ierr);
349 ierr = PetscOptionsGetString(NULL,NULL,"-euler_subdir",simCtx->euler_subdir,sizeof(simCtx->euler_subdir),NULL);CHKERRQ(ierr);
350 ierr = PetscOptionsGetString(NULL,NULL,"-particle_subdir",simCtx->particle_subdir,sizeof(simCtx->particle_subdir),NULL);CHKERRQ(ierr);
351 ierr = PetscOptionsGetBool(NULL, NULL, "-walltime_guard_enabled", &simCtx->walltimeGuardEnabled, NULL); CHKERRQ(ierr);
352 ierr = PetscOptionsGetInt(NULL, NULL, "-walltime_guard_warmup_steps", &simCtx->walltimeGuardWarmupSteps, NULL); CHKERRQ(ierr);
353 ierr = PetscOptionsGetReal(NULL, NULL, "-walltime_guard_multiplier", &simCtx->walltimeGuardMultiplier, NULL); CHKERRQ(ierr);
354 ierr = PetscOptionsGetReal(NULL, NULL, "-walltime_guard_min_seconds", &simCtx->walltimeGuardMinSeconds, NULL); CHKERRQ(ierr);
355 ierr = PetscOptionsGetReal(NULL, NULL, "-walltime_guard_estimator_alpha", &simCtx->walltimeGuardEstimatorAlpha, NULL); CHKERRQ(ierr);
356
357 if (simCtx->walltimeGuardWarmupSteps <= 0) {
358 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Invalid value for -walltime_guard_warmup_steps: %d. Must be > 0.", simCtx->walltimeGuardWarmupSteps);
359 }
360 if (simCtx->walltimeGuardMultiplier <= 0.0 || simCtx->walltimeGuardMultiplier > 5.0) {
361 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Invalid value for -walltime_guard_multiplier: %.6f. Must be in (0, 5].", (double)simCtx->walltimeGuardMultiplier);
362 }
363 if (simCtx->walltimeGuardMinSeconds <= 0.0) {
364 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Invalid value for -walltime_guard_min_seconds: %.6f. Must be > 0.", (double)simCtx->walltimeGuardMinSeconds);
365 }
366 if (simCtx->walltimeGuardEstimatorAlpha <= 0.0 || simCtx->walltimeGuardEstimatorAlpha > 1.0) {
367 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Invalid value for -walltime_guard_estimator_alpha: %.6f. Must be in (0, 1].", (double)simCtx->walltimeGuardEstimatorAlpha);
368 }
369
370 if(strcmp(simCtx->eulerianSource,"solve")!= 0 && strcmp(simCtx->eulerianSource,"load") != 0 && strcmp(simCtx->eulerianSource,"analytical")!=0){
371 SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Invalid value for -euler_field_source. Must be 'load','analytical' or 'solve'. You provided '%s'.",simCtx->eulerianSource);
372 }
373
374 if (simCtx->walltimeGuardEnabled) {
375 const char *job_start_env = getenv("PICURV_JOB_START_EPOCH");
376 const char *limit_env = getenv("PICURV_WALLTIME_LIMIT_SECONDS");
377 PetscBool job_start_ok = RuntimeWalltimeGuardParsePositiveSeconds(job_start_env, &simCtx->walltimeGuardJobStartEpochSeconds);
378 PetscBool limit_ok = RuntimeWalltimeGuardParsePositiveSeconds(limit_env, &simCtx->walltimeGuardLimitSeconds);
379
380 if (!job_start_ok || !limit_ok) {
381 simCtx->walltimeGuardActive = PETSC_FALSE;
383 simCtx->walltimeGuardLimitSeconds = 0.0;
384 LOG_ALLOW(
385 GLOBAL,
387 "Runtime walltime guard enabled but %s/%s are missing or invalid. Falling back to external shutdown signals only.\n",
388 "PICURV_JOB_START_EPOCH",
389 "PICURV_WALLTIME_LIMIT_SECONDS"
390 );
391 } else {
392 simCtx->walltimeGuardActive = PETSC_TRUE;
393 }
394 }
395
396 // --- Group 3
397 LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 3: High-Level Physics & Model Selection Flags\n");
398 ierr = PetscOptionsGetInt(NULL, NULL, "-imm", &simCtx->immersed, NULL); CHKERRQ(ierr);
399 ierr = PetscOptionsGetInt(NULL, NULL, "-fsi", &simCtx->movefsi, NULL); CHKERRQ(ierr);
400 ierr = PetscOptionsGetInt(NULL, NULL, "-rfsi", &simCtx->rotatefsi, NULL); CHKERRQ(ierr);
401 ierr = PetscOptionsGetInt(NULL, NULL, "-sediment", &simCtx->sediment, NULL); CHKERRQ(ierr);
402 ierr = PetscOptionsGetInt(NULL, NULL, "-rheology", &simCtx->rheology, NULL); CHKERRQ(ierr);
403 ierr = PetscOptionsGetInt(NULL, NULL, "-inv", &simCtx->invicid, NULL); CHKERRQ(ierr);
404 ierr = PetscOptionsGetInt(NULL, NULL, "-TwoD", &simCtx->TwoD, NULL); CHKERRQ(ierr);
405 ierr = PetscOptionsGetInt(NULL, NULL, "-thin", &simCtx->thin, NULL); CHKERRQ(ierr);
406 ierr = PetscOptionsGetInt(NULL, NULL, "-mframe", &simCtx->moveframe, NULL); CHKERRQ(ierr);
407 ierr = PetscOptionsGetInt(NULL, NULL, "-rframe", &simCtx->rotateframe, NULL); CHKERRQ(ierr);
408 ierr = PetscOptionsGetInt(NULL, NULL, "-blk", &simCtx->blank, NULL); CHKERRQ(ierr);
409 ierr = PetscOptionsGetInt(NULL, NULL, "-dgf_z", &simCtx->dgf_z, NULL); CHKERRQ(ierr);
410 ierr = PetscOptionsGetInt(NULL, NULL, "-dgf_y", &simCtx->dgf_y, NULL); CHKERRQ(ierr);
411 ierr = PetscOptionsGetInt(NULL, NULL, "-dgf_x", &simCtx->dgf_x, NULL); CHKERRQ(ierr);
412 ierr = PetscOptionsGetInt(NULL, NULL, "-dgf_az", &simCtx->dgf_az, NULL); CHKERRQ(ierr);
413 ierr = PetscOptionsGetInt(NULL, NULL, "-dgf_ay", &simCtx->dgf_ay, NULL); CHKERRQ(ierr);
414 ierr = PetscOptionsGetInt(NULL, NULL, "-dgf_ax", &simCtx->dgf_ax, NULL); CHKERRQ(ierr);
415 ierr = PetscOptionsGetString(NULL,NULL,"-analytical_type",simCtx->AnalyticalSolutionType,sizeof(simCtx->AnalyticalSolutionType),NULL);CHKERRQ(ierr);
416
417 // --- Group 4
418 LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 4: Specific Simulation Case Flags \n");
419 ierr = PetscOptionsGetInt(NULL, NULL, "-cop", &simCtx->cop, NULL); CHKERRQ(ierr);
420 ierr = PetscOptionsGetInt(NULL, NULL, "-fish", &simCtx->fish, NULL); CHKERRQ(ierr);
421 ierr = PetscOptionsGetInt(NULL, NULL, "-pizza", &simCtx->pizza, NULL); CHKERRQ(ierr);
422 ierr = PetscOptionsGetInt(NULL, NULL, "-turbine", &simCtx->turbine, NULL); CHKERRQ(ierr);
423 ierr = PetscOptionsGetInt(NULL, NULL, "-fishcyl", &simCtx->fishcyl, NULL); CHKERRQ(ierr);
424 ierr = PetscOptionsGetInt(NULL, NULL, "-eel", &simCtx->eel, NULL); CHKERRQ(ierr);
425 ierr = PetscOptionsGetInt(NULL, NULL, "-cstart", &simCtx->fish_c, NULL); CHKERRQ(ierr);
426 ierr = PetscOptionsGetInt(NULL, NULL, "-wing", &simCtx->wing, NULL); CHKERRQ(ierr);
427 ierr = PetscOptionsGetInt(NULL, NULL, "-mhv", &simCtx->MHV, NULL); CHKERRQ(ierr);
428 ierr = PetscOptionsGetInt(NULL, NULL, "-hydro", &simCtx->hydro, NULL); CHKERRQ(ierr);
429 ierr = PetscOptionsGetInt(NULL, NULL, "-lv", &simCtx->LV, NULL); CHKERRQ(ierr);
430 ierr = PetscOptionsGetInt(NULL, NULL, "-Pipe", &simCtx->Pipe, NULL); CHKERRQ(ierr);
431 ierr = PetscOptionsGetInt(NULL, NULL, "-Turbulent_Channel_z", &simCtx->channelz, NULL); CHKERRQ(ierr);
432 ierr = PetscOptionsGetReal(NULL,NULL,"-Turbulent_Channel_z_Driving_Force",&simCtx->drivingForceMagnitude,NULL);CHKERRQ(ierr);
433 ierr = PetscOptionsGetReal(NULL,NULL,"-Turbulent_Channel_z_Scaling_Factor",&simCtx->forceScalingFactor,NULL);CHKERRQ(ierr);
434 // --- Group 5
435 LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 5: Solver & Numerics Parameters \n");
436 char mom_solver_type_char[PETSC_MAX_PATH_LEN];
437 PetscBool mom_solver_type_flg = PETSC_FALSE;
438 ierr = PetscOptionsGetString(NULL, NULL, "-mom_solver_type", mom_solver_type_char, sizeof(mom_solver_type_char), &mom_solver_type_flg); CHKERRQ(ierr);
439 ierr = PetscOptionsGetInt(NULL, NULL, "-mom_max_pseudo_steps", &simCtx->mom_max_pseudo_steps, NULL); CHKERRQ(ierr);
440 ierr = PetscOptionsGetReal(NULL, NULL, "-mom_atol", &simCtx->mom_atol, NULL); CHKERRQ(ierr);
441 ierr = PetscOptionsGetReal(NULL, NULL, "-mom_rtol", &simCtx->mom_rtol, NULL); CHKERRQ(ierr);
442 ierr = PetscOptionsGetReal(NULL, NULL, "-imp_stol", &simCtx->imp_stol, NULL); CHKERRQ(ierr);
443 ierr = PetscOptionsGetInt(NULL, NULL, "-central", &simCtx->central, NULL); CHKERRQ(ierr);
444
445 // Keep parser acceptance aligned with the enum and FlowSolver dispatch.
446 if (mom_solver_type_flg) {
447 if(strcmp(mom_solver_type_char, "DUALTIME_PICARD_RK4") == 0) {
449 } else if (strcmp(mom_solver_type_char, "EXPLICIT_RK") == 0) {
451 } else {
452 LOG(GLOBAL, LOG_ERROR, "Invalid value for -mom_solver_type: '%s'. Valid options are: 'DUALTIME_PICARD_RK4', 'EXPLICIT_RK'.\n", mom_solver_type_char);
453 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Invalid value for -mom_solver_type: '%s'.", mom_solver_type_char);
454 }
455 }
456
457 // --- Multigrid Options ---
458 ierr = PetscOptionsGetInt(NULL, NULL, "-mg_level", &simCtx->mglevels, NULL); CHKERRQ(ierr);
459 ierr = PetscOptionsGetInt(NULL, NULL, "-mg_max_it", &simCtx->mg_MAX_IT, NULL); CHKERRQ(ierr);
460 ierr = PetscOptionsGetInt(NULL, NULL, "-mg_idx", &simCtx->mg_idx, NULL); CHKERRQ(ierr);
461 ierr = PetscOptionsGetInt(NULL, NULL, "-mg_pre_it", &simCtx->mg_preItr, NULL); CHKERRQ(ierr);
462 ierr = PetscOptionsGetInt(NULL, NULL, "-mg_post_it", &simCtx->mg_poItr, NULL); CHKERRQ(ierr);
463
464 // --- Other Solver Options ---
465 ierr = PetscOptionsGetInt(NULL, NULL, "-poisson", &simCtx->poisson, NULL); CHKERRQ(ierr);
466 ierr = PetscOptionsGetReal(NULL, NULL, "-poisson_tol", &simCtx->poisson_tol, NULL); CHKERRQ(ierr);
467 ierr = PetscOptionsGetInt(NULL, NULL, "-str", &simCtx->STRONG_COUPLING, NULL); CHKERRQ(ierr);
468 ierr = PetscOptionsGetReal(NULL, NULL, "-ren", &simCtx->ren, NULL); CHKERRQ(ierr);
469 ierr = PetscOptionsGetReal(NULL, NULL, "-pseudo_cfl", &simCtx->pseudo_cfl, NULL); CHKERRQ(ierr);
470 ierr = PetscOptionsGetReal(NULL, NULL, "-max_pseudo_cfl", &simCtx->max_pseudo_cfl, NULL); CHKERRQ(ierr);
471 ierr = PetscOptionsGetReal(NULL, NULL, "-min_pseudo_cfl", &simCtx->min_pseudo_cfl, NULL); CHKERRQ(ierr);
472 ierr = PetscOptionsGetReal(NULL, NULL, "-pseudo_cfl_reduction_factor", &simCtx->pseudo_cfl_reduction_factor, NULL); CHKERRQ(ierr);
473 ierr = PetscOptionsGetReal(NULL, NULL, "-pseudo_cfl_growth_factor", &simCtx->pseudo_cfl_growth_factor, NULL); CHKERRQ(ierr);
474 ierr = PetscOptionsGetReal(NULL,NULL, "-mom_dt_rk4_residual_norm_noise_allowance_factor",&simCtx->mom_dt_rk4_residual_norm_noise_allowance_factor,NULL);CHKERRQ(ierr);
475 ierr = PetscOptionsHasName(NULL, NULL, "-ps_ksp_pic_monitor_true_residual", &simCtx->ps_ksp_pic_monitor_true_residual); CHKERRQ(ierr);
476 ierr = PetscOptionsGetInt(NULL, NULL, "-finit", &simCtx->FieldInitialization, NULL); CHKERRQ(ierr);
477 ierr = PetscOptionsGetReal(NULL, NULL, "-ucont_x", &simCtx->InitialConstantContra.x, NULL); CHKERRQ(ierr);
478 ierr = PetscOptionsGetReal(NULL, NULL, "-ucont_y", &simCtx->InitialConstantContra.y, NULL); CHKERRQ(ierr);
479 ierr = PetscOptionsGetReal(NULL, NULL, "-ucont_z", &simCtx->InitialConstantContra.z, NULL); CHKERRQ(ierr);
480 ierr = PetscOptionsGetReal(NULL, NULL, "-analytical_uniform_u", &simCtx->AnalyticalUniformVelocity.x, NULL); CHKERRQ(ierr);
481 ierr = PetscOptionsGetReal(NULL, NULL, "-analytical_uniform_v", &simCtx->AnalyticalUniformVelocity.y, NULL); CHKERRQ(ierr);
482 ierr = PetscOptionsGetReal(NULL, NULL, "-analytical_uniform_w", &simCtx->AnalyticalUniformVelocity.z, NULL); CHKERRQ(ierr);
483 PetscBool verification_scalar_value_set = PETSC_FALSE;
484 PetscBool verification_scalar_phi0_set = PETSC_FALSE;
485 PetscBool verification_scalar_slope_x_set = PETSC_FALSE;
486 PetscBool verification_scalar_amplitude_set = PETSC_FALSE;
487 PetscBool verification_scalar_kx_set = PETSC_FALSE;
488 PetscBool verification_scalar_ky_set = PETSC_FALSE;
489 PetscBool verification_scalar_kz_set = PETSC_FALSE;
490 ierr = PetscOptionsGetString(NULL, NULL, "-verification_diffusivity_mode",
492 sizeof(simCtx->verificationDiffusivity.mode), NULL); CHKERRQ(ierr);
493 ierr = PetscOptionsGetString(NULL, NULL, "-verification_diffusivity_profile",
495 sizeof(simCtx->verificationDiffusivity.profile), NULL); CHKERRQ(ierr);
496 ierr = PetscOptionsGetReal(NULL, NULL, "-verification_diffusivity_gamma0",
497 &simCtx->verificationDiffusivity.gamma0, NULL); CHKERRQ(ierr);
498 ierr = PetscOptionsGetReal(NULL, NULL, "-verification_diffusivity_slope_x",
499 &simCtx->verificationDiffusivity.slope_x, NULL); CHKERRQ(ierr);
500 ierr = PetscOptionsGetString(NULL, NULL, "-verification_scalar_mode",
501 simCtx->verificationScalar.mode,
502 sizeof(simCtx->verificationScalar.mode), NULL); CHKERRQ(ierr);
503 ierr = PetscOptionsGetString(NULL, NULL, "-verification_scalar_profile",
505 sizeof(simCtx->verificationScalar.profile), NULL); CHKERRQ(ierr);
506 ierr = PetscOptionsGetReal(NULL, NULL, "-verification_scalar_value",
507 &simCtx->verificationScalar.value, &verification_scalar_value_set); CHKERRQ(ierr);
508 ierr = PetscOptionsGetReal(NULL, NULL, "-verification_scalar_phi0",
509 &simCtx->verificationScalar.phi0, &verification_scalar_phi0_set); CHKERRQ(ierr);
510 ierr = PetscOptionsGetReal(NULL, NULL, "-verification_scalar_slope_x",
511 &simCtx->verificationScalar.slope_x, &verification_scalar_slope_x_set); CHKERRQ(ierr);
512 ierr = PetscOptionsGetReal(NULL, NULL, "-verification_scalar_amplitude",
513 &simCtx->verificationScalar.amplitude, &verification_scalar_amplitude_set); CHKERRQ(ierr);
514 ierr = PetscOptionsGetReal(NULL, NULL, "-verification_scalar_kx",
515 &simCtx->verificationScalar.kx, &verification_scalar_kx_set); CHKERRQ(ierr);
516 ierr = PetscOptionsGetReal(NULL, NULL, "-verification_scalar_ky",
517 &simCtx->verificationScalar.ky, &verification_scalar_ky_set); CHKERRQ(ierr);
518 ierr = PetscOptionsGetReal(NULL, NULL, "-verification_scalar_kz",
519 &simCtx->verificationScalar.kz, &verification_scalar_kz_set); CHKERRQ(ierr);
521 (PetscBool)(simCtx->verificationDiffusivity.mode[0] != '\0' ||
522 simCtx->verificationDiffusivity.profile[0] != '\0');
524 (PetscBool)(simCtx->verificationScalar.mode[0] != '\0' ||
525 simCtx->verificationScalar.profile[0] != '\0');
526 if (simCtx->verificationDiffusivity.enabled) {
527 if (strcmp(simCtx->eulerianSource, "analytical") != 0) {
528 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE,
529 "verification diffusivity overrides require -euler_field_source \"analytical\".");
530 }
531 if (strcmp(simCtx->verificationDiffusivity.mode, "analytical") != 0) {
532 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
533 "Unsupported -verification_diffusivity_mode '%s'. Only 'analytical' is supported.",
535 }
536 if (strcmp(simCtx->verificationDiffusivity.profile, "LINEAR_X") != 0) {
537 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
538 "Unsupported -verification_diffusivity_profile '%s'. Only 'LINEAR_X' is supported.",
540 }
541 }
542 if (simCtx->verificationScalar.enabled) {
543 if (strcmp(simCtx->eulerianSource, "analytical") != 0) {
544 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE,
545 "verification scalar overrides require -euler_field_source \"analytical\".");
546 }
547 if (strcmp(simCtx->verificationScalar.mode, "analytical") != 0) {
548 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
549 "Unsupported -verification_scalar_mode '%s'. Only 'analytical' is supported.",
550 simCtx->verificationScalar.mode);
551 }
552 if (strcmp(simCtx->verificationScalar.profile, "CONSTANT") == 0) {
553 if (!verification_scalar_value_set) {
554 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
555 "verification scalar profile CONSTANT requires -verification_scalar_value.");
556 }
557 } else if (strcmp(simCtx->verificationScalar.profile, "LINEAR_X") == 0) {
558 if (!verification_scalar_phi0_set || !verification_scalar_slope_x_set) {
559 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
560 "verification scalar profile LINEAR_X requires -verification_scalar_phi0 and -verification_scalar_slope_x.");
561 }
562 } else if (strcmp(simCtx->verificationScalar.profile, "SIN_PRODUCT") == 0) {
563 if (!verification_scalar_amplitude_set || !verification_scalar_kx_set ||
564 !verification_scalar_ky_set || !verification_scalar_kz_set) {
565 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
566 "verification scalar profile SIN_PRODUCT requires -verification_scalar_amplitude, -verification_scalar_kx, -verification_scalar_ky, and -verification_scalar_kz.");
567 }
568 } else {
569 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
570 "Unsupported -verification_scalar_profile '%s'. Supported profiles: CONSTANT, LINEAR_X, SIN_PRODUCT.",
572 }
573 }
574 // NOTE: cdisx,cdisy,cdisz haven't been parsed, add if necessary.
575
576 // --- Group 6
577 LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 6: Physical & Geometric Parameters \n");
578 ierr = PetscOptionsGetReal(NULL,NULL,"-schmidt_number",&simCtx->schmidt_number,NULL);CHKERRQ(ierr);
579 ierr = PetscOptionsGetReal(NULL,NULL,"-turb_schmidt_number",&simCtx->Turbulent_schmidt_number,NULL);CHKERRQ(ierr);
580 ierr = PetscOptionsGetInt(NULL, NULL, "-no_of_bodies", &simCtx->NumberOfBodies, NULL); CHKERRQ(ierr);
581 ierr = PetscOptionsGetReal(NULL,NULL,"-wall_roughness",&simCtx->wall_roughness_height,NULL);CHKERRQ(ierr);
582 // NOTE: angle is not parsed in the original code, it set programmatically. We will follow that.
583 // NOTE: max_angle is calculated based on other flags (like MHV) in the legacy code.
584 // We will defer that logic to a later setup stage and not parse them directly.
585 // The Scaling Information is calculated here
586 ierr = ParseScalingInformation(simCtx); CHKERRQ(ierr);
587
588 // --- Group 7
589 LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 7: Grid, Domain, and Boundary Condition Settings \n");
590 ierr = PetscOptionsGetInt(NULL, NULL, "-nblk", &simCtx->block_number, NULL); CHKERRQ(ierr); // This is also a modern option
591 ierr = PetscOptionsGetInt(NULL, NULL, "-inlet", &simCtx->inletprofile, NULL); CHKERRQ(ierr);
592 ierr = PetscOptionsGetInt(NULL, NULL, "-Ogrid", &simCtx->Ogrid, NULL); CHKERRQ(ierr);
593 // NOTE: channelz was not parsed, likely set programmatically. We will omit its parsing call.
594 ierr = PetscOptionsGetInt(NULL, NULL, "-grid1d", &simCtx->grid1d, NULL); CHKERRQ(ierr);
595 ierr = PetscOptionsGetBool(NULL, NULL, "-grid", &simCtx->generate_grid, NULL); CHKERRQ(ierr);
596 ierr = PetscOptionsGetString(NULL, NULL, "-grid_file", simCtx->grid_file, PETSC_MAX_PATH_LEN, NULL); CHKERRQ(ierr);
597 ierr = PetscOptionsGetInt(NULL, NULL, "-da_processors_x", &simCtx->da_procs_x, NULL); CHKERRQ(ierr);
598 ierr = PetscOptionsGetInt(NULL, NULL, "-da_processors_y", &simCtx->da_procs_y, NULL); CHKERRQ(ierr);
599 ierr = PetscOptionsGetInt(NULL, NULL, "-da_processors_z", &simCtx->da_procs_z, NULL); CHKERRQ(ierr);
600 ierr = PetscOptionsGetInt(NULL, NULL, "-i_periodic", &simCtx->i_periodic, NULL); CHKERRQ(ierr);
601 ierr = PetscOptionsGetInt(NULL, NULL, "-j_periodic", &simCtx->j_periodic, NULL); CHKERRQ(ierr);
602 ierr = PetscOptionsGetInt(NULL, NULL, "-k_periodic", &simCtx->k_periodic, NULL); CHKERRQ(ierr);
603 ierr = PetscOptionsGetInt(NULL, NULL, "-pbc_domain", &simCtx->blkpbc, NULL); CHKERRQ(ierr);
604 // NOTE: pseudo_periodic was not parsed. We will omit its parsing call.
605 ierr = PetscOptionsGetReal(NULL, NULL, "-grid_rotation_angle", &simCtx->grid_rotation_angle, NULL); CHKERRQ(ierr);
606 ierr = PetscOptionsGetReal(NULL, NULL, "-Croty", &simCtx->Croty, NULL); CHKERRQ(ierr);
607 ierr = PetscOptionsGetReal(NULL, NULL, "-Crotz", &simCtx->Crotz, NULL); CHKERRQ(ierr);
608 PetscBool bcs_flg;
609 char file_list_str[PETSC_MAX_PATH_LEN * 10]; // Buffer for comma-separated list
610
611 ierr = PetscOptionsGetString(NULL, NULL, "-bcs_files", file_list_str, sizeof(file_list_str), &bcs_flg); CHKERRQ(ierr);
612 ierr = PetscOptionsGetReal(NULL, NULL, "-U_bc", &simCtx->U_bc, NULL); CHKERRQ(ierr);
613
614 if (bcs_flg) {
615 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Found -bcs_files option, overriding default.\n");
616
617 // A. Clean up the default memory we allocated in Phase 1.
618 ierr = PetscFree(simCtx->bcs_files[0]); CHKERRQ(ierr);
619 ierr = PetscFree(simCtx->bcs_files); CHKERRQ(ierr);
620 simCtx->num_bcs_files = 0;
621 simCtx->bcs_files = NULL;
622
623 // B. Parse the user-provided comma-separated list.
624 char *token;
625 char *str_copy;
626 ierr = PetscStrallocpy(file_list_str, &str_copy); CHKERRQ(ierr);
627
628 // First pass: count the number of files.
629 token = strtok(str_copy, ",");
630 while (token) {
631 simCtx->num_bcs_files++;
632 token = strtok(NULL, ",");
633 }
634 ierr = PetscFree(str_copy); CHKERRQ(ierr);
635
636 // Second pass: allocate memory and store the filenames.
637 ierr = PetscMalloc1(simCtx->num_bcs_files, &simCtx->bcs_files); CHKERRQ(ierr);
638 ierr = PetscStrallocpy(file_list_str, &str_copy); CHKERRQ(ierr);
639 token = strtok(str_copy, ",");
640 for (PetscInt i = 0; i < simCtx->num_bcs_files; i++) {
641 ierr = PetscStrallocpy(token, &simCtx->bcs_files[i]); CHKERRQ(ierr);
642 token = strtok(NULL, ",");
643 }
644 ierr = PetscFree(str_copy); CHKERRQ(ierr);
645 }
646
647
648 // --- Group 8
649 LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 8: Turbulence Modeling (LES/RANS) \n");
650 PetscInt temp_les_model;
651 ierr = PetscOptionsGetInt(NULL, NULL, "-les", &temp_les_model, NULL); CHKERRQ(ierr);
652 simCtx->les = (LESModelType)temp_les_model;
653 ierr = PetscOptionsGetInt(NULL, NULL, "-rans", &simCtx->rans, NULL); CHKERRQ(ierr);
654 ierr = PetscOptionsGetInt(NULL, NULL, "-wallfunction", &simCtx->wallfunction, NULL); CHKERRQ(ierr);
655 ierr = PetscOptionsGetInt(NULL, NULL, "-mixed", &simCtx->mixed, NULL); CHKERRQ(ierr);
656 ierr = PetscOptionsGetInt(NULL, NULL, "-clark", &simCtx->clark, NULL); CHKERRQ(ierr);
657 ierr = PetscOptionsGetInt(NULL, NULL, "-dynamic_freq", &simCtx->dynamic_freq, NULL); CHKERRQ(ierr);
658 ierr = PetscOptionsGetReal(NULL, NULL, "-max_cs", &simCtx->max_cs, NULL); CHKERRQ(ierr);
659 ierr = PetscOptionsGetReal(NULL, NULL, "-const_cs", &simCtx->Const_CS, NULL); CHKERRQ(ierr);
660 ierr = PetscOptionsGetInt(NULL, NULL, "-testfilter_ik", &simCtx->testfilter_ik, NULL); CHKERRQ(ierr);
661 ierr = PetscOptionsGetInt(NULL, NULL, "-testfilter_1d", &simCtx->testfilter_1d, NULL); CHKERRQ(ierr);
662 ierr = PetscOptionsGetInt(NULL, NULL, "-i_homo_filter", &simCtx->i_homo_filter, NULL); CHKERRQ(ierr);
663 ierr = PetscOptionsGetInt(NULL, NULL, "-j_homo_filter", &simCtx->j_homo_filter, NULL); CHKERRQ(ierr);
664 ierr = PetscOptionsGetInt(NULL, NULL, "-k_homo_filter", &simCtx->k_homo_filter, NULL); CHKERRQ(ierr);
665 ierr = PetscOptionsGetBool(NULL, NULL, "-averaging", &simCtx->averaging, NULL); CHKERRQ(ierr);
666
667 // --- Group 9
668 LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 9: Particle / DMSwarm Data & Settings \n");
669 ierr = PetscOptionsGetInt(NULL, NULL, "-numParticles", &simCtx->np, NULL); CHKERRQ(ierr);
670 ierr = PetscOptionsGetBool(NULL, NULL, "-read_fields", &simCtx->readFields, NULL); CHKERRQ(ierr);
671 PetscInt temp_pinit = (PetscInt)PARTICLE_INIT_SURFACE_RANDOM;
672 ierr = PetscOptionsGetInt(NULL, NULL, "-pinit", &temp_pinit, NULL); CHKERRQ(ierr);
674 PetscInt temp_interp = (PetscInt)INTERP_TRILINEAR;
675 ierr = PetscOptionsGetInt(NULL, NULL, "-interpolation_method", &temp_interp, NULL); CHKERRQ(ierr);
676 simCtx->interpolationMethod = (InterpolationMethod)temp_interp;
677 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Interpolation method: %s\n",
678 simCtx->interpolationMethod == INTERP_TRILINEAR ? "Trilinear (direct cell-center)" : "CornerAveraged (legacy)");
679 ierr = PetscOptionsGetReal(NULL, NULL, "-psrc_x", &simCtx->psrc_x, NULL); CHKERRQ(ierr);
680 ierr = PetscOptionsGetReal(NULL, NULL, "-psrc_y", &simCtx->psrc_y, NULL); CHKERRQ(ierr);
681 ierr = PetscOptionsGetReal(NULL, NULL, "-psrc_z", &simCtx->psrc_z, NULL); CHKERRQ(ierr);
682 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Particle initialization mode: %s. Point source: (%.6f, %.6f, %.6f)\n",
684 simCtx->psrc_x, simCtx->psrc_y, simCtx->psrc_z);
685 ierr = PetscOptionsGetString(NULL,NULL,"-particle_restart_mode",simCtx->particleRestartMode,sizeof(simCtx->particleRestartMode),NULL); CHKERRQ(ierr);
686 // Validation for Particle Restart Mode
687 if (strcmp(simCtx->particleRestartMode, "load") != 0 && strcmp(simCtx->particleRestartMode, "init") != 0) {
688 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Invalid value for -particle_restart_mode. Must be 'load' or 'init'. You provided '%s'.", simCtx->particleRestartMode);
689 }
690 ierr = InitializeBrownianRNG(simCtx); CHKERRQ(ierr);
691 // --- Group 10
692 LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 10: Immersed Boundary & FSI Data Object Pointers \n");
693 ierr = PetscOptionsGetBool(NULL, NULL, "-rs_fsi", &simCtx->rstart_fsi, NULL); CHKERRQ(ierr);
694 ierr = PetscOptionsGetInt(NULL, NULL, "-duplicate", &simCtx->duplicate, NULL); CHKERRQ(ierr);
695
696 // --- Group 11
697 LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 11: Top-Level Managers & Custom Configuration \n");
698 ierr = PetscOptionsGetInt(NULL, NULL, "-logfreq", &simCtx->LoggingFrequency, NULL); CHKERRQ(ierr);
699
700 if (simCtx->num_bcs_files != simCtx->block_number) {
701 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);
702 }
703
704 // --- Group 12
705 LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 12: Post-Processing Information.\n");
706 // This logic determines the Post Processing configuration and STORES it in simCtx for later reference and cleanup.
707 ierr = PetscOptionsGetString(NULL,NULL,"-postprocessing_config_file",simCtx->PostprocessingControlFile,PETSC_MAX_PATH_LEN,NULL); CHKERRQ(ierr);
708 /* Parse post settings for both solver and post-processor binaries using the single pre-allocated pps object. */
709 ierr = ParsePostProcessingSettings(simCtx);
710
711 // === 5. Dependent Parameter Calculations ================================
712 // Some parameters depend on others, so we calculate them here.
713 simCtx->StartTime = (PetscReal)simCtx->StartStep*simCtx->dt;
714 simCtx->ti = simCtx->StartTime;
715 simCtx->step = simCtx->StartStep;
716
717 // === 5. Log Summary and Finalize Setup ==================================
718 LOG_ALLOW(GLOBAL, LOG_DEBUG, "-- Console Output Functions [Total : %d] : --\n", simCtx->nAllowed);
719 for (PetscInt i = 0; i < simCtx->nAllowed; ++i) {
720 LOG_ALLOW(GLOBAL, LOG_DEBUG, " [%2d] «%s»\n", i, simCtx->allowedFuncs[i]);
721 }
722
723 LOG_ALLOW(GLOBAL, LOG_INFO, "Configuration complete. Key parameters:\n");
724 LOG_ALLOW(GLOBAL, LOG_INFO, " - Run mode: %s\n", simCtx->OnlySetup ? "SETUP ONLY" : "Full Simulation");
725 LOG_ALLOW(GLOBAL, LOG_INFO, " - Time steps: %d (from %d to %d)\n", simCtx->StepsToRun, simCtx->StartStep, simCtx->StartStep + simCtx->StepsToRun);
726 LOG_ALLOW(GLOBAL, LOG_INFO, " - Time step size (dt): %g\n", simCtx->dt);
727 if (simCtx->tiout > 0) {
728 LOG_ALLOW(GLOBAL, LOG_INFO, " - Field/restart output cadence: every %d step(s)\n", simCtx->tiout);
729 } else {
730 LOG_ALLOW(GLOBAL, LOG_INFO, " - Field/restart output cadence: DISABLED\n");
731 }
732 LOG_ALLOW(GLOBAL, LOG_INFO, " - Immersed Boundary: %s\n", simCtx->immersed ? "ENABLED" : "DISABLED");
733 LOG_ALLOW(GLOBAL, LOG_INFO, " - Particles: %d\n", simCtx->np);
734 if (simCtx->np > 0) {
735 if (simCtx->particleConsoleOutputFreq > 0) {
736 LOG_ALLOW(GLOBAL, LOG_INFO, " - Particle console cadence: every %d step(s)\n", simCtx->particleConsoleOutputFreq);
737 } else {
738 LOG_ALLOW(GLOBAL, LOG_INFO, " - Particle console cadence: DISABLED\n");
739 }
740 LOG_ALLOW(GLOBAL, LOG_INFO, " - Particle console row subsampling: every %d particle(s)\n", simCtx->LoggingFrequency);
741 }
742 if (simCtx->StartStep > 0 && simCtx->np > 0) {
743 LOG_ALLOW(GLOBAL, LOG_INFO, " - Particle Restart Mode: %s\n", simCtx->particleRestartMode);
744 }
745
746 // --- Initialize PETSc's internal performance logging stage ---
747 ierr = PetscLogDefaultBegin(); CHKERRQ(ierr); // REDUNDANT but safe.
748
749 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Finished CreateSimulationContext successfully on rank %d.\n", simCtx->rank);
750
752 PetscFunctionReturn(0);
753}
PetscErrorCode ParsePostProcessingSettings(SimCtx *simCtx)
Initializes post-processing settings from a config file and command-line overrides.
Definition io.c:2209
PetscErrorCode ParseScalingInformation(SimCtx *simCtx)
Parses physical scaling parameters from command-line options.
Definition io.c:2357
void set_allowed_functions(const char **functionList, int count)
Sets the global list of function names that are allowed to log.
Definition logging.c:152
#define LOG_ALLOW_SYNC(scope, level, fmt,...)
Synchronized logging macro that checks both the log level and whether the calling function is in the ...
Definition logging.h:252
#define LOCAL
Logging scope definitions for controlling message output.
Definition logging.h:44
#define GLOBAL
Scope for global logging across all processes.
Definition logging.h:45
#define LOG_ALLOW(scope, level, fmt,...)
Logging macro that checks both the log level and whether the calling function is in the allowed-funct...
Definition logging.h:199
PetscErrorCode print_log_level(void)
Prints the current logging level to the console.
Definition logging.c:116
#define PROFILE_FUNCTION_END
Marks the end of a profiled code block.
Definition logging.h:779
#define LOG(scope, level, fmt,...)
Logging macro for PETSc-based applications with scope control.
Definition logging.h:83
PetscErrorCode LoadAllowedFunctionsFromFile(const char filename[], char ***funcsOut, PetscInt *nOut)
Load function names from a text file.
Definition logging.c:596
PetscErrorCode ProfilingInitialize(SimCtx *simCtx)
Initializes the custom profiling system using configuration from SimCtx.
Definition logging.c:1044
@ LOG_ERROR
Critical errors that may halt the program.
Definition logging.h:28
@ LOG_INFO
Informational messages about program execution.
Definition logging.h:30
@ LOG_WARNING
Non-critical issues that warrant attention.
Definition logging.h:29
@ LOG_DEBUG
Detailed debugging information.
Definition logging.h:31
#define PROFILE_FUNCTION_BEGIN
Marks the beginning of a profiled code block (typically a function).
Definition logging.h:770
const char * ParticleInitializationToString(ParticleInitializationType ParticleInitialization)
Helper function to convert ParticleInitialization to a string representation.
Definition logging.c:703
PetscErrorCode InitializeBrownianRNG(SimCtx *simCtx)
Internal helper implementation: InitializeBrownianRNG().
Definition setup.c:2769
PetscBool RuntimeWalltimeGuardParsePositiveSeconds(const char *text, PetscReal *seconds_out)
Implementation of RuntimeWalltimeGuardParsePositiveSeconds().
Definition setup.c:18
LESModelType
Identifies the six logical faces of a structured computational block.
Definition variables.h:488
@ NO_LES_MODEL
Definition variables.h:489
PetscInt MHV
Definition variables.h:679
PetscInt turbine
Definition variables.h:679
PetscInt fishcyl
Definition variables.h:679
PetscInt clark
Definition variables.h:733
PetscInt movefsi
Definition variables.h:673
PetscBool continueMode
Definition variables.h:660
PetscInt moveframe
Definition variables.h:674
PetscInt TwoD
Definition variables.h:674
PetscInt pseudo_periodic
Definition variables.h:713
PetscInt fish_c
Definition variables.h:679
PetscInt dgf_z
Definition variables.h:675
PetscReal poisson_tol
Definition variables.h:688
PetscBool profilingFinalSummary
Definition variables.h:780
PetscReal schmidt_number
Definition variables.h:709
PetscMPIInt rank
Definition variables.h:646
char profilingTimestepFile[PETSC_MAX_PATH_LEN]
Definition variables.h:779
PetscInt fish
Definition variables.h:679
PetscInt LV
Definition variables.h:679
PetscReal angle
Definition variables.h:704
PetscReal Turbulent_schmidt_number
Definition variables.h:709
PetscInt64 searchLocatedCount
Definition variables.h:224
PetscInt thin
Definition variables.h:674
PetscInt grid1d
Definition variables.h:712
PetscInt block_number
Definition variables.h:712
PetscReal mom_rtol
Definition variables.h:685
PetscInt64 searchLostCount
Definition variables.h:225
PetscInt da_procs_z
Definition variables.h:718
PetscInt blkpbc
Definition variables.h:713
PetscInt sediment
Definition variables.h:673
PetscReal targetVolumetricFlux
Definition variables.h:724
PetscBool walltimeGuardActive
Definition variables.h:782
PetscInt channelz
Definition variables.h:680
char euler_subdir[PETSC_MAX_PATH_LEN]
Definition variables.h:666
PetscReal walltimeGuardWarmupTotalSeconds
Definition variables.h:790
PetscReal forceScalingFactor
Definition variables.h:723
PetscReal pseudo_cfl_reduction_factor
Definition variables.h:692
PetscInt rans
Definition variables.h:732
ParticleInitializationType
Enumerator to identify the particle initialization strategy.
Definition variables.h:508
@ PARTICLE_INIT_SURFACE_RANDOM
Random placement on the inlet face.
Definition variables.h:509
PetscReal StartTime
Definition variables.h:657
PetscInt dgf_az
Definition variables.h:675
PetscReal FluxOutSum
Definition variables.h:721
PetscBool walltimeGuardHasEWMA
Definition variables.h:792
PetscReal CMy_c
Definition variables.h:705
char ** bcs_files
Definition variables.h:720
PetscReal max_angle
Definition variables.h:704
PetscReal min_pseudo_cfl
Definition variables.h:693
PetscInt64 boundaryClampCount
Definition variables.h:231
PetscInt particlesLostLastStep
Definition variables.h:746
PetscInt duplicate
Definition variables.h:761
PetscInt tiout
Definition variables.h:655
PetscReal walltimeGuardMinSeconds
Definition variables.h:785
char allowedFile[PETSC_MAX_PATH_LEN]
Definition variables.h:765
PetscInt da_procs_y
Definition variables.h:718
PetscInt64 traversalStepsSum
Definition variables.h:226
PetscInt testfilter_1d
Definition variables.h:735
PetscReal psrc_x
Definition variables.h:706
PetscInt FieldInitialization
Definition variables.h:696
PetscReal ren
Definition variables.h:691
PetscReal Crotz
Definition variables.h:716
PetscInt mixed
Definition variables.h:733
IBMVNodes * ibmv
Definition variables.h:758
PetscInt64 searchPopulation
Definition variables.h:223
char output_dir[PETSC_MAX_PATH_LEN]
Definition variables.h:665
PetscReal walltimeGuardLatestStepSeconds
Definition variables.h:794
PetscReal dt
Definition variables.h:658
PetscInt occupiedCellCount
Definition variables.h:750
PetscInt StepsToRun
Definition variables.h:654
char profilingTimestepMode[32]
Definition variables.h:778
PetscInt k_periodic
Definition variables.h:713
PetscInt inletprofile
Definition variables.h:712
PetscReal cdisy
Definition variables.h:691
PetscReal mom_atol
Definition variables.h:685
PetscBool rstart_fsi
Definition variables.h:760
PetscInt currentSettlementPass
Definition variables.h:235
PetscInt np
Definition variables.h:739
PetscBool averaging
Definition variables.h:736
PetscReal C_IEM
Definition variables.h:754
PetscInt ccc
Definition variables.h:728
PetscReal ratio
Definition variables.h:729
PetscInt mg_idx
Definition variables.h:686
PetscInt StartStep
Definition variables.h:653
PetscInt mg_MAX_IT
Definition variables.h:686
PetscBool OnlySetup
Definition variables.h:659
PetscInt rotatefsi
Definition variables.h:673
@ MOMENTUM_SOLVER_EXPLICIT_RK
Definition variables.h:503
@ MOMENTUM_SOLVER_DUALTIME_PICARD_RK4
Definition variables.h:504
PetscReal cdisz
Definition variables.h:691
PetscScalar x
Definition variables.h:101
PetscInt64 reSearchCount
Definition variables.h:227
PetscInt dgf_x
Definition variables.h:675
char * current_io_directory
Definition variables.h:670
PetscInt pizza
Definition variables.h:679
PetscReal MaxDiv
Definition variables.h:771
char grid_file[PETSC_MAX_PATH_LEN]
Definition variables.h:717
PetscReal max_cs
Definition variables.h:734
PetscInt invicid
Definition variables.h:674
char ** allowedFuncs
Definition variables.h:767
PetscInt64 bboxGuessFallbackCount
Definition variables.h:233
InterpolationMethod interpolationMethod
Definition variables.h:744
PetscReal psrc_z
Point source location for PARTICLE_INIT_POINT_SOURCE.
Definition variables.h:706
PetscInt mg_poItr
Definition variables.h:686
PetscInt STRONG_COUPLING
Definition variables.h:689
VerificationScalarConfig verificationScalar
Definition variables.h:700
PetscReal max_pseudo_cfl
Definition variables.h:693
PetscInt MaxDivx
Definition variables.h:772
PetscInt poisson
Definition variables.h:687
PetscInt k_homo_filter
Definition variables.h:735
char profilingSelectedFuncsFile[PETSC_MAX_PATH_LEN]
Definition variables.h:774
PetscInt MaxDivy
Definition variables.h:772
PetscInt NumberOfBodies
Definition variables.h:703
char particleRestartMode[16]
Definition variables.h:745
PetscInt Ogrid
Definition variables.h:712
PetscInt64 bboxGuessSuccessCount
Definition variables.h:232
char particle_subdir[PETSC_MAX_PATH_LEN]
Definition variables.h:667
PetscInt MaxDivz
Definition variables.h:772
BoundingBox * bboxlist
Definition variables.h:742
PetscInt j_homo_filter
Definition variables.h:735
PetscInt eel
Definition variables.h:679
char log_dir[PETSC_MAX_PATH_LEN]
Definition variables.h:668
PetscInt MaxDivFlatArg
Definition variables.h:772
PetscReal FluxInSum
Definition variables.h:721
PetscInt walltimeGuardCompletedSteps
Definition variables.h:789
PetscInt64 maxParticlePassDepth
Definition variables.h:234
PetscReal CMz_c
Definition variables.h:705
PetscInt64 maxTraversalSteps
Definition variables.h:228
PetscBool generate_grid
Definition variables.h:714
Cmpnts AnalyticalUniformVelocity
Definition variables.h:698
char eulerianSource[PETSC_MAX_PATH_LEN]
Definition variables.h:663
PetscReal imp_stol
Definition variables.h:685
PetscInt nAllowed
Definition variables.h:768
PetscBool walltimeGuardEnabled
Definition variables.h:781
PetscReal wall_roughness_height
Definition variables.h:708
PetscBool useProfilingSelectedFuncsCfg
Definition variables.h:775
PetscInt walltimeGuardWarmupSteps
Definition variables.h:783
ParticleInitializationType ParticleInitialization
Definition variables.h:743
PetscScalar z
Definition variables.h:101
InterpolationMethod
Selects the grid-to-particle interpolation method.
Definition variables.h:521
@ INTERP_TRILINEAR
Definition variables.h:522
PetscReal Const_CS
Definition variables.h:734
PetscInt i_homo_filter
Definition variables.h:735
PetscInt wallfunction
Definition variables.h:733
PetscInt rheology
Definition variables.h:673
PetscReal Flux_in
Definition variables.h:704
char ** profilingSelectedFuncs
Definition variables.h:776
PetscReal cdisx
Definition variables.h:691
PetscInt dgf_ax
Definition variables.h:675
PetscInt mglevels
Definition variables.h:686
PetscInt num_bcs_files
Definition variables.h:719
DM dm_swarm
Definition variables.h:741
PetscBool useCfg
Definition variables.h:766
PetscReal psrc_y
Definition variables.h:706
PetscBool readFields
Definition variables.h:740
PetscInt central
Definition variables.h:689
PetscReal Fluxsum
Definition variables.h:721
PetscReal pseudo_cfl_growth_factor
Definition variables.h:692
PetscReal Croty
Definition variables.h:716
PetscInt particlesLostCumulative
Definition variables.h:747
PetscInt nProfilingSelectedFuncs
Definition variables.h:777
PetscReal mom_dt_rk4_residual_norm_noise_allowance_factor
Definition variables.h:694
PetscInt particlesMigratedLastStep
Definition variables.h:749
PetscReal grid_rotation_angle
Definition variables.h:715
PetscInt dynamic_freq
Definition variables.h:733
char AnalyticalSolutionType[PETSC_MAX_PATH_LEN]
Definition variables.h:676
PetscInt da_procs_x
Definition variables.h:718
PetscReal U_bc
Definition variables.h:727
PetscReal walltimeGuardWarmupAverageSeconds
Definition variables.h:791
PetscInt particleConsoleOutputFreq
Definition variables.h:656
Cmpnts InitialConstantContra
Definition variables.h:697
SearchMetricsState searchMetrics
Definition variables.h:752
PetscInt i_periodic
Definition variables.h:713
PetscInt step
Definition variables.h:651
PetscReal walltimeGuardEWMASeconds
Definition variables.h:793
PetscReal AreaOutSum
Definition variables.h:726
PetscInt dgf_ay
Definition variables.h:675
PetscInt mom_max_pseudo_steps
Definition variables.h:684
PetscRandom BrownianMotionRNG
Definition variables.h:753
PetscInt testfilter_ik
Definition variables.h:735
PetscInt hydro
Definition variables.h:679
PostProcessParams * pps
Definition variables.h:798
PetscInt migrationPassesLastStep
Definition variables.h:748
PetscScalar y
Definition variables.h:101
PetscMPIInt size
Definition variables.h:647
char _io_context_buffer[PETSC_MAX_PATH_LEN]
Definition variables.h:669
PetscReal walltimeGuardLimitSeconds
Definition variables.h:788
PetscBool ps_ksp_pic_monitor_true_residual
Definition variables.h:695
PetscReal walltimeGuardEstimatorAlpha
Definition variables.h:786
PetscInt les
Definition variables.h:732
PetscInt mg_preItr
Definition variables.h:686
PetscViewer logviewer
Definition variables.h:661
PetscInt64 searchAttempts
Definition variables.h:222
PetscInt64 tieBreakCount
Definition variables.h:230
PetscInt cop
Definition variables.h:679
PetscReal ti
Definition variables.h:652
PetscReal walltimeGuardMultiplier
Definition variables.h:784
PetscInt Pipe
Definition variables.h:679
PetscInt rotateframe
Definition variables.h:674
IBMNodes * ibm
Definition variables.h:757
PetscReal AreaInSum
Definition variables.h:726
MomentumSolverType mom_solver_type
Definition variables.h:683
PetscReal summationRHS
Definition variables.h:770
PetscInt immersed
Definition variables.h:673
PetscInt64 maxTraversalFailCount
Definition variables.h:229
char PostprocessingControlFile[PETSC_MAX_PATH_LEN]
Definition variables.h:797
char restart_dir[PETSC_MAX_PATH_LEN]
Definition variables.h:664
VerificationDiffusivityConfig verificationDiffusivity
Definition variables.h:699
PetscInt blank
Definition variables.h:674
PetscInt dgf_y
Definition variables.h:675
PetscReal walltimeGuardJobStartEpochSeconds
Definition variables.h:787
PetscReal pseudo_cfl
Definition variables.h:691
PetscInt LoggingFrequency
Definition variables.h:769
PetscReal CMx_c
Definition variables.h:705
PetscReal drivingForceMagnitude
Definition variables.h:723
PetscReal particleLoadImbalance
Definition variables.h:751
PetscInt j_periodic
Definition variables.h:713
PetscInt wing
Definition variables.h:679
FSInfo * fsi
Definition variables.h:759
The master context for the entire simulation.
Definition variables.h:643
Here is the call graph for this function:
Here is the caller graph for this function:

◆ PetscMkdirRecursive()

static PetscErrorCode PetscMkdirRecursive ( const char *  path)
static

Internal helper implementation: PetscMkdirRecursive().

Local to this translation unit.

Definition at line 761 of file setup.c.

762{
763 PetscErrorCode ierr;
764 char tmp_path[PETSC_MAX_PATH_LEN];
765 char *p = NULL;
766 size_t len;
767 PetscBool exists;
768
769 PetscFunctionBeginUser;
770
771 // Create a mutable copy of the path
772 len = strlen(path);
773 if (len >= sizeof(tmp_path)) {
774 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Path is too long to process: %s", path);
775 }
776 strcpy(tmp_path, path);
777
778 // If the path ends with a separator, remove it
779 if (tmp_path[len - 1] == '/') {
780 tmp_path[len - 1] = 0;
781 }
782
783 // Iterate through the path, creating each directory level
784 for (p = tmp_path + 1; *p; p++) {
785 if (*p == '/') {
786 *p = 0; // Temporarily terminate the string
787
788 // Check if this directory level exists
789 ierr = PetscTestDirectory(tmp_path, 'r', &exists); CHKERRQ(ierr);
790 if (!exists) {
791 ierr = PetscMkdir(tmp_path); CHKERRQ(ierr);
792 }
793
794 *p = '/'; // Restore the separator
795 }
796 }
797
798 // Create the final, full directory path
799 ierr = PetscTestDirectory(tmp_path, 'r', &exists); CHKERRQ(ierr);
800 if (!exists) {
801 ierr = PetscMkdir(tmp_path); CHKERRQ(ierr);
802 }
803
804 PetscFunctionReturn(0);
805}
Here is the caller graph for this function:

◆ SetupSimulationEnvironment()

PetscErrorCode SetupSimulationEnvironment ( SimCtx simCtx)

Internal helper implementation: SetupSimulationEnvironment().

Verifies and prepares the complete I/O environment for a simulation run.

Local to this translation unit.

Definition at line 813 of file setup.c.

814{
815 PetscErrorCode ierr;
816 PetscMPIInt rank;
817 PetscBool exists;
818
819 PetscFunctionBeginUser;
820 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
821
822 LOG_ALLOW(GLOBAL, LOG_INFO, "--- Setting up simulation environment ---\n");
823
824 /* =====================================================================
825 * Phase 1: Check for all required and optional INPUT files.
826 * ===================================================================== */
827 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Phase 1: Verifying input files...\n");
828
829 // --- Mandatory Inputs ---
830 if (!simCtx->generate_grid) {
831 ierr = VerifyPathExistence(simCtx->grid_file, PETSC_FALSE, PETSC_FALSE, "Grid file", &exists); CHKERRQ(ierr);
832 }
833 for (PetscInt i = 0; i < simCtx->num_bcs_files; i++) {
834 char desc[128];
835 ierr = PetscSNPrintf(desc, sizeof(desc), "BCS file #%d", i + 1); CHKERRQ(ierr);
836 ierr = VerifyPathExistence(simCtx->bcs_files[i], PETSC_FALSE, PETSC_FALSE, desc, &exists); CHKERRQ(ierr);
837 }
838
839 // --- Optional Inputs (these produce warnings if missing) ---
840 if (simCtx->useCfg) {
841 ierr = VerifyPathExistence(simCtx->allowedFile, PETSC_FALSE, PETSC_TRUE, "Whitelist config file", &exists); CHKERRQ(ierr);
842 }
843 if (simCtx->useProfilingSelectedFuncsCfg) {
844 ierr = VerifyPathExistence(simCtx->profilingSelectedFuncsFile, PETSC_FALSE, PETSC_TRUE, "Profiling config file", &exists); CHKERRQ(ierr);
845 }
846 if (simCtx->exec_mode == EXEC_MODE_POSTPROCESSOR) {
847 ierr = VerifyPathExistence(simCtx->PostprocessingControlFile, PETSC_FALSE, PETSC_TRUE, "Post-processing control file", &exists); CHKERRQ(ierr);
848 }
849
850
851 /* =====================================================================
852 * Phase 2: Validate directories specific to the execution mode.
853 * ===================================================================== */
854 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Phase 2: Verifying execution mode directories...\n");
855 // The data source directory must exist if we intend to load any data from it.
856 // This is true if:
857 // 1. We are restarting from a previous time step (StartStep > 0), which implies
858 // loading Eulerian fields and/or particle fields.
859 // 2. We are starting from t=0 but are explicitly told to load the initial
860 // Eulerian fields from a file (eulerianSource == "load").
861 if (simCtx->StartStep > 0 || strcmp(simCtx->eulerianSource,"load")== 0){ // If this is a restart run
862 ierr = VerifyPathExistence(simCtx->restart_dir, PETSC_TRUE, PETSC_FALSE, "Restart source directory", &exists); CHKERRQ(ierr);
863 }
864 if (simCtx->exec_mode == EXEC_MODE_POSTPROCESSOR) {
865 ierr = VerifyPathExistence(simCtx->pps->source_dir, PETSC_TRUE, PETSC_FALSE, "Post-processing source directory", &exists); CHKERRQ(ierr);
866 }
867
868 /* =====================================================================
869 * Phase 3: Create and prepare all OUTPUT directories.
870 * ===================================================================== */
871 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Phase 3: Preparing output directories...\n");
872
873 if (rank == 0){
874 if(simCtx->exec_mode == EXEC_MODE_SOLVER){
875 // --- Prepare Log Directory ---
876 if (!simCtx->continueMode) {
877 // Only wipe logs on fresh runs; continue mode appends to existing logs.
878 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Creating/cleaning log directory: %s\n", simCtx->log_dir);
879 ierr = PetscRMTree(simCtx->log_dir); // Wipes the directory and its contents
880 if (ierr) { /* Ignore file-not-found error, but fail on others */
881 PetscError(PETSC_COMM_SELF, __LINE__, __FUNCT__, __FILE__, ierr, PETSC_ERROR_INITIAL, "Could not remove existing log directory '%s'. Check permissions.", simCtx->log_dir);
882 }
883 ierr = PetscMkdir(simCtx->log_dir); CHKERRQ(ierr);
884 } else {
885 // In continue mode, ensure log directory exists but don't wipe it.
886 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Continue mode: preserving existing log directory: %s\n", simCtx->log_dir);
887 ierr = PetscMkdir(simCtx->log_dir); CHKERRQ(ierr);
888 }
889
890 // --- Prepare Output Directories ---
891 char path_buffer[PETSC_MAX_PATH_LEN];
892
893 // 1. Check/Create the main output directory
894 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Verifying main output directory: %s\n", simCtx->output_dir);
895 ierr = PetscTestDirectory(simCtx->output_dir, 'r', &exists); CHKERRQ(ierr);
896 if (!exists) {
897 LOG_ALLOW(GLOBAL, LOG_INFO, "Output directory not found. Creating: %s\n", simCtx->output_dir);
898 ierr = PetscMkdir(simCtx->output_dir); CHKERRQ(ierr);
899 }
900
901 // 2. Check/Create the Eulerian subdirectory
902 ierr = PetscSNPrintf(path_buffer, sizeof(path_buffer), "%s/%s", simCtx->output_dir, simCtx->euler_subdir); CHKERRQ(ierr);
903 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Verifying Eulerian subdirectory: %s\n", path_buffer);
904 ierr = PetscTestDirectory(path_buffer, 'r', &exists); CHKERRQ(ierr);
905 if (!exists) {
906 LOG_ALLOW(GLOBAL, LOG_INFO, "Eulerian subdirectory not found. Creating: %s\n", path_buffer);
907 ierr = PetscMkdir(path_buffer); CHKERRQ(ierr);
908 }
909
910 // 3. Check/Create the Particle subdirectory if needed
911 if (simCtx->np > 0) {
912 ierr = PetscSNPrintf(path_buffer, sizeof(path_buffer), "%s/%s", simCtx->output_dir, simCtx->particle_subdir); CHKERRQ(ierr);
913 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Verifying Particle subdirectory: %s\n", path_buffer);
914 ierr = PetscTestDirectory(path_buffer, 'r', &exists); CHKERRQ(ierr);
915 if (!exists) {
916 LOG_ALLOW(GLOBAL, LOG_INFO, "Particle subdirectory not found. Creating: %s\n", path_buffer);
917 ierr = PetscMkdir(path_buffer); CHKERRQ(ierr);
918 }
919 }
920 } else if(simCtx->exec_mode == EXEC_MODE_POSTPROCESSOR){
921 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Preparing post-processing output directories ...\n");
922
923 PostProcessParams *pps = simCtx->pps;
924 char path_buffer[PETSC_MAX_PATH_LEN];
925
926 const char *last_slash_euler = strrchr(pps->output_prefix, '/');
927 if(last_slash_euler){
928 size_t dir_len = last_slash_euler - pps->output_prefix;
929 if(dir_len > 0){
930 if(dir_len >= sizeof(path_buffer)) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Post-processing output prefix path is too long.");
931 strncpy(path_buffer, pps->output_prefix, dir_len);
932 path_buffer[dir_len] = '\0';
933
934 ierr = PetscTestDirectory(path_buffer, 'r', &exists); CHKERRQ(ierr);
935 if (!exists){
936 LOG_ALLOW(GLOBAL, LOG_INFO, "Creating post-processing Eulerian output directory: %s\n", path_buffer);
937 ierr = PetscMkdirRecursive(path_buffer); CHKERRQ(ierr);
938 }
939 }
940 }
941
942 // Particle output directory
943 if(pps->outputParticles){
944 const char *last_slash_particle = strrchr(pps->particle_output_prefix, '/');
945 if(last_slash_particle){
946 size_t dir_len = last_slash_particle - pps->particle_output_prefix;
947 if(dir_len > 0){
948 if(dir_len > sizeof(path_buffer)) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Post-processing particle output prefix path is too long.");
949 strncpy(path_buffer, pps->particle_output_prefix, dir_len);
950 path_buffer[dir_len] = '\0';
951
952 ierr = PetscTestDirectory(path_buffer, 'r', &exists); CHKERRQ(ierr);
953
954 if (!exists){
955 LOG_ALLOW(GLOBAL, LOG_INFO, "Creating post-processing Particle output directory: %s\n", path_buffer);
956 ierr = PetscMkdirRecursive(path_buffer); CHKERRQ(ierr);
957 }
958 }
959 }
960 }
961
962 // Statistics output directory
963 if(pps->statistics_pipeline[0] != '\0'){
964 const char *last_slash_stats = strrchr(pps->statistics_output_prefix, '/');
965 if(last_slash_stats){
966 size_t dir_len = last_slash_stats - pps->statistics_output_prefix;
967 if(dir_len > 0){
968 if(dir_len >= sizeof(path_buffer)) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Post-processing statistics output prefix path is too long.");
969 strncpy(path_buffer, pps->statistics_output_prefix, dir_len);
970 path_buffer[dir_len] = '\0';
971
972 ierr = PetscTestDirectory(path_buffer, 'r', &exists); CHKERRQ(ierr);
973 if (!exists){
974 LOG_ALLOW(GLOBAL, LOG_INFO, "Creating post-processing Statistics output directory: %s\n", path_buffer);
975 ierr = PetscMkdirRecursive(path_buffer); CHKERRQ(ierr);
976 }
977 }
978 }
979 }
980 }
981 }
982
983 // Synchronize all processes before proceeding
984 ierr = MPI_Barrier(PETSC_COMM_WORLD); CHKERRMPI(ierr);
985
986 LOG_ALLOW(GLOBAL, LOG_INFO, "--- Environment setup complete ---\n");
987
988 PetscFunctionReturn(0);
989}
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:738
static PetscErrorCode PetscMkdirRecursive(const char *path)
Internal helper implementation: PetscMkdirRecursive().
Definition setup.c:761
#define __FUNCT__
Definition setup.c:42
char statistics_output_prefix[256]
basename for CSV output, e.g.
Definition variables.h:575
char particle_output_prefix[256]
Definition variables.h:570
char output_prefix[256]
Definition variables.h:567
char statistics_pipeline[1024]
e.g.
Definition variables.h:574
char source_dir[PETSC_MAX_PATH_LEN]
Definition variables.h:555
PetscBool outputParticles
Definition variables.h:561
@ EXEC_MODE_SOLVER
Definition variables.h:616
@ EXEC_MODE_POSTPROCESSOR
Definition variables.h:617
ExecutionMode exec_mode
Definition variables.h:662
Holds all configuration parameters for a post-processing run.
Definition variables.h:553
Here is the call graph for this function:
Here is the caller graph for this function:

◆ AllocateContextHierarchy()

static PetscErrorCode AllocateContextHierarchy ( SimCtx simCtx)
static

Internal helper implementation: AllocateContextHierarchy().

Local to this translation unit.

Definition at line 997 of file setup.c.

998{
999 PetscErrorCode ierr;
1000 UserMG *usermg = &simCtx->usermg;
1001 MGCtx *mgctx;
1002 PetscInt nblk = simCtx->block_number;
1003 PetscBool found;
1004 PetscFunctionBeginUser;
1006
1007 LOG_ALLOW(GLOBAL, LOG_INFO, "Allocating context hierarchy for %d levels and %d blocks...\n", simCtx->mglevels, nblk);
1008
1009 // Store the number of levels in the UserMG struct itself
1010 usermg->mglevels = simCtx->mglevels;
1011
1012 // --- 1. Allocate the array of MGCtx structs ---
1013 ierr = PetscMalloc(usermg->mglevels * sizeof(MGCtx), &usermg->mgctx); CHKERRQ(ierr);
1014 // Zero-initialize to ensure all pointers (especially packer) are NULL
1015 ierr = PetscMemzero(usermg->mgctx, usermg->mglevels * sizeof(MGCtx)); CHKERRQ(ierr);
1016 mgctx = usermg->mgctx;
1017 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Allocated MGCtx array of size %d.\n", simCtx->rank, usermg->mglevels);
1018
1019 // --- 2. Parse semi-coarsening options (logic from MG_Initial) ---
1020 // These flags determine if a dimension is coarsened in the multigrid hierarchy.
1021 PetscInt *isc, *jsc, *ksc;
1022 ierr = PetscMalloc3(nblk, &isc, nblk, &jsc, nblk, &ksc); CHKERRQ(ierr);
1023 // Set defaults to FALSE (full coarsening)
1024 for (PetscInt i = 0; i < nblk; ++i) {
1025 isc[i] = 0; jsc[i] = 0; ksc[i] = 0;
1026 }
1027
1028// Use a temporary variable for the 'count' argument to the parsing function.
1029 // This protects the original 'nblk' which is needed for the loop bounds.
1030 PetscInt n_opts_found = nblk;
1031 ierr = PetscOptionsGetIntArray(NULL, NULL, "-mg_i_semi", isc, &n_opts_found, &found); CHKERRQ(ierr);
1032
1033 n_opts_found = nblk; // Reset the temp variable before the next call
1034 ierr = PetscOptionsGetIntArray(NULL, NULL, "-mg_j_semi", jsc, &n_opts_found, &found); CHKERRQ(ierr);
1035
1036 n_opts_found = nblk; // Reset the temp variable before the next call
1037 ierr = PetscOptionsGetIntArray(NULL, NULL, "-mg_k_semi", ksc, &n_opts_found, &found); CHKERRQ(ierr);
1038
1039 // --- 3. Loop over levels and blocks to allocate UserCtx arrays ---
1040 for (PetscInt level = 0; level < simCtx->mglevels; level++) {
1041
1042 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Setting up MG Level %d...\n", simCtx->rank, level);
1043 // Allocate the array of UserCtx structs for this level
1044 ierr = PetscMalloc(nblk * sizeof(UserCtx), &mgctx[level].user); CHKERRQ(ierr);
1045 // It's good practice to zero out the memory to avoid uninitialized values
1046 ierr = PetscMemzero(mgctx[level].user, nblk * sizeof(UserCtx)); CHKERRQ(ierr);
1047 mgctx[level].thislevel = level;
1048
1049 for (PetscInt bi = 0; bi < nblk; bi++) {
1050 UserCtx *currentUser = &mgctx[level].user[bi];
1051 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Initializing UserCtx for Level %d, Block %d.\n", simCtx->rank, level, bi);
1052
1053 // --- CRITICAL STEP: Set the back-pointer to the master context ---
1054 currentUser->simCtx = simCtx;
1055
1056 // Initialize other per-context values
1057 currentUser->thislevel = level;
1058 currentUser->_this = bi; //
1059 currentUser->mglevels = usermg->mglevels;
1060
1061 // Assign semi-coarsening flags
1062 currentUser->isc = isc[bi];
1063 currentUser->jsc = jsc[bi];
1064 currentUser->ksc = ksc[bi];
1065
1066 // Link to finer/coarser contexts for multigrid operations
1067 if (level > 0) {
1068 currentUser->user_c = &mgctx[level-1].user[bi];
1069 LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG, "Rank %d: -> Linked to coarser context (user_c).\n", simCtx->rank);
1070 }
1071 if (level < usermg->mglevels - 1) {
1072 currentUser->user_f = &mgctx[level+1].user[bi];
1073 LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG, "Rank %d: -> Linked to finer context (user_f).\n", simCtx->rank);
1074 }
1075 }
1076 }
1077
1078 // Log a summary of the parsed flags on each rank.
1079 if (get_log_level() >= LOG_DEBUG && nblk > 0) {
1080 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Final semi-coarsening configuration view:\n", simCtx->rank);
1081 for (PetscInt bi = 0; bi < nblk; ++bi) {
1082 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]);
1083 }
1084 }
1085
1086 // Clean up temporary arrays
1087 ierr = PetscFree3(isc, jsc, ksc); CHKERRQ(ierr);
1088
1089 LOG_ALLOW(GLOBAL, LOG_INFO, "Context hierarchy allocation complete.\n");
1091 PetscFunctionReturn(0);
1092}
LogLevel get_log_level()
Retrieves the current logging level from the environment variable LOG_LEVEL.
Definition logging.c:84
PetscInt isc
Definition variables.h:824
UserCtx * user
Definition variables.h:528
PetscInt mglevels
Definition variables.h:874
UserCtx * user_f
Definition variables.h:875
SimCtx * simCtx
Back-pointer to the master simulation context.
Definition variables.h:814
PetscInt ksc
Definition variables.h:824
UserMG usermg
Definition variables.h:764
PetscInt _this
Definition variables.h:824
PetscInt jsc
Definition variables.h:824
PetscInt thislevel
Definition variables.h:529
UserCtx * user_c
Definition variables.h:875
PetscInt thislevel
Definition variables.h:874
PetscInt mglevels
Definition variables.h:535
MGCtx * mgctx
Definition variables.h:538
Context for Multigrid operations.
Definition variables.h:527
User-defined context containing data specific to a single computational grid level.
Definition variables.h:811
User-level context for managing the entire multigrid hierarchy.
Definition variables.h:534
Here is the call graph for this function:
Here is the caller graph for this function:

◆ SetupSolverParameters()

static PetscErrorCode SetupSolverParameters ( SimCtx simCtx)
static

Internal helper implementation: SetupSolverParameters().

Local to this translation unit.

Definition at line 1100 of file setup.c.

1100 {
1101
1102 PetscFunctionBeginUser;
1104
1105 LOG_ALLOW(GLOBAL,LOG_INFO, " -- Setting up solver parameters -- .\n");
1106
1107 UserMG *usermg = &simCtx->usermg;
1108 MGCtx *mgctx = usermg->mgctx;
1109 PetscInt nblk = simCtx->block_number;
1110
1111 for (PetscInt level = usermg->mglevels-1; level >=0; level--) {
1112 for (PetscInt bi = 0; bi < nblk; bi++) {
1113 UserCtx *user = &mgctx[level].user[bi];
1114 LOG_ALLOW_SYNC(LOCAL, LOG_DEBUG, "Rank %d: Setting up parameters for level %d, block %d\n", simCtx->rank, level, bi);
1115
1116 user->assignedA = PETSC_FALSE;
1117 user->multinullspace = PETSC_FALSE;
1118 }
1119 }
1121 PetscFunctionReturn(0);
1122}
PetscBool assignedA
Definition variables.h:854
PetscBool multinullspace
Definition variables.h:851
Here is the caller graph for this function:

◆ SetupGridAndSolvers()

PetscErrorCode SetupGridAndSolvers ( SimCtx simCtx)

Implementation of SetupGridAndSolvers().

The main orchestrator for setting up all grid-related components.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/setup.h.

See also
SetupGridAndSolvers()

Definition at line 1132 of file setup.c.

1133{
1134 PetscErrorCode ierr;
1135 PetscFunctionBeginUser;
1136
1138
1139 LOG_ALLOW(GLOBAL, LOG_INFO, "--- Starting Grid and Solvers Setup ---\n");
1140
1141 // Phase 1: Allocate the UserMG and UserCtx hierarchy
1142 ierr = AllocateContextHierarchy(simCtx); CHKERRQ(ierr);
1143
1144 ierr = DefineAllGridDimensions(simCtx); CHKERRQ(ierr);
1145 ierr = InitializeAllGridDMs(simCtx); CHKERRQ(ierr);
1146 ierr = AssignAllGridCoordinates(simCtx);
1147 ierr = CreateAndInitializeAllVectors(simCtx); CHKERRQ(ierr);
1148 ierr = SetupSolverParameters(simCtx); CHKERRQ(ierr);
1149
1150 // NOTE: CalculateAllGridMetrics is now called inside SetupBoundaryConditions (not here) to ensure:
1151 // 1. Boundary condition configuration data (boundary_faces) is available for periodic BC corrections
1152 // 2. Computed metrics are available for inlet/outlet area calculations
1153 // This resolves the circular dependency between BC setup and metric calculations.
1154
1155 LOG_ALLOW(GLOBAL, LOG_INFO, "--- Grid and Solvers Setup Complete ---\n");
1156
1158 PetscFunctionReturn(0);
1159}
PetscErrorCode DefineAllGridDimensions(SimCtx *simCtx)
Orchestrates the parsing and setting of grid dimensions for all blocks.
Definition grid.c:57
PetscErrorCode InitializeAllGridDMs(SimCtx *simCtx)
Orchestrates the creation of DMDA objects for every block and multigrid level.
Definition grid.c:235
PetscErrorCode AssignAllGridCoordinates(SimCtx *simCtx)
Orchestrates the assignment of physical coordinates to all DMDA objects.
Definition grid.c:317
static PetscErrorCode SetupSolverParameters(SimCtx *simCtx)
Internal helper implementation: SetupSolverParameters().
Definition setup.c:1100
PetscErrorCode CreateAndInitializeAllVectors(SimCtx *simCtx)
Internal helper implementation: CreateAndInitializeAllVectors().
Definition setup.c:1168
static PetscErrorCode AllocateContextHierarchy(SimCtx *simCtx)
Internal helper implementation: AllocateContextHierarchy().
Definition setup.c:997
Here is the call graph for this function:
Here is the caller graph for this function:

◆ CreateAndInitializeAllVectors()

PetscErrorCode CreateAndInitializeAllVectors ( SimCtx simCtx)

Internal helper implementation: CreateAndInitializeAllVectors().

Creates and initializes all PETSc Vec objects for all fields.

Local to this translation unit.

Definition at line 1168 of file setup.c.

1169{
1170 PetscErrorCode ierr;
1171 UserMG *usermg = &simCtx->usermg;
1172 MGCtx *mgctx = usermg->mgctx;
1173 PetscInt nblk = simCtx->block_number;
1174
1175 PetscFunctionBeginUser;
1176
1178
1179 LOG_ALLOW(GLOBAL, LOG_INFO, "Creating and initializing all simulation vectors...\n");
1180
1181 for (PetscInt level = usermg->mglevels-1; level >=0; level--) {
1182 for (PetscInt bi = 0; bi < nblk; bi++) {
1183 UserCtx *user = &mgctx[level].user[bi];
1184
1185 if(!user->da || !user->fda) {
1186 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "DMs not properly initialized in UserCtx before vector creation.");
1187 }
1188
1189 LOG_ALLOW_SYNC(LOCAL, LOG_DEBUG, "Rank %d: Creating vectors for level %d, block %d\n", simCtx->rank, level, bi);
1190
1191 // --- Group A: Primary Flow Fields (Global and Local) ---
1192 // These are the core solution variables.
1193 ierr = DMCreateGlobalVector(user->fda, &user->Ucont); CHKERRQ(ierr); ierr = VecSet(user->Ucont, 0.0); CHKERRQ(ierr);
1194 ierr = DMCreateGlobalVector(user->fda, &user->Ucat); CHKERRQ(ierr); ierr = VecSet(user->Ucat, 0.0); CHKERRQ(ierr);
1195 ierr = DMCreateGlobalVector(user->da, &user->P); CHKERRQ(ierr); ierr = VecSet(user->P, 0.0); CHKERRQ(ierr);
1196 ierr = DMCreateGlobalVector(user->da, &user->Nvert); CHKERRQ(ierr); ierr = VecSet(user->Nvert, 0.0); CHKERRQ(ierr);
1197
1198 ierr = DMCreateLocalVector(user->fda, &user->lUcont); CHKERRQ(ierr); ierr = VecSet(user->lUcont, 0.0); CHKERRQ(ierr);
1199 ierr = DMCreateLocalVector(user->fda, &user->lUcat); CHKERRQ(ierr); ierr = VecSet(user->lUcat, 0.0); CHKERRQ(ierr);
1200 ierr = DMCreateLocalVector(user->da, &user->lP); CHKERRQ(ierr); ierr = VecSet(user->lP, 0.0); CHKERRQ(ierr);
1201 ierr = DMCreateLocalVector(user->da, &user->lNvert); CHKERRQ(ierr); ierr = VecSet(user->lNvert, 0.0); CHKERRQ(ierr);
1202
1203 // -- Group A2: Derived Flow Fields (Global and Local) ---
1204 ierr = VecDuplicate(user->P,&user->Diffusivity); CHKERRQ(ierr); ierr = VecSet(user->Diffusivity, 0.0); CHKERRQ(ierr);
1205 ierr = VecDuplicate(user->lP,&user->lDiffusivity); CHKERRQ(ierr); ierr = VecSet(user->lDiffusivity, 0.0); CHKERRQ(ierr);
1206 ierr = VecDuplicate(user->Ucat,&user->DiffusivityGradient); CHKERRQ(ierr); ierr = VecSet(user->DiffusivityGradient, 0.0); CHKERRQ(ierr);
1207 ierr = VecDuplicate(user->lUcat,&user->lDiffusivityGradient); CHKERRQ(ierr); ierr = VecSet(user->lDiffusivityGradient, 0.0); CHKERRQ(ierr);
1208
1209 // -- Group B: Solver Work Vectors (Global and Local) ---
1210 ierr = VecDuplicate(user->P, &user->Phi); CHKERRQ(ierr); ierr = VecSet(user->Phi, 0.0); CHKERRQ(ierr);
1211 ierr = VecDuplicate(user->lP, &user->lPhi); CHKERRQ(ierr); ierr = VecSet(user->lPhi, 0.0); CHKERRQ(ierr);
1212
1213 // --- Group C: Time-Stepping & Workspace Fields (Finest Level Only) ---
1214 if (level == usermg->mglevels - 1) {
1215 ierr = VecDuplicate(user->Ucont, &user->Ucont_o); CHKERRQ(ierr); ierr = VecSet(user->Ucont_o, 0.0); CHKERRQ(ierr);
1216 ierr = VecDuplicate(user->Ucont, &user->Ucont_rm1); CHKERRQ(ierr); ierr = VecSet(user->Ucont_rm1, 0.0); CHKERRQ(ierr);
1217 ierr = VecDuplicate(user->Ucat, &user->Ucat_o); CHKERRQ(ierr); ierr = VecSet(user->Ucat_o, 0.0); CHKERRQ(ierr);
1218 ierr = VecDuplicate(user->P, &user->P_o); CHKERRQ(ierr); ierr = VecSet(user->P_o, 0.0); CHKERRQ(ierr);
1219 ierr = VecDuplicate(user->lUcont, &user->lUcont_o); CHKERRQ(ierr); ierr = VecSet(user->lUcont_o, 0.0); CHKERRQ(ierr);
1220 ierr = VecDuplicate(user->lUcont, &user->lUcont_rm1); CHKERRQ(ierr); ierr = VecSet(user->lUcont_rm1, 0.0); CHKERRQ(ierr);
1221 ierr = DMCreateLocalVector(user->da, &user->lNvert_o); CHKERRQ(ierr); ierr = VecSet(user->lNvert_o, 0.0); CHKERRQ(ierr);
1222 ierr = VecDuplicate(user->Nvert, &user->Nvert_o); CHKERRQ(ierr); ierr = VecSet(user->Nvert_o, 0.0); CHKERRQ(ierr);
1223 }
1224
1225 // --- Group D: Grid Metrics (Face-Centered) ---
1226 ierr = DMCreateGlobalVector(user->fda, &user->Csi); CHKERRQ(ierr); ierr = VecSet(user->Csi, 0.0); CHKERRQ(ierr);
1227 ierr = VecDuplicate(user->Csi, &user->Eta); CHKERRQ(ierr); ierr = VecSet(user->Eta, 0.0); CHKERRQ(ierr);
1228 ierr = VecDuplicate(user->Csi, &user->Zet); CHKERRQ(ierr); ierr = VecSet(user->Zet, 0.0); CHKERRQ(ierr);
1229 ierr = DMCreateGlobalVector(user->da, &user->Aj); CHKERRQ(ierr); ierr = VecSet(user->Aj, 0.0); CHKERRQ(ierr);
1230
1231 ierr = DMCreateLocalVector(user->fda, &user->lCsi); CHKERRQ(ierr); ierr = VecSet(user->lCsi, 0.0); CHKERRQ(ierr);
1232 ierr = VecDuplicate(user->lCsi, &user->lEta); CHKERRQ(ierr); ierr = VecSet(user->lEta, 0.0); CHKERRQ(ierr);
1233 ierr = VecDuplicate(user->lCsi, &user->lZet); CHKERRQ(ierr); ierr = VecSet(user->lZet, 0.0); CHKERRQ(ierr);
1234 ierr = DMCreateLocalVector(user->da, &user->lAj); CHKERRQ(ierr); ierr = VecSet(user->lAj, 0.0); CHKERRQ(ierr);
1235
1236
1237 // --- Group E: Grid Metrics (Face-Centered) ---
1238 // Vector metrics are duplicated from Csi (DOF=3, fda-based)
1239 ierr = VecDuplicate(user->Csi, &user->ICsi); CHKERRQ(ierr); ierr = VecSet(user->ICsi, 0.0); CHKERRQ(ierr);
1240 ierr = VecDuplicate(user->Csi, &user->IEta); CHKERRQ(ierr); ierr = VecSet(user->IEta, 0.0); CHKERRQ(ierr);
1241 ierr = VecDuplicate(user->Csi, &user->IZet); CHKERRQ(ierr); ierr = VecSet(user->IZet, 0.0); CHKERRQ(ierr);
1242 ierr = VecDuplicate(user->Csi, &user->JCsi); CHKERRQ(ierr); ierr = VecSet(user->JCsi, 0.0); CHKERRQ(ierr);
1243 ierr = VecDuplicate(user->Csi, &user->JEta); CHKERRQ(ierr); ierr = VecSet(user->JEta, 0.0); CHKERRQ(ierr);
1244 ierr = VecDuplicate(user->Csi, &user->JZet); CHKERRQ(ierr); ierr = VecSet(user->JZet, 0.0); CHKERRQ(ierr);
1245 ierr = VecDuplicate(user->Csi, &user->KCsi); CHKERRQ(ierr); ierr = VecSet(user->KCsi, 0.0); CHKERRQ(ierr);
1246 ierr = VecDuplicate(user->Csi, &user->KEta); CHKERRQ(ierr); ierr = VecSet(user->KEta, 0.0); CHKERRQ(ierr);
1247 ierr = VecDuplicate(user->Csi, &user->KZet); CHKERRQ(ierr); ierr = VecSet(user->KZet, 0.0); CHKERRQ(ierr);
1248 // Scalar metrics are duplicated from Aj (DOF=1, da-based)
1249 ierr = VecDuplicate(user->Aj, &user->IAj); CHKERRQ(ierr); ierr = VecSet(user->IAj, 0.0); CHKERRQ(ierr);
1250 ierr = VecDuplicate(user->Aj, &user->JAj); CHKERRQ(ierr); ierr = VecSet(user->JAj, 0.0); CHKERRQ(ierr);
1251 ierr = VecDuplicate(user->Aj, &user->KAj); CHKERRQ(ierr); ierr = VecSet(user->KAj, 0.0); CHKERRQ(ierr);
1252
1253 ierr = VecDuplicate(user->lCsi, &user->lICsi); CHKERRQ(ierr); ierr = VecSet(user->lICsi, 0.0); CHKERRQ(ierr);
1254 ierr = VecDuplicate(user->lCsi, &user->lIEta); CHKERRQ(ierr); ierr = VecSet(user->lIEta, 0.0); CHKERRQ(ierr);
1255 ierr = VecDuplicate(user->lCsi, &user->lIZet); CHKERRQ(ierr); ierr = VecSet(user->lIZet, 0.0); CHKERRQ(ierr);
1256 ierr = VecDuplicate(user->lCsi, &user->lJCsi); CHKERRQ(ierr); ierr = VecSet(user->lJCsi, 0.0); CHKERRQ(ierr);
1257 ierr = VecDuplicate(user->lCsi, &user->lJEta); CHKERRQ(ierr); ierr = VecSet(user->lJEta, 0.0); CHKERRQ(ierr);
1258 ierr = VecDuplicate(user->lCsi, &user->lJZet); CHKERRQ(ierr); ierr = VecSet(user->lJZet, 0.0); CHKERRQ(ierr);
1259 ierr = VecDuplicate(user->lCsi, &user->lKCsi); CHKERRQ(ierr); ierr = VecSet(user->lKCsi, 0.0); CHKERRQ(ierr);
1260 ierr = VecDuplicate(user->lCsi, &user->lKEta); CHKERRQ(ierr); ierr = VecSet(user->lKEta, 0.0); CHKERRQ(ierr);
1261 ierr = VecDuplicate(user->lCsi, &user->lKZet); CHKERRQ(ierr); ierr = VecSet(user->lKZet, 0.0); CHKERRQ(ierr);
1262
1263 ierr = VecDuplicate(user->lAj, &user->lIAj); CHKERRQ(ierr); ierr = VecSet(user->lIAj, 0.0); CHKERRQ(ierr);
1264 ierr = VecDuplicate(user->lAj, &user->lJAj); CHKERRQ(ierr); ierr = VecSet(user->lJAj, 0.0); CHKERRQ(ierr);
1265 ierr = VecDuplicate(user->lAj, &user->lKAj); CHKERRQ(ierr); ierr = VecSet(user->lKAj, 0.0); CHKERRQ(ierr);
1266
1267 // --- Group F: Cell/Face Center Coordinates and Grid Spacing ---
1268 ierr = DMCreateGlobalVector(user->fda, &user->Cent); CHKERRQ(ierr); ierr = VecSet(user->Cent, 0.0); CHKERRQ(ierr);
1269 ierr = DMCreateLocalVector(user->fda, &user->lCent); CHKERRQ(ierr); ierr = VecSet(user->lCent, 0.0); CHKERRQ(ierr);
1270
1271 ierr = VecDuplicate(user->Cent, &user->GridSpace); CHKERRQ(ierr); ierr = VecSet(user->GridSpace, 0.0); CHKERRQ(ierr);
1272 ierr = VecDuplicate(user->lCent, &user->lGridSpace); CHKERRQ(ierr); ierr = VecSet(user->lGridSpace, 0.0); CHKERRQ(ierr);
1273
1274 // Face-center coordinate vectors are LOCAL to hold calculated values before scattering
1275 ierr = VecDuplicate(user->lCent, &user->Centx); CHKERRQ(ierr); ierr = VecSet(user->Centx, 0.0); CHKERRQ(ierr);
1276 ierr = VecDuplicate(user->lCent, &user->Centy); CHKERRQ(ierr); ierr = VecSet(user->Centy, 0.0); CHKERRQ(ierr);
1277 ierr = VecDuplicate(user->lCent, &user->Centz); CHKERRQ(ierr); ierr = VecSet(user->Centz, 0.0); CHKERRQ(ierr);
1278
1279 if(level == usermg->mglevels -1){
1280 // --- Group G: Turbulence Models (Finest Level Only) ---
1281 if (simCtx->les || simCtx->rans) {
1282 ierr = DMCreateGlobalVector(user->da, &user->Nu_t); CHKERRQ(ierr); ierr = VecSet(user->Nu_t, 0.0); CHKERRQ(ierr);
1283 ierr = DMCreateLocalVector(user->da, &user->lNu_t); CHKERRQ(ierr); ierr = VecSet(user->lNu_t, 0.0); CHKERRQ(ierr);
1284 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Turbulence viscosity (Nu_t) vectors created for LES/RANS model.\n");
1285 if(simCtx->les){
1286 ierr = DMCreateGlobalVector(user->da,&user->CS); CHKERRQ(ierr); ierr = VecSet(user->CS,0.0); CHKERRQ(ierr);
1287 ierr = DMCreateLocalVector(user->da,&user->lCs); CHKERRQ(ierr); ierr = VecSet(user->lCs,0.0); CHKERRQ(ierr);
1288 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Smagorinsky constant (CS) vectors created for LES model.\n");
1289 }
1290
1291 if(simCtx->wallfunction){
1292 ierr = DMCreateLocalVector(user->fda,&user->lFriction_Velocity); CHKERRQ(ierr); ierr = VecSet(user->lFriction_Velocity,0.0);
1293 }
1294 // Add K_Omega etc. here as needed
1295
1296 // Note: Add any other vectors from the legacy MG_Initial here as needed.
1297 // For example: Rhs, Forcing, turbulence Vecs (K_Omega, Nu_t)...
1298
1299 }
1300 // --- Group H: Particle Methods
1301 if(simCtx->np>0){
1302 ierr = DMCreateGlobalVector(user->da,&user->ParticleCount); CHKERRQ(ierr); ierr = VecSet(user->ParticleCount,0.0); CHKERRQ(ierr);
1303 ierr = DMCreateLocalVector(user->da,&user->lParticleCount); CHKERRQ(ierr); ierr = VecSet(user->lParticleCount,0.0); CHKERRQ(ierr);
1304 // Scalar field to hold particle scalar property (e.g., temperature, concentration)
1305 ierr = DMCreateGlobalVector(user->da,&user->Psi); CHKERRQ(ierr); ierr = VecSet(user->Psi,0.0); CHKERRQ(ierr);
1306 ierr = DMCreateLocalVector(user->da,&user->lPsi); CHKERRQ(ierr); ierr = VecSet(user->lPsi,0.0); CHKERRQ(ierr);
1307 LOG_ALLOW(GLOBAL,LOG_DEBUG,"ParticleCount & Scalar(Psi) created for %d particles.\n",simCtx->np);
1308 }
1309 }
1310 // --- Group I: Boundary Condition vectors ---
1311 ierr = DMCreateGlobalVector(user->fda, &user->Bcs.Ubcs); CHKERRQ(ierr);
1312 ierr = VecSet(user->Bcs.Ubcs, 0.0); CHKERRQ(ierr);
1313 ierr = DMCreateGlobalVector(user->fda, &user->Bcs.Uch); CHKERRQ(ierr);
1314 ierr = VecSet(user->Bcs.Uch, 0.0); CHKERRQ(ierr);
1315
1316 if(level == usermg->mglevels - 1){
1317 if(simCtx->exec_mode == EXEC_MODE_POSTPROCESSOR){
1318 LOG_ALLOW(LOCAL, LOG_DEBUG, "Post-processor mode detected. Allocating derived field vectors.\n");
1319
1320 ierr = VecDuplicate(user->P, &user->P_nodal); CHKERRQ(ierr);
1321 ierr = VecSet(user->P_nodal, 0.0); CHKERRQ(ierr);
1322
1323 ierr = VecDuplicate(user->Ucat, &user->Ucat_nodal); CHKERRQ(ierr);
1324 ierr = VecSet(user->Ucat_nodal, 0.0); CHKERRQ(ierr);
1325
1326 ierr = VecDuplicate(user->P, &user->Qcrit); CHKERRQ(ierr);
1327 ierr = VecSet(user->Qcrit, 0.0); CHKERRQ(ierr);
1328
1329 LOG_ALLOW(LOCAL, LOG_DEBUG, "Derived field vectors P_nodal, Ucat_nodal, and Qcrit created.\n");
1330
1331 if(simCtx->np>0){
1332 ierr = VecDuplicate(user->Psi, &user->Psi_nodal); CHKERRQ(ierr);
1333 ierr = VecSet(user->Psi_nodal, 0.0); CHKERRQ(ierr);
1334
1335 LOG_ALLOW(LOCAL, LOG_DEBUG, "Derived field vector Psi_nodal created for particle scalar property.\n");
1336
1337 }
1338 }else{
1339 user->P_nodal = NULL;
1340 user->Ucat_nodal = NULL;
1341 user->Qcrit = NULL;
1342 user->Psi_nodal = NULL;
1343 }
1344 }
1345
1346 }
1347}
1348
1349 LOG_ALLOW(GLOBAL, LOG_INFO, "All simulation vectors created and initialized.\n");
1350
1352 PetscFunctionReturn(0);
1353}
Vec lFriction_Velocity
Definition variables.h:833
Vec lDiffusivityGradient
Definition variables.h:841
Vec lCent
Definition variables.h:858
Vec GridSpace
Definition variables.h:858
Vec P_nodal
Definition variables.h:887
Vec JCsi
Definition variables.h:861
Vec KAj
Definition variables.h:862
Vec JEta
Definition variables.h:861
Vec Zet
Definition variables.h:858
Vec lIEta
Definition variables.h:860
Vec lIZet
Definition variables.h:860
Vec lNvert
Definition variables.h:837
Vec Phi
Definition variables.h:837
Vec IZet
Definition variables.h:860
Vec Centz
Definition variables.h:859
Vec IEta
Definition variables.h:860
Vec lZet
Definition variables.h:858
Vec Csi
Definition variables.h:858
Vec lUcont_rm1
Definition variables.h:845
Vec lIAj
Definition variables.h:860
Vec lKEta
Definition variables.h:862
Vec Ucat_nodal
Definition variables.h:888
Vec lPsi
Definition variables.h:883
Vec DiffusivityGradient
Definition variables.h:841
Vec lJCsi
Definition variables.h:861
Vec lCs
Definition variables.h:865
Vec Ucont
Definition variables.h:837
Vec Ubcs
Physical Cartesian velocity at boundary faces. Full 3D array but only boundary-face entries are meani...
Definition variables.h:121
Vec Qcrit
Definition variables.h:889
Vec JZet
Definition variables.h:861
Vec Centx
Definition variables.h:859
BCS Bcs
Definition variables.h:832
Vec lPhi
Definition variables.h:837
Vec lParticleCount
Definition variables.h:882
Vec lUcont_o
Definition variables.h:844
Vec Ucat_o
Definition variables.h:844
Vec lKZet
Definition variables.h:862
Vec Eta
Definition variables.h:858
Vec lNu_t
Definition variables.h:865
Vec Nu_t
Definition variables.h:865
Vec lJEta
Definition variables.h:861
Vec lCsi
Definition variables.h:858
Vec lGridSpace
Definition variables.h:858
Vec ICsi
Definition variables.h:860
Vec lKCsi
Definition variables.h:862
Vec Ucat
Definition variables.h:837
Vec ParticleCount
Definition variables.h:882
Vec Ucont_o
Definition variables.h:844
Vec lJZet
Definition variables.h:861
Vec Nvert_o
Definition variables.h:844
Vec IAj
Definition variables.h:860
Vec Psi_nodal
Definition variables.h:890
Vec JAj
Definition variables.h:861
Vec KEta
Definition variables.h:862
Vec Ucont_rm1
Definition variables.h:845
Vec lUcont
Definition variables.h:837
Vec Diffusivity
Definition variables.h:840
Vec lAj
Definition variables.h:858
Vec lICsi
Definition variables.h:860
Vec lUcat
Definition variables.h:837
Vec lEta
Definition variables.h:858
Vec KZet
Definition variables.h:862
Vec Cent
Definition variables.h:858
Vec Nvert
Definition variables.h:837
Vec KCsi
Definition variables.h:862
Vec lDiffusivity
Definition variables.h:840
Vec lNvert_o
Definition variables.h:844
Vec Centy
Definition variables.h:859
Vec lJAj
Definition variables.h:861
Vec lKAj
Definition variables.h:862
Vec Psi
Definition variables.h:883
Vec P_o
Definition variables.h:844
Vec Uch
Characteristic velocity for boundary conditions.
Definition variables.h:122
Here is the caller graph for this function:

◆ UpdateLocalGhosts()

PetscErrorCode UpdateLocalGhosts ( UserCtx user,
const char *  fieldName 
)

Internal helper implementation: UpdateLocalGhosts().

Updates the local vector (including ghost points) from its corresponding global vector.

Local to this translation unit.

Definition at line 1361 of file setup.c.

1362{
1363 PetscErrorCode ierr;
1364 PetscMPIInt rank;
1365 Vec globalVec = NULL;
1366 Vec localVec = NULL;
1367 DM dm = NULL; // The DM associated with this field pair
1368
1369 PetscFunctionBeginUser; // Use User version for application code
1371 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
1372 LOG_ALLOW(GLOBAL, LOG_INFO, "Rank %d: Starting ghost update for field '%s'.\n", rank, fieldName);
1373
1374 // --- 1. Identify the correct Vectors and DM ---
1375 if (strcmp(fieldName, "Ucat") == 0) {
1376 globalVec = user->Ucat;
1377 localVec = user->lUcat;
1378 dm = user->fda;
1379 } else if (strcmp(fieldName, "Ucont") == 0) {
1380 globalVec = user->Ucont;
1381 localVec = user->lUcont;
1382 dm = user->fda;
1383 } else if (strcmp(fieldName, "P") == 0) {
1384 globalVec = user->P;
1385 localVec = user->lP;
1386 dm = user->da;
1387 } else if (strcmp(fieldName, "Diffusivity") == 0) {
1388 globalVec = user->Diffusivity;
1389 localVec = user->lDiffusivity;
1390 dm = user->da;
1391 } else if (strcmp(fieldName, "DiffusivityGradient") == 0) {
1392 globalVec = user->DiffusivityGradient;
1393 localVec = user->lDiffusivityGradient;
1394 dm = user->fda;
1395 } else if (strcmp(fieldName, "Csi") == 0) {
1396 globalVec = user->Csi;
1397 localVec = user->lCsi;
1398 dm = user->fda;
1399 } else if (strcmp(fieldName, "Eta") == 0) {
1400 globalVec = user->Eta;
1401 localVec = user->lEta;
1402 dm = user->fda;
1403 } else if (strcmp(fieldName, "Zet") == 0) {
1404 globalVec = user->Zet;
1405 localVec = user->lZet;
1406 dm = user->fda;
1407 }else if (strcmp(fieldName, "Nvert") == 0) {
1408 globalVec = user->Nvert;
1409 localVec = user->lNvert;
1410 dm = user->da;
1411 // Add other fields as needed
1412 } else if (strcmp(fieldName, "Aj") == 0) {
1413 globalVec = user->Aj;
1414 localVec = user->lAj;
1415 dm = user->da;
1416 } else if (strcmp(fieldName, "Cent") == 0) {
1417 globalVec = user->Cent;
1418 localVec = user->lCent;
1419 dm = user->fda;
1420 }else if (strcmp(fieldName, "GridSpace") == 0) {
1421 globalVec = user->GridSpace;
1422 localVec = user->lGridSpace;
1423 dm = user->fda;
1424 }else if (strcmp(fieldName,"ICsi") == 0){
1425 globalVec = user->ICsi;
1426 localVec = user->lICsi;
1427 dm = user->fda;
1428 }else if (strcmp(fieldName,"IEta") == 0){
1429 globalVec = user->IEta;
1430 localVec = user->lIEta;
1431 dm = user->fda;
1432 }else if (strcmp(fieldName,"IZet") == 0){
1433 globalVec = user->IZet;
1434 localVec = user->lIZet;
1435 dm = user->fda;
1436 }else if (strcmp(fieldName,"JCsi") == 0){
1437 globalVec = user->JCsi;
1438 localVec = user->lJCsi;
1439 dm = user->fda;
1440 }else if (strcmp(fieldName,"JEta") == 0){
1441 globalVec = user->JEta;
1442 localVec = user->lJEta;
1443 dm = user->fda;
1444 }else if (strcmp(fieldName,"JZet") == 0){
1445 globalVec = user->JZet;
1446 localVec = user->lJZet;
1447 dm = user->fda;
1448 }else if (strcmp(fieldName,"KCsi") == 0){
1449 globalVec = user->KCsi;
1450 localVec = user->lKCsi;
1451 dm = user->fda;
1452 }else if (strcmp(fieldName,"KEta") == 0){
1453 globalVec = user->KEta;
1454 localVec = user->lKEta;
1455 dm = user->fda;
1456 }else if (strcmp(fieldName,"KZet") == 0){
1457 globalVec = user->KZet;
1458 localVec = user->lKZet;
1459 dm = user->fda;
1460 }else if (strcmp(fieldName,"IAj") == 0){
1461 globalVec = user->IAj;
1462 localVec = user->lIAj;
1463 dm = user->da;
1464 }else if (strcmp(fieldName,"JAj") == 0){
1465 globalVec = user->JAj;
1466 localVec = user->lJAj;
1467 dm = user->da;
1468 }else if (strcmp(fieldName,"KAj") == 0){
1469 globalVec = user->KAj;
1470 localVec = user->lKAj;
1471 dm = user->da;
1472 }else if (strcmp(fieldName,"Phi") == 0){ // Pressure correction term.
1473 globalVec = user->Phi;
1474 localVec = user->lPhi;
1475 dm = user->da;
1476 }else if (strcmp(fieldName,"Psi") == 0){ // Particle scalar property.
1477 globalVec = user->Psi;
1478 localVec = user->lPsi;
1479 dm = user->da;
1480 }else {
1481 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Field '%s' not recognized for ghost update.", fieldName);
1482 }
1483
1484 // --- 2. Check if components were found ---
1485 if (!globalVec) {
1486 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Global vector for field '%s' is NULL.", fieldName);
1487 }
1488 if (!localVec) {
1489 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Local vector for field '%s' is NULL.", fieldName);
1490 }
1491 if (!dm) {
1492 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "DM for field '%s' is NULL.", fieldName);
1493 }
1494
1495 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Identified components for '%s': DM=%p, GlobalVec=%p, LocalVec=%p.\n",
1496 rank, fieldName, (void*)dm, (void*)globalVec, (void*)localVec);
1497
1498 // --- 3. Optional Debugging: Norm Before Update ---
1499 // Use your logging convention check
1500 // if (get_log_level() >= LOG_LEVEL_DEBUG && is_function_allowed("UpdateLocalGhosts")) { // Example check
1501 if(get_log_level() == LOG_DEBUG && is_function_allowed(__func__)){
1502 PetscReal norm_global_before;
1503 ierr = VecNorm(globalVec, NORM_INFINITY, &norm_global_before); CHKERRQ(ierr);
1504 LOG_ALLOW(GLOBAL, LOG_INFO,"Max norm '%s' (Global) BEFORE Ghost Update: %g\n", fieldName, norm_global_before);
1505 // Optional: Norm of local vector before update (might contain old ghost values)
1506 // PetscReal norm_local_before;
1507 // ierr = VecNorm(localVec, NORM_INFINITY, &norm_local_before); CHKERRQ(ierr);
1508 // LOG_ALLOW(GLOBAL, LOG_DEBUG,"Max norm '%s' (Local) BEFORE Ghost Update: %g\n", fieldName, norm_local_before);
1509 }
1510
1511 // --- 4. Perform the Global-to-Local Transfer (Ghost Update) ---
1512 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Calling DMGlobalToLocalBegin/End for '%s'.\n", rank, fieldName);
1513 ierr = DMGlobalToLocalBegin(dm, globalVec, INSERT_VALUES, localVec); CHKERRQ(ierr);
1514 ierr = DMGlobalToLocalEnd(dm, globalVec, INSERT_VALUES, localVec); CHKERRQ(ierr);
1515 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Completed DMGlobalToLocalBegin/End for '%s'.\n", rank, fieldName);
1516
1517 // --- 5. Optional Debugging: Norm After Update ---
1518 // Use your logging convention check
1519 // if (get_log_level() >= LOG_LEVEL_DEBUG && is_function_allowed("UpdateLocalGhosts")) { // Example check
1520 if(get_log_level() == LOG_DEBUG && is_function_allowed(__func__)){ // Using your specific check
1521 PetscReal norm_local_after;
1522 ierr = VecNorm(localVec, NORM_INFINITY, &norm_local_after); CHKERRQ(ierr);
1523 LOG_ALLOW(GLOBAL, LOG_INFO,"Max norm '%s' (Local) AFTER Ghost Update: %g\n", fieldName, norm_local_after);
1524
1525 // --- 6. Optional Debugging: Specific Point Checks (Example for Ucat on Rank 0/1) ---
1526 // (Keep this conditional if it's only for specific debug scenarios)
1527 if (strcmp(fieldName, "Ucat") == 0) { // Only do detailed checks for Ucat for now
1528 PetscMPIInt rank_test;
1529 MPI_Comm_rank(PETSC_COMM_WORLD, &rank_test);
1530
1531 // Get Local Info needed for indexing checks
1532 DMDALocalInfo info_check;
1533 ierr = DMDAGetLocalInfo(dm, &info_check); CHKERRQ(ierr); // Use the correct dm
1534
1535 // Buffer for array pointer
1536 Cmpnts ***lUcat_arr_test = NULL;
1537 PetscErrorCode ierr_test = 0;
1538
1539 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Testing '%s' access immediately after ghost update...\n", rank_test, fieldName);
1540 ierr_test = DMDAVecGetArrayDOFRead(dm, localVec, &lUcat_arr_test); // Use correct dm and localVec
1541
1542 if (ierr_test) {
1543 LOG_ALLOW(LOCAL, LOG_ERROR, "Rank %d: ERROR %d getting '%s' array after ghost update!\n", rank_test, ierr_test, fieldName);
1544 } else if (!lUcat_arr_test) {
1545 LOG_ALLOW(LOCAL, LOG_ERROR, "Rank %d: ERROR NULL pointer getting '%s' array after ghost update!\n", rank_test, fieldName);
1546 }
1547 else {
1548 // Check owned interior point (e.g., first interior point)
1549 PetscInt k_int = info_check.zs + (info_check.zm > 1 ? 1 : 0); // Global k index (at least zs+1 if possible)
1550 PetscInt j_int = info_check.ys + (info_check.ym > 1 ? 1 : 0); // Global j index
1551 PetscInt i_int = info_check.xs + (info_check.xm > 1 ? 1 : 0); // Global i index
1552 // Ensure indices are within global bounds if domain is very small
1553 //if (k_int >= info_check.mz-1) k_int = info_check.mz-2; if (k_int < 1) k_int = 1;
1554 //if (j_int >= info_check.my-1) j_int = info_check.my-2; if (j_int < 1) j_int = 1;
1555 // if (i_int >= info_check.mx-1) i_int = info_check.mx-2; if (i_int < 1) i_int = 1;
1556 // clamp k_int to [1 .. mz-2]
1557 if (k_int >= info_check.mz - 1) {
1558 k_int = info_check.mz - 2;
1559 }
1560 if (k_int < 1) {
1561 k_int = 1;
1562 }
1563
1564 // clamp j_int to [1 .. my-2]
1565 if (j_int >= info_check.my - 1) {
1566 j_int = info_check.my - 2;
1567 }
1568 if (j_int < 1) {
1569 j_int = 1;
1570 }
1571
1572 // clamp i_int to [1 .. mx-2]
1573 if (i_int >= info_check.mx - 1) {
1574 i_int = info_check.mx - 2;
1575 }
1576 if (i_int < 1) {
1577 i_int = 1;
1578 }
1579
1580 // Only attempt read if indices are actually owned (relevant for multi-rank)
1581 if (k_int >= info_check.zs && k_int < info_check.zs + info_check.zm &&
1582 j_int >= info_check.ys && j_int < info_check.ys + info_check.ym &&
1583 i_int >= info_check.xs && i_int < info_check.xs + info_check.xm)
1584 {
1585 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);
1586 Cmpnts test_val_owned_interior = lUcat_arr_test[k_int][j_int][i_int];
1587 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: SUCCESS reading owned interior: x=%g\n", rank_test, test_val_owned_interior.x);
1588 } else {
1589 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);
1590 }
1591
1592
1593 // Check owned boundary point (e.g., first owned point)
1594 PetscInt k_bnd = info_check.zs; // Global k index
1595 PetscInt j_bnd = info_check.ys; // Global j index
1596 PetscInt i_bnd = info_check.xs; // Global i index
1597 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);
1598 Cmpnts test_val_owned_boundary = lUcat_arr_test[k_bnd][j_bnd][i_bnd];
1599 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: SUCCESS reading owned boundary: x=%g\n", rank_test, test_val_owned_boundary.x);
1600
1601
1602 // Check ghost point (e.g., one layer below in k, if applicable)
1603 if (info_check.zs > 0) { // Only if there's a rank below
1604 PetscInt k_ghost = info_check.zs - 1;
1605 PetscInt j_ghost = info_check.ys; // Use start of owned y, simple example
1606 PetscInt i_ghost = info_check.xs; // Use start of owned x, simple example
1607 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Attempting test read GHOST [%d][%d][%d] (Global)\n", rank_test, k_ghost, j_ghost, i_ghost);
1608 Cmpnts test_val_ghost = lUcat_arr_test[k_ghost][j_ghost][i_ghost];
1609 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: SUCCESS reading ghost: x=%g\n", rank_test, test_val_ghost.x);
1610 } else {
1611 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Skipping ghost test read (zs=0).\n", rank_test);
1612 }
1613
1614 // Restore the array
1615 ierr_test = DMDAVecRestoreArrayDOFRead(dm, localVec, &lUcat_arr_test);
1616 if(ierr_test){ LOG_ALLOW(LOCAL, LOG_ERROR, "Rank %d: ERROR %d restoring '%s' array after test read!\n", rank_test, ierr_test, fieldName); }
1617 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Finished testing '%s' access.\n", rank_test, fieldName);
1618 }
1619 } // end if Ucat
1620 } // end debug logging check
1621
1622 LOG_ALLOW(GLOBAL, LOG_INFO, "Rank %d: Completed ghost update for field '%s'.\n", rank, fieldName);
1624 PetscFunctionReturn(0);
1625}
PetscBool is_function_allowed(const char *functionName)
Checks if a given function is in the allow-list.
Definition logging.c:183
A 3D point or vector with PetscScalar components.
Definition variables.h:100
Here is the call graph for this function:
Here is the caller graph for this function:

◆ SetupBoundaryConditions()

PetscErrorCode SetupBoundaryConditions ( SimCtx simCtx)

Internal helper implementation: SetupBoundaryConditions().

(Orchestrator) Sets up all boundary conditions for the simulation.

Local to this translation unit.

Definition at line 1633 of file setup.c.

1634{
1635 PetscErrorCode ierr;
1636 PetscFunctionBeginUser;
1637
1639
1640 LOG_ALLOW(GLOBAL,LOG_INFO, "--- Setting up Boundary Conditions ---\n");
1641 // --- Phase 1: Parse and initialize BC configuration for all blocks ---
1642 LOG_ALLOW(GLOBAL,LOG_INFO,"Parsing BC configuration files and initializing boundary condition data structures.\n");
1643 UserCtx *user_finest = simCtx->usermg.mgctx[simCtx->usermg.mglevels-1].user;
1644 for (PetscInt bi = 0; bi < simCtx->block_number; bi++) {
1645 LOG_ALLOW(GLOBAL,LOG_DEBUG, " -> Processing Block %d:\n", bi);
1646
1647 // --- Generate the filename for the current block ---
1648 const char *current_bc_filename = simCtx->bcs_files[bi];
1649 LOG_ALLOW(GLOBAL,LOG_DEBUG," -> Processing Block %d using config file '%s'\n", bi, current_bc_filename);
1650 // This will populate user_finest[bi].boundary_faces
1651
1652 //ierr = ParseAllBoundaryConditions(&user_finest[bi],current_bc_filename); CHKERRQ(ierr);
1653
1654 ierr = BoundarySystem_Initialize(&user_finest[bi], current_bc_filename); CHKERRQ(ierr);
1655 }
1656
1657 // Propogate BC Configuration to coarser levels.
1658 ierr = PropagateBoundaryConfigToCoarserLevels(simCtx); CHKERRQ(ierr);
1659
1660 // --- Calculate Grid Metrics (requires BC configuration) ---
1661 // NOTE: This MUST be called here (after BC initialization but before inlet/outlet calculations) because:
1662 // 1. Periodic BC corrections in metric calculations need boundary_faces data to be populated
1663 // 2. Inlet/Outlet area calculations (below) require computed metrics (Csi, Eta, Zet) to be available
1664 // Previously this was in SetupGridAndSolvers, but that caused metrics to be computed without BC info.
1665 LOG_ALLOW(GLOBAL,LOG_INFO,"Computing grid metrics with boundary condition information.\n");
1666 ierr = CalculateAllGridMetrics(simCtx); CHKERRQ(ierr);
1667
1668 // --- Phase 2: Calculate inlet/outlet properties (requires computed metrics) ---
1669 LOG_ALLOW(GLOBAL,LOG_INFO,"Calculating inlet and outlet face properties.\n");
1670 for (PetscInt bi = 0; bi < simCtx->block_number; bi++) {
1671 // Call the function to calculate the center of the inlet face & the inlet area, which may be used to calculate Boundary values.
1672 ierr = CalculateInletProperties(&user_finest[bi]); CHKERRQ(ierr);
1673
1674 // Call the function to calculate the center of the outlet face & the outlet area, which may be used to calculate Boundary values.
1675 ierr = CalculateOutletProperties(&user_finest[bi]); CHKERRQ(ierr);
1676 }
1677
1678 LOG_ALLOW(GLOBAL,LOG_INFO, "--- Boundary Conditions setup complete ---\n");
1679
1680
1682 PetscFunctionReturn(0);
1683}
PetscErrorCode BoundarySystem_Initialize(UserCtx *user, const char *bcs_filename)
Initializes the entire boundary system.
Definition Boundaries.c:857
PetscErrorCode PropagateBoundaryConfigToCoarserLevels(SimCtx *simCtx)
Propagates boundary condition configuration from finest to all coarser multigrid levels.
Definition Boundaries.c:982
PetscErrorCode CalculateAllGridMetrics(SimCtx *simCtx)
Orchestrates the calculation of all grid metrics.
Definition Metric.c:2168
PetscErrorCode CalculateOutletProperties(UserCtx *user)
Calculates the center and area of the primary OUTLET face.
Definition grid.c:986
PetscErrorCode CalculateInletProperties(UserCtx *user)
Calculates the center and area of the primary INLET face.
Definition grid.c:933
Here is the call graph for this function:
Here is the caller graph for this function:

◆ Allocate3DArrayScalar()

PetscErrorCode Allocate3DArrayScalar ( PetscReal ****  array,
PetscInt  nz,
PetscInt  ny,
PetscInt  nx 
)

Internal helper implementation: Allocate3DArrayScalar().

Allocates a 3D array of PetscReal values using PetscCalloc.

Local to this translation unit.

Definition at line 1689 of file setup.c.

1690{
1691 PetscErrorCode ierr;
1692 PetscReal ***data;
1693 PetscReal *dataContiguous;
1694 PetscInt k, j;
1695
1696 PetscFunctionBegin;
1697 /* Step 1: Allocate memory for an array of nz layer pointers (zero-initialized) */
1698 ierr = PetscCalloc1(nz, &data); CHKERRQ(ierr);
1699
1700 /* Step 2: Allocate memory for all row pointers (nz * ny pointers) */
1701 ierr = PetscCalloc1(nz * ny, &data[0]); CHKERRQ(ierr);
1702 for (k = 1; k < nz; k++) {
1703 data[k] = data[0] + k * ny;
1704 }
1705
1706 /* Step 3: Allocate one contiguous block for all data elements (nz*ny*nx) */
1707 ierr = PetscCalloc1(nz * ny * nx, &dataContiguous); CHKERRQ(ierr);
1708
1709 /* Build the 3D pointer structure: each row pointer gets the correct segment of data */
1710 for (k = 0; k < nz; k++) {
1711 for (j = 0; j < ny; j++) {
1712 data[k][j] = dataContiguous + (k * ny + j) * nx;
1713 /* Memory is already zeroed by PetscCalloc1, so no manual initialization is needed */
1714 }
1715 }
1716 *array = data;
1717 PetscFunctionReturn(0);
1718}

◆ Deallocate3DArrayScalar()

PetscErrorCode Deallocate3DArrayScalar ( PetscReal ***  array,
PetscInt  nz,
PetscInt  ny 
)

Internal helper implementation: Deallocate3DArrayScalar().

Deallocates a 3D array of PetscReal values allocated by Allocate3DArrayScalar.

Local to this translation unit.

Definition at line 1724 of file setup.c.

1725{
1726 PetscErrorCode ierr;
1727 (void)nz;
1728 (void)ny;
1729
1730 PetscFunctionBegin;
1731 if (!array || !array[0] || !array[0][0] ) { // Added more robust check
1732 LOG_ALLOW(GLOBAL, LOG_WARNING, "Deallocate3DArrayScalar called with potentially unallocated or NULL array.\n");
1733 if (array) {
1734 if (array[0]) { // Check if row pointers might exist
1735 // Cannot safely access array[0][0] if array[0] might be invalid/freed
1736 // Standard deallocation below assumes valid pointers.
1737 ierr = PetscFree(array[0]); CHKERRQ(ierr); // Free row pointers if they exist
1738 }
1739 ierr = PetscFree(array); CHKERRQ(ierr); // Free layer pointers if they exist
1740 }
1741 PetscFunctionReturn(0);
1742 }
1743
1744 // --- Standard Deallocation (assuming valid allocation) ---
1745
1746 /* 1. Free the contiguous block of PetscReal values.
1747 The starting address was stored in array[0][0]. */
1748 ierr = PetscFree(array[0][0]); CHKERRQ(ierr); // Free the ACTUAL DATA
1749
1750 /* 2. Free the contiguous block of row pointers.
1751 The starting address was stored in array[0]. */
1752 ierr = PetscFree(array[0]); CHKERRQ(ierr); // Free the ROW POINTERS
1753
1754 /* 3. Free the layer pointer array.
1755 The starting address is 'array' itself. */
1756 ierr = PetscFree(array); CHKERRQ(ierr); // Free the LAYER POINTERS
1757
1758 PetscFunctionReturn(0);
1759}

◆ Allocate3DArrayVector()

PetscErrorCode Allocate3DArrayVector ( Cmpnts ****  array,
PetscInt  nz,
PetscInt  ny,
PetscInt  nx 
)

Implementation of Allocate3DArrayVector().

Allocates a contiguous 3D array of Cmpnts values.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/setup.h.

See also
Allocate3DArrayVector()

Definition at line 1767 of file setup.c.

1768{
1769 PetscErrorCode ierr;
1770 Cmpnts ***data;
1771 Cmpnts *dataContiguous;
1772 PetscInt k, j;
1773 PetscMPIInt rank;
1774
1775 PetscFunctionBegin;
1776
1777 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
1778
1779 /* Step 1: Allocate memory for nz layer pointers (zeroed) */
1780 ierr = PetscCalloc1(nz, &data); CHKERRQ(ierr);
1781
1782 LOG_ALLOW(LOCAL,LOG_DEBUG," [Rank %d] memory allocated for outermost layer (%d k-layer pointers).\n",rank,nz);
1783
1784 /* Step 2: Allocate memory for all row pointers (nz * ny pointers) */
1785 ierr = PetscCalloc1(nz * ny, &data[0]); CHKERRQ(ierr);
1786 for (k = 1; k < nz; k++) {
1787 data[k] = data[0] + k * ny;
1788 }
1789
1790 LOG_ALLOW(LOCAL,LOG_DEBUG,"[Rank %d] memory allocated for %dx%d row pointers.\n",rank,nz,ny);
1791
1792 /* Step 3: Allocate one contiguous block for nz*ny*nx Cmpnts structures (zeroed) */
1793 ierr = PetscCalloc1(nz * ny * nx, &dataContiguous); CHKERRQ(ierr);
1794
1795 LOG_ALLOW(GLOBAL,LOG_DEBUG,"[Rank %d] memory allocated for contigous block of %dx%dx%d Cmpnts structures).\n",rank,nz,ny,nx);
1796
1797 /* Build the 3D pointer structure for vector data */
1798 for (k = 0; k < nz; k++) {
1799 for (j = 0; j < ny; j++) {
1800 data[k][j] = dataContiguous + (k * ny + j) * nx;
1801 /* The PetscCalloc1 call has already initialized each Cmpnts to zero. */
1802 }
1803 }
1804
1805 LOG_ALLOW(GLOBAL,LOG_DEBUG,"[Rank %d] 3D pointer structure for vector data created. \n",rank);
1806
1807 *array = data;
1808 PetscFunctionReturn(0);
1809}

◆ Deallocate3DArrayVector()

PetscErrorCode Deallocate3DArrayVector ( Cmpnts ***  array,
PetscInt  nz,
PetscInt  ny 
)

Implementation of Deallocate3DArrayVector().

Deallocates a 3D array of Cmpnts structures allocated by Allocate3DArrayVector.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/setup.h.

See also
Deallocate3DArrayVector()

Definition at line 1817 of file setup.c.

1818{
1819 PetscErrorCode ierr;
1820 (void)nz;
1821 (void)ny;
1822
1823 PetscFunctionBegin;
1824 // If array is NULL or hasn't been allocated properly, just return.
1825 if (!array || !array[0] || !array[0][0] ) {
1826 LOG_ALLOW(GLOBAL, LOG_WARNING, "Deallocate3DArrayVector called with potentially unallocated or NULL array.\n");
1827 // Attempt to free what might exist, but be cautious
1828 if (array) {
1829 if (array[0]) { // Check if row pointers were allocated
1830 // We don't have a direct pointer to the contiguous data block
1831 // saved separately in this allocation scheme. The allocation relies
1832 // on array[0][0] pointing to it. If array[0] was freed first,
1833 // accessing array[0][0] is unsafe.
1834 // The allocation scheme where the contiguous data block is not
1835 // stored separately makes safe deallocation tricky if freeing
1836 // happens out of order or if parts are NULL.
1837
1838 // A SAFER ALLOCATION/DEALLOCATION would store the data pointer separately.
1839 // Given the current allocation scheme, the order MUST be:
1840 // 1. Free the data block (pointed to by array[0][0])
1841 // 2. Free the row pointer block (pointed to by array[0])
1842 // 3. Free the layer pointer block (pointed to by array)
1843
1844 // Let's assume the allocation was successful and pointers are valid.
1845 // Get pointer to the contiguous data block *before* freeing row pointers
1846 Cmpnts *dataContiguous = array[0][0];
1847 ierr = PetscFree(dataContiguous); CHKERRQ(ierr); // Free data block
1848
1849 // Now free the row pointers block
1850 ierr = PetscFree(array[0]); CHKERRQ(ierr); // Free row pointers
1851
1852 }
1853 // Finally, free the array of layer pointers
1854 ierr = PetscFree(array); CHKERRQ(ierr);
1855 }
1856 PetscFunctionReturn(0); // Return gracefully if input was NULL initially
1857 }
1858
1859
1860 // --- Standard Deallocation (assuming valid allocation) ---
1861
1862 /* 1. Free the contiguous block of Cmpnts structures.
1863 The starting address was stored in array[0][0] by Allocate3DArrayVector. */
1864 ierr = PetscFree(array[0][0]); CHKERRQ(ierr); // Free the ACTUAL DATA
1865
1866 /* 2. Free the contiguous block of row pointers.
1867 The starting address was stored in array[0]. */
1868 ierr = PetscFree(array[0]); CHKERRQ(ierr); // Free the ROW POINTERS
1869
1870 /* 3. Free the layer pointer array.
1871 The starting address is 'array' itself. */
1872 ierr = PetscFree(array); CHKERRQ(ierr); // Free the LAYER POINTERS
1873
1874 PetscFunctionReturn(0);
1875}

◆ GetOwnedCellRange()

PetscErrorCode GetOwnedCellRange ( const DMDALocalInfo *  info_nodes,
PetscInt  dim,
PetscInt *  xs_cell_global_out,
PetscInt *  xm_cell_local_out 
)

Internal helper implementation: GetOwnedCellRange().

Determines the global starting index and number of CELLS owned by the current processor in a specified dimension.

Local to this translation unit.

Definition at line 1883 of file setup.c.

1887{
1888 PetscErrorCode ierr = 0; // Standard PETSc error code, not explicitly set here but good practice.
1889 PetscInt xs_node_global_rank; // Global index of the first node owned by this rank in the specified dimension.
1890 PetscInt num_nodes_owned_rank; // Number of nodes owned by this rank in this dimension (local count, excluding ghosts).
1891 PetscInt GlobalNodesInDim_from_info; // Total number of DA points in this dimension, from DMDALocalInfo.
1892
1893 PetscFunctionBeginUser;
1894
1895 // --- 1. Input Validation ---
1896 if (!info_nodes || !xs_cell_global_out || !xm_cell_local_out) {
1897 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Null pointer passed to GetOwnedCellRange.");
1898 }
1899
1900 // --- 2. Extract Node Ownership and Global Dimension Information from DMDALocalInfo ---
1901 if (dim == 0) { // I-direction
1902 xs_node_global_rank = info_nodes->xs;
1903 num_nodes_owned_rank = info_nodes->xm;
1904 GlobalNodesInDim_from_info = info_nodes->mx;
1905 } else if (dim == 1) { // J-direction
1906 xs_node_global_rank = info_nodes->ys;
1907 num_nodes_owned_rank = info_nodes->ym;
1908 GlobalNodesInDim_from_info = info_nodes->my;
1909 } else if (dim == 2) { // K-direction
1910 xs_node_global_rank = info_nodes->zs;
1911 num_nodes_owned_rank = info_nodes->zm;
1912 GlobalNodesInDim_from_info = info_nodes->mz;
1913 } else {
1914 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d in GetOwnedCellRange. Must be 0, 1, or 2.", dim);
1915 }
1916
1917 // --- 3. Correct for User-Defined Ghost Node ---
1918 // Per the function's contract (@warning), the DA size includes an extra, non-physical
1919 // node. We subtract 1 to get the true number of physical nodes for cell calculations.
1920 const PetscInt physical_nodes_in_dim = GlobalNodesInDim_from_info - 1;
1921
1922 // --- 4. Handle Edge Cases for Physical Domain Size ---
1923 // If the physical domain has 0 or 1 node, no cells can be formed.
1924 if (physical_nodes_in_dim <= 1) {
1925 *xs_cell_global_out = xs_node_global_rank; // Still report the rank's starting node
1926 *xm_cell_local_out = 0; // But 0 cells
1927 PetscFunctionReturn(0);
1928 }
1929
1930 // --- 5. Determine Cell Ownership Based on Corrected Node Ownership ---
1931 // The first cell this rank *could* define has its origin at the first node this rank owns.
1932 *xs_cell_global_out = xs_node_global_rank;
1933
1934 // If the rank owns no nodes in this dimension, it can't form any cell origins.
1935 if (num_nodes_owned_rank == 0) {
1936 *xm_cell_local_out = 0;
1937 } else {
1938 // --- BUG FIX APPLIED HERE ---
1939 // The previous logic incorrectly assumed a cell's end node (N_{k+1}) must be on the
1940 // same rank as its origin node (N_k). The correct logic is to find the intersection
1941 // between the nodes this rank owns and the nodes that are valid origins globally.
1942
1943 // The first node owned by the rank is its first potential origin.
1944 PetscInt first_owned_origin = xs_node_global_rank;
1945
1946 // The absolute last node owned by this rank. Any node up to and including this one
1947 // is a potential cell origin from this rank's perspective.
1948 PetscInt last_node_owned_by_rank = xs_node_global_rank + num_nodes_owned_rank - 1;
1949
1950 // The absolute last node in the entire PHYSICAL domain that can serve as a cell origin.
1951 // If there are `N` physical nodes (0 to N-1), this index is `N-2`.
1952 PetscInt last_possible_origin_global_idx = physical_nodes_in_dim - 2;
1953
1954 // The actual last origin this rank can provide is the *minimum* of what it owns
1955 // and what is globally possible. This correctly handles both ranks in the middle of
1956 // the domain and the very last rank.
1957 PetscInt actual_last_origin_this_rank_can_form = PetscMin(last_node_owned_by_rank, last_possible_origin_global_idx);
1958
1959 // If the first potential origin this rank owns is already beyond the actual last
1960 // origin it can form, then this rank forms no valid cell origins. This happens if
1961 // the rank only owns the very last physical node.
1962 if (first_owned_origin > actual_last_origin_this_rank_can_form) {
1963 *xm_cell_local_out = 0;
1964 } else {
1965 // The number of cells is the count of valid origins this rank owns.
1966 // (Count = Last Index - First Index + 1)
1967 *xm_cell_local_out = actual_last_origin_this_rank_can_form - first_owned_origin + 1;
1968 }
1969 }
1970
1971 PetscFunctionReturn(ierr);
1972}
Here is the caller graph for this function:

◆ ComputeAndStoreNeighborRanks()

PetscErrorCode ComputeAndStoreNeighborRanks ( UserCtx user)

Internal helper implementation: ComputeAndStoreNeighborRanks().

Computes and stores the Cartesian neighbor ranks for the DMDA decomposition.

Local to this translation unit.

Definition at line 1980 of file setup.c.

1981{
1982 PetscErrorCode ierr;
1983 PetscMPIInt rank;
1984 PetscMPIInt size; // MPI communicator size
1985 const PetscMPIInt *neighbor_ranks_ptr; // Pointer to raw neighbor data from PETSc
1986
1987 PetscFunctionBeginUser;
1989 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
1990 ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size); CHKERRQ(ierr); // Get MPI size for validation
1991
1992 LOG_ALLOW(GLOBAL, LOG_INFO, "Rank %d: Computing DMDA neighbor ranks.\n", rank);
1993
1994 if (!user || !user->da) {
1995 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "UserCtx or user->da is NULL in ComputeAndStoreNeighborRanks.");
1996 }
1997
1998 // Get the neighbor information from the DMDA
1999 // neighbor_ranks_ptr will point to an internal PETSc array of 27 ranks.
2000 ierr = DMDAGetNeighbors(user->da, &neighbor_ranks_ptr); CHKERRQ(ierr);
2001
2002 // Log the raw values from DMDAGetNeighbors for boundary-relevant directions for debugging
2003 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",
2004 rank,
2005 neighbor_ranks_ptr[12], neighbor_ranks_ptr[14],
2006 neighbor_ranks_ptr[10], neighbor_ranks_ptr[16],
2007 neighbor_ranks_ptr[4], neighbor_ranks_ptr[22],
2008 (int)MPI_PROC_NULL);
2009
2010 // PETSc standard indices for 3D face neighbors from the 27-point stencil:
2011 // Index = k_offset*9 + j_offset*3 + i_offset (where offsets -1,0,1 map to 0,1,2)
2012 // Center: (i_off=1, j_off=1, k_off=1) => 1*9 + 1*3 + 1 = 13
2013 // X-min: (i_off=0, j_off=1, k_off=1) => 1*9 + 1*3 + 0 = 12
2014 // X-plus: (i_off=2, j_off=1, k_off=1) => 1*9 + 1*3 + 2 = 14
2015 // Y-min: (i_off=1, j_off=0, k_off=1) => 1*9 + 0*3 + 1 = 10
2016 // Y-plus: (i_off=1, j_off=2, k_off=1) => 1*9 + 2*3 + 1 = 16
2017 // Z-min: (i_off=1, j_off=1, k_off=0) => 0*9 + 1*3 + 1 = 4
2018 // Z-plus: (i_off=1, j_off=1, k_off=2) => 2*9 + 1*3 + 1 = 22
2019
2020 if (neighbor_ranks_ptr[13] != rank) {
2021 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",
2022 rank, neighbor_ranks_ptr[13], rank);
2023 // This warning is important. If the center isn't the current rank, the offsets are likely wrong.
2024 // However, PETSc should ensure this unless the DM is too small for a 3x3x3 stencil.
2025 }
2026
2027 // Assign and sanitize each neighbor rank
2028 PetscMPIInt temp_neighbor;
2029
2030 temp_neighbor = neighbor_ranks_ptr[12]; // xm
2031 if (temp_neighbor < 0 || temp_neighbor >= size) {
2032 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);
2033 user->neighbors.rank_xm = MPI_PROC_NULL;
2034 } else {
2035 user->neighbors.rank_xm = temp_neighbor;
2036 }
2037
2038 temp_neighbor = neighbor_ranks_ptr[14]; // xp
2039 if (temp_neighbor < 0 || temp_neighbor >= size) {
2040 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);
2041 user->neighbors.rank_xp = MPI_PROC_NULL;
2042 } else {
2043 user->neighbors.rank_xp = temp_neighbor;
2044 }
2045
2046 temp_neighbor = neighbor_ranks_ptr[10]; // ym
2047 if (temp_neighbor < 0 || temp_neighbor >= size) {
2048 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);
2049 user->neighbors.rank_ym = MPI_PROC_NULL;
2050 } else {
2051 user->neighbors.rank_ym = temp_neighbor;
2052 }
2053
2054 temp_neighbor = neighbor_ranks_ptr[16]; // yp
2055 if (temp_neighbor < 0 || temp_neighbor >= size) {
2056 // The log for index 16 was "zm" in your output, should be yp
2057 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);
2058 user->neighbors.rank_yp = MPI_PROC_NULL;
2059 } else {
2060 user->neighbors.rank_yp = temp_neighbor;
2061 }
2062
2063 temp_neighbor = neighbor_ranks_ptr[4]; // zm
2064 if (temp_neighbor < 0 || temp_neighbor >= size) {
2065 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);
2066 user->neighbors.rank_zm = MPI_PROC_NULL;
2067 } else {
2068 user->neighbors.rank_zm = temp_neighbor;
2069 }
2070
2071 temp_neighbor = neighbor_ranks_ptr[22]; // zp
2072 if (temp_neighbor < 0 || temp_neighbor >= size) {
2073 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);
2074 user->neighbors.rank_zp = MPI_PROC_NULL;
2075 } else {
2076 user->neighbors.rank_zp = temp_neighbor;
2077 }
2078
2079 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,
2080 user->neighbors.rank_xm, user->neighbors.rank_xp,
2081 user->neighbors.rank_ym, user->neighbors.rank_yp,
2082 user->neighbors.rank_zm, user->neighbors.rank_zp);
2083 PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT); // Ensure logs are flushed
2084
2085 // Note: neighbor_ranks_ptr memory is managed by PETSc, do not free it.
2087 PetscFunctionReturn(0);
2088}
PetscMPIInt rank_zm
Definition variables.h:182
PetscMPIInt rank_yp
Definition variables.h:181
PetscMPIInt rank_ym
Definition variables.h:181
PetscMPIInt rank_xp
Definition variables.h:180
RankNeighbors neighbors
Definition variables.h:823
PetscMPIInt rank_xm
Definition variables.h:180
PetscMPIInt rank_zp
Definition variables.h:182
Here is the caller graph for this function:

◆ SetDMDAProcLayout()

PetscErrorCode SetDMDAProcLayout ( DM  dm,
UserCtx user 
)

Internal helper implementation: SetDMDAProcLayout().

Sets the processor layout for a given DMDA based on PETSc options.

Local to this translation unit.

Definition at line 2096 of file setup.c.

2097{
2098 PetscErrorCode ierr;
2099 PetscMPIInt size, rank;
2100 PetscInt px = PETSC_DECIDE, py = PETSC_DECIDE, pz = PETSC_DECIDE;
2101 PetscBool px_set = PETSC_FALSE, py_set = PETSC_FALSE, pz_set = PETSC_FALSE;
2102 SimCtx *simCtx = user->simCtx;
2103
2104 // Set no.of processors in direction 1
2105 if(simCtx->da_procs_x) {
2106 px_set = PETSC_TRUE;
2107 px = simCtx->da_procs_x;
2108 }
2109 // Set no.of processors in direction 2
2110 if(simCtx->da_procs_y) {
2111 py_set = PETSC_TRUE;
2112 py = simCtx->da_procs_y;
2113 }
2114 // Set no.of processors in direction 1
2115 if(simCtx->da_procs_z) {
2116 pz_set = PETSC_TRUE;
2117 pz = simCtx->da_procs_z;
2118 }
2119
2120 PetscFunctionBeginUser;
2122 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size); CHKERRQ(ierr);
2123 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank); CHKERRQ(ierr);
2124 LOG_ALLOW(GLOBAL, LOG_INFO, "Rank %d: Configuring DMDA processor layout for %d total processes.\n", rank, size);
2125
2126 // --- Validate User Input (Optional but Recommended) ---
2127 // Check if specified processor counts multiply to the total MPI size
2128 if (px_set && py_set && pz_set) {
2129 if (px * py * pz != size) {
2130 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_INCOMP,
2131 "Specified processor layout %d x %d x %d = %d does not match MPI size %d",
2132 px, py, pz, px * py * pz, size);
2133 }
2134 LOG_ALLOW(GLOBAL, LOG_INFO, "Using specified processor layout: %d x %d x %d\n", px, py, pz);
2135 } else if (px_set || py_set || pz_set) {
2136 // If only some are set, PETSC_DECIDE will be used for others
2137 LOG_ALLOW(GLOBAL, LOG_INFO, "Using partially specified processor layout: %d x %d x %d (PETSC_DECIDE for unspecified)\n", px, py, pz);
2138 } else {
2139 LOG_ALLOW(GLOBAL, LOG_INFO, "Using fully automatic processor layout (PETSC_DECIDE x PETSC_DECIDE x PETSC_DECIDE)\n");
2140 }
2141 // Additional checks: Ensure px, py, pz are positive if set
2142 if ((px_set && px <= 0) || (py_set && py <= 0) || (pz_set && pz <= 0)) {
2143 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Specified processor counts must be positive.");
2144 }
2145
2146
2147 // --- Apply the layout to the DMDA ---
2148 ierr = DMDASetNumProcs(dm, px, py, pz); CHKERRQ(ierr);
2149 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Rank %d: DMDASetNumProcs called with px=%d, py=%d, pz=%d.\n", rank, px, py, pz);
2150
2151 // --- Store the values in UserCtx (Optional) ---
2152 // Note: If PETSC_DECIDE was used, PETSc calculates the actual values during DMSetUp.
2153 // We store the *requested* values here. To get the *actual* values used,
2154 // you would need to call DMDAGetInfo after DMSetUp.
2155 /*
2156 if (user) {
2157 user->procs_x = px;
2158 user->procs_y = py;
2159 user->procs_z = pz;
2160 }
2161 */
2163 PetscFunctionReturn(0);
2164}

◆ SetupDomainRankInfo()

PetscErrorCode SetupDomainRankInfo ( SimCtx simCtx)

Implementation of SetupDomainRankInfo().

Sets up the full rank communication infrastructure, including neighbor ranks and bounding box exchange.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/setup.h.

See also
SetupDomainRankInfo()

Definition at line 2174 of file setup.c.

2175{
2176 PetscErrorCode ierr;
2177 PetscInt nblk = simCtx->block_number;
2178 PetscInt size = simCtx->size;
2179 BoundingBox *final_bboxlist = NULL;
2180
2181 PetscFunctionBeginUser;
2183
2184 LOG_ALLOW(GLOBAL, LOG_INFO, "Starting full rank communication setup for %d block(s).\n", nblk);
2185
2186 UserCtx *user_finest = simCtx->usermg.mgctx[simCtx->usermg.mglevels - 1].user;
2187
2188 // --- Step 1: Compute neighbor ranks (unchanged) ---
2189 for (int bi = 0; bi < nblk; bi++) {
2190 ierr = ComputeAndStoreNeighborRanks(&user_finest[bi]); CHKERRQ(ierr);
2191 }
2192 LOG_ALLOW(GLOBAL, LOG_INFO, "Neighbor ranks computed and stored for all blocks.\n");
2193
2194 // --- Step 2: Allocate the final, unified list on ALL ranks ---
2195 // Every rank will build this list in parallel.
2196 ierr = PetscMalloc1(size * nblk, &final_bboxlist); CHKERRQ(ierr);
2197
2198 // --- Step 3: Loop through each block, gather then broadcast its bbox list ---
2199 for (int bi = 0; bi < nblk; bi++) {
2200 // This is a temporary pointer for the current block's list.
2201 BoundingBox *block_bboxlist = NULL;
2202
2203 LOG_ALLOW(GLOBAL, LOG_INFO, "Processing bounding boxes for block %d...\n", bi);
2204
2205 // A) GATHER: On rank 0, block_bboxlist is allocated and filled. On others, it's NULL.
2206 ierr = GatherAllBoundingBoxes(&user_finest[bi], &block_bboxlist); CHKERRQ(ierr);
2207 LOG_ALLOW(GLOBAL, LOG_DEBUG, " -> Gather complete for block %d.\n", bi);
2208
2209 // B) BROADCAST: On non-root ranks, block_bboxlist is allocated. Then, the data
2210 // from rank 0 is broadcast to all ranks. After this call, ALL ranks have
2211 // an identical, complete copy of the bounding boxes for the current block.
2212 ierr = BroadcastAllBoundingBoxes(&user_finest[bi], &block_bboxlist); CHKERRQ(ierr);
2213 LOG_ALLOW(GLOBAL, LOG_DEBUG, " -> Broadcast complete for block %d.\n", bi);
2214
2215 // C) ASSEMBLE: Every rank now copies the data for this block into the
2216 // correct segment of its final, unified list.
2217 for (int r = 0; r < size; r++) {
2218 // The layout is [r0b0, r1b0, ..., r(size-1)b0, r0b1, r1b1, ...]
2219 final_bboxlist[bi * size + r] = block_bboxlist[r];
2220 }
2221 LOG_ALLOW(GLOBAL, LOG_DEBUG, " -> Assembly into final list complete for block %d.\n", bi);
2222
2223 // D) CLEANUP: Free the temporary list for this block on ALL ranks before the next iteration.
2224 // Your helper functions use malloc, so we must use free.
2225 free(block_bboxlist);
2226 }
2227
2228 // --- Step 4: Assign the final pointer and run the last setup step ---
2229 simCtx->bboxlist = final_bboxlist;
2230 LOG_ALLOW(GLOBAL, LOG_INFO, "Final unified bboxlist created on all ranks and stored in SimCtx.\n");
2231
2232 ierr = SetupDomainCellDecompositionMap(&user_finest[0]); CHKERRQ(ierr);
2233 LOG_ALLOW(GLOBAL, LOG_INFO, "Domain Cell Composition set and broadcasted.\n");
2234
2235 LOG_ALLOW(GLOBAL, LOG_INFO, "SetupDomainRankInfo: Completed successfully.\n");
2236
2238 PetscFunctionReturn(0);
2239}
PetscErrorCode BroadcastAllBoundingBoxes(UserCtx *user, BoundingBox **bboxlist)
Broadcasts the bounding box information collected on rank 0 to all other ranks.
Definition grid.c:883
PetscErrorCode GatherAllBoundingBoxes(UserCtx *user, BoundingBox **allBBoxes)
Gathers local bounding boxes from all MPI processes to rank 0.
Definition grid.c:821
PetscErrorCode ComputeAndStoreNeighborRanks(UserCtx *user)
Internal helper implementation: ComputeAndStoreNeighborRanks().
Definition setup.c:1980
PetscErrorCode SetupDomainCellDecompositionMap(UserCtx *user)
Internal helper implementation: SetupDomainCellDecompositionMap().
Definition setup.c:2383
Defines a 3D axis-aligned bounding box.
Definition variables.h:154
Here is the call graph for this function:
Here is the caller graph for this function:

◆ Contra2Cart()

PetscErrorCode Contra2Cart ( UserCtx user)

Internal helper implementation: Contra2Cart().

Reconstructs Cartesian velocity (Ucat) at cell centers from contravariant velocity (Ucont) defined on cell faces.

Local to this translation unit.

Definition at line 2247 of file setup.c.

2248{
2249 PetscErrorCode ierr;
2250 DMDALocalInfo info;
2251 Cmpnts ***lcsi_arr, ***leta_arr, ***lzet_arr; // Local metric arrays
2252 Cmpnts ***lucont_arr; // Local contravariant velocity array
2253 Cmpnts ***gucat_arr; // Global Cartesian velocity array
2254 PetscReal ***lnvert_arr; // Local Nvert array
2255 PetscReal ***laj_arr; // Local Jacobian Determinant inverse array
2256
2257 PetscFunctionBeginUser;
2259 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Starting Contravariant-to-Cartesian velocity transformation.\n");
2260
2261 // --- 1. Get DMDA Info and Check for Valid Inputs ---
2262 // All inputs (lUcont, lCsi, etc.) and outputs (Ucat) are on DMs from the UserCtx.
2263 // We get local info from fda, which governs the layout of most arrays here.
2264 ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
2265 if (!user->lUcont || !user->lCsi || !user->lEta || !user->lZet || !user->lNvert || !user->Ucat) {
2266 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Contra2Cart requires lUcont, lCsi/Eta/Zet, lNvert, and Ucat to be non-NULL.");
2267 }
2268
2269
2270 // --- 2. Get Read-Only Array Access to Local Input Vectors (with ghosts) ---
2271 ierr = DMDAVecGetArrayRead(user->fda, user->lUcont, &lucont_arr); CHKERRQ(ierr);
2272 ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, &lcsi_arr); CHKERRQ(ierr);
2273 ierr = DMDAVecGetArrayRead(user->fda, user->lEta, &leta_arr); CHKERRQ(ierr);
2274 ierr = DMDAVecGetArrayRead(user->fda, user->lZet, &lzet_arr); CHKERRQ(ierr);
2275 ierr = DMDAVecGetArrayRead(user->da, user->lNvert, &lnvert_arr); CHKERRQ(ierr);
2276 ierr = DMDAVecGetArrayRead(user->da, user->lAj, &laj_arr); CHKERRQ(ierr);
2277
2278 // --- 3. Get Write-Only Array Access to the Global Output Vector ---
2279 // We compute for local owned cells and write into the global vector.
2280 // PETSc handles mapping the global indices to the correct local memory locations.
2281 ierr = DMDAVecGetArray(user->fda, user->Ucat, &gucat_arr); CHKERRQ(ierr);
2282
2283
2284 // --- 4. Define Loop Bounds for INTERIOR Cells ---
2285 // We use adjusted bounds to avoid calculating Ucat on the physical domain boundaries,
2286 // as these are typically set explicitly by boundary condition functions.
2287 // The stencils use indices like i-1, j-1, k-1, so we must start loops at least at index 1.
2288 PetscInt i_start = (info.xs == 0) ? info.xs + 1 : info.xs;
2289 PetscInt i_end = (info.xs + info.xm == info.mx) ? info.xs + info.xm - 1 : info.xs + info.xm;
2290
2291 PetscInt j_start = (info.ys == 0) ? info.ys + 1 : info.ys;
2292 PetscInt j_end = (info.ys + info.ym == info.my) ? info.ys + info.ym - 1 : info.ys + info.ym;
2293
2294 PetscInt k_start = (info.zs == 0) ? info.zs + 1 : info.zs;
2295 PetscInt k_end = (info.zs + info.zm == info.mz) ? info.zs + info.zm - 1 : info.zs + info.zm;
2296
2297 // --- 5. Main Computation Loop ---
2298 // Loops over the GLOBAL indices of interior cells owned by this rank.
2299 for (PetscInt k_cell = k_start; k_cell < k_end; ++k_cell) {
2300 for (PetscInt j_cell = j_start; j_cell < j_end; ++j_cell) {
2301 for (PetscInt i_cell = i_start; i_cell < i_end; ++i_cell) {
2302
2303 // Check if the cell is a fluid cell (not solid/blanked)
2304 // if (lnvert_arr[k_cell][j_cell][i_cell] > 0.1) continue; // Skip solid/blanked cells
2305
2306 // Transformation matrix [mat] is the metric tensor at the cell center,
2307 // estimated by averaging metrics from adjacent faces.
2308 PetscReal mat[3][3];
2309
2310 // PetscReal aj_center = laj_arr[k_cell+1][j_cell+1][i_cell+1];
2311
2312 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;
2313 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;
2314 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;
2315
2316 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;
2317 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;
2318 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;
2319
2320 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;
2321 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;
2322 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;
2323
2324 // Contravariant velocity vector `q` at the cell center,
2325 // estimated by averaging face-based contravariant velocities.
2326 PetscReal q[3];
2327 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
2328 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
2329 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
2330
2331 // Solve the 3x3 system `mat * ucat = q` using Cramer's rule.
2332 PetscReal det = mat[0][0] * (mat[1][1] * mat[2][2] - mat[1][2] * mat[2][1]) -
2333 mat[0][1] * (mat[1][0] * mat[2][2] - mat[1][2] * mat[2][0]) +
2334 mat[0][2] * (mat[1][0] * mat[2][1] - mat[1][1] * mat[2][0]);
2335
2336 if (PetscAbsReal(det) < 1.0e-18) {
2337 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);
2338 }
2339
2340 PetscReal det_inv = 1.0 / det;
2341
2342 PetscReal det0 = q[0] * (mat[1][1] * mat[2][2] - mat[1][2] * mat[2][1]) -
2343 q[1] * (mat[0][1] * mat[2][2] - mat[0][2] * mat[2][1]) +
2344 q[2] * (mat[0][1] * mat[1][2] - mat[0][2] * mat[1][1]);
2345
2346 PetscReal det1 = -q[0] * (mat[1][0] * mat[2][2] - mat[1][2] * mat[2][0]) +
2347 q[1] * (mat[0][0] * mat[2][2] - mat[0][2] * mat[2][0]) -
2348 q[2] * (mat[0][0] * mat[1][2] - mat[0][2] * mat[1][0]);
2349
2350 PetscReal det2 = q[0] * (mat[1][0] * mat[2][1] - mat[1][1] * mat[2][0]) -
2351 q[1] * (mat[0][0] * mat[2][1] - mat[0][1] * mat[2][0]) +
2352 q[2] * (mat[0][0] * mat[1][1] - mat[0][1] * mat[1][0]);
2353
2354 // Store computed Cartesian velocity in the GLOBAL Ucat array at the
2355 // array index corresponding to the cell's origin node.
2356 gucat_arr[k_cell][j_cell][i_cell].x = det0 * det_inv;
2357 gucat_arr[k_cell][j_cell][i_cell].y = det1 * det_inv;
2358 gucat_arr[k_cell][j_cell][i_cell].z = det2 * det_inv;
2359 }
2360 }
2361 }
2362
2363 // --- 6. Restore Array Access ---
2364 ierr = DMDAVecRestoreArrayRead(user->fda, user->lUcont, &lucont_arr); CHKERRQ(ierr);
2365 ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, &lcsi_arr); CHKERRQ(ierr);
2366 ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, &leta_arr); CHKERRQ(ierr);
2367 ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, &lzet_arr); CHKERRQ(ierr);
2368 ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, &lnvert_arr); CHKERRQ(ierr);
2369 ierr = DMDAVecRestoreArrayRead(user->da, user->lAj, &laj_arr); CHKERRQ(ierr);
2370 ierr = DMDAVecRestoreArray(user->fda, user->Ucat, &gucat_arr); CHKERRQ(ierr);
2371
2372 LOG_ALLOW(GLOBAL, LOG_INFO, "Completed Contravariant-to-Cartesian velocity transformation. \n");
2374 PetscFunctionReturn(0);
2375}
Here is the caller graph for this function:

◆ SetupDomainCellDecompositionMap()

PetscErrorCode SetupDomainCellDecompositionMap ( UserCtx user)

Internal helper implementation: SetupDomainCellDecompositionMap().

Creates and distributes a map of the domain's cell decomposition to all ranks.

Local to this translation unit.

Definition at line 2383 of file setup.c.

2384{
2385 PetscErrorCode ierr;
2386 DMDALocalInfo local_node_info;
2387 RankCellInfo my_cell_info;
2388 PetscMPIInt rank, size;
2389
2390 PetscFunctionBeginUser;
2392
2393 // --- 1. Input Validation and MPI Info ---
2394 if (!user) {
2395 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "UserCtx pointer is NULL in SetupDomainCellDecompositionMap.");
2396 }
2397 if (!user->da) {
2398 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "user->da is not initialized in SetupDomainCellDecompositionMap.");
2399 }
2400
2401 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
2402 ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size); CHKERRQ(ierr);
2403
2404 LOG_ALLOW(GLOBAL, LOG_INFO, "Setting up domain cell decomposition map for %d ranks.\n", size);
2405
2406 // --- 2. Determine Local Cell Ownership ---
2407 // Get the local node ownership information from the primary DMDA.
2408 ierr = DMDAGetLocalInfo(user->da, &local_node_info); CHKERRQ(ierr);
2409
2410 // Use the robust helper function to convert node ownership to cell ownership.
2411 // A cell's index is defined by its origin node.
2412
2413 ierr = GetOwnedCellRange(&local_node_info, 0, &my_cell_info.xs_cell, &my_cell_info.xm_cell); CHKERRQ(ierr);
2414 ierr = GetOwnedCellRange(&local_node_info, 1, &my_cell_info.ys_cell, &my_cell_info.ym_cell); CHKERRQ(ierr);
2415 ierr = GetOwnedCellRange(&local_node_info, 2, &my_cell_info.zs_cell, &my_cell_info.zm_cell); CHKERRQ(ierr);
2416
2417 // Log the calculated local ownership for debugging purposes.
2418 LOG_ALLOW(LOCAL, LOG_DEBUG, "[Rank %d] Owns cells: i[%d, %d), j[%d, %d), k[%d, %d)\n",
2419 rank, my_cell_info.xs_cell, my_cell_info.xs_cell + my_cell_info.xm_cell,
2420 my_cell_info.ys_cell, my_cell_info.ys_cell + my_cell_info.ym_cell,
2421 my_cell_info.zs_cell, my_cell_info.zs_cell + my_cell_info.zm_cell);
2422
2423 // --- 3. Allocate and Distribute the Global Map ---
2424 // Allocate memory for the global map that will hold information from all ranks.
2425 ierr = PetscMalloc1(size, &user->RankCellInfoMap); CHKERRQ(ierr);
2426
2427 // Perform the collective communication to gather the `RankCellInfo` struct from every rank.
2428 // Each rank sends its `my_cell_info` and receives the complete array in `user->RankCellInfoMap`.
2429 // We use MPI_BYTE to ensure portability across different systems and struct padding.
2430 ierr = MPI_Allgather(&my_cell_info, sizeof(RankCellInfo), MPI_BYTE,
2431 user->RankCellInfoMap, sizeof(RankCellInfo), MPI_BYTE,
2432 PETSC_COMM_WORLD); CHKERRQ(ierr);
2433
2434 LOG_ALLOW(GLOBAL, LOG_INFO, "Domain cell decomposition map created and distributed successfully.\n");
2435
2437 PetscFunctionReturn(0);
2438}
PetscErrorCode GetOwnedCellRange(const DMDALocalInfo *info_nodes, PetscInt dim, PetscInt *xs_cell_global_out, PetscInt *xm_cell_local_out)
Internal helper implementation: GetOwnedCellRange().
Definition setup.c:1883
PetscInt ys_cell
Definition variables.h:187
PetscInt xs_cell
Definition variables.h:187
PetscInt zm_cell
Definition variables.h:188
PetscInt zs_cell
Definition variables.h:187
PetscInt xm_cell
Definition variables.h:188
RankCellInfo * RankCellInfoMap
Definition variables.h:881
PetscInt ym_cell
Definition variables.h:188
A lean struct to hold the global cell ownership range for a single MPI rank.
Definition variables.h:186
Here is the call graph for this function:
Here is the caller graph for this function:

◆ BinarySearchInt64()

PetscErrorCode BinarySearchInt64 ( PetscInt  n,
const PetscInt64  arr[],
PetscInt64  key,
PetscBool *  found 
)

Implementation of BinarySearchInt64().

Performs a binary search for a key in a sorted array of PetscInt64.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/setup.h.

See also
BinarySearchInt64()

Definition at line 2448 of file setup.c.

2449{
2450 PetscInt low = 0, high = n - 1;
2451
2452 PetscFunctionBeginUser;
2454
2455 // --- 1. Input Validation ---
2456 if (!found) {
2457 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Output pointer 'found' is NULL in PetscBinarySearchInt64.");
2458 }
2459 if (n > 0 && !arr) {
2460 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Input array 'arr' is NULL for n > 0.");
2461 }
2462
2463 // Initialize output
2464 *found = PETSC_FALSE;
2465
2466 // --- 2. Binary Search Algorithm ---
2467 while (low <= high) {
2468 // Use this form to prevent potential integer overflow on very large arrays
2469 PetscInt mid = low + (high - low) / 2;
2470
2471 if (arr[mid] == key) {
2472 *found = PETSC_TRUE; // Key found!
2473 break; // Exit the loop
2474 }
2475
2476 if (arr[mid] < key) {
2477 low = mid + 1; // Search in the right half
2478 } else {
2479 high = mid - 1; // Search in the left half
2480 }
2481 }
2482
2484 PetscFunctionReturn(0);
2485}
Here is the caller graph for this function:

◆ Gidx()

static PetscInt Gidx ( PetscInt  i,
PetscInt  j,
PetscInt  k,
UserCtx user 
)
static

Internal helper implementation: Gidx().

Local to this translation unit.

Definition at line 2492 of file setup.c.

2493{
2494 PetscInt nidx;
2495 DMDALocalInfo info = user->info;
2496
2497 PetscInt mx = info.mx, my = info.my;
2498
2499 AO ao;
2500 DMDAGetAO(user->da, &ao);
2501 nidx=i+j*mx+k*mx*my;
2502
2503 AOApplicationToPetsc(ao,1,&nidx);
2504
2505 return (nidx);
2506}
DMDALocalInfo info
Definition variables.h:818
Here is the caller graph for this function:

◆ ComputeDivergence()

PetscErrorCode ComputeDivergence ( UserCtx user)

Implementation of ComputeDivergence().

Computes the discrete divergence of the contravariant velocity field.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/setup.h.

See also
ComputeDivergence()

Definition at line 2518 of file setup.c.

2519{
2520 DM da = user->da, fda = user->fda;
2521 DMDALocalInfo info = user->info;
2522
2523 PetscInt ti = user->simCtx->step;
2524
2525 PetscInt xs = info.xs, xe = info.xs + info.xm;
2526 PetscInt ys = info.ys, ye = info.ys + info.ym;
2527 PetscInt zs = info.zs, ze = info.zs + info.zm;
2528 PetscInt mx = info.mx, my = info.my, mz = info.mz;
2529
2530 PetscInt lxs, lys, lzs, lxe, lye, lze;
2531 PetscInt i, j, k;
2532
2533 Vec Div;
2534 PetscReal ***div, ***aj, ***nvert,***p;
2535 Cmpnts ***ucont;
2536 PetscReal maxdiv;
2537
2538 lxs = xs; lxe = xe;
2539 lys = ys; lye = ye;
2540 lzs = zs; lze = ze;
2541
2542 if (xs==0) lxs = xs+1;
2543 if (ys==0) lys = ys+1;
2544 if (zs==0) lzs = zs+1;
2545
2546 if (xe==mx) lxe = xe-1;
2547 if (ye==my) lye = ye-1;
2548 if (ze==mz) lze = ze-1;
2549
2550 PetscFunctionBeginUser;
2552
2553 DMDAVecGetArray(fda,user->lUcont, &ucont);
2554 DMDAVecGetArray(da, user->lAj, &aj);
2555 VecDuplicate(user->P, &Div);
2556 DMDAVecGetArray(da, Div, &div);
2557 DMDAVecGetArray(da, user->lNvert, &nvert);
2558 DMDAVecGetArray(da, user->P, &p);
2559 for (k=lzs; k<lze; k++) {
2560 for (j=lys; j<lye; j++){
2561 for (i=lxs; i<lxe; i++) {
2562 if (k==10 && j==10 && i==1){
2563 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]);
2564 }
2565
2566 if (k==10 && j==10 && i==mx-3)
2567 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]);
2568 }
2569 }
2570 }
2571 DMDAVecRestoreArray(da, user->P, &p);
2572
2573
2574 for (k=lzs; k<lze; k++) {
2575 for (j=lys; j<lye; j++) {
2576 for (i=lxs; i<lxe; i++) {
2577 maxdiv = fabs((ucont[k][j][i].x - ucont[k][j][i-1].x +
2578 ucont[k][j][i].y - ucont[k][j-1][i].y +
2579 ucont[k][j][i].z - ucont[k-1][j][i].z)*aj[k][j][i]);
2580 if (nvert[k][j][i] + nvert[k+1][j][i] + nvert[k-1][j][i] +
2581 nvert[k][j+1][i] + nvert[k][j-1][i] +
2582 nvert[k][j][i+1] + nvert[k][j][i-1] > 0.1) maxdiv = 0.;
2583 div[k][j][i] = maxdiv;
2584
2585 }
2586 }
2587 }
2588
2589 if (zs==0) {
2590 k=0;
2591 for (j=ys; j<ye; j++) {
2592 for (i=xs; i<xe; i++) {
2593 div[k][j][i] = 0.;
2594 }
2595 }
2596 }
2597
2598 if (ze == mz) {
2599 k=mz-1;
2600 for (j=ys; j<ye; j++) {
2601 for (i=xs; i<xe; i++) {
2602 div[k][j][i] = 0.;
2603 }
2604 }
2605 }
2606
2607 if (xs==0) {
2608 i=0;
2609 for (k=zs; k<ze; k++) {
2610 for (j=ys; j<ye; j++) {
2611 div[k][j][i] = 0.;
2612 }
2613 }
2614 }
2615
2616 if (xe==mx) {
2617 i=mx-1;
2618 for (k=zs; k<ze; k++) {
2619 for (j=ys; j<ye; j++) {
2620 div[k][j][i] = 0;
2621 }
2622 }
2623 }
2624
2625 if (ys==0) {
2626 j=0;
2627 for (k=zs; k<ze; k++) {
2628 for (i=xs; i<xe; i++) {
2629 div[k][j][i] = 0.;
2630 }
2631 }
2632 }
2633
2634 if (ye==my) {
2635 j=my-1;
2636 for (k=zs; k<ze; k++) {
2637 for (i=xs; i<xe; i++) {
2638 div[k][j][i] = 0.;
2639 }
2640 }
2641 }
2642 DMDAVecRestoreArray(da, Div, &div);
2643 PetscInt MaxFlatIndex;
2644
2645 VecMax(Div, &MaxFlatIndex, &maxdiv);
2646
2647 LOG_ALLOW(GLOBAL,LOG_INFO,"[Step %d]] The Maximum Divergence is %e at flat index %d.\n",ti,maxdiv,MaxFlatIndex);
2648
2649 user->simCtx->MaxDivFlatArg = MaxFlatIndex;
2650 user->simCtx->MaxDiv = maxdiv;
2651
2652 for (k=zs; k<ze; k++) {
2653 for (j=ys; j<ye; j++) {
2654 for (i=xs; i<xe; i++) {
2655 if (Gidx(i,j,k,user) == MaxFlatIndex) {
2656 LOG_ALLOW(GLOBAL,LOG_INFO,"[Step %d] The Maximum Divergence(%e) is at location [%d][%d][%d]. \n", ti, maxdiv,k,j,i);
2657 user->simCtx->MaxDivz = k;
2658 user->simCtx->MaxDivy = j;
2659 user->simCtx->MaxDivx = i;
2660 }
2661 }
2662 }
2663 }
2664
2665
2666 DMDAVecRestoreArray(da, user->lNvert, &nvert);
2667 DMDAVecRestoreArray(fda, user->lUcont, &ucont);
2668 DMDAVecRestoreArray(da, user->lAj, &aj);
2669 VecDestroy(&Div);
2670
2672 PetscFunctionReturn(0);
2673}
static PetscInt Gidx(PetscInt i, PetscInt j, PetscInt k, UserCtx *user)
Internal helper implementation: Gidx().
Definition setup.c:2492
Here is the call graph for this function:
Here is the caller graph for this function:

◆ InitializeRandomGenerators()

PetscErrorCode InitializeRandomGenerators ( UserCtx user,
PetscRandom *  randx,
PetscRandom *  randy,
PetscRandom *  randz 
)

Implementation of InitializeRandomGenerators().

Initializes random number generators for assigning particle properties.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/setup.h.

See also
InitializeRandomGenerators()

Definition at line 2684 of file setup.c.

2684 {
2685 PetscErrorCode ierr; // Error code for PETSc functions
2686 PetscMPIInt rank;
2687 PetscFunctionBeginUser;
2689 MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
2690
2691 // Initialize RNG for x-coordinate
2692 ierr = PetscRandomCreate(PETSC_COMM_SELF, randx); CHKERRQ(ierr);
2693 ierr = PetscRandomSetType((*randx), PETSCRAND48); CHKERRQ(ierr);
2694 ierr = PetscRandomSetInterval(*randx, user->bbox.min_coords.x, user->bbox.max_coords.x); CHKERRQ(ierr);
2695 ierr = PetscRandomSetSeed(*randx, rank + 12345); CHKERRQ(ierr); // Unique seed per rank
2696 ierr = PetscRandomSeed(*randx); CHKERRQ(ierr);
2697 LOG_ALLOW_SYNC(LOCAL,LOG_VERBOSE, "[Rank %d]Initialized RNG for X-axis.\n",rank);
2698
2699 // Initialize RNG for y-coordinate
2700 ierr = PetscRandomCreate(PETSC_COMM_SELF, randy); CHKERRQ(ierr);
2701 ierr = PetscRandomSetType((*randy), PETSCRAND48); CHKERRQ(ierr);
2702 ierr = PetscRandomSetInterval(*randy, user->bbox.min_coords.y, user->bbox.max_coords.y); CHKERRQ(ierr);
2703 ierr = PetscRandomSetSeed(*randy, rank + 67890); CHKERRQ(ierr); // Unique seed per rank
2704 ierr = PetscRandomSeed(*randy); CHKERRQ(ierr);
2705 LOG_ALLOW_SYNC(LOCAL,LOG_VERBOSE, "[Rank %d]Initialized RNG for Y-axis.\n",rank);
2706
2707 // Initialize RNG for z-coordinate
2708 ierr = PetscRandomCreate(PETSC_COMM_SELF, randz); CHKERRQ(ierr);
2709 ierr = PetscRandomSetType((*randz), PETSCRAND48); CHKERRQ(ierr);
2710 ierr = PetscRandomSetInterval(*randz, user->bbox.min_coords.z, user->bbox.max_coords.z); CHKERRQ(ierr);
2711 ierr = PetscRandomSetSeed(*randz, rank + 54321); CHKERRQ(ierr); // Unique seed per rank
2712 ierr = PetscRandomSeed(*randz); CHKERRQ(ierr);
2713 LOG_ALLOW_SYNC(LOCAL,LOG_VERBOSE, "[Rank %d]Initialized RNG for Z-axis.\n",rank);
2714
2716 PetscFunctionReturn(0);
2717}
@ LOG_VERBOSE
Extremely detailed logs, typically for development use only.
Definition logging.h:33
Cmpnts max_coords
Maximum x, y, z coordinates of the bounding box.
Definition variables.h:156
Cmpnts min_coords
Minimum x, y, z coordinates of the bounding box.
Definition variables.h:155
BoundingBox bbox
Definition variables.h:822
Here is the caller graph for this function:

◆ InitializeLogicalSpaceRNGs()

PetscErrorCode InitializeLogicalSpaceRNGs ( PetscRandom *  rand_logic_i,
PetscRandom *  rand_logic_j,
PetscRandom *  rand_logic_k 
)

Internal helper implementation: InitializeLogicalSpaceRNGs().

Initializes random number generators for logical space operations [0.0, 1.0).

Local to this translation unit.

Definition at line 2725 of file setup.c.

2725 {
2726 PetscErrorCode ierr;
2727 PetscMPIInt rank;
2728 PetscFunctionBeginUser;
2729
2731
2732 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
2733
2734 // --- RNG for i-logical dimension ---
2735 ierr = PetscRandomCreate(PETSC_COMM_SELF, rand_logic_i); CHKERRQ(ierr);
2736 ierr = PetscRandomSetType((*rand_logic_i), PETSCRAND48); CHKERRQ(ierr);
2737 ierr = PetscRandomSetInterval(*rand_logic_i, 0.0, 1.0); CHKERRQ(ierr); // Key change: [0,1)
2738 ierr = PetscRandomSetSeed(*rand_logic_i, rank + 202401); CHKERRQ(ierr); // Unique seed
2739 ierr = PetscRandomSeed(*rand_logic_i); CHKERRQ(ierr);
2740 LOG_ALLOW(LOCAL,LOG_VERBOSE, "[Rank %d] Initialized RNG for i-logical dimension [0,1).\n",rank);
2741
2742 // --- RNG for j-logical dimension ---
2743 ierr = PetscRandomCreate(PETSC_COMM_SELF, rand_logic_j); CHKERRQ(ierr);
2744 ierr = PetscRandomSetType((*rand_logic_j), PETSCRAND48); CHKERRQ(ierr);
2745 ierr = PetscRandomSetInterval(*rand_logic_j, 0.0, 1.0); CHKERRQ(ierr); // Key change: [0,1)
2746 ierr = PetscRandomSetSeed(*rand_logic_j, rank + 202402); CHKERRQ(ierr);
2747 ierr = PetscRandomSeed(*rand_logic_j); CHKERRQ(ierr);
2748 LOG_ALLOW(LOCAL,LOG_VERBOSE, "[Rank %d] Initialized RNG for j-logical dimension [0,1).\n",rank);
2749
2750 // --- RNG for k-logical dimension ---
2751 ierr = PetscRandomCreate(PETSC_COMM_SELF, rand_logic_k); CHKERRQ(ierr);
2752 ierr = PetscRandomSetType((*rand_logic_k), PETSCRAND48); CHKERRQ(ierr);
2753 ierr = PetscRandomSetInterval(*rand_logic_k, 0.0, 1.0); CHKERRQ(ierr); // Key change: [0,1)
2754 ierr = PetscRandomSetSeed(*rand_logic_k, rank + 202403); CHKERRQ(ierr);
2755 ierr = PetscRandomSeed(*rand_logic_k); CHKERRQ(ierr);
2756 LOG_ALLOW(LOCAL,LOG_VERBOSE, "[Rank %d] Initialized RNG for k-logical dimension [0,1).\n",rank);
2757
2758
2760 PetscFunctionReturn(0);
2761}
Here is the caller graph for this function:

◆ InitializeBrownianRNG()

PetscErrorCode InitializeBrownianRNG ( SimCtx simCtx)

Internal helper implementation: InitializeBrownianRNG().

Initializes a single master RNG for time-stepping physics (Brownian motion).

Local to this translation unit.

Definition at line 2769 of file setup.c.

2769 {
2770 PetscErrorCode ierr;
2771 PetscMPIInt rank;
2772
2773 PetscFunctionBeginUser;
2775
2776 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
2777
2778 // 1. Create the generator (stored in SimCtx, not UserCtx, as it is global physics)
2779 ierr = PetscRandomCreate(PETSC_COMM_WORLD, &simCtx->BrownianMotionRNG); CHKERRQ(ierr);
2780 ierr = PetscRandomSetType(simCtx->BrownianMotionRNG, PETSCRAND48); CHKERRQ(ierr);
2781
2782 // 2. CRITICAL: Set interval to [0, 1).
2783 // This is required for the Gaussian math to work.
2784 ierr = PetscRandomSetInterval(simCtx->BrownianMotionRNG, 0.0, 1.0); CHKERRQ(ierr);
2785
2786 // 3. Seed based on Rank to ensure spatial randomness
2787 // Multiplying by a large prime helps separate the streams significantly
2788 unsigned long seed = (unsigned long)rank * 987654321 + (unsigned long)time(NULL);
2789 ierr = PetscRandomSetSeed(simCtx->BrownianMotionRNG, seed); CHKERRQ(ierr);
2790 ierr = PetscRandomSeed(simCtx->BrownianMotionRNG); CHKERRQ(ierr);
2791
2792 LOG_ALLOW(LOCAL, LOG_VERBOSE, "[Rank %d] Initialized Brownian Physics RNG.\n", rank);
2793
2795 PetscFunctionReturn(0);
2796}
Here is the caller graph for this function:

◆ TransformScalarDerivativesToPhysical()

void TransformScalarDerivativesToPhysical ( PetscReal  jacobian,
Cmpnts  csi_metrics,
Cmpnts  eta_metrics,
Cmpnts  zet_metrics,
PetscReal  dPhi_dcsi,
PetscReal  dPhi_deta,
PetscReal  dPhi_dzet,
Cmpnts gradPhi 
)

Implementation of TransformScalarDerivativesToPhysical().

Transforms scalar derivatives from computational space to physical space.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/setup.h.

See also
TransformScalarDerivativesToPhysical()

Definition at line 2808 of file setup.c.

2816{
2817 // Gradient X component
2818 gradPhi->x = jacobian * (dPhi_dcsi * csi_metrics.x + dPhi_deta * eta_metrics.x + dPhi_dzet * zet_metrics.x);
2819
2820 // Gradient Y component
2821 gradPhi->y = jacobian * (dPhi_dcsi * csi_metrics.y + dPhi_deta * eta_metrics.y + dPhi_dzet * zet_metrics.y);
2822
2823 // Gradient Z component
2824 gradPhi->z = jacobian * (dPhi_dcsi * csi_metrics.z + dPhi_deta * eta_metrics.z + dPhi_dzet * zet_metrics.z);
2825}
Here is the caller graph for this function:

◆ 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 
)
static

Internal helper implementation: TransformDerivativesToPhysical().

Local to this translation unit.

Definition at line 2833 of file setup.c.

2836{
2837 // Derivatives of the first component (u)
2838 dudx->x = jacobian * (deriv_csi.x * csi_metrics.x + deriv_eta.x * eta_metrics.x + deriv_zet.x * zet_metrics.x);
2839 dudx->y = jacobian * (deriv_csi.x * csi_metrics.y + deriv_eta.x * eta_metrics.y + deriv_zet.x * zet_metrics.y);
2840 dudx->z = jacobian * (deriv_csi.x * csi_metrics.z + deriv_eta.x * eta_metrics.z + deriv_zet.x * zet_metrics.z);
2841 // Derivatives of the second component (v)
2842 dvdx->x = jacobian * (deriv_csi.y * csi_metrics.x + deriv_eta.y * eta_metrics.x + deriv_zet.y * zet_metrics.x);
2843 dvdx->y = jacobian * (deriv_csi.y * csi_metrics.y + deriv_eta.y * eta_metrics.y + deriv_zet.y * zet_metrics.y);
2844 dvdx->z = jacobian * (deriv_csi.y * csi_metrics.z + deriv_eta.y * eta_metrics.z + deriv_zet.y * zet_metrics.z);
2845 // Derivatives of the third component (w)
2846 dwdx->x = jacobian * (deriv_csi.z * csi_metrics.x + deriv_eta.z * eta_metrics.x + deriv_zet.z * zet_metrics.x);
2847 dwdx->y = jacobian * (deriv_csi.z * csi_metrics.y + deriv_eta.z * eta_metrics.y + deriv_zet.z * zet_metrics.y);
2848 dwdx->z = jacobian * (deriv_csi.z * csi_metrics.z + deriv_eta.z * eta_metrics.z + deriv_zet.z * zet_metrics.z);
2849}
Here is the caller graph for this function:

◆ ComputeScalarFieldDerivatives()

PetscErrorCode ComputeScalarFieldDerivatives ( UserCtx user,
PetscInt  i,
PetscInt  j,
PetscInt  k,
PetscReal ***  field_data,
Cmpnts grad 
)

Internal helper implementation: ComputeScalarFieldDerivatives().

Computes the gradient of a cell-centered SCALAR field at a specific grid point.

Local to this translation unit.

Definition at line 2857 of file setup.c.

2859{
2860 PetscErrorCode ierr;
2861 Cmpnts ***csi, ***eta, ***zet;
2862 PetscReal ***jac;
2863 PetscReal d_csi, d_eta, d_zet;
2864
2865 PetscFunctionBeginUser;
2866
2867 // 1. Get read-only access to metrics
2868 ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, &csi); CHKERRQ(ierr);
2869 ierr = DMDAVecGetArrayRead(user->fda, user->lEta, &eta); CHKERRQ(ierr);
2870 ierr = DMDAVecGetArrayRead(user->fda, user->lZet, &zet); CHKERRQ(ierr);
2871 ierr = DMDAVecGetArrayRead(user->da, user->lAj, &jac); CHKERRQ(ierr);
2872
2873 // 2. Compute derivatives in computational space (Central Difference)
2874 // Assumes ghosts are available at i+/-1
2875 d_csi = 0.5 * (field_data[k][j][i+1] - field_data[k][j][i-1]);
2876 d_eta = 0.5 * (field_data[k][j+1][i] - field_data[k][j-1][i]);
2877 d_zet = 0.5 * (field_data[k+1][j][i] - field_data[k-1][j][i]);
2878
2879 // 3. Transform to physical space
2881 csi[k][j][i], eta[k][j][i], zet[k][j][i],
2882 d_csi, d_eta, d_zet,
2883 grad);
2884
2885 // 4. Restore arrays
2886 ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, &csi); CHKERRQ(ierr);
2887 ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, &eta); CHKERRQ(ierr);
2888 ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, &zet); CHKERRQ(ierr);
2889 ierr = DMDAVecRestoreArrayRead(user->da, user->lAj, &jac); CHKERRQ(ierr);
2890
2891 PetscFunctionReturn(0);
2892}
void TransformScalarDerivativesToPhysical(PetscReal jacobian, Cmpnts csi_metrics, Cmpnts eta_metrics, Cmpnts zet_metrics, PetscReal dPhi_dcsi, PetscReal dPhi_deta, PetscReal dPhi_dzet, Cmpnts *gradPhi)
Implementation of TransformScalarDerivativesToPhysical().
Definition setup.c:2808
Here is the call graph for this function:

◆ ComputeVectorFieldDerivatives()

PetscErrorCode ComputeVectorFieldDerivatives ( UserCtx user,
PetscInt  i,
PetscInt  j,
PetscInt  k,
Cmpnts ***  field_data,
Cmpnts dudx,
Cmpnts dvdx,
Cmpnts dwdx 
)

Internal helper implementation: ComputeVectorFieldDerivatives().

Computes the derivatives of a cell-centered vector field at a specific grid point.

Local to this translation unit.

Definition at line 2900 of file setup.c.

2902{
2903 PetscErrorCode ierr;
2904 Cmpnts ***csi, ***eta, ***zet;
2905 PetscReal ***jac;
2906 PetscFunctionBeginUser;
2907
2908 // 1. Get read-only access to the necessary metric data arrays
2909 ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, &csi); CHKERRQ(ierr);
2910 ierr = DMDAVecGetArrayRead(user->fda, user->lEta, &eta); CHKERRQ(ierr);
2911 ierr = DMDAVecGetArrayRead(user->fda, user->lZet, &zet); CHKERRQ(ierr);
2912 ierr = DMDAVecGetArrayRead(user->da, user->lAj, &jac); CHKERRQ(ierr);
2913
2914 // 2. Calculate derivatives in computational space using central differencing
2915 Cmpnts deriv_csi, deriv_eta, deriv_zet;
2916 deriv_csi.x = (field_data[k][j][i+1].x - field_data[k][j][i-1].x) * 0.5;
2917 deriv_csi.y = (field_data[k][j][i+1].y - field_data[k][j][i-1].y) * 0.5;
2918 deriv_csi.z = (field_data[k][j][i+1].z - field_data[k][j][i-1].z) * 0.5;
2919
2920 deriv_eta.x = (field_data[k][j+1][i].x - field_data[k][j-1][i].x) * 0.5;
2921 deriv_eta.y = (field_data[k][j+1][i].y - field_data[k][j-1][i].y) * 0.5;
2922 deriv_eta.z = (field_data[k][j+1][i].z - field_data[k][j-1][i].z) * 0.5;
2923
2924 deriv_zet.x = (field_data[k+1][j][i].x - field_data[k-1][j][i].x) * 0.5;
2925 deriv_zet.y = (field_data[k+1][j][i].y - field_data[k-1][j][i].y) * 0.5;
2926 deriv_zet.z = (field_data[k+1][j][i].z - field_data[k-1][j][i].z) * 0.5;
2927
2928 // 3. Transform derivatives to physical space
2929 TransformDerivativesToPhysical(jac[k][j][i], csi[k][j][i], eta[k][j][i], zet[k][j][i],
2930 deriv_csi, deriv_eta, deriv_zet,
2931 dudx, dvdx, dwdx);
2932
2933 // 4. Restore access to the PETSc data arrays
2934 ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, &csi); CHKERRQ(ierr);
2935 ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, &eta); CHKERRQ(ierr);
2936 ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, &zet); CHKERRQ(ierr);
2937 ierr = DMDAVecRestoreArrayRead(user->da, user->lAj, &jac); CHKERRQ(ierr);
2938
2939 PetscFunctionReturn(0);
2940}
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)
Internal helper implementation: TransformDerivativesToPhysical().
Definition setup.c:2833
Here is the call graph for this function:
Here is the caller graph for this function:

◆ DestroyUserVectors()

PetscErrorCode DestroyUserVectors ( UserCtx user)

Internal helper implementation: DestroyUserVectors().

Destroys all PETSc Vec objects within a single UserCtx structure.

Local to this translation unit.

Definition at line 2954 of file setup.c.

2955{
2956 PetscErrorCode ierr;
2957 PetscFunctionBeginUser;
2958
2959 // --- Group A: Primary Flow Fields (Always allocated at all levels) ---
2960 if (user->Ucont) { ierr = VecDestroy(&user->Ucont); CHKERRQ(ierr); }
2961 if (user->lUcont) { ierr = VecDestroy(&user->lUcont); CHKERRQ(ierr); }
2962 if (user->Ucat) { ierr = VecDestroy(&user->Ucat); CHKERRQ(ierr); }
2963 if (user->lUcat) { ierr = VecDestroy(&user->lUcat); CHKERRQ(ierr); }
2964 if (user->P) { ierr = VecDestroy(&user->P); CHKERRQ(ierr); }
2965 if (user->lP) { ierr = VecDestroy(&user->lP); CHKERRQ(ierr); }
2966 if (user->Nvert) { ierr = VecDestroy(&user->Nvert); CHKERRQ(ierr); }
2967 if (user->lNvert) { ierr = VecDestroy(&user->lNvert); CHKERRQ(ierr); }
2968
2969 // --- Group A2: Derived Flow Fields (Conditional) ---
2970 if(user->Diffusivity) {ierr = VecDestroy(&user->Diffusivity); CHKERRQ(ierr);}
2971 if(user->lDiffusivity){ierr = VecDestroy(&user->lDiffusivity); CHKERRQ(ierr);}
2972 if(user->DiffusivityGradient){ierr = VecDestroy(&user->DiffusivityGradient); CHKERRQ(ierr);}
2973 if(user->lDiffusivityGradient){ierr = VecDestroy(&user->lDiffusivityGradient); CHKERRQ(ierr);}
2974
2975 // --- Group B: Solver Work Vectors (All levels) ---
2976 if (user->Phi) { ierr = VecDestroy(&user->Phi); CHKERRQ(ierr); }
2977 if (user->lPhi) { ierr = VecDestroy(&user->lPhi); CHKERRQ(ierr); }
2978
2979 // --- Group C: Time-Stepping Vectors (Finest level only) ---
2980 if (user->Ucont_o) { ierr = VecDestroy(&user->Ucont_o); CHKERRQ(ierr); }
2981 if (user->Ucont_rm1) { ierr = VecDestroy(&user->Ucont_rm1); CHKERRQ(ierr); }
2982 if (user->Ucat_o) { ierr = VecDestroy(&user->Ucat_o); CHKERRQ(ierr); }
2983 if (user->P_o) { ierr = VecDestroy(&user->P_o); CHKERRQ(ierr); }
2984 if (user->Nvert_o) { ierr = VecDestroy(&user->Nvert_o); CHKERRQ(ierr); }
2985 if (user->lUcont_o) { ierr = VecDestroy(&user->lUcont_o); CHKERRQ(ierr); }
2986 if (user->lUcont_rm1) { ierr = VecDestroy(&user->lUcont_rm1); CHKERRQ(ierr); }
2987 if (user->lNvert_o) { ierr = VecDestroy(&user->lNvert_o); CHKERRQ(ierr); }
2988
2989 // --- Group D: Grid Metrics - Face Centered (All levels) ---
2990 if (user->Csi) { ierr = VecDestroy(&user->Csi); CHKERRQ(ierr); }
2991 if (user->Eta) { ierr = VecDestroy(&user->Eta); CHKERRQ(ierr); }
2992 if (user->Zet) { ierr = VecDestroy(&user->Zet); CHKERRQ(ierr); }
2993 if (user->Aj) { ierr = VecDestroy(&user->Aj); CHKERRQ(ierr); }
2994 if (user->lCsi) { ierr = VecDestroy(&user->lCsi); CHKERRQ(ierr); }
2995 if (user->lEta) { ierr = VecDestroy(&user->lEta); CHKERRQ(ierr); }
2996 if (user->lZet) { ierr = VecDestroy(&user->lZet); CHKERRQ(ierr); }
2997 if (user->lAj) { ierr = VecDestroy(&user->lAj); CHKERRQ(ierr); }
2998
2999 // --- Group E: Grid Metrics - Face Centered (All levels) ---
3000 if (user->ICsi) { ierr = VecDestroy(&user->ICsi); CHKERRQ(ierr); }
3001 if (user->IEta) { ierr = VecDestroy(&user->IEta); CHKERRQ(ierr); }
3002 if (user->IZet) { ierr = VecDestroy(&user->IZet); CHKERRQ(ierr); }
3003 if (user->JCsi) { ierr = VecDestroy(&user->JCsi); CHKERRQ(ierr); }
3004 if (user->JEta) { ierr = VecDestroy(&user->JEta); CHKERRQ(ierr); }
3005 if (user->JZet) { ierr = VecDestroy(&user->JZet); CHKERRQ(ierr); }
3006 if (user->KCsi) { ierr = VecDestroy(&user->KCsi); CHKERRQ(ierr); }
3007 if (user->KEta) { ierr = VecDestroy(&user->KEta); CHKERRQ(ierr); }
3008 if (user->KZet) { ierr = VecDestroy(&user->KZet); CHKERRQ(ierr); }
3009 if (user->IAj) { ierr = VecDestroy(&user->IAj); CHKERRQ(ierr); }
3010 if (user->JAj) { ierr = VecDestroy(&user->JAj); CHKERRQ(ierr); }
3011 if (user->KAj) { ierr = VecDestroy(&user->KAj); CHKERRQ(ierr); }
3012 if (user->lICsi) { ierr = VecDestroy(&user->lICsi); CHKERRQ(ierr); }
3013 if (user->lIEta) { ierr = VecDestroy(&user->lIEta); CHKERRQ(ierr); }
3014 if (user->lIZet) { ierr = VecDestroy(&user->lIZet); CHKERRQ(ierr); }
3015 if (user->lJCsi) { ierr = VecDestroy(&user->lJCsi); CHKERRQ(ierr); }
3016 if (user->lJEta) { ierr = VecDestroy(&user->lJEta); CHKERRQ(ierr); }
3017 if (user->lJZet) { ierr = VecDestroy(&user->lJZet); CHKERRQ(ierr); }
3018 if (user->lKCsi) { ierr = VecDestroy(&user->lKCsi); CHKERRQ(ierr); }
3019 if (user->lKEta) { ierr = VecDestroy(&user->lKEta); CHKERRQ(ierr); }
3020 if (user->lKZet) { ierr = VecDestroy(&user->lKZet); CHKERRQ(ierr); }
3021 if (user->lIAj) { ierr = VecDestroy(&user->lIAj); CHKERRQ(ierr); }
3022 if (user->lJAj) { ierr = VecDestroy(&user->lJAj); CHKERRQ(ierr); }
3023 if (user->lKAj) { ierr = VecDestroy(&user->lKAj); CHKERRQ(ierr); }
3024
3025 // --- Group F: Cell/Face Coordinates and Grid Spacing (All levels) ---
3026 if (user->Cent) { ierr = VecDestroy(&user->Cent); CHKERRQ(ierr); }
3027 if (user->lCent) { ierr = VecDestroy(&user->lCent); CHKERRQ(ierr); }
3028 if (user->GridSpace) { ierr = VecDestroy(&user->GridSpace); CHKERRQ(ierr); }
3029 if (user->lGridSpace) { ierr = VecDestroy(&user->lGridSpace); CHKERRQ(ierr); }
3030 if (user->Centx) { ierr = VecDestroy(&user->Centx); CHKERRQ(ierr); }
3031 if (user->Centy) { ierr = VecDestroy(&user->Centy); CHKERRQ(ierr); }
3032 if (user->Centz) { ierr = VecDestroy(&user->Centz); CHKERRQ(ierr); }
3033
3034 // --- Group G: Turbulence Model Vectors (Finest level, conditional on les/rans) ---
3035 if (user->Nu_t) { ierr = VecDestroy(&user->Nu_t); CHKERRQ(ierr); }
3036 if (user->lNu_t) { ierr = VecDestroy(&user->lNu_t); CHKERRQ(ierr); }
3037 if (user->CS) { ierr = VecDestroy(&user->CS); CHKERRQ(ierr); }
3038 if (user->lCs) { ierr = VecDestroy(&user->lCs); CHKERRQ(ierr); }
3039 if (user->lFriction_Velocity) { ierr = VecDestroy(&user->lFriction_Velocity); CHKERRQ(ierr); }
3040 if (user->K_Omega) { ierr = VecDestroy(&user->K_Omega); CHKERRQ(ierr); }
3041 if (user->lK_Omega) { ierr = VecDestroy(&user->lK_Omega); CHKERRQ(ierr); }
3042 if (user->K_Omega_o) { ierr = VecDestroy(&user->K_Omega_o); CHKERRQ(ierr); }
3043 if (user->lK_Omega_o) { ierr = VecDestroy(&user->lK_Omega_o); CHKERRQ(ierr); }
3044
3045 // --- Group H: Particle Vectors (Finest level, conditional on np > 0) ---
3046 if (user->ParticleCount) { ierr = VecDestroy(&user->ParticleCount); CHKERRQ(ierr); }
3047 if (user->lParticleCount) { ierr = VecDestroy(&user->lParticleCount); CHKERRQ(ierr); }
3048 if (user->Psi) { ierr = VecDestroy(&user->Psi); CHKERRQ(ierr); }
3049 if (user->lPsi) { ierr = VecDestroy(&user->lPsi); CHKERRQ(ierr); }
3050
3051 // --- Group I: Boundary Condition Vectors (All levels) ---
3052 if (user->Bcs.Ubcs) { ierr = VecDestroy(&user->Bcs.Ubcs); CHKERRQ(ierr); }
3053 if (user->Bcs.Uch) { ierr = VecDestroy(&user->Bcs.Uch); CHKERRQ(ierr); }
3054
3055 // --- Group J: Post-Processing Vectors (Finest level, postprocessor mode) ---
3056 if (user->P_nodal) { ierr = VecDestroy(&user->P_nodal); CHKERRQ(ierr); }
3057 if (user->Ucat_nodal) { ierr = VecDestroy(&user->Ucat_nodal); CHKERRQ(ierr); }
3058 if (user->Qcrit) { ierr = VecDestroy(&user->Qcrit); CHKERRQ(ierr); }
3059 if (user->Psi_nodal) { ierr = VecDestroy(&user->Psi_nodal); CHKERRQ(ierr); }
3060
3061 // --- Group K: Interpolation Vectors (Lazy allocation) ---
3062 if (user->CellFieldAtCorner) { ierr = VecDestroy(&user->CellFieldAtCorner); CHKERRQ(ierr); }
3063 if (user->lCellFieldAtCorner) { ierr = VecDestroy(&user->lCellFieldAtCorner); CHKERRQ(ierr); }
3064
3065 // --- Group L: Statistical Averaging Vectors (If allocated) ---
3066 if (user->Ucat_sum) { ierr = VecDestroy(&user->Ucat_sum); CHKERRQ(ierr); }
3067 if (user->Ucat_cross_sum) { ierr = VecDestroy(&user->Ucat_cross_sum); CHKERRQ(ierr); }
3068 if (user->Ucat_square_sum) { ierr = VecDestroy(&user->Ucat_square_sum); CHKERRQ(ierr); }
3069 if (user->P_sum) { ierr = VecDestroy(&user->P_sum); CHKERRQ(ierr); }
3070
3071 // --- Group M: Implicit Solver Temporary Vectors (Destroyed after use, but check anyway) ---
3072 if (user->Rhs) { ierr = VecDestroy(&user->Rhs); CHKERRQ(ierr); }
3073 if (user->dUcont) { ierr = VecDestroy(&user->dUcont); CHKERRQ(ierr); }
3074 if (user->pUcont) { ierr = VecDestroy(&user->pUcont); CHKERRQ(ierr); }
3075
3076 // --- Group N: Poisson Solver Vectors (Destroyed after solve, but check anyway) ---
3077 if (user->B) { ierr = VecDestroy(&user->B); CHKERRQ(ierr); }
3078 if (user->R) { ierr = VecDestroy(&user->R); CHKERRQ(ierr); }
3079
3080 LOG_ALLOW(LOCAL, LOG_DEBUG, "All vectors destroyed for UserCtx.\n");
3081 PetscFunctionReturn(0);
3082}
Vec Rhs
Definition variables.h:845
Vec K_Omega_o
Definition variables.h:865
Vec K_Omega
Definition variables.h:865
Vec lCellFieldAtCorner
Definition variables.h:846
Vec Ucat_square_sum
Definition variables.h:868
Vec P_sum
Definition variables.h:868
Vec pUcont
Definition variables.h:845
Vec lK_Omega_o
Definition variables.h:865
Vec lK_Omega
Definition variables.h:865
Vec CellFieldAtCorner
Definition variables.h:846
Vec Ucat_sum
Definition variables.h:868
Vec dUcont
Definition variables.h:845
Vec Ucat_cross_sum
Definition variables.h:868
Here is the caller graph for this function:

◆ DestroyUserContext()

PetscErrorCode DestroyUserContext ( UserCtx user)

Internal helper implementation: DestroyUserContext().

Destroys all resources allocated within a single UserCtx structure.

Local to this translation unit.

Definition at line 3089 of file setup.c.

3090{
3091 PetscErrorCode ierr;
3092 PetscFunctionBeginUser;
3093
3094 if (!user) {
3095 LOG_ALLOW(LOCAL, LOG_WARNING, "DestroyUserContext called with NULL user pointer.\n");
3096 PetscFunctionReturn(0);
3097 }
3098
3099 LOG_ALLOW(LOCAL, LOG_INFO, "Destroying UserCtx at level %d...\n", user->thislevel);
3100
3101 // --- Step 1: Destroy Boundary Condition System ---
3102 // This handles all BC handlers and their private data.
3103 ierr = BoundarySystem_Destroy(user); CHKERRQ(ierr);
3104 LOG_ALLOW(LOCAL, LOG_DEBUG, " Boundary system destroyed.\n");
3105
3106 // --- Step 2: Destroy All Vectors ---
3107 // Handles ~74 Vec objects with proper NULL checking.
3108 ierr = DestroyUserVectors(user); CHKERRQ(ierr);
3109 LOG_ALLOW(LOCAL, LOG_DEBUG, " All vectors destroyed.\n");
3110
3111 // --- Step 3: Destroy Matrix and Solver Objects ---
3112 // Destroy pressure-Poisson matrices and solver.
3113 if (user->A) {
3114 ierr = MatDestroy(&user->A); CHKERRQ(ierr);
3115 LOG_ALLOW(LOCAL, LOG_DEBUG, " Matrix A destroyed.\n");
3116 }
3117 if (user->C) {
3118 ierr = MatDestroy(&user->C); CHKERRQ(ierr);
3119 LOG_ALLOW(LOCAL, LOG_DEBUG, " Matrix C destroyed.\n");
3120 }
3121 if (user->MR) {
3122 ierr = MatDestroy(&user->MR); CHKERRQ(ierr);
3123 LOG_ALLOW(LOCAL, LOG_DEBUG, " Matrix MR destroyed.\n");
3124 }
3125 if (user->MP) {
3126 ierr = MatDestroy(&user->MP); CHKERRQ(ierr);
3127 LOG_ALLOW(LOCAL, LOG_DEBUG, " Matrix MP destroyed.\n");
3128 }
3129 if (user->ksp) {
3130 ierr = KSPDestroy(&user->ksp); CHKERRQ(ierr);
3131 LOG_ALLOW(LOCAL, LOG_DEBUG, " KSP solver destroyed.\n");
3132 }
3133 if (user->nullsp) {
3134 ierr = MatNullSpaceDestroy(&user->nullsp); CHKERRQ(ierr);
3135 LOG_ALLOW(LOCAL, LOG_DEBUG, " MatNullSpace destroyed.\n");
3136 }
3137
3138 // --- Step 4: Destroy Application Ordering ---
3139 if (user->ao) {
3140 ierr = AODestroy(&user->ao); CHKERRQ(ierr);
3141 LOG_ALLOW(LOCAL, LOG_DEBUG, " AO destroyed.\n");
3142 }
3143
3144 // --- Step 5: Destroy DM Objects ---
3145 // Destroy in reverse order of dependency: post_swarm, swarm, fda2, fda, da
3146 if (user->post_swarm) {
3147 ierr = DMDestroy(&user->post_swarm); CHKERRQ(ierr);
3148 LOG_ALLOW(LOCAL, LOG_DEBUG, " post_swarm DM destroyed.\n");
3149 }
3150 if (user->swarm) {
3151 ierr = DMDestroy(&user->swarm); CHKERRQ(ierr);
3152 LOG_ALLOW(LOCAL, LOG_DEBUG, " swarm DM destroyed.\n");
3153 }
3154 if (user->fda2) {
3155 ierr = DMDestroy(&user->fda2); CHKERRQ(ierr);
3156 LOG_ALLOW(LOCAL, LOG_DEBUG, " fda2 DM destroyed.\n");
3157 }
3158 if (user->da) {
3159 ierr = DMDestroy(&user->da); CHKERRQ(ierr);
3160 LOG_ALLOW(LOCAL, LOG_DEBUG, " da DM destroyed.\n");
3161 }
3162
3163 // --- Step 6: Free PetscMalloc'd Arrays ---
3164 // Free arrays allocated with PetscMalloc1
3165 if (user->RankCellInfoMap) {
3166 ierr = PetscFree(user->RankCellInfoMap); CHKERRQ(ierr);
3167 user->RankCellInfoMap = NULL;
3168 LOG_ALLOW(LOCAL, LOG_DEBUG, " RankCellInfoMap freed.\n");
3169 }
3170 if (user->KSKE) {
3171 ierr = PetscFree(user->KSKE); CHKERRQ(ierr);
3172 user->KSKE = NULL;
3173 LOG_ALLOW(LOCAL, LOG_DEBUG, " KSKE array freed.\n");
3174 }
3175
3176 LOG_ALLOW(LOCAL, LOG_INFO, "UserCtx at level %d fully destroyed.\n", user->thislevel);
3177 PetscFunctionReturn(0);
3178}
PetscErrorCode BoundarySystem_Destroy(UserCtx *user)
Cleans up and destroys all boundary system resources.
PetscErrorCode DestroyUserVectors(UserCtx *user)
Internal helper implementation: DestroyUserVectors().
Definition setup.c:2954
MatNullSpace nullsp
Definition variables.h:849
PetscInt * KSKE
Definition variables.h:850
DM post_swarm
Definition variables.h:886
KSP ksp
Definition variables.h:849
Here is the call graph for this function:
Here is the caller graph for this function:

◆ FinalizeSimulation()

PetscErrorCode FinalizeSimulation ( SimCtx simCtx)

Implementation of FinalizeSimulation().

Main cleanup function for the entire simulation context.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/setup.h.

See also
FinalizeSimulation()

Definition at line 3188 of file setup.c.

3189{
3190 PetscErrorCode ierr;
3191 PetscFunctionBeginUser;
3192
3193 if (!simCtx) {
3194 LOG_ALLOW(GLOBAL, LOG_WARNING, "FinalizeSimulation called with NULL SimCtx pointer.\n");
3195 PetscFunctionReturn(0);
3196 }
3197
3198 LOG_ALLOW(GLOBAL, LOG_INFO, "========================================\n");
3199 LOG_ALLOW(GLOBAL, LOG_INFO, "Beginning simulation memory cleanup...\n");
3200 LOG_ALLOW(GLOBAL, LOG_INFO, "========================================\n");
3201
3202 // ============================================================================
3203 // PHASE 1: DESTROY MULTIGRID HIERARCHY (All UserCtx structures)
3204 // ============================================================================
3205
3206 if (simCtx->usermg.mgctx) {
3207 LOG_ALLOW(GLOBAL, LOG_INFO, "Destroying multigrid hierarchy (%d levels)...\n",
3208 simCtx->usermg.mglevels);
3209
3210 // Destroy each UserCtx from finest to coarsest (reverse order is safer)
3211 for (PetscInt level = simCtx->usermg.mglevels - 1; level >= 0; level--) {
3212 UserCtx *user = simCtx->usermg.mgctx[level].user;
3213 if (user) {
3214 LOG_ALLOW(LOCAL, LOG_INFO, " Destroying level %d of %d...\n",
3215 level, simCtx->usermg.mglevels - 1);
3216 ierr = DestroyUserContext(user); CHKERRQ(ierr);
3217
3218 // Free the UserCtx structure itself
3219 ierr = PetscFree(user); CHKERRQ(ierr);
3220 simCtx->usermg.mgctx[level].user = NULL;
3221 }
3222
3223 // Destroy the MGCtx-level packer DM
3224 if (simCtx->usermg.mgctx[level].packer) {
3225 ierr = DMDestroy(&simCtx->usermg.mgctx[level].packer); CHKERRQ(ierr);
3226 LOG_ALLOW(LOCAL, LOG_DEBUG, " MGCtx[%d].packer destroyed.\n", level);
3227 }
3228 }
3229
3230 // Free the MGCtx array itself
3231 ierr = PetscFree(simCtx->usermg.mgctx); CHKERRQ(ierr);
3232 simCtx->usermg.mgctx = NULL;
3233 LOG_ALLOW(GLOBAL, LOG_INFO, "All multigrid levels destroyed.\n");
3234 }
3235
3236 // ============================================================================
3237 // PHASE 2: DESTROY USERMG-LEVEL OBJECTS
3238 // ============================================================================
3239
3240 if (simCtx->usermg.packer) {
3241 ierr = DMDestroy(&simCtx->usermg.packer); CHKERRQ(ierr);
3242 LOG_ALLOW(LOCAL, LOG_DEBUG, "UserMG.packer DM destroyed.\n");
3243 }
3244
3245 if (simCtx->usermg.snespacker) {
3246 ierr = SNESDestroy(&simCtx->usermg.snespacker); CHKERRQ(ierr);
3247 LOG_ALLOW(LOCAL, LOG_DEBUG, "UserMG.snespacker SNES destroyed.\n");
3248 }
3249
3250 // ============================================================================
3251 // PHASE 3: DESTROY SIMCTX-LEVEL OBJECTS
3252 // ============================================================================
3253
3254 LOG_ALLOW(GLOBAL, LOG_INFO, "Destroying SimCtx-level objects...\n");
3255
3256 // --- PetscViewer for logging ---
3257 if (simCtx->logviewer) {
3258 ierr = PetscViewerDestroy(&simCtx->logviewer); CHKERRQ(ierr);
3259 LOG_ALLOW(LOCAL, LOG_DEBUG, " logviewer destroyed.\n");
3260 }
3261
3262 // --- Particle System DM ---
3263 if (simCtx->dm_swarm) {
3264 ierr = DMDestroy(&simCtx->dm_swarm); CHKERRQ(ierr);
3265 LOG_ALLOW(LOCAL, LOG_DEBUG, " dm_swarm destroyed.\n");
3266 }
3267
3268 // --- BoundingBox List (Array of BoundingBox structs) ---
3269 if (simCtx->bboxlist) {
3270 ierr = PetscFree(simCtx->bboxlist); CHKERRQ(ierr);
3271 simCtx->bboxlist = NULL;
3272 LOG_ALLOW(LOCAL, LOG_DEBUG, " bboxlist freed.\n");
3273 }
3274
3275 // --- Boundary Condition Files (Array of strings) ---
3276 if (simCtx->bcs_files) {
3277 for (PetscInt i = 0; i < simCtx->num_bcs_files; i++) {
3278 if (simCtx->bcs_files[i]) {
3279 ierr = PetscFree(simCtx->bcs_files[i]); CHKERRQ(ierr);
3280 }
3281 }
3282 ierr = PetscFree(simCtx->bcs_files); CHKERRQ(ierr);
3283 simCtx->bcs_files = NULL;
3284 LOG_ALLOW(LOCAL, LOG_DEBUG, " bcs_files array freed (%d files).\n", simCtx->num_bcs_files);
3285 }
3286
3287 // --- Brownian Motion RNG ---
3288 if (simCtx->BrownianMotionRNG) {
3289 ierr = PetscRandomDestroy(&simCtx->BrownianMotionRNG); CHKERRQ(ierr);
3290 LOG_ALLOW(LOCAL, LOG_DEBUG, " BrownianMotionRNG destroyed.\n");
3291 }
3292 // --- Post-Processing Parameters ---
3293 // pps is allocated with PetscNew and contains only static char arrays and basic types.
3294 // No internal dynamic allocations need to be freed.
3295 if (simCtx->pps) {
3296 ierr = PetscFree(simCtx->pps); CHKERRQ(ierr);
3297 simCtx->pps = NULL;
3298 LOG_ALLOW(LOCAL, LOG_DEBUG, " PostProcessParams freed.\n");
3299 }
3300
3301 // --- IBM/FSI Objects ---
3302 // Note: These are initialized to NULL and currently have no dedicated destroy functions.
3303 // If these modules are extended with cleanup routines, call them here.
3304 if (simCtx->ibm != NULL) {
3305 LOG_ALLOW(GLOBAL, LOG_WARNING, " WARNING: simCtx->ibm is non-NULL but no destroy function exists. Potential memory leak.\n");
3306 }
3307 if (simCtx->ibmv != NULL) {
3308 LOG_ALLOW(GLOBAL, LOG_WARNING, " WARNING: simCtx->ibmv is non-NULL but no destroy function exists. Potential memory leak.\n");
3309 }
3310 if (simCtx->fsi != NULL) {
3311 LOG_ALLOW(GLOBAL, LOG_WARNING, " WARNING: simCtx->fsi is non-NULL but no destroy function exists. Potential memory leak.\n");
3312 }
3313
3314 // --- Logging Allowed Functions (Array of strings) ---
3315 // Note: The logging system maintains its own copy via set_allowed_functions(),
3316 // so freeing simCtx->allowedFuncs will NOT affect LOG_ALLOW functionality.
3317 if (simCtx->allowedFuncs) {
3318 for (PetscInt i = 0; i < simCtx->nAllowed; i++) {
3319 if (simCtx->allowedFuncs[i]) {
3320 ierr = PetscFree(simCtx->allowedFuncs[i]); CHKERRQ(ierr);
3321 }
3322 }
3323 ierr = PetscFree(simCtx->allowedFuncs); CHKERRQ(ierr);
3324 simCtx->allowedFuncs = NULL;
3325 LOG_ALLOW(LOCAL, LOG_DEBUG, " allowedFuncs array freed (%d functions).\n", simCtx->nAllowed);
3326 }
3327
3328 // --- Profiling Critical Functions (Array of strings) ---
3329 if (simCtx->profilingSelectedFuncs) {
3330 for (PetscInt i = 0; i < simCtx->nProfilingSelectedFuncs; i++) {
3331 if (simCtx->profilingSelectedFuncs[i]) {
3332 ierr = PetscFree(simCtx->profilingSelectedFuncs[i]); CHKERRQ(ierr);
3333 }
3334 }
3335 ierr = PetscFree(simCtx->profilingSelectedFuncs); CHKERRQ(ierr);
3336 simCtx->profilingSelectedFuncs = NULL;
3337 LOG_ALLOW(LOCAL, LOG_DEBUG, " profilingSelectedFuncs array freed (%d functions).\n", simCtx->nProfilingSelectedFuncs);
3338 }
3339
3340 // ============================================================================
3341 // PHASE 4: FINAL SUMMARY
3342 // ============================================================================
3343
3344 LOG_ALLOW(GLOBAL, LOG_INFO, "========================================\n");
3345 LOG_ALLOW(GLOBAL, LOG_INFO, "Simulation cleanup completed successfully.\n");
3346 LOG_ALLOW(GLOBAL, LOG_INFO, "All PETSc objects have been destroyed.\n");
3347 LOG_ALLOW(GLOBAL, LOG_INFO, "========================================\n");
3348
3349 PetscFunctionReturn(0);
3350}
PetscErrorCode DestroyUserContext(UserCtx *user)
Internal helper implementation: DestroyUserContext().
Definition setup.c:3089
DM packer
Definition variables.h:539
SNES snespacker
Definition variables.h:540
DM packer
Definition variables.h:530
Here is the call graph for this function:
Here is the caller graph for this function: