setup.h File Reference

#include <petscpf.h>
#include <petscdmswarm.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <petsctime.h>
#include <petscsys.h>
#include <petscdmcomposite.h>
#include <petscsystypes.h>
#include "variables.h"
#include "field_catalog.h"
#include "ParticleSwarm.h"
#include "walkingsearch.h"
#include "grid.h"
#include "logging.h"
#include "io.h"
#include "interpolation.h"
#include "ParticleMotion.h"
#include "Boundaries.h"

Include dependency graph for setup.h:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Macros
#define	Allocate3DArray(array, nz, ny, nx)
#define	Deallocate3DArray(array, nz, ny)

Functions
int	PicurvHandleVersionArgument (int argc, char **argv, const char *executable_name)
	Print the shared native build identity when a version flag is present.
PetscErrorCode	CreateSimulationContext (int argc, char **argv, SimCtx **p_simCtx)
	Allocates and populates the master SimulationContext object.
PetscBool	RuntimeWalltimeGuardParsePositiveSeconds (const char *text, PetscReal *seconds_out)
	Parse a positive floating-point seconds value from runtime metadata.
PetscErrorCode	SetupSimulationEnvironment (SimCtx *simCtx)
	Verifies and prepares the complete I/O environment for a simulation run.
PetscErrorCode	SetupGridAndSolvers (SimCtx *simCtx)
	The main orchestrator for setting up all grid-related components.
PetscErrorCode	CreateAndInitializeAllVectors (SimCtx *simCtx)
	Creates and initializes all PETSc Vec objects for all fields.
PetscErrorCode	UpdateLocalGhosts (UserCtx *user, FieldId field_id)
	Updates the local vector (including ghost points) from its corresponding global vector.
PetscErrorCode	SetupBoundaryConditions (SimCtx *simCtx)
	(Orchestrator) Sets up all boundary conditions for the simulation.
PetscErrorCode	InitializeSolutionConvergenceState (SimCtx *simCtx)
	Allocates any runtime storage required by solution-convergence logging.
PetscErrorCode	DestroySolutionConvergenceState (SimCtx *simCtx)
	Frees any runtime storage allocated for solution-convergence logging.
PetscErrorCode	Allocate3DArrayScalar (PetscReal ****array, PetscInt nz, PetscInt ny, PetscInt nx)
	Allocates a 3D array of PetscReal values using PetscCalloc.
PetscErrorCode	Deallocate3DArrayScalar (PetscReal ***array, PetscInt nz, PetscInt ny)
	Deallocates a 3D array of PetscReal values allocated by Allocate3DArrayScalar.
PetscErrorCode	Allocate3DArrayVector (Cmpnts ****array, PetscInt nz, PetscInt ny, PetscInt nx)
	Allocates a contiguous 3D array of Cmpnts values.
PetscErrorCode	Deallocate3DArrayVector (Cmpnts ***array, PetscInt nz, PetscInt ny)
	Deallocates a 3D array of Cmpnts structures allocated by Allocate3DArrayVector.
PetscErrorCode	GetOwnedCellRange (const DMDALocalInfo *info_nodes, PetscInt dim, PetscInt *xs_cell_global_out, PetscInt *xm_cell_local_out)
	Determines the global starting index and number of CELLS owned by the current processor in a specified dimension.
PetscErrorCode	ComputeAndStoreNeighborRanks (UserCtx *user)
	Computes and stores the Cartesian neighbor ranks for the DMDA decomposition.
PetscErrorCode	SetDMDAProcLayout (DM dm, UserCtx *user)
	Sets the processor layout for a given DMDA based on PETSc options.
PetscErrorCode	SetupDomainRankInfo (SimCtx *simCtx)
	Sets up the full rank communication infrastructure, including neighbor ranks and bounding box exchange.
PetscErrorCode	Contra2Cart (UserCtx *user)
	Reconstructs Cartesian velocity (Ucat) at cell centers from contravariant velocity (Ucont) defined on cell faces.
PetscErrorCode	Cart2Contra (UserCtx *user)
	Convert the ghosted Cartesian velocity field to contravariant face fluxes.
PetscErrorCode	UniformCart2Contra (UserCtx *user, PetscReal u, PetscReal v, PetscReal w)
	Populate contravariant fluxes from one uniform Cartesian velocity.
PetscErrorCode	SetupDomainCellDecompositionMap (UserCtx *user)
	Creates and distributes a map of the domain's cell decomposition to all ranks.
PetscErrorCode	BinarySearchInt64 (PetscInt n, const PetscInt64 arr[], PetscInt64 key, PetscBool *found)
	Performs a binary search for a key in a sorted array of PetscInt64.
PetscErrorCode	ComputeDivergence (UserCtx *user)
	Computes the discrete divergence of the contravariant velocity field.
PetscErrorCode	InitializeRandomGenerators (UserCtx *user, PetscRandom *randx, PetscRandom *randy, PetscRandom *randz)
	Initializes random number generators for assigning particle properties.
PetscErrorCode	InitializeLogicalSpaceRNGs (PetscInt base_seed, PetscRandom *rand_logic_i, PetscRandom *rand_logic_j, PetscRandom *rand_logic_k)
	Initializes random number generators for logical space operations [0.0, 1.0).
PetscErrorCode	InitializeBrownianRNG (SimCtx *simCtx)
	Initializes a single master RNG for time-stepping physics (Brownian motion).
void	TransformScalarDerivativesToPhysical (PetscReal jacobian, Cmpnts csi_metrics, Cmpnts eta_metrics, Cmpnts zet_metrics, PetscReal dPhi_dcsi, PetscReal dPhi_deta, PetscReal dPhi_dzet, Cmpnts *gradPhi)
	Transforms scalar derivatives from computational space to physical space.
PetscErrorCode	ComputeScalarFieldDerivatives (UserCtx *user, PetscInt i, PetscInt j, PetscInt k, PetscReal ***field_data, Cmpnts *grad)
	Computes the gradient of a cell-centered SCALAR field at a specific grid point.
PetscErrorCode	ComputeVectorFieldDerivatives (UserCtx *user, PetscInt i, PetscInt j, PetscInt k, Cmpnts ***field_data, Cmpnts *dudx, Cmpnts *dvdx, Cmpnts *dwdx)
	Computes the derivatives of a cell-centered vector field at a specific grid point.
PetscErrorCode	LESConfigSetDefaults (LESConfig *config)
	Installs the shipped defaults for every LES closure parameter.
PetscErrorCode	DestroyUserVectors (UserCtx *user)
	Destroys all PETSc Vec objects within a single UserCtx structure.
PetscErrorCode	DestroyUserContext (UserCtx *user)
	Destroys all resources allocated within a single UserCtx structure.
PetscErrorCode	FinalizeSimulation (SimCtx *simCtx)
	Main cleanup function for the entire simulation context.

Macro Definition Documentation

Allocate3DArray

#define Allocate3DArray	(		array,
			nz,
			ny,
			nx
	)

Value:

  _Generic((array),                          \
    PetscReal ****: Allocate3DArrayScalar,       \
    Cmpnts    ****: Allocate3DArrayVector            \
  )(array, nz, ny, nx)

Definition at line 41 of file setup.h.

                  : Allocate3DArrayScalar,       \
    Cmpnts    ****: Allocate3DArrayVector            \
  )(array, nz, ny, nx)

Deallocate3DArray

#define Deallocate3DArray	(		array,
			nz,
			ny
	)

Value:

  _Generic((array),                        \
    PetscReal ***: Deallocate3DArrayScalar,     \
    Cmpnts    ***: Deallocate3DArrayVector          \
  )(array, nz, ny)

Definition at line 48 of file setup.h.

                 : Deallocate3DArrayScalar,     \
    Cmpnts    ***: Deallocate3DArrayVector          \
  )(array, nz, ny)

Function Documentation

PicurvHandleVersionArgument()

int PicurvHandleVersionArgument	(	int	argc,
		char **	argv,
		const char *	executable_name
	)

Print the shared native build identity when a version flag is present.

Prints two lines: <name> <release>+g<commit>[.dirty], then the PETSc this binary was compiled against, petsc <version> <debug|optimized> <arch> <dir>. The second line is also kept in the binary behind a PICURV_PETSC_BUILD: marker, so the conductor can read it from a binary that cannot start.

Parameters

[in]	argc	Native executable argument count.
[in]	argv	Native executable argument vector.
[in]	executable_name	User-facing executable name printed with the identity.

Returns

Nonzero when a version flag was handled and execution should stop.

Print the shared native build identity when a version flag is present.

The public argument and return contract is documented in include/setup.h.

See also

PicurvHandleVersionArgument()

Definition at line 47 of file setup.c.

{
    for (int index = 1; index < argc; ++index) {
        if (!strcmp(argv[index], "--version") || !strcmp(argv[index], "-version")) {
            printf("%s %s+g%.12s%s\n", executable_name, PICURV_RELEASE_VERSION,
                   PICURV_GIT_COMMIT,
                   !strcmp(PICURV_BUILD_DIRTY, "true") ? ".dirty" : "");
            printf("%s\n", picurv_petsc_build_stamp + sizeof(PICURV_PETSC_STAMP_MARKER) - 1);
            return 1;
        }
    }
    return 0;
}

Here is the caller graph for this function:

CreateSimulationContext()

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

Allocates and populates the master SimulationContext object.

This function serves as the single, authoritative entry point for all simulation configuration. It merges the setup logic from both the legacy FSI/IBM solver and the modern particle solver into a unified, robust process.

The function follows a strict sequence:

1. Allocate Context & Set Defaults: It first allocates the SimulationContext and populates every field with a sane, hardcoded default value. This ensures the simulation starts from a known, predictable state.
2. Configure Logging System: It configures the custom logging framework. It parses the -func_config_file option to load a list of function names allowed to produce log output. This configuration (the file path and the list of function names) is stored within the SimulationContext for later reference and cleanup.
3. Parse All Options: It performs a comprehensive pass of PetscOptionsGet... calls for every possible command-line flag, overriding the default values set in step 1.
4. Log Summary: After all options are parsed, it uses the now-active logging system to print a summary of the key simulation parameters.

Parameters

[in]	argc	Argument count passed from main.
[in]	argv	Argument vector passed from main.
[out]	p_simCtx	On success, this will point to the newly created and fully configured SimulationContext pointer. The caller is responsible for eventually destroying this object by calling FinalizeSimulation().

Returns

PetscErrorCode Returns 0 on success, or a non-zero PETSc error code on failure.

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 371 of file setup.c.

{
    PetscErrorCode ierr;
    (void)argc;
    (void)argv;
    SimCtx *simCtx;
    char control_filename[PETSC_MAX_PATH_LEN] = ""; // Temporary placeholder for control file name.
    PetscBool control_flg; // Temporary placeholder for control file tag existence check flag.
    PetscBool particle_console_output_freq_flg = PETSC_FALSE;
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    // === 1. Allocate the Context Struct and Set ALL Defaults ==================
    ierr = PetscNew(p_simCtx); CHKERRQ(ierr);
    simCtx = *p_simCtx;
 
    // --- Group 1: Parallelism & MPI ---
    simCtx->rank = 0; simCtx->size = 1;
 
    // --- Group 2: Simulation Control, Time, and I/O ---
    simCtx->step = 0; simCtx->ti = 0.0; simCtx->StartStep = 0; simCtx->StepsToRun = 10;
    simCtx->tiout = 10; simCtx->particleConsoleOutputFreq = simCtx->tiout;
    simCtx->StartTime = 0.0; simCtx->dt = 0.001;
    simCtx->OnlySetup = PETSC_FALSE;
    simCtx->continueMode = PETSC_FALSE;
    simCtx->logviewer = NULL;
    strcpy(simCtx->eulerianSource,"solve");
    strcpy(simCtx->restart_dir,"restart");
    strcpy(simCtx->output_dir,"output");
    strcpy(simCtx->log_dir,"logs");
    strcpy(simCtx->analysis_dir,"output/analysis/metrics");
    simCtx->_io_context_buffer[0] = '\0';
    simCtx->current_io_directory = NULL;
    simCtx->checkpointGeometrySHA256[0] = '\0';
    simCtx->checkpointGeometryHashReady = PETSC_FALSE;
    simCtx->restartHistoryAvailable = PETSC_FALSE;
 
    // --- Group 3: High-Level Physics & Model Selection Flags ---
    simCtx->immersed = 0; simCtx->movefsi = 0; simCtx->rotatefsi = 0;
    simCtx->invicid = 0;
    simCtx->TwoD = 0; simCtx->moveframe = 0; simCtx->rotateframe = 0;
    strcpy(simCtx->AnalyticalSolutionType,"TGV3D");
 
    // --- Group 4: Immersed-body flux corrections (refused; see Group 3 parsing) ---
    simCtx->MHV=0; simCtx->LV=0;
 
    // --- Group 5: Solver & Numerics Parameters ---
    simCtx->mom_solver_type = MOMENTUM_SOLVER_DUALTIME_PICARD_JAMESON_RK; simCtx->mom_max_pseudo_steps = 50;
    simCtx->mom_dt_jameson_residual_norm_noise_allowance_factor = 1.1; // raised from 1.05; less aggressive rejection
    simCtx->mom_atol = 1e-7; simCtx->mom_rtol = 1e-4;
    /* Residual-based convergence is the default. Leaving both at 0.0 selected the
     * update-only branch, where |dU| <= mom_atol can pass purely because dtau collapsed
     * (|dU| ~ dtau*|R|), converging on a state that does not satisfy the equations.
     * Shipped configs all set these explicitly; these defaults protect a minimal
     * user-written config. Set both non-positive to opt back out deliberately. */
    simCtx->mom_resid_atol = 1e-8; simCtx->mom_resid_rtol = 1e-3;
    simCtx->imp_stol = 1.e-8;
    simCtx->mglevels = 3;
    simCtx->mg_preItr = 1; simCtx->mg_poItr = 1;
    simCtx->poisson = 0;
    simCtx->central=0;
    /* pseudo_cfl and its bounds are now dimensionless Courant numbers: CFL = dtau * lambda_max,
       where lambda_max is the global spectral radius computed at each physical timestep.
       Stable range for 4-stage Jameson RK: ~0–2.83. Initial 0.5 gives a comfortable margin. */
    simCtx->ren = 100.0; simCtx->pseudo_cfl = 0.5;
    simCtx->max_pseudo_cfl = 2.0; simCtx->min_pseudo_cfl = 0.001;
    simCtx->pseudo_cfl_reduction_factor = 0.75;
    simCtx->pseudo_cfl_growth_factor = 1.1; // raised from 1.0; controller can now increase CFL
    simCtx->no_pseudo_cfl_backtrack = PETSC_FALSE;
    simCtx->mom_ratio_ema_alpha = 0.3;   /* moderate smoothing; set to 1.0 to recover original raw-ratio behavior */
    simCtx->mom_last_converged  = PETSC_TRUE;
    simCtx->mom_last_lambda_max = 0.0;   /* populated after first momentum solve */
    simCtx->mom_nk_monitor_history = PETSC_FALSE;
    simCtx->ps_ksp_pic_monitor_true_residual = PETSC_FALSE;
    simCtx->cdisx = 0.0; simCtx->cdisy = 0.0; simCtx->cdisz = 0.0;
    simCtx->initialConditionMode = IC_MODE_ZERO;
    simCtx->initialConditionField = IC_FIELD_UCAT;
    strcpy(simCtx->initialConditionDirectory, "config/initial_condition");
    simCtx->InitialConstantContra.x = 0.0;
    simCtx->InitialConstantContra.y = 0.0;
    simCtx->InitialConstantContra.z = 0.0;
    simCtx->flowDirection = FLOW_DIR_UNSET;
    simCtx->icVelocityPhysical = 0.0;
    simCtx->AnalyticalUniformVelocity.x = 0.0;
    simCtx->AnalyticalUniformVelocity.y = 0.0;
    simCtx->AnalyticalUniformVelocity.z = 0.0;
    simCtx->fieldStatisticsEnabled = PETSC_FALSE;
    simCtx->fieldStatisticsWindowCount = 0;
    simCtx->fieldStatisticsWindows = NULL;
    simCtx->statisticsConsoleOutputFreq = 0;
    simCtx->fieldStatisticsContinue = PETSC_FALSE;
    simCtx->solutionConvergenceEnabled = PETSC_TRUE;
    simCtx->solutionConvergenceMode = SOLUTION_CONVERGENCE_STEADY_DETERMINISTIC;
    simCtx->solutionConvergencePeriodSteps = 0;
    simCtx->solutionConvergenceWindowSteps = 0;
    simCtx->solutionConvergenceSamplesRecorded = 0;
    simCtx->solutionConvergenceMeanSpeedHistory = NULL;
    simCtx->solutionConvergenceMeanKEHistory = NULL;
    simCtx->verificationDiffusivity.enabled = PETSC_FALSE;
    strcpy(simCtx->verificationDiffusivity.mode, "");
    strcpy(simCtx->verificationDiffusivity.profile, "");
    simCtx->verificationDiffusivity.gamma0 = 0.0;
    simCtx->verificationDiffusivity.slope_x = 0.0;
 
    // --- Group 6: Physical & Geometric Parameters ---
    simCtx->Flux_in = 1.0; simCtx->angle = 0.0;
    simCtx->max_angle = -54. * 3.1415926 / 180.;
    simCtx->CMx_c=0.0; simCtx->CMy_c=0.0; simCtx->CMz_c=0.0;
    simCtx->wall_roughness_height = 1e-16;
    simCtx->schmidt_number = 1.0; simCtx->Turbulent_schmidt_number = 0.7;
    simCtx->iem_constant = 2.0;
 
    // --- Group 7: Grid, Domain, and Boundary Condition Settings ---
    simCtx->block_number = 1; simCtx->inletprofile = 1;
    simCtx->i_periodic = 0; simCtx->j_periodic = 0; simCtx->k_periodic = 0;
    simCtx->pseudo_periodic = 0;
    strcpy(simCtx->grid_file, "config/grid.run");
    simCtx->generate_grid = PETSC_FALSE;
    simCtx->da_procs_x = PETSC_DECIDE;
    simCtx->da_procs_y = PETSC_DECIDE;
    simCtx->da_procs_z = PETSC_DECIDE;
    simCtx->num_bcs_files = 1;
    ierr = PetscMalloc1(1, &simCtx->bcs_files); CHKERRQ(ierr);
    ierr = PetscStrallocpy("config/bcs.run", &simCtx->bcs_files[0]); CHKERRQ(ierr);
    simCtx->FluxInSum = 0.0; simCtx->FluxOutSum = 0.0; simCtx->Fluxsum = 0.0;
    simCtx->drivingForceMagnitude = 0.0, simCtx->forceScalingFactor = 1.8;
    simCtx->drivingForceStep = -1;
    simCtx->drivenFluxMeasured = 0.0; simCtx->drivenFluxArea = 0.0;
    simCtx->targetVolumetricFlux  = 0.0;
    simCtx->drivenFluxTargetLatched = PETSC_FALSE;
    simCtx->bulkVelocityCorrection = 0.0;
    simCtx->boundaryVelocityCorrection = 0.0;
    simCtx->AreaInSum = 0.0; simCtx->AreaOutSum = 0.0;
    simCtx->ccc = 0;
    simCtx->ratio = 0.0;
    
    
    // --- Group 8: Turbulence Modeling (LES) ---
    simCtx->les = NO_LES_MODEL;
    simCtx->wallfunction = 0; simCtx->les_gradient_model = 0;
    ierr = LESConfigSetDefaults(&simCtx->les_config); CHKERRQ(ierr);
 
    // --- Group 9: Particle / DMSwarm Data & Settings ---
    simCtx->np = 0;
    simCtx->particleRandomSeed = 12345;
    simCtx->dm_swarm = NULL; simCtx->bboxlist = NULL;
    simCtx->ParticleInitialization = PARTICLE_INIT_SURFACE_RANDOM;
    simCtx->interpolationMethod = INTERP_TRILINEAR;
    strcpy(simCtx->particleRestartMode,"load");
    simCtx->particlesLostLastStep = 0;
    simCtx->particlesLostCumulative = 0;
    simCtx->particlesLostScalarLastStep = 0.0;
    simCtx->particlesMigratedLastStep = 0;
    simCtx->occupiedCellCount = 0;
    simCtx->particleLoadImbalance = 0.0;
    simCtx->migrationPassesLastStep = 0;
    simCtx->searchMetrics.searchAttempts = 0;
    simCtx->searchMetrics.searchPopulation = 0;
    simCtx->searchMetrics.searchLocatedCount = 0;
    simCtx->searchMetrics.searchLostCount = 0;
    simCtx->searchMetrics.traversalStepsSum = 0;
    simCtx->searchMetrics.reSearchCount = 0;
    simCtx->searchMetrics.maxTraversalSteps = 0;
    simCtx->searchMetrics.maxTraversalFailCount = 0;
    simCtx->searchMetrics.tieBreakCount = 0;
    simCtx->searchMetrics.boundaryClampCount = 0;
    simCtx->searchMetrics.bboxGuessSuccessCount = 0;
    simCtx->searchMetrics.bboxGuessFallbackCount = 0;
    simCtx->searchMetrics.maxParticlePassDepth = 0;
    simCtx->searchMetrics.currentSettlementPass = 0;
    simCtx->BrownianMotionRNG = NULL;
    simCtx->C_IEM = 2.0;
 
    // --- Group 10: Immersed Boundary & FSI Data Object Pointers ---
    simCtx->ibm = NULL; simCtx->ibmv = NULL; simCtx->fsi = NULL;
 
    // --- Group 11: Logging and Custom Configuration ---
    strcpy(simCtx->allowedFile, "config/whitelist.run");
    simCtx->useCfg = PETSC_FALSE;
    simCtx->allowedFuncs = NULL;
    simCtx->nAllowed = 0;
    simCtx->LoggingFrequency = 10;
    simCtx->poissonSourceImbalance = 0.0;
    simCtx->MaxDiv = 0.0;
    simCtx->MaxDivFlatArg = 0; simCtx->MaxDivx = 0; simCtx->MaxDivy = 0; simCtx->MaxDivz = 0;
    strcpy(simCtx->profilingSelectedFuncsFile, "config/profile.run");
    simCtx->useProfilingSelectedFuncsCfg =  PETSC_FALSE;
    simCtx->profilingSelectedFuncs = NULL;
    simCtx->nProfilingSelectedFuncs = 0;
    strcpy(simCtx->profilingTimestepMode, "selected");
    strcpy(simCtx->profilingTimestepFile, "Profiling_Timestep_Summary.csv");
    simCtx->profilingFinalSummary = PETSC_TRUE;
    simCtx->walltimeGuardEnabled = PETSC_FALSE;
    simCtx->walltimeGuardActive = PETSC_FALSE;
    simCtx->walltimeGuardWarmupSteps = 10;
    simCtx->walltimeGuardMultiplier = 2.0;
    simCtx->walltimeGuardMinSeconds = 60.0;
    simCtx->walltimeGuardEstimatorAlpha = 0.35;
    simCtx->walltimeGuardJobStartEpochSeconds = 0.0;
    simCtx->walltimeGuardLimitSeconds = 0.0;
    simCtx->walltimeGuardCompletedSteps = 0;
    simCtx->walltimeGuardWarmupTotalSeconds = 0.0;
    simCtx->walltimeGuardWarmupAverageSeconds = 0.0;
    simCtx->walltimeGuardHasEWMA = PETSC_FALSE;
    simCtx->walltimeGuardEWMASeconds = 0.0;
    simCtx->walltimeGuardLatestStepSeconds = 0.0;
    simCtx->runtimeMemoryLogEnabled = PETSC_TRUE;
    strcpy(simCtx->runtimeMemoryLogFile, "Runtime_Memory.log");
    simCtx->runtimeMemoryLogStarted = PETSC_FALSE;
    simCtx->runtimeMemoryLogHasPrevious = PETSC_FALSE;
    simCtx->runtimeMemoryLogPreviousProcessMB = 0.0;
    // --- Group 11: Post-Processing Information ---
    strcpy(simCtx->PostprocessingControlFile, "config/post.run");
    ierr = PetscNew(&simCtx->pps); CHKERRQ(ierr);
 
    // === 2. Get MPI Info and Handle Config File =============================
    // -- Group 1:  Parallelism & MPI Information
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &simCtx->rank); CHKERRQ(ierr);
    ierr = MPI_Comm_size(PETSC_COMM_WORLD, &simCtx->size); CHKERRQ(ierr);
    
 // First, check if the -control_file argument was provided by the user/script.
    ierr = PetscOptionsGetString(NULL, NULL, "-control_file", control_filename, sizeof(control_filename), &control_flg); CHKERRQ(ierr);
 
    // If the flag is NOT present or the filename is empty, abort with a helpful error.
    if (!control_flg || strlen(control_filename) == 0) {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
                "\n\n*** MANDATORY ARGUMENT MISSING ***\n"
                "The -control_file argument was not provided.\n"
                "This program must be launched with a configuration file.\n"
                "Example: mpiexec -n 4 ./simulator -control_file /path/to/your/config.control\n"
                "This is typically handled automatically by the 'picurv' script.\n");
    }
 
    // At this point, we have a valid filename. Attempt to load it.
    LOG(GLOBAL, LOG_INFO, "Loading mandatory configuration from: %s\n", control_filename);
    ierr = PetscOptionsInsertFile(PETSC_COMM_WORLD, NULL, control_filename, PETSC_FALSE);
    if (ierr == PETSC_ERR_FILE_OPEN) {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_FILE_OPEN, "The specified control file was not found or could not be opened: %s", control_filename);
    }
    CHKERRQ(ierr);
 
    {
        PetscBool legacy_averaging = PETSC_FALSE;
        ierr = PetscOptionsHasName(NULL, NULL, "-averaging", &legacy_averaging); CHKERRQ(ierr);
        PetscCheck(!legacy_averaging, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
                   "Legacy -averaging was removed. Use instantaneous output and offline "
                   "postprocessing until the replacement field-statistics pipeline is available.");
    }
 
    // === 3. A Configure Logging System ========================================
    // This logic determines the logging configuration and STORES it in simCtx for
    // later reference and cleanup.
    ierr = PetscOptionsGetString(NULL, NULL, "-whitelist_config_file", simCtx->allowedFile, PETSC_MAX_PATH_LEN, &simCtx->useCfg); CHKERRQ(ierr);
    
    if (simCtx->useCfg) {
        ierr = LoadAllowedFunctionsFromFile(simCtx->allowedFile, &simCtx->allowedFuncs, &simCtx->nAllowed);
        if (ierr) {
            // Use direct PetscPrintf as logging system isn't fully active yet.
            PetscPrintf(PETSC_COMM_SELF, "[%s] WARNING: Failed to load allowed functions from '%s'. Falling back to default list.\n", __func__, simCtx->allowedFile);
            simCtx->useCfg = PETSC_FALSE; // Mark as failed.
            ierr = 0; // Clear the error to allow fallback.
        } else if (simCtx->nAllowed == 0) {
            SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
                    "Whitelist config file '%s' is empty. Omit -whitelist_config_file to use the default allow-list, or list at least one function.",
                    simCtx->allowedFile);
        }
    }
    if (!simCtx->useCfg) {
        // Fallback to default logging functions if no file was used or if loading failed.
        simCtx->nAllowed = 2;
        ierr = PetscMalloc1(simCtx->nAllowed, &simCtx->allowedFuncs); CHKERRQ(ierr);
        ierr = PetscStrallocpy("main", &simCtx->allowedFuncs[0]); CHKERRQ(ierr);
        ierr = PetscStrallocpy("CreateSimulationContext", &simCtx->allowedFuncs[1]); CHKERRQ(ierr);
    }
    
    // Activate the configuration by passing it to the logging module's setup function.
    set_allowed_functions((const char**)simCtx->allowedFuncs, (size_t)simCtx->nAllowed);
    
    // Now that the logger is configured, we can use it.
    LOG_ALLOW_SYNC(LOCAL, LOG_INFO, "Context created. Initializing on rank %d of %d.\n", simCtx->rank, simCtx->size);
    print_log_level(); // This will now correctly reflect the LOG_LEVEL environment variable.
   
    // === 3.B Configure Profiling System ========================================
    ierr = PetscOptionsGetString(NULL, NULL, "-profiling_timestep_mode", simCtx->profilingTimestepMode, sizeof(simCtx->profilingTimestepMode), NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetString(NULL, NULL, "-profiling_timestep_file", simCtx->profilingTimestepFile, PETSC_MAX_PATH_LEN, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetBool(NULL, NULL, "-profiling_final_summary", &simCtx->profilingFinalSummary, NULL); CHKERRQ(ierr);
    if (strcmp(simCtx->profilingTimestepMode, "off") != 0 &&
        strcmp(simCtx->profilingTimestepMode, "selected") != 0 &&
        strcmp(simCtx->profilingTimestepMode, "all") != 0) {
        PetscPrintf(PETSC_COMM_SELF, "[%s] WARNING: Unknown profiling timestep mode '%s'. Falling back to 'selected'.\n", __func__, simCtx->profilingTimestepMode);
        strcpy(simCtx->profilingTimestepMode, "selected");
    }
 
    if (strcmp(simCtx->profilingTimestepMode, "selected") == 0) {
        ierr = PetscOptionsGetString(NULL, NULL, "-profile_config_file", simCtx->profilingSelectedFuncsFile, PETSC_MAX_PATH_LEN, &simCtx->useProfilingSelectedFuncsCfg); CHKERRQ(ierr);
        if (simCtx->useProfilingSelectedFuncsCfg) {
            ierr = LoadAllowedFunctionsFromFile(simCtx->profilingSelectedFuncsFile, &simCtx->profilingSelectedFuncs, &simCtx->nProfilingSelectedFuncs);
            if (ierr) {
                PetscPrintf(PETSC_COMM_SELF, "[%s] WARNING: Failed to load selected profiling functions from '%s'. Falling back to default list.\n", __func__, simCtx->profilingSelectedFuncsFile);
                simCtx->useProfilingSelectedFuncsCfg = PETSC_FALSE;
                ierr = 0;
            }
        }
        if (!simCtx->useProfilingSelectedFuncsCfg) {
            // Fallback to a hardcoded default list if no file was provided or loading failed.
            simCtx->nProfilingSelectedFuncs = 4;
            ierr = PetscMalloc1(simCtx->nProfilingSelectedFuncs, &simCtx->profilingSelectedFuncs); CHKERRQ(ierr);
            ierr = PetscStrallocpy("FlowSolver", &simCtx->profilingSelectedFuncs[0]); CHKERRQ(ierr);
            ierr = PetscStrallocpy("AdvanceSimulation", &simCtx->profilingSelectedFuncs[1]); CHKERRQ(ierr);
            ierr = PetscStrallocpy("LocateAllParticlesInGrid", &simCtx->profilingSelectedFuncs[2]); CHKERRQ(ierr);
            ierr = PetscStrallocpy("InterpolateAllFieldsToSwarm", &simCtx->profilingSelectedFuncs[3]); CHKERRQ(ierr);
        }
    }
 
    // Initialize the profiling system with the current updated simulation context.
    ierr = ProfilingInitialize(simCtx); CHKERRQ(ierr);
 
    // === 4. Parse All Command Line Options ==================================
    LOG_ALLOW(GLOBAL, LOG_INFO, "Parsing command-line options...\n");
 
    // --- Group 2
    LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 2: Simulation Control,Time and I/O.\n");
    // Read the physical time to start from.
    // The default is already 0.0, so this will only be non-zero if the user provides it.
    ierr = PetscOptionsGetInt(NULL, NULL, "-start_step", &simCtx->StartStep, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL,NULL, "-totalsteps", &simCtx->StepsToRun, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetBool(NULL, NULL, "-only_setup", &simCtx->OnlySetup, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetBool(NULL, NULL, "-continue_mode", &simCtx->continueMode, NULL); CHKERRQ(ierr);
    /* Resuming in the same run directory continues accumulated statistics by
     * default; branching with --restart-from must opt in, because the branch may
     * follow a different physical trajectory than the samples already collected. */
    simCtx->fieldStatisticsContinue = simCtx->continueMode;
    ierr = PetscOptionsGetBool(NULL, NULL, "-field_statistics_continue",
                               &simCtx->fieldStatisticsContinue, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-dt", &simCtx->dt, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-tio", &simCtx->tiout, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-particle_console_output_freq", &simCtx->particleConsoleOutputFreq, &particle_console_output_freq_flg); CHKERRQ(ierr);
    if (!particle_console_output_freq_flg) {
        simCtx->particleConsoleOutputFreq = simCtx->tiout;
    }
    ierr = PetscOptionsGetString(NULL,NULL,"-euler_field_source",simCtx->eulerianSource,sizeof(simCtx->eulerianSource),NULL);CHKERRQ(ierr);
    ierr = PetscOptionsGetString(NULL,NULL,"-output_dir",simCtx->output_dir,sizeof(simCtx->output_dir),NULL);CHKERRQ(ierr);
    ierr = PetscOptionsGetString(NULL,NULL,"-restart_dir",simCtx->restart_dir,sizeof(simCtx->restart_dir),NULL);CHKERRQ(ierr);
    ierr = PetscOptionsGetString(NULL,NULL,"-log_dir",simCtx->log_dir,sizeof(simCtx->log_dir),NULL);CHKERRQ(ierr);
    ierr = PetscOptionsGetString(NULL,NULL,"-analysis_dir",simCtx->analysis_dir,sizeof(simCtx->analysis_dir),NULL);CHKERRQ(ierr);
    ierr = PetscOptionsGetBool(NULL, NULL, "-walltime_guard_enabled", &simCtx->walltimeGuardEnabled, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-walltime_guard_warmup_steps", &simCtx->walltimeGuardWarmupSteps, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-walltime_guard_multiplier", &simCtx->walltimeGuardMultiplier, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetBool(NULL, NULL, "-runtime_memory_log_enabled", &simCtx->runtimeMemoryLogEnabled, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetString(NULL, NULL, "-runtime_memory_log_file", simCtx->runtimeMemoryLogFile, PETSC_MAX_PATH_LEN, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-walltime_guard_min_seconds", &simCtx->walltimeGuardMinSeconds, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-walltime_guard_estimator_alpha", &simCtx->walltimeGuardEstimatorAlpha, NULL); CHKERRQ(ierr);
 
    if (simCtx->walltimeGuardWarmupSteps <= 0) {
      SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Invalid value for -walltime_guard_warmup_steps: %d. Must be > 0.", simCtx->walltimeGuardWarmupSteps);
    }
    if (simCtx->walltimeGuardMultiplier <= 0.0 || simCtx->walltimeGuardMultiplier > 5.0) {
      SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Invalid value for -walltime_guard_multiplier: %.6f. Must be in (0, 5].", (double)simCtx->walltimeGuardMultiplier);
    }
    if (simCtx->walltimeGuardMinSeconds <= 0.0) {
      SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Invalid value for -walltime_guard_min_seconds: %.6f. Must be > 0.", (double)simCtx->walltimeGuardMinSeconds);
    }
    if (simCtx->walltimeGuardEstimatorAlpha <= 0.0 || simCtx->walltimeGuardEstimatorAlpha > 1.0) {
      SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Invalid value for -walltime_guard_estimator_alpha: %.6f. Must be in (0, 1].", (double)simCtx->walltimeGuardEstimatorAlpha);
    }
 
    if(strcmp(simCtx->eulerianSource,"solve")!= 0 && strcmp(simCtx->eulerianSource,"load") != 0 && strcmp(simCtx->eulerianSource,"analytical")!=0){
      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);
    }
    if (simCtx->walltimeGuardEnabled) {
      const char *job_start_env = getenv("PICURV_JOB_START_EPOCH");
      const char *limit_env = getenv("PICURV_WALLTIME_LIMIT_SECONDS");
      PetscBool job_start_ok = RuntimeWalltimeGuardParsePositiveSeconds(job_start_env, &simCtx->walltimeGuardJobStartEpochSeconds);
      PetscBool limit_ok = RuntimeWalltimeGuardParsePositiveSeconds(limit_env, &simCtx->walltimeGuardLimitSeconds);
 
      if (!job_start_ok || !limit_ok) {
        simCtx->walltimeGuardActive = PETSC_FALSE;
        simCtx->walltimeGuardJobStartEpochSeconds = 0.0;
        simCtx->walltimeGuardLimitSeconds = 0.0;
        LOG_ALLOW(
          GLOBAL,
          LOG_WARNING,
          "Runtime walltime guard enabled but %s/%s are missing or invalid. Falling back to external shutdown signals only.\n",
          "PICURV_JOB_START_EPOCH",
          "PICURV_WALLTIME_LIMIT_SECONDS"
        );
      } else {
        simCtx->walltimeGuardActive = PETSC_TRUE;
      }
    }
 
    //  --- Group 3
    LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 3: High-Level Physics & Model Selection Flags\n");    
    ierr = PetscOptionsGetInt(NULL, NULL, "-imm", &simCtx->immersed, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-fsi", &simCtx->movefsi, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-rfsi", &simCtx->rotatefsi, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-inv", &simCtx->invicid, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-TwoD", &simCtx->TwoD, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-mframe", &simCtx->moveframe, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-rframe", &simCtx->rotateframe, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetString(NULL,NULL,"-analytical_type",simCtx->AnalyticalSolutionType,sizeof(simCtx->AnalyticalSolutionType),NULL);CHKERRQ(ierr);
 
    /* Immersed boundaries, moving bodies and moving reference frames are planned, not
       implemented. No body geometry is ever loaded, and the IBM interpolation the
       momentum solvers would call is commented out and defined nowhere, so -imm only
       reconfigures the Poisson solve around a body that is not there. The moving-frame
       convection branch in ComputeRHS is commented out as well, so -mframe/-rframe would
       drop the convective term. Each would run a different problem without saying so;
       the conductor refuses the YAML switches, and this refuses the same flags arriving
       through a PETSc passthrough. */
    PetscCheck(!simCtx->immersed && !simCtx->movefsi && !simCtx->rotatefsi, PETSC_COMM_WORLD, PETSC_ERR_SUP,
               "Immersed boundaries and moving bodies (-imm, -fsi, -rfsi) are planned, not implemented.");
    PetscCheck(!simCtx->moveframe && !simCtx->rotateframe, PETSC_COMM_WORLD, PETSC_ERR_SUP,
               "Moving and rotating reference frames (-mframe, -rframe) are planned, not implemented: "
               "their convection branch is absent, so the convective term would be dropped.");
 
    //  --- Group 4
    LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 4: Immersed-Body Flux Corrections and Driven-Flow Controls\n");
    ierr = PetscOptionsGetInt(NULL, NULL, "-mhv", &simCtx->MHV, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-lv", &simCtx->LV, NULL); CHKERRQ(ierr);
    /* The heart-valve volume-flux corrections balance flux across an immersed body, which
       cannot exist while immersed boundaries are unimplemented (see the Group 3 check). */
    PetscCheck(!simCtx->MHV && !simCtx->LV, PETSC_COMM_WORLD, PETSC_ERR_SUP,
               "The -mhv/-lv immersed-body flux corrections need immersed boundaries, which are planned, not implemented.");
    ierr = PetscOptionsGetReal(NULL,NULL,"-driven_flow_initial_force",&simCtx->drivingForceMagnitude,NULL);CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL,NULL,"-driven_flow_scaling_factor",&simCtx->forceScalingFactor,NULL);CHKERRQ(ierr);
    //  --- Group 5
    LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 5: Solver & Numerics Parameters \n");
    char mom_solver_type_char[PETSC_MAX_PATH_LEN];
    char solution_convergence_mode_char[PETSC_MAX_PATH_LEN];
    PetscBool mom_solver_type_flg = PETSC_FALSE;
    PetscBool solution_convergence_mode_flg = PETSC_FALSE;
    ierr = PetscOptionsGetString(NULL, NULL, "-mom_solver_type", mom_solver_type_char, sizeof(mom_solver_type_char), &mom_solver_type_flg); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-mom_max_pseudo_steps", &simCtx->mom_max_pseudo_steps, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-mom_atol", &simCtx->mom_atol, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-mom_rtol", &simCtx->mom_rtol, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-mom_resid_atol", &simCtx->mom_resid_atol, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-mom_resid_rtol", &simCtx->mom_resid_rtol, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-imp_stol", &simCtx->imp_stol, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-central", &simCtx->central, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetString(NULL, NULL, "-solution_convergence_mode",
                                 solution_convergence_mode_char, sizeof(solution_convergence_mode_char),
                                 &solution_convergence_mode_flg); CHKERRQ(ierr);
    ierr = PetscOptionsGetBool(NULL, NULL, "-solution_convergence_enabled", &simCtx->solutionConvergenceEnabled, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-solution_convergence_period_steps", &simCtx->solutionConvergencePeriodSteps, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-solution_convergence_window_steps", &simCtx->solutionConvergenceWindowSteps, NULL); CHKERRQ(ierr);
    /* Field statistics resolve here rather than inline, because a variable-arity
     * window list needs its own parse site; see statistics_config.c. It must
     * precede CreateAndInitializeAllVectors, which sizes the per-window
     * accumulators from the window count resolved here. */
    ierr = ParseFieldStatisticsConfig(simCtx); CHKERRQ(ierr);
 
    // Keep parser acceptance aligned with the enum and FlowSolver dispatch.
    if (mom_solver_type_flg) {
        if(strcmp(mom_solver_type_char, "DUALTIME_PICARD_JAMESON_RK") == 0 ||
           strcmp(mom_solver_type_char, "DUALTIME_PICARD_RK4") == 0) {
            simCtx->mom_solver_type = MOMENTUM_SOLVER_DUALTIME_PICARD_JAMESON_RK;
        } else if (strcmp(mom_solver_type_char, "EXPLICIT_RK") == 0) {
            simCtx->mom_solver_type = MOMENTUM_SOLVER_EXPLICIT_RK;
        } else if (strcmp(mom_solver_type_char, "newton_krylov") == 0) {
            simCtx->mom_solver_type = MOMENTUM_SOLVER_NEWTON_KRYLOV;
        } else {
            LOG(GLOBAL, LOG_ERROR, "Invalid value for -mom_solver_type: '%s'. Valid options are: 'DUALTIME_PICARD_JAMESON_RK', 'EXPLICIT_RK', 'newton_krylov'.\n", mom_solver_type_char);
            SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Invalid value for -mom_solver_type: '%s'.", mom_solver_type_char);
        }
    }
 
    if (solution_convergence_mode_flg) {
        if (strcmp(solution_convergence_mode_char, "STEADY_DETERMINISTIC") == 0) {
            simCtx->solutionConvergenceMode = SOLUTION_CONVERGENCE_STEADY_DETERMINISTIC;
        } else if (strcmp(solution_convergence_mode_char, "PERIODIC_DETERMINISTIC") == 0) {
            simCtx->solutionConvergenceMode = SOLUTION_CONVERGENCE_PERIODIC_DETERMINISTIC;
        } else if (strcmp(solution_convergence_mode_char, "STATISTICAL_STEADY") == 0) {
            simCtx->solutionConvergenceMode = SOLUTION_CONVERGENCE_STATISTICAL_STEADY;
        } else if (strcmp(solution_convergence_mode_char, "TRANSIENT") == 0) {
            simCtx->solutionConvergenceMode = SOLUTION_CONVERGENCE_TRANSIENT;
        } else {
            SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
                    "Invalid value for -solution_convergence_mode: '%s'.", solution_convergence_mode_char);
        }
    }
 
    if (simCtx->solutionConvergenceMode == SOLUTION_CONVERGENCE_PERIODIC_DETERMINISTIC &&
        simCtx->solutionConvergencePeriodSteps <= 0) {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
                "solution convergence mode PERIODIC_DETERMINISTIC requires -solution_convergence_period_steps > 0.");
    }
    if (simCtx->solutionConvergenceMode == SOLUTION_CONVERGENCE_STATISTICAL_STEADY &&
        simCtx->solutionConvergenceWindowSteps <= 0) {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
                "solution convergence mode STATISTICAL_STEADY requires -solution_convergence_window_steps > 0.");
    }
 
    // --- Multigrid Options ---
    ierr = PetscOptionsGetInt(NULL, NULL, "-mg_level", &simCtx->mglevels, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-mg_pre_it", &simCtx->mg_preItr, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-mg_post_it", &simCtx->mg_poItr, NULL); CHKERRQ(ierr);
 
    // --- Other Solver Options ---
    ierr = PetscOptionsGetInt(NULL, NULL, "-poisson", &simCtx->poisson, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-ren", &simCtx->ren, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-pseudo_cfl", &simCtx->pseudo_cfl, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-max_pseudo_cfl", &simCtx->max_pseudo_cfl, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-min_pseudo_cfl", &simCtx->min_pseudo_cfl, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-pseudo_cfl_reduction_factor", &simCtx->pseudo_cfl_reduction_factor, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-pseudo_cfl_growth_factor", &simCtx->pseudo_cfl_growth_factor, NULL); CHKERRQ(ierr);
    // Read the deprecated RK4 spelling first so the canonical Jameson option wins if both are present.
    ierr = PetscOptionsGetReal(NULL,NULL, "-mom_dt_rk4_residual_norm_noise_allowance_factor",&simCtx->mom_dt_jameson_residual_norm_noise_allowance_factor,NULL);CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL,NULL, "-mom_dt_jameson_residual_norm_noise_allowance_factor",&simCtx->mom_dt_jameson_residual_norm_noise_allowance_factor,NULL);CHKERRQ(ierr);
    ierr = PetscOptionsGetBool(NULL, NULL, "-no_pseudo_cfl_backtrack", &simCtx->no_pseudo_cfl_backtrack, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-mom_ratio_ema_alpha", &simCtx->mom_ratio_ema_alpha, NULL); CHKERRQ(ierr);
    if (simCtx->min_pseudo_cfl <= 0.0 ||
        simCtx->pseudo_cfl < simCtx->min_pseudo_cfl ||
        simCtx->pseudo_cfl > simCtx->max_pseudo_cfl) {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE,
                "Pseudo-CFL controls require 0 < minimum <= initial <= maximum.");
    }
    if (simCtx->pseudo_cfl_growth_factor < 1.0 ||
        simCtx->pseudo_cfl_reduction_factor <= 0.0 ||
        simCtx->pseudo_cfl_reduction_factor >= 1.0 ||
        simCtx->mom_dt_jameson_residual_norm_noise_allowance_factor < 1.0) {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE,
                "Pseudo-CFL controls require growth_factor >= 1, 0 < reduction_factor < 1, and noise allowance >= 1.");
    }
    if (simCtx->mom_ratio_ema_alpha < 0.0 || simCtx->mom_ratio_ema_alpha > 1.0) {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE,
                "-mom_ratio_ema_alpha must be in [0, 1].");
    }
    ierr = PetscOptionsHasName(NULL, NULL, "-ps_ksp_pic_monitor_true_residual", &simCtx->ps_ksp_pic_monitor_true_residual); CHKERRQ(ierr);
    ierr = PetscOptionsGetBool(NULL, NULL, "-mom_nk_pic_monitor", &simCtx->mom_nk_monitor_history, NULL); CHKERRQ(ierr);
    {
        PetscInt ic_mode = (PetscInt)simCtx->initialConditionMode;
        PetscInt ic_field = (PetscInt)simCtx->initialConditionField;
        ierr = PetscOptionsGetInt(NULL, NULL, "-finit", &ic_mode, NULL); CHKERRQ(ierr);
        ierr = PetscOptionsGetInt(NULL, NULL, "-ic_field", &ic_field, NULL); CHKERRQ(ierr);
        simCtx->initialConditionMode = (InitialConditionMode)ic_mode;
        simCtx->initialConditionField = (InitialConditionField)ic_field;
    }
    ierr = PetscOptionsGetString(NULL, NULL, "-ic_dir", simCtx->initialConditionDirectory,
                                 sizeof(simCtx->initialConditionDirectory), NULL); CHKERRQ(ierr);
    if (simCtx->initialConditionMode < IC_MODE_ZERO || simCtx->initialConditionMode > IC_MODE_FILE) {
      SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE,
              "Invalid value for -finit. Expected an initial-condition mode in [0,4], got %d.",
              simCtx->initialConditionMode);
    }
    if (simCtx->initialConditionField < IC_FIELD_UCAT || simCtx->initialConditionField > IC_FIELD_UCONT) {
      SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE,
              "Invalid value for -ic_field. Expected 0 (Ucat) or 1 (Ucont), got %d.",
              simCtx->initialConditionField);
    }
    ierr = PetscOptionsGetReal(NULL, NULL, "-ucont_x", &simCtx->InitialConstantContra.x, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-ucont_y", &simCtx->InitialConstantContra.y, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-ucont_z", &simCtx->InitialConstantContra.z, NULL); CHKERRQ(ierr);
    {
        PetscInt  fd_int = (PetscInt)FLOW_DIR_UNSET;
        PetscBool fd_set = PETSC_FALSE;
        ierr = PetscOptionsGetInt(NULL, NULL, "-flow_direction", &fd_int, &fd_set); CHKERRQ(ierr);
        if (fd_set) simCtx->flowDirection = (FlowDirection)fd_int;
    }
    ierr = PetscOptionsGetReal(NULL, NULL, "-ic_velocity_physical", &simCtx->icVelocityPhysical, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-analytical_uniform_u", &simCtx->AnalyticalUniformVelocity.x, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-analytical_uniform_v", &simCtx->AnalyticalUniformVelocity.y, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-analytical_uniform_w", &simCtx->AnalyticalUniformVelocity.z, NULL); CHKERRQ(ierr);
    PetscBool verification_scalar_value_set = PETSC_FALSE;
    PetscBool verification_scalar_phi0_set = PETSC_FALSE;
    PetscBool verification_scalar_slope_x_set = PETSC_FALSE;
    PetscBool verification_scalar_amplitude_set = PETSC_FALSE;
    PetscBool verification_scalar_kx_set = PETSC_FALSE;
    PetscBool verification_scalar_ky_set = PETSC_FALSE;
    PetscBool verification_scalar_kz_set = PETSC_FALSE;
    ierr = PetscOptionsGetString(NULL, NULL, "-verification_diffusivity_mode",
                                 simCtx->verificationDiffusivity.mode,
                                 sizeof(simCtx->verificationDiffusivity.mode), NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetString(NULL, NULL, "-verification_diffusivity_profile",
                                 simCtx->verificationDiffusivity.profile,
                                 sizeof(simCtx->verificationDiffusivity.profile), NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-verification_diffusivity_gamma0",
                               &simCtx->verificationDiffusivity.gamma0, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-verification_diffusivity_slope_x",
                               &simCtx->verificationDiffusivity.slope_x, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetString(NULL, NULL, "-verification_scalar_mode",
                                 simCtx->verificationScalar.mode,
                                 sizeof(simCtx->verificationScalar.mode), NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetString(NULL, NULL, "-verification_scalar_profile",
                                 simCtx->verificationScalar.profile,
                                 sizeof(simCtx->verificationScalar.profile), NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-verification_scalar_value",
                               &simCtx->verificationScalar.value, &verification_scalar_value_set); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-verification_scalar_phi0",
                               &simCtx->verificationScalar.phi0, &verification_scalar_phi0_set); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-verification_scalar_slope_x",
                               &simCtx->verificationScalar.slope_x, &verification_scalar_slope_x_set); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-verification_scalar_amplitude",
                               &simCtx->verificationScalar.amplitude, &verification_scalar_amplitude_set); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-verification_scalar_kx",
                               &simCtx->verificationScalar.kx, &verification_scalar_kx_set); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-verification_scalar_ky",
                               &simCtx->verificationScalar.ky, &verification_scalar_ky_set); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-verification_scalar_kz",
                               &simCtx->verificationScalar.kz, &verification_scalar_kz_set); CHKERRQ(ierr);
    simCtx->verificationDiffusivity.enabled =
      (PetscBool)(simCtx->verificationDiffusivity.mode[0] != '\0' ||
                  simCtx->verificationDiffusivity.profile[0] != '\0');
    simCtx->verificationScalar.enabled =
      (PetscBool)(simCtx->verificationScalar.mode[0] != '\0' ||
                  simCtx->verificationScalar.profile[0] != '\0');
    if (simCtx->verificationDiffusivity.enabled) {
      if (strcmp(simCtx->eulerianSource, "analytical") != 0) {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE,
                "verification diffusivity overrides require -euler_field_source \"analytical\".");
      }
      if (strcmp(simCtx->verificationDiffusivity.mode, "analytical") != 0) {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
                "Unsupported -verification_diffusivity_mode '%s'. Only 'analytical' is supported.",
                simCtx->verificationDiffusivity.mode);
      }
      if (strcmp(simCtx->verificationDiffusivity.profile, "LINEAR_X") != 0) {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
                "Unsupported -verification_diffusivity_profile '%s'. Only 'LINEAR_X' is supported.",
                simCtx->verificationDiffusivity.profile);
      }
    }
    if (simCtx->verificationScalar.enabled) {
      if (strcmp(simCtx->eulerianSource, "analytical") != 0) {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE,
                "verification scalar overrides require -euler_field_source \"analytical\".");
      }
      if (strcmp(simCtx->verificationScalar.mode, "analytical") != 0) {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
                "Unsupported -verification_scalar_mode '%s'. Only 'analytical' is supported.",
                simCtx->verificationScalar.mode);
      }
      if (strcmp(simCtx->verificationScalar.profile, "CONSTANT") == 0) {
        if (!verification_scalar_value_set) {
          SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
                  "verification scalar profile CONSTANT requires -verification_scalar_value.");
        }
      } else if (strcmp(simCtx->verificationScalar.profile, "LINEAR_X") == 0) {
        if (!verification_scalar_phi0_set || !verification_scalar_slope_x_set) {
          SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
                  "verification scalar profile LINEAR_X requires -verification_scalar_phi0 and -verification_scalar_slope_x.");
        }
      } else if (strcmp(simCtx->verificationScalar.profile, "SIN_PRODUCT") == 0) {
        if (!verification_scalar_amplitude_set || !verification_scalar_kx_set ||
            !verification_scalar_ky_set || !verification_scalar_kz_set) {
          SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
                  "verification scalar profile SIN_PRODUCT requires -verification_scalar_amplitude, -verification_scalar_kx, -verification_scalar_ky, and -verification_scalar_kz.");
        }
      } else {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
                "Unsupported -verification_scalar_profile '%s'. Supported profiles: CONSTANT, LINEAR_X, SIN_PRODUCT.",
                simCtx->verificationScalar.profile);
      }
    }
    // NOTE: cdisx,cdisy,cdisz haven't been parsed, add if necessary.
 
     //  --- Group 6
    LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 6: Physical & Geometric Parameters \n");   
    ierr = PetscOptionsGetReal(NULL,NULL,"-schmidt_number",&simCtx->schmidt_number,NULL);CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL,NULL,"-turb_schmidt_number",&simCtx->Turbulent_schmidt_number,NULL);CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL,NULL,"-iem_constant",&simCtx->iem_constant,NULL);CHKERRQ(ierr);
    /* Zero is valid: it switches micromixing off, so a particle scalar is carried unchanged. */
    PetscCheck(simCtx->iem_constant >= 0.0 && !PetscIsInfOrNanReal(simCtx->iem_constant),
               PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE,
               "-iem_constant must be a non-negative finite number (got %g).", (double)simCtx->iem_constant);
    ierr = PetscOptionsGetReal(NULL,NULL,"-wall_roughness",&simCtx->wall_roughness_height,NULL);CHKERRQ(ierr);
    // NOTE: angle is not parsed in the original code, it set programmatically. We will follow that.
    // NOTE: max_angle is calculated based on other flags (like MHV) in the legacy code.
    // We will defer that logic to a later setup stage and not parse them directly.
    // The Scaling Information is calculated here
    ierr = ParseScalingInformation(simCtx); CHKERRQ(ierr);
 
     //  --- Group 7
    LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 7: Grid, Domain, and Boundary Condition Settings \n");
    ierr = PetscOptionsGetInt(NULL, NULL, "-nblk", &simCtx->block_number, NULL); CHKERRQ(ierr); // This is also a modern option
    /* Every inter-block exchange (Block_Interface_U) is commented out in the momentum
       solvers and the initial condition, so blocks would be solved as isolated domains.
       Multi-block coupling is planned, not implemented. */
    PetscCheck(simCtx->block_number == 1, PETSC_COMM_WORLD, PETSC_ERR_SUP,
               "Multi-block domains are planned, not implemented: blocks are never coupled (got -nblk %" PetscInt_FMT ").",
               simCtx->block_number);
    ierr = PetscOptionsGetInt(NULL, NULL, "-inlet", &simCtx->inletprofile, NULL); CHKERRQ(ierr);
    // NOTE: channelz was not parsed, likely set programmatically. We will omit its parsing call.
    ierr = PetscOptionsGetBool(NULL, NULL, "-grid", &simCtx->generate_grid, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetString(NULL, NULL, "-grid_file", simCtx->grid_file, PETSC_MAX_PATH_LEN, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-da_processors_x", &simCtx->da_procs_x, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-da_processors_y", &simCtx->da_procs_y, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-da_processors_z", &simCtx->da_procs_z, NULL); CHKERRQ(ierr);
    // NOTE: pseudo_periodic was not parsed. We will omit its parsing call. 
    PetscBool bcs_flg;
    char      file_list_str[PETSC_MAX_PATH_LEN * 10]; // Buffer for comma-separated list
 
    ierr = PetscOptionsGetString(NULL, NULL, "-bcs_files", file_list_str, sizeof(file_list_str), &bcs_flg); CHKERRQ(ierr);
 
    if (bcs_flg) {
      LOG_ALLOW(GLOBAL, LOG_DEBUG, "Found -bcs_files option, overriding default.\n");
    
      // A. Clean up the default memory we allocated in Phase 1.
      ierr = PetscFree(simCtx->bcs_files[0]); CHKERRQ(ierr);
      ierr = PetscFree(simCtx->bcs_files); CHKERRQ(ierr);
      simCtx->num_bcs_files = 0;
      simCtx->bcs_files = NULL;
 
      // B. Parse the user-provided comma-separated list.
      char *token;
      char *str_copy;
      ierr = PetscStrallocpy(file_list_str, &str_copy); CHKERRQ(ierr);
 
      // First pass: count the number of files.
      token = strtok(str_copy, ",");
      while (token) {
        simCtx->num_bcs_files++;
        token = strtok(NULL, ",");
      }
      ierr = PetscFree(str_copy); CHKERRQ(ierr);
    
      // Second pass: allocate memory and store the filenames.
      ierr = PetscMalloc1(simCtx->num_bcs_files, &simCtx->bcs_files); CHKERRQ(ierr);
      ierr = PetscStrallocpy(file_list_str, &str_copy); CHKERRQ(ierr);
      token = strtok(str_copy, ",");
      for (PetscInt i = 0; i < simCtx->num_bcs_files; i++) {
        ierr = PetscStrallocpy(token, &simCtx->bcs_files[i]); CHKERRQ(ierr);
        token = strtok(NULL, ",");
      }
      ierr = PetscFree(str_copy); CHKERRQ(ierr);
    }
    
 
     //  --- Group 8
    LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 8: Turbulence Modeling (LES) \n");
    // Seeded from the default already in simCtx: PetscOptionsGetInt leaves the target
    // untouched when -les is absent, and an uninitialized value would become the model.
    PetscInt temp_les_model = (PetscInt)simCtx->les;
    ierr = PetscOptionsGetInt(NULL, NULL, "-les", &temp_les_model, NULL); CHKERRQ(ierr);
    PetscCheck(temp_les_model >= NO_LES_MODEL && temp_les_model <= WALE, PETSC_COMM_WORLD,
               PETSC_ERR_ARG_OUTOFRANGE,
               "-les must be 0 (none), 1 (constant_smagorinsky), 2 (dynamic_smagorinsky), "
               "3 (vreman), or 4 (wale); received %" PetscInt_FMT ".", temp_les_model);
    simCtx->les = (LESModelType)temp_les_model;
    {
        /* RANS is planned, not implemented: the k-omega transport update was never
           written and its fields were never allocated. The selector and its dead
           hooks were removed on 2026-09-22 (see src/guide.md); this refuses the flag
           arriving through a PETSc passthrough. */
        PetscInt requested_rans = 0;
        ierr = PetscOptionsGetInt(NULL, NULL, "-rans", &requested_rans, NULL); CHKERRQ(ierr);
        PetscCheck(requested_rans == 0, PETSC_COMM_WORLD, PETSC_ERR_SUP,
                   "RANS closures (-rans) are planned, not implemented: no transport equation is "
                   "solved for them. Use an LES closure, or run laminar.");
    }
    ierr = PetscOptionsGetInt(NULL, NULL, "-wallfunction", &simCtx->wallfunction, NULL); CHKERRQ(ierr);
    PetscCheck(simCtx->wallfunction >= WALL_FUNCTION_NONE && simCtx->wallfunction <= WALL_FUNCTION_CABOT,
               PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE,
               "-wallfunction must be 0 (none), 1 (log_law), 2 (werner), or 3 (cabot); received %" PetscInt_FMT ".",
               simCtx->wallfunction);
    ierr = PetscOptionsGetInt(NULL, NULL, "-les_gradient_model", &simCtx->les_gradient_model, NULL); CHKERRQ(ierr);
    ierr = ParseLESConfiguration(simCtx); CHKERRQ(ierr);
 
     //  --- Group 9
    LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 9:  Particle / DMSwarm Data & Settings \n");
    ierr = PetscOptionsGetInt(NULL, NULL, "-numParticles", &simCtx->np, NULL); CHKERRQ(ierr);
    PetscInt temp_pinit = (PetscInt)PARTICLE_INIT_SURFACE_RANDOM;
    ierr = PetscOptionsGetInt(NULL, NULL, "-pinit", &temp_pinit, NULL); CHKERRQ(ierr);
    simCtx->ParticleInitialization = (ParticleInitializationType)temp_pinit;
    PetscInt temp_interp = (PetscInt)INTERP_TRILINEAR;
    ierr = PetscOptionsGetInt(NULL, NULL, "-interpolation_method", &temp_interp, NULL); CHKERRQ(ierr);
    simCtx->interpolationMethod = (InterpolationMethod)temp_interp;
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Interpolation method: %s\n",
              simCtx->interpolationMethod == INTERP_TRILINEAR ? "Trilinear (direct cell-center)" : "CornerAveraged (legacy)");
    ierr = PetscOptionsGetReal(NULL, NULL, "-psrc_x", &simCtx->psrc_x, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-psrc_y", &simCtx->psrc_y, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-psrc_z", &simCtx->psrc_z, NULL); CHKERRQ(ierr);
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Particle initialization mode: %s. Point source: (%.6f, %.6f, %.6f)\n",
              ParticleInitializationToString(simCtx->ParticleInitialization),
              simCtx->psrc_x, simCtx->psrc_y, simCtx->psrc_z);
    ierr = PetscOptionsGetString(NULL,NULL,"-particle_restart_mode",simCtx->particleRestartMode,sizeof(simCtx->particleRestartMode),NULL); CHKERRQ(ierr);
    // Validation for Particle Restart Mode
    if (strcmp(simCtx->particleRestartMode, "load") != 0 && strcmp(simCtx->particleRestartMode, "init") != 0) {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Invalid value for -particle_restart_mode. Must be 'load' or 'init'. You provided '%s'.", simCtx->particleRestartMode);
    }
    ierr = PetscOptionsGetInt(NULL, NULL, "-particle_random_seed", &simCtx->particleRandomSeed, NULL); CHKERRQ(ierr);
    PetscCheck(simCtx->particleRandomSeed >= 0, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE,
               "-particle_random_seed must be non-negative (got %" PetscInt_FMT ").", simCtx->particleRandomSeed);
    ierr = InitializeBrownianRNG(simCtx); CHKERRQ(ierr);
    ierr = ParticleFieldPlanCreate(&simCtx->particleFieldPlan); CHKERRQ(ierr);
    PetscCheck(!simCtx->particleFieldPlan || simCtx->np > 0, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
               "Particle field initial values are configured but the run has no particles.");
    PetscCheck(!simCtx->particleFieldPlan || !VerificationScalarOverrideActive(simCtx), PETSC_COMM_WORLD,
               PETSC_ERR_ARG_WRONG,
               "Particle field initial values cannot be combined with the verification scalar source, "
               "which prescribes Psi at every step.");
    // --- Group 10
    LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 10: Immersed Boundary & FSI Data Object Pointers \n");
 
    // --- Group 11
    LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 11: Top-Level Managers & Custom Configuration \n");
    ierr = PetscOptionsGetInt(NULL, NULL, "-logfreq", &simCtx->LoggingFrequency, NULL); CHKERRQ(ierr);
 
    if (simCtx->num_bcs_files != simCtx->block_number) {
      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);
    }
    
    // --- Group 12
    LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 12: Post-Processing Information.\n");
    // This logic determines the Post Processing configuration and STORES it in simCtx for later reference and cleanup.
    ierr = PetscOptionsGetString(NULL,NULL,"-postprocessing_config_file",simCtx->PostprocessingControlFile,PETSC_MAX_PATH_LEN,NULL); CHKERRQ(ierr);
    /* Parse post settings for both solver and post-processor binaries using the single pre-allocated pps object. */
    ierr = ParsePostProcessingSettings(simCtx); CHKERRQ(ierr);
 
    // === 5. Dependent Parameter Calculations ================================
    // Some parameters depend on others, so we calculate them here.
    simCtx->StartTime = (PetscReal)simCtx->StartStep*simCtx->dt;
    simCtx->ti = simCtx->StartTime;
    simCtx->step = simCtx->StartStep;
 
    // === 5. Log Summary and Finalize Setup ==================================
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "-- Console Output Functions [Total : %d] : --\n", simCtx->nAllowed);
    for (PetscInt i = 0; i < simCtx->nAllowed; ++i) {
      LOG_ALLOW(GLOBAL, LOG_DEBUG, "   [%2d] «%s»\n", i, simCtx->allowedFuncs[i]);
    }
    
    LOG_ALLOW(GLOBAL, LOG_INFO, "Configuration complete. Key parameters:\n");
    LOG_ALLOW(GLOBAL, LOG_INFO, "  - Run mode: %s\n", simCtx->OnlySetup ? "SETUP ONLY" : "Full Simulation");
    LOG_ALLOW(GLOBAL, LOG_INFO, "  - Time steps: %d (from %d to %d)\n", simCtx->StepsToRun, simCtx->StartStep, simCtx->StartStep + simCtx->StepsToRun);
    LOG_ALLOW(GLOBAL, LOG_INFO, "  - Time step size (dt): %g\n", simCtx->dt);
    if (simCtx->tiout > 0) {
      LOG_ALLOW(GLOBAL, LOG_INFO, "  - Field/restart output cadence: every %d step(s)\n", simCtx->tiout);
    } else {
      LOG_ALLOW(GLOBAL, LOG_INFO, "  - Field/restart output cadence: DISABLED\n");
    }
    LOG_ALLOW(GLOBAL, LOG_INFO, "  - Immersed Boundary: %s\n", simCtx->immersed ? "ENABLED" : "DISABLED");
    LOG_ALLOW(GLOBAL, LOG_INFO, "  - Particles: %d\n", simCtx->np);
    if (simCtx->np > 0) {
      if (simCtx->particleConsoleOutputFreq > 0) {
        LOG_ALLOW(GLOBAL, LOG_INFO, "  - Particle console cadence: every %d step(s)\n", simCtx->particleConsoleOutputFreq);
      } else {
        LOG_ALLOW(GLOBAL, LOG_INFO, "  - Particle console cadence: DISABLED\n");
      }
      LOG_ALLOW(GLOBAL, LOG_INFO, "  - Particle console row subsampling: every %d particle(s)\n", simCtx->LoggingFrequency);
    }
    if (simCtx->StartStep > 0 && simCtx->np > 0) {
      LOG_ALLOW(GLOBAL, LOG_INFO, "    - Particle Restart Mode: %s\n", simCtx->particleRestartMode);
    }
    
    // --- Initialize PETSc's internal performance logging stage ---
    ierr = PetscLogDefaultBegin(); CHKERRQ(ierr); // REDUNDANT but safe.
    ierr = PetscMemorySetGetMaximumUsage(); CHKERRQ(ierr);
    
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Finished CreateSimulationContext successfully on rank %d.\n", simCtx->rank);
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

RuntimeWalltimeGuardParsePositiveSeconds()

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

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

This helper is used for shell-exported walltime metadata such as PICURV_JOB_START_EPOCH and PICURV_WALLTIME_LIMIT_SECONDS.

Parameters

[in]	text	String to parse.
[out]	seconds_out	Parsed positive seconds value when successful.

Returns

PetscBool PETSC_TRUE when parsing succeeds, else PETSC_FALSE.

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 67 of file setup.c.

{
    char   *endptr = NULL;
    double parsed_value;
 
    if (seconds_out) *seconds_out = 0.0;
    if (!text || text[0] == '\0') return PETSC_FALSE;
 
    errno        = 0;
    parsed_value = strtod(text, &endptr);
    if (endptr == text || errno == ERANGE || !isfinite(parsed_value) || parsed_value <= 0.0) {
        return PETSC_FALSE;
    }
 
    while (*endptr != '\0' && isspace((unsigned char)*endptr)) {
        endptr++;
    }
    if (*endptr != '\0') return PETSC_FALSE;
 
    if (seconds_out) *seconds_out = (PetscReal)parsed_value;
    return PETSC_TRUE;
}

Here is the caller graph for this function:

SetupSimulationEnvironment()

PetscErrorCode SetupSimulationEnvironment

(

SimCtx *

simCtx

)

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

This function performs a comprehensive series of checks and setup actions to ensure a valid and clean environment. It is parallel-safe; all filesystem operations and checks are performed by Rank 0, with collective error handling.

The function's responsibilities include:

1. Checking Mandatory Inputs: Verifies existence of grid and BCs files.
2. Checking Optional Inputs: Warns if optional config files (whitelist, profile) are missing.
3. Validating Run Mode Paths: Ensures restart_dir or post-processing source directories exist when needed.
4. Preparing Log Directory: Creates the log directory and cleans it of previous logs.
5. Preparing Output Directories: Creates the main output directory and its required subdirectories.

Parameters

[in]

simCtx

The fully configured SimulationContext object.

Returns

PetscErrorCode Returns 0 on success, or a non-zero error code if a mandatory file/directory is missing or a critical operation fails.

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

Local to this translation unit.

Definition at line 1685 of file setup.c.

{
    PetscErrorCode ierr;
    PetscMPIInt    rank;
    PetscBool      exists;
 
    PetscFunctionBeginUser;
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "--- Setting up simulation environment ---\n");
 
    /* =====================================================================
     *  Phase 1: Check for all required and optional INPUT files.
     * ===================================================================== */
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Phase 1: Verifying input files...\n");
 
    // --- Mandatory Inputs ---
    if (!simCtx->generate_grid) {
        ierr = VerifyPathExistence(simCtx->grid_file, PETSC_FALSE, PETSC_FALSE, "Grid file", &exists); CHKERRQ(ierr);
    }
    for (PetscInt i = 0; i < simCtx->num_bcs_files; i++) {
        char desc[128];
        ierr = PetscSNPrintf(desc, sizeof(desc), "BCS file #%d", i + 1); CHKERRQ(ierr);
        ierr = VerifyPathExistence(simCtx->bcs_files[i], PETSC_FALSE, PETSC_FALSE, desc, &exists); CHKERRQ(ierr);
    }
 
    // --- Optional Inputs (these produce warnings if missing) ---
    if (simCtx->useCfg) {
        ierr = VerifyPathExistence(simCtx->allowedFile, PETSC_FALSE, PETSC_TRUE, "Whitelist config file", &exists); CHKERRQ(ierr);
    }
    if (simCtx->useProfilingSelectedFuncsCfg) {
        ierr = VerifyPathExistence(simCtx->profilingSelectedFuncsFile, PETSC_FALSE, PETSC_TRUE, "Profiling config file", &exists); CHKERRQ(ierr);
    }
    if (simCtx->exec_mode == EXEC_MODE_POSTPROCESSOR) {
        ierr = VerifyPathExistence(simCtx->PostprocessingControlFile, PETSC_FALSE, PETSC_TRUE, "Post-processing control file", &exists); CHKERRQ(ierr);
    }
 
 
    /* =====================================================================
     *  Phase 2: Validate directories specific to the execution mode.
     * ===================================================================== */
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Phase 2: Verifying execution mode directories...\n");
    // The data source directory must exist if we intend to load any data from it.
    // This is true if:
    //  1. We are restarting from a previous time step (StartStep > 0), which implies
    //     loading Eulerian fields and/or particle fields.
    //  2. We are starting from t=0 but are explicitly told to load the initial
    //     Eulerian fields from a file (eulerianSource == "load").
    if (simCtx->StartStep > 0 || strcmp(simCtx->eulerianSource,"load")== 0){ // If this is a restart run
        ierr = VerifyPathExistence(simCtx->restart_dir, PETSC_TRUE, PETSC_FALSE, "Restart source directory", &exists); CHKERRQ(ierr);
    }
    if (simCtx->StartStep == 0 && strcmp(simCtx->eulerianSource, "solve") == 0 &&
        simCtx->initialConditionMode == IC_MODE_FILE) {
        ierr = VerifyPathExistence(simCtx->initialConditionDirectory, PETSC_TRUE, PETSC_FALSE,
                                   "Initial-condition source directory", &exists); CHKERRQ(ierr);
    }
    if (simCtx->exec_mode == EXEC_MODE_POSTPROCESSOR) {
        ierr = VerifyPathExistence(simCtx->pps->source_dir, PETSC_TRUE, PETSC_FALSE, "Post-processing source directory", &exists); CHKERRQ(ierr);
    }
 
    /* =====================================================================
     *  Phase 3: Create and prepare all OUTPUT directories.
     * ===================================================================== */
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Phase 3: Preparing output directories...\n");
 
    if (rank == 0){
      if(simCtx->exec_mode == EXEC_MODE_SOLVER){
        // --- Prepare Log Directory ---
        if (!simCtx->continueMode) {
            // Only wipe logs on fresh runs; continue mode appends to existing logs.
            // Final guard: this delete is recursive and irreversible, so re-check
            // containment here rather than trusting that the launcher validated it.
            {
                PetscBool   unsafe_authorized = PETSC_FALSE;
                const char *refusal = NULL;
                ierr = PetscOptionsGetBool(NULL, NULL, "-allow_unsafe_log_dir",
                                           &unsafe_authorized, NULL); CHKERRQ(ierr);
                if (!LogDirectoryIsSafeToWipe(simCtx->log_dir, simCtx->output_dir,
                                              unsafe_authorized, &refusal)) {
                    SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE,
                            "Refusing to delete log directory '%s': it %s. Configure run "
                            "directories through monitor.io.directories.",
                            simCtx->log_dir, refusal);
                }
                if (unsafe_authorized) {
                    LOG_ALLOW(GLOBAL, LOG_WARNING,
                              "Deleting log directory '%s' outside the working directory, "
                              "authorized by -allow_unsafe_log_dir.\n", simCtx->log_dir);
                }
            }
            LOG_ALLOW(GLOBAL, LOG_DEBUG, "Creating/cleaning log directory: %s\n", simCtx->log_dir);
            ierr = PetscRMTree(simCtx->log_dir); // Wipes the directory and its contents
            if (ierr) { /* Ignore file-not-found error, but fail on others */
                PetscError(PETSC_COMM_SELF, __LINE__, __FUNCT__, __FILE__, ierr, PETSC_ERROR_INITIAL, "Could not remove existing log directory '%s'. Check permissions.", simCtx->log_dir);
            }
            ierr = PetscMkdir(simCtx->log_dir); CHKERRQ(ierr);
        } else {
            // In continue mode, ensure log directory exists but don't wipe it.
            LOG_ALLOW(GLOBAL, LOG_DEBUG, "Continue mode: preserving existing log directory: %s\n", simCtx->log_dir);
            ierr = PetscMkdir(simCtx->log_dir); CHKERRQ(ierr);
        }
 
        // --- Prepare Output Directory ---
        // The checkpoint coordinator creates its fixed internal hierarchy
        // transactionally when a checkpoint is due.
        LOG_ALLOW(GLOBAL, LOG_DEBUG, "Verifying main output directory: %s\n", simCtx->output_dir);
        ierr = PetscTestDirectory(simCtx->output_dir, 'r', &exists); CHKERRQ(ierr);
        if (!exists) {
            LOG_ALLOW(GLOBAL, LOG_INFO, "Output directory not found. Creating: %s\n", simCtx->output_dir);
            ierr = PetscMkdir(simCtx->output_dir); CHKERRQ(ierr);
        }
        // The analysis home is nested (output/analysis/metrics by default), so it
        // needs the recursive create the postprocessor branch already uses; PetscMkdir
        // fails on a path whose parents do not exist yet.
        ierr = PetscMkdirRecursive(simCtx->analysis_dir); CHKERRQ(ierr);
      } else if(simCtx->exec_mode == EXEC_MODE_POSTPROCESSOR){
        LOG_ALLOW(GLOBAL, LOG_DEBUG, "Preparing post-processing output directories ...\n");
 
        PostProcessParams *pps = simCtx->pps;
        char path_buffer[PETSC_MAX_PATH_LEN];
 
        const char *last_slash_euler = strrchr(pps->output_prefix, '/');
        if(last_slash_euler){
          size_t dir_len = last_slash_euler - pps->output_prefix;
          if(dir_len > 0){
            if(dir_len >= sizeof(path_buffer)) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Post-processing output prefix path is too long.");
            strncpy(path_buffer, pps->output_prefix, dir_len);
            path_buffer[dir_len] = '\0';
 
            ierr = PetscTestDirectory(path_buffer, 'r', &exists); CHKERRQ(ierr);
            if (!exists){
              LOG_ALLOW(GLOBAL, LOG_INFO, "Creating post-processing Eulerian output directory: %s\n", path_buffer);
              ierr = PetscMkdirRecursive(path_buffer); CHKERRQ(ierr);
            }
          }
        }
 
        // Particle output directory
        if(pps->outputParticles){
          const char  *last_slash_particle = strrchr(pps->particle_output_prefix, '/');
          if(last_slash_particle){
            size_t dir_len = last_slash_particle - pps->particle_output_prefix;
            if(dir_len >  0){
              if(dir_len > sizeof(path_buffer)) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Post-processing particle output prefix path is too long.");
              strncpy(path_buffer, pps->particle_output_prefix, dir_len);
              path_buffer[dir_len] = '\0';
              
              ierr = PetscTestDirectory(path_buffer, 'r', &exists); CHKERRQ(ierr);
 
              if (!exists){
                LOG_ALLOW(GLOBAL, LOG_INFO, "Creating post-processing Particle output directory: %s\n", path_buffer);
                ierr = PetscMkdirRecursive(path_buffer); CHKERRQ(ierr);
              }
            }
          }
        }
 
        // Statistics output directory
        if(pps->statistics_pipeline[0] != '\0'){
          const char *last_slash_stats = strrchr(pps->statistics_output_prefix, '/');
          if(last_slash_stats){
            size_t dir_len = last_slash_stats - pps->statistics_output_prefix;
            if(dir_len > 0){
              if(dir_len >= sizeof(path_buffer)) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Post-processing statistics output prefix path is too long.");
              strncpy(path_buffer, pps->statistics_output_prefix, dir_len);
              path_buffer[dir_len] = '\0';
 
              ierr = PetscTestDirectory(path_buffer, 'r', &exists); CHKERRQ(ierr);
              if (!exists){
                LOG_ALLOW(GLOBAL, LOG_INFO, "Creating post-processing Statistics output directory: %s\n", path_buffer);
                ierr = PetscMkdirRecursive(path_buffer); CHKERRQ(ierr);
              }
            }
          }
        }
      }
    }
    
    // Synchronize all processes before proceeding
    ierr = MPI_Barrier(PETSC_COMM_WORLD); CHKERRMPI(ierr);
 
    /* PETSc opens a -info file during PetscInitialize, before the fresh-run wipe above
       removes the log directory it usually lives in. Every later record would then go
       to an unlinked file and the run would end with no info log at all, so each rank
       reopens its file once the directory exists again. Records from initialization
       itself are the only ones lost. */
    if (simCtx->exec_mode == EXEC_MODE_SOLVER && !simCtx->continueMode) {
        char *info_name = NULL;
        FILE *info_file = NULL;
        ierr = PetscInfoGetFile(&info_name, &info_file); CHKERRQ(ierr);
        if (info_name && info_name[0] != '\0') {
            char reopen_name[PETSC_MAX_PATH_LEN];
            ierr = PetscStrncpy(reopen_name, info_name, sizeof(reopen_name)); CHKERRQ(ierr);
            if (info_file && info_file != PETSC_STDOUT) {
                ierr = PetscFClose(PETSC_COMM_SELF, info_file); CHKERRQ(ierr);
            }
            ierr = PetscInfoSetFile(reopen_name, "w"); CHKERRQ(ierr);
        }
        /* PetscInfoGetFile hands back its own copy of the name. */
        ierr = PetscFree(info_name); CHKERRQ(ierr);
    }
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "--- Environment setup complete ---\n");
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

SetupGridAndSolvers()

PetscErrorCode SetupGridAndSolvers

(

SimCtx *

simCtx

)

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

This function is the high-level driver for creating the entire computational domain, including the multigrid hierarchy, PETSc DMDA and Vec objects, and calculating all necessary grid metrics.

Parameters

simCtx

The fully configured SimulationContext.

Returns

PetscErrorCode

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 2001 of file setup.c.

{
    PetscErrorCode ierr;
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "--- Starting Grid and Solvers Setup ---\n");
 
    // Phase 1: Allocate the UserMG and UserCtx hierarchy
    ierr = AllocateContextHierarchy(simCtx); CHKERRQ(ierr);
 
    ierr = DefineAllGridDimensions(simCtx); CHKERRQ(ierr);
    ierr = InitializeAllGridDMs(simCtx); CHKERRQ(ierr);
    ierr = AssignAllGridCoordinates(simCtx); CHKERRQ(ierr);
    ierr = CreateAndInitializeAllVectors(simCtx); CHKERRQ(ierr);
 
    // NOTE: CalculateAllGridMetrics is now called inside SetupBoundaryConditions (not here) to ensure:
    // 1. Boundary condition configuration data (boundary_faces) is available for periodic BC corrections
    // 2. Computed metrics are available for inlet/outlet area calculations
    // This resolves the circular dependency between BC setup and metric calculations.
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "--- Grid and Solvers Setup Complete ---\n");
    
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

CreateAndInitializeAllVectors()

PetscErrorCode CreateAndInitializeAllVectors

(

SimCtx *

simCtx

)

Creates and initializes all PETSc Vec objects for all fields.

This function iterates through every UserCtx in the multigrid and multi-block hierarchy. For each context, it creates the comprehensive set of global and local PETSc Vecs required by the flow solver (e.g., Ucont, P, Nvert, metrics, turbulence fields, etc.). Each vector is initialized to zero.

Parameters

simCtx

The master SimCtx, containing the configured UserCtx hierarchy.

Returns

PetscErrorCode

Creates and initializes all PETSc Vec objects for all fields.

Local to this translation unit.

Definition at line 2036 of file setup.c.

{
    PetscErrorCode ierr;
    UserMG         *usermg = &simCtx->usermg;
    MGCtx          *mgctx = usermg->mgctx;
    PetscInt       nblk = simCtx->block_number;
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
    
    LOG_ALLOW(GLOBAL, LOG_INFO, "Creating and initializing all simulation vectors...\n");
        
    for (PetscInt level  = usermg->mglevels-1; level >=0; level--) {
        for (PetscInt bi = 0; bi < nblk; bi++) {
            UserCtx *user = &mgctx[level].user[bi];
 
            if(!user->da || !user->fda) {
              SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "DMs not properly initialized in UserCtx before vector creation.");
            }
 
            LOG_ALLOW_SYNC(LOCAL, LOG_DEBUG, "Rank %d: Creating vectors for level %d, block %d\n", simCtx->rank, level, bi);
 
        // --- Group A: Primary Flow Fields (Global and Local) ---
            // These are the core solution variables.
            ierr = DMCreateGlobalVector(user->fda, &user->Ucont); CHKERRQ(ierr); ierr = VecSet(user->Ucont, 0.0); CHKERRQ(ierr);
            ierr = DMCreateGlobalVector(user->fda, &user->Ucat);  CHKERRQ(ierr); ierr = VecSet(user->Ucat, 0.0); CHKERRQ(ierr);
            ierr = DMCreateGlobalVector(user->da,  &user->P);     CHKERRQ(ierr); ierr = VecSet(user->P, 0.0); CHKERRQ(ierr);
            ierr = DMCreateGlobalVector(user->da,  &user->Nvert); CHKERRQ(ierr); ierr = VecSet(user->Nvert, 0.0); CHKERRQ(ierr);
 
            ierr = DMCreateLocalVector(user->fda, &user->lUcont); CHKERRQ(ierr); ierr = VecSet(user->lUcont, 0.0); CHKERRQ(ierr);
            ierr = DMCreateLocalVector(user->fda, &user->lUcat);  CHKERRQ(ierr); ierr = VecSet(user->lUcat, 0.0); CHKERRQ(ierr);
            ierr = DMCreateLocalVector(user->da,  &user->lP);     CHKERRQ(ierr); ierr = VecSet(user->lP, 0.0); CHKERRQ(ierr);
            ierr = DMCreateLocalVector(user->da,  &user->lNvert); CHKERRQ(ierr); ierr = VecSet(user->lNvert, 0.0); CHKERRQ(ierr);
 
            // -- Group A2: Derived Flow Fields (Global and Local) ---
            ierr = VecDuplicate(user->P,&user->Diffusivity); CHKERRQ(ierr); ierr = VecSet(user->Diffusivity, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->lP,&user->lDiffusivity); CHKERRQ(ierr); ierr = VecSet(user->lDiffusivity, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->Ucat,&user->DiffusivityGradient); CHKERRQ(ierr); ierr = VecSet(user->DiffusivityGradient, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->lUcat,&user->lDiffusivityGradient); CHKERRQ(ierr); ierr = VecSet(user->lDiffusivityGradient, 0.0); CHKERRQ(ierr);
 
            // -- Group B: Solver Work Vectors (Global and Local) ---
            ierr = VecDuplicate(user->P, &user->Phi);       CHKERRQ(ierr); ierr = VecSet(user->Phi, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->lP, &user->lPhi);       CHKERRQ(ierr); ierr = VecSet(user->lPhi, 0.0); CHKERRQ(ierr);
 
            // --- Group C: Time-Stepping & Workspace Fields (Finest Level Only) ---
            if (level == usermg->mglevels - 1) {
                ierr = VecDuplicate(user->Ucont, &user->Ucont_o);   CHKERRQ(ierr); ierr = VecSet(user->Ucont_o, 0.0); CHKERRQ(ierr);
                ierr = VecDuplicate(user->Ucont, &user->Ucont_rm1); CHKERRQ(ierr); ierr = VecSet(user->Ucont_rm1, 0.0); CHKERRQ(ierr);
                ierr = VecDuplicate(user->Ucat,  &user->Ucat_o);    CHKERRQ(ierr); ierr = VecSet(user->Ucat_o, 0.0); CHKERRQ(ierr);
                ierr = VecDuplicate(user->P,     &user->P_o);       CHKERRQ(ierr); ierr = VecSet(user->P_o, 0.0); CHKERRQ(ierr);
                  ierr = VecDuplicate(user->lUcont, &user->lUcont_o);   CHKERRQ(ierr); ierr = VecSet(user->lUcont_o, 0.0); CHKERRQ(ierr);
                ierr = VecDuplicate(user->lUcont, &user->lUcont_rm1); CHKERRQ(ierr); ierr = VecSet(user->lUcont_rm1, 0.0); CHKERRQ(ierr);
                ierr = DMCreateLocalVector(user->da, &user->lNvert_o); CHKERRQ(ierr); ierr = VecSet(user->lNvert_o, 0.0); CHKERRQ(ierr);
                    ierr = VecDuplicate(user->Nvert, &user->Nvert_o); CHKERRQ(ierr); ierr = VecSet(user->Nvert_o, 0.0); CHKERRQ(ierr);
              }
 
        // --- Group D: Grid Metrics (Face-Centered) ---
            ierr = DMCreateGlobalVector(user->fda, &user->Csi); CHKERRQ(ierr); ierr = VecSet(user->Csi, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->Csi, &user->Eta);         CHKERRQ(ierr); ierr = VecSet(user->Eta, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->Csi, &user->Zet);         CHKERRQ(ierr); ierr = VecSet(user->Zet, 0.0); CHKERRQ(ierr);
            ierr = DMCreateGlobalVector(user->da,  &user->Aj);  CHKERRQ(ierr); ierr = VecSet(user->Aj, 0.0); CHKERRQ(ierr);
            
            ierr = DMCreateLocalVector(user->fda, &user->lCsi); CHKERRQ(ierr); ierr = VecSet(user->lCsi, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->lCsi, &user->lEta);       CHKERRQ(ierr); ierr = VecSet(user->lEta, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->lCsi, &user->lZet);       CHKERRQ(ierr); ierr = VecSet(user->lZet, 0.0); CHKERRQ(ierr);
            ierr = DMCreateLocalVector(user->da,  &user->lAj);  CHKERRQ(ierr); ierr = VecSet(user->lAj, 0.0); CHKERRQ(ierr);
         
        
       // --- Group E: Grid Metrics (Face-Centered) ---
            // Vector metrics are duplicated from Csi (DOF=3, fda-based)
            ierr = VecDuplicate(user->Csi, &user->ICsi); CHKERRQ(ierr); ierr = VecSet(user->ICsi, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->Csi, &user->IEta); CHKERRQ(ierr); ierr = VecSet(user->IEta, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->Csi, &user->IZet); CHKERRQ(ierr); ierr = VecSet(user->IZet, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->Csi, &user->JCsi); CHKERRQ(ierr); ierr = VecSet(user->JCsi, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->Csi, &user->JEta); CHKERRQ(ierr); ierr = VecSet(user->JEta, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->Csi, &user->JZet); CHKERRQ(ierr); ierr = VecSet(user->JZet, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->Csi, &user->KCsi); CHKERRQ(ierr); ierr = VecSet(user->KCsi, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->Csi, &user->KEta); CHKERRQ(ierr); ierr = VecSet(user->KEta, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->Csi, &user->KZet); CHKERRQ(ierr); ierr = VecSet(user->KZet, 0.0); CHKERRQ(ierr);
            // Scalar metrics are duplicated from Aj (DOF=1, da-based)
            ierr = VecDuplicate(user->Aj, &user->IAj); CHKERRQ(ierr); ierr = VecSet(user->IAj, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->Aj, &user->JAj); CHKERRQ(ierr); ierr = VecSet(user->JAj, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->Aj, &user->KAj); CHKERRQ(ierr); ierr = VecSet(user->KAj, 0.0); CHKERRQ(ierr);
 
            ierr = VecDuplicate(user->lCsi, &user->lICsi); CHKERRQ(ierr); ierr = VecSet(user->lICsi, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->lCsi, &user->lIEta); CHKERRQ(ierr); ierr = VecSet(user->lIEta, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->lCsi, &user->lIZet); CHKERRQ(ierr); ierr = VecSet(user->lIZet, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->lCsi, &user->lJCsi); CHKERRQ(ierr); ierr = VecSet(user->lJCsi, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->lCsi, &user->lJEta); CHKERRQ(ierr); ierr = VecSet(user->lJEta, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->lCsi, &user->lJZet); CHKERRQ(ierr); ierr = VecSet(user->lJZet, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->lCsi, &user->lKCsi); CHKERRQ(ierr); ierr = VecSet(user->lKCsi, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->lCsi, &user->lKEta); CHKERRQ(ierr); ierr = VecSet(user->lKEta, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->lCsi, &user->lKZet); CHKERRQ(ierr); ierr = VecSet(user->lKZet, 0.0); CHKERRQ(ierr);
 
            ierr = VecDuplicate(user->lAj, &user->lIAj); CHKERRQ(ierr); ierr = VecSet(user->lIAj, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->lAj, &user->lJAj); CHKERRQ(ierr); ierr = VecSet(user->lJAj, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->lAj, &user->lKAj); CHKERRQ(ierr); ierr = VecSet(user->lKAj, 0.0); CHKERRQ(ierr);
        
        // --- Group F: Cell/Face Center Coordinates and Grid Spacing ---
            ierr = DMCreateGlobalVector(user->fda, &user->Cent); CHKERRQ(ierr); ierr = VecSet(user->Cent, 0.0); CHKERRQ(ierr);
            ierr = DMCreateLocalVector(user->fda, &user->lCent); CHKERRQ(ierr); ierr = VecSet(user->lCent, 0.0); CHKERRQ(ierr);
            
            ierr = VecDuplicate(user->Cent, &user->GridSpace); CHKERRQ(ierr); ierr = VecSet(user->GridSpace, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->lCent, &user->lGridSpace); CHKERRQ(ierr); ierr = VecSet(user->lGridSpace, 0.0); CHKERRQ(ierr);
 
            ierr = VecDuplicate(user->Cent, &user->Centx); CHKERRQ(ierr); ierr = VecSet(user->Centx, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->Cent, &user->Centy); CHKERRQ(ierr); ierr = VecSet(user->Centy, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->Cent, &user->Centz); CHKERRQ(ierr); ierr = VecSet(user->Centz, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->lCent, &user->lCentx); CHKERRQ(ierr); ierr = VecSet(user->lCentx, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->lCent, &user->lCenty); CHKERRQ(ierr); ierr = VecSet(user->lCenty, 0.0); CHKERRQ(ierr);
            ierr = VecDuplicate(user->lCent, &user->lCentz); CHKERRQ(ierr); ierr = VecSet(user->lCentz, 0.0); CHKERRQ(ierr);
 
        if(level == usermg->mglevels -1){
        // --- Group G: Turbulence Models (Finest Level Only) ---
            if (simCtx->les) {
                ierr = DMCreateGlobalVector(user->da, &user->Nu_t); CHKERRQ(ierr); ierr = VecSet(user->Nu_t, 0.0); CHKERRQ(ierr);
                ierr = DMCreateLocalVector(user->da, &user->lNu_t); CHKERRQ(ierr); ierr = VecSet(user->lNu_t, 0.0); CHKERRQ(ierr);
                LOG_ALLOW(GLOBAL, LOG_DEBUG, "Turbulence viscosity (Nu_t) vectors created for the LES model.\n");
                // Only the dynamic model carries a coefficient field. The constant model
                // reads its coefficient straight from configuration, so allocating,
                // synchronizing, and checkpointing a field of one repeated number would
                // buy nothing.
                if(simCtx->les == DYNAMIC_SMAGORINSKY){
                ierr = DMCreateGlobalVector(user->da,&user->CS); CHKERRQ(ierr); ierr = VecSet(user->CS,0.0); CHKERRQ(ierr);
                ierr = DMCreateLocalVector(user->da,&user->lCs); CHKERRQ(ierr); ierr = VecSet(user->lCs,0.0); CHKERRQ(ierr);
                LOG_ALLOW(GLOBAL, LOG_DEBUG, "Dynamic Smagorinsky coefficient (CS) vectors created.\n");
                }
 
                  // Add K_Omega etc. here as needed
 
                // Note: Add any other vectors from the legacy MG_Initial here as needed.
                // For example: Rhs, Forcing, turbulence Vecs (K_Omega, Nu_t)...
        
            }
 
            /* A wall model is configured independently of LES and RANS - `wall_function`
               is their sibling in the schema, not their child - so its storage cannot sit
               inside their block. It did, which left `ApplyWallFunction` dereferencing a
               null vector on any case that enabled a wall model without one of them. */
            if (simCtx->wallfunction) {
                /* Scalar per cell, so it lives on `da`. It was previously created from
                   `fda`, whose dof is 3, while every reader opens it through `da` - a
                   pairing PETSc rejects outright. */
                ierr = DMCreateGlobalVector(user->da, &user->Friction_Velocity); CHKERRQ(ierr);
                ierr = VecSet(user->Friction_Velocity, 0.0); CHKERRQ(ierr);
                ierr = DMCreateLocalVector(user->da, &user->lFriction_Velocity); CHKERRQ(ierr);
                ierr = VecSet(user->lFriction_Velocity, 0.0); CHKERRQ(ierr);
                ierr = DMCreateGlobalVector(user->da, &user->Nu_Wall); CHKERRQ(ierr);
                ierr = VecSet(user->Nu_Wall, 0.0); CHKERRQ(ierr);
                ierr = DMCreateLocalVector(user->da, &user->lNu_Wall); CHKERRQ(ierr);
                ierr = VecSet(user->lNu_Wall, 0.0); CHKERRQ(ierr);
            }
        // --- Group H: Particle Methods    
        if(simCtx->np>0){
          ierr = DMCreateGlobalVector(user->da,&user->ParticleCount); CHKERRQ(ierr); ierr = VecSet(user->ParticleCount,0.0); CHKERRQ(ierr);
        ierr = DMCreateLocalVector(user->da,&user->lParticleCount); CHKERRQ(ierr); ierr = VecSet(user->lParticleCount,0.0); CHKERRQ(ierr);
        // Scalar field to hold particle scalar property (e.g., temperature, concentration)
          ierr = DMCreateGlobalVector(user->da,&user->Psi); CHKERRQ(ierr); ierr = VecSet(user->Psi,0.0); CHKERRQ(ierr);
        ierr = DMCreateLocalVector(user->da,&user->lPsi); CHKERRQ(ierr); ierr = VecSet(user->lPsi,0.0); CHKERRQ(ierr);
          LOG_ALLOW(GLOBAL,LOG_DEBUG,"ParticleCount & Scalar(Psi) created for %d particles.\n",simCtx->np);
          }
        }
        // --- Group I: Boundary Condition vectors ---
        ierr = DMCreateGlobalVector(user->fda, &user->Bcs.Ubcs); CHKERRQ(ierr);
        ierr = VecSet(user->Bcs.Ubcs, 0.0); CHKERRQ(ierr);
        ierr = DMCreateGlobalVector(user->fda, &user->Bcs.Uch); CHKERRQ(ierr);
        ierr = VecSet(user->Bcs.Uch, 0.0); CHKERRQ(ierr);
        
        // --- Group L: Field-Statistics Accumulators (Finest Level Only) ---
        // Config-counted, like the convergence-state reference fields: the window
        // count is resolved before this factory runs, and each accumulator is
        // duplicated from a vector created above so it inherits DM and layout.
        if (level == usermg->mglevels - 1 && FieldStatisticsIsActive(simCtx)) {
            ierr = PetscCalloc1((size_t)simCtx->fieldStatisticsWindowCount,
                                &user->fieldStatisticsStorage); CHKERRQ(ierr);
            for (PetscInt w = 0; w < simCtx->fieldStatisticsWindowCount; ++w) {
                ierr = PicurvWindowStorageCreate(user,
                           &simCtx->fieldStatisticsWindows[w].definition,
                           &user->fieldStatisticsStorage[w]); CHKERRQ(ierr);
            }
            LOG_ALLOW(LOCAL, LOG_INFO,
                      "Allocated accumulators for %d statistics window(s).\n",
                      simCtx->fieldStatisticsWindowCount);
        }
 
        // --- Group K: Corner-Staging Workspace (Finest Level Only) ---
        // One pair per block size. Previously these were created lazily inside
        // the interpolation routine and rebuilt whenever the block size changed;
        // allocating both here removes that churn and the cached-size check.
        if (level == usermg->mglevels - 1) {
            ierr = DMCreateGlobalVector(user->da,  &user->CellScalarAtCorner);  CHKERRQ(ierr);
            ierr = VecSet(user->CellScalarAtCorner, 0.0); CHKERRQ(ierr);
            ierr = DMCreateLocalVector(user->da,   &user->lCellScalarAtCorner); CHKERRQ(ierr);
            ierr = VecSet(user->lCellScalarAtCorner, 0.0); CHKERRQ(ierr);
            ierr = DMCreateGlobalVector(user->fda, &user->CellVectorAtCorner);  CHKERRQ(ierr);
            ierr = VecSet(user->CellVectorAtCorner, 0.0); CHKERRQ(ierr);
            ierr = DMCreateLocalVector(user->fda,  &user->lCellVectorAtCorner); CHKERRQ(ierr);
            ierr = VecSet(user->lCellVectorAtCorner, 0.0); CHKERRQ(ierr);
        }
 
        /* --- Group N: Statistics Output Staging (Finest Level Only) ---
         * The same kind of object as the corner workspace: a named buffer the
         * shared kernels can address, rather than simulation state. It exists only
         * to carry derived statistics to the nodal and VTK paths, so a run without
         * statistics allocates none of it. */
        if (level == usermg->mglevels - 1 && FieldStatisticsIsActive(simCtx)) {
            ierr = DMCreateGlobalVector(user->da,  &user->PostScalar);       CHKERRQ(ierr);
            ierr = VecSet(user->PostScalar, 0.0); CHKERRQ(ierr);
            ierr = DMCreateLocalVector(user->da,   &user->lPostScalar);      CHKERRQ(ierr);
            ierr = VecSet(user->lPostScalar, 0.0); CHKERRQ(ierr);
            ierr = DMCreateGlobalVector(user->da,  &user->PostScalarNodal);  CHKERRQ(ierr);
            ierr = VecSet(user->PostScalarNodal, 0.0); CHKERRQ(ierr);
            ierr = DMCreateGlobalVector(user->fda, &user->PostVector);       CHKERRQ(ierr);
            ierr = VecSet(user->PostVector, 0.0); CHKERRQ(ierr);
            ierr = DMCreateLocalVector(user->fda,  &user->lPostVector);      CHKERRQ(ierr);
            ierr = VecSet(user->lPostVector, 0.0); CHKERRQ(ierr);
            ierr = DMCreateGlobalVector(user->fda, &user->PostVectorNodal);  CHKERRQ(ierr);
            ierr = VecSet(user->PostVectorNodal, 0.0); CHKERRQ(ierr);
        }
 
      if(level == usermg->mglevels - 1){
        if(simCtx->exec_mode == EXEC_MODE_POSTPROCESSOR){
                LOG_ALLOW(LOCAL, LOG_DEBUG, "Post-processor mode detected. Allocating derived field vectors.\n");
 
                ierr = VecDuplicate(user->P, &user->P_nodal); CHKERRQ(ierr);
                ierr = VecSet(user->P_nodal, 0.0); CHKERRQ(ierr);
 
                ierr = VecDuplicate(user->Ucat, &user->Ucat_nodal); CHKERRQ(ierr);
                ierr = VecSet(user->Ucat_nodal, 0.0); CHKERRQ(ierr);
                
                ierr = VecDuplicate(user->P, &user->Qcrit); CHKERRQ(ierr);
                ierr = VecSet(user->Qcrit, 0.0); CHKERRQ(ierr);
                ierr = DMCreateLocalVector(user->da, &user->lQcrit); CHKERRQ(ierr);
                ierr = VecSet(user->lQcrit, 0.0); CHKERRQ(ierr);
                ierr = VecDuplicate(user->P, &user->Qcrit_nodal); CHKERRQ(ierr);
                ierr = VecSet(user->Qcrit_nodal, 0.0); CHKERRQ(ierr);
 
                LOG_ALLOW(LOCAL, LOG_DEBUG, "Derived field vectors P_nodal, Ucat_nodal, and Qcrit created.\n");
 
                if(simCtx->np>0){
                ierr = VecDuplicate(user->Psi, &user->Psi_nodal); CHKERRQ(ierr);
                ierr = VecSet(user->Psi_nodal, 0.0); CHKERRQ(ierr);
 
                LOG_ALLOW(LOCAL, LOG_DEBUG, "Derived field vector Psi_nodal created for particle scalar property.\n");
 
                }
        }else{
                user->P_nodal = NULL;
                user->Ucat_nodal = NULL;
                user->Qcrit = NULL;
                user->lQcrit = NULL;
                user->Qcrit_nodal = NULL;
                user->Psi_nodal = NULL;
        }
      }
    
  }
}
 
    /* Config-counted vectors belong to this factory too. The convergence state
     * allocates a Vec array whose length comes from configuration, so it is
     * invoked here rather than by the caller, keeping every vector this run
     * owns created in one place. It depends only on option-parsed configuration
     * and on the vectors created above. */
    ierr = InitializeSolutionConvergenceState(simCtx); CHKERRQ(ierr);
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "All simulation vectors created and initialized.\n");
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

UpdateLocalGhosts()

PetscErrorCode UpdateLocalGhosts	(	UserCtx *	user,
		FieldId	field_id
	)

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

This function resolves the correct global vector, local vector, and DM from the authoritative field catalog and performs the standard PETSc DMGlobalToLocalBegin/End sequence. For registered single-face-family and component-staggered fields, it then repairs the adjacent periodic ghost in each normal direction so central differences see the correct neighboring face. Cell-centered fields and tangential face directions retain PETSc's native periodic wraparound. Includes optional debugging output (max norms before/after).

Parameters

user	The UserCtx structure containing the vectors and DMs.
field_id	Typed identity of the catalogued field to update.

Returns

PetscErrorCode 0 on success, non-zero on failure.

Note

This function assumes the global vector associated with field_id has already been populated with the desired data (including any boundary conditions and periodic duplicate planes). It scatters that current global state into the local vector; it does not repair stale global periodic duplicate planes.

Wider QUICK-style normal ghost layers require a separate dedicated exchange; this function repairs the adjacent layer used by central stencils.

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

Resolves the typed field view, performs the existing PETSc scatter, and applies the established periodic normal-face repair when needed.

Definition at line 2489 of file setup.c.

{
    PetscErrorCode ierr;
    PetscMPIInt    rank;
    FieldView      field_view;
    const char    *field_name;
    Vec            globalVec;
    Vec            localVec;
    DM             dm;
    PetscInt       dof;
    char           face_direction = '\0';
    PetscBool      component_staggered = PETSC_FALSE;
 
    PetscFunctionBeginUser; // Use User version for application code
    PROFILE_FUNCTION_BEGIN;
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
    ierr = FieldGetView(user, field_id, &field_view); CHKERRQ(ierr);
    PetscCheck((field_view.descriptor->capabilities & FIELD_CAPABILITY_GHOST_UPDATE) != 0u,
               PETSC_COMM_SELF, PETSC_ERR_SUP,
               "Field '%s' does not support ghost updates.",
               field_view.descriptor->canonical_name);
    field_name = field_view.descriptor->canonical_name;
    globalVec = field_view.global_vec;
    localVec = field_view.local_vec;
    dm = field_view.dm;
    dof = field_view.descriptor->dof;
 
    switch (field_view.descriptor->sync_class) {
        case FIELD_SYNC_I_FACE:
            face_direction = 'i';
            break;
        case FIELD_SYNC_J_FACE:
            face_direction = 'j';
            break;
        case FIELD_SYNC_K_FACE:
            face_direction = 'k';
            break;
        case FIELD_SYNC_COMPONENT_STAGGERED:
            component_staggered = PETSC_TRUE;
            break;
        case FIELD_SYNC_STANDARD:
            break;
        default:
            SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB,
                    "Field '%s' has an invalid ghost synchronization class.", field_name);
    }
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Rank %d: Starting ghost update for field '%s'.\n", rank, field_name);
 
    LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Identified components for '%s': DM=%p, GlobalVec=%p, LocalVec=%p.\n",
              rank, field_name, (void*)dm, (void*)globalVec, (void*)localVec);
 
    // --- 3. Optional Debugging: Norm Before Update ---
    // Use your logging convention check
    // if (get_log_level() >= LOG_LEVEL_DEBUG && is_function_allowed("UpdateLocalGhosts")) { // Example check
    if(get_log_level() == LOG_DEBUG && is_function_allowed(__func__)){
        PetscReal norm_global_before;
        ierr = VecNorm(globalVec, NORM_INFINITY, &norm_global_before); CHKERRQ(ierr);
        LOG_ALLOW(GLOBAL, LOG_INFO,"Max norm '%s' (Global) BEFORE Ghost Update: %g\n", field_name, norm_global_before);
        // Optional: Norm of local vector before update (might contain old ghost values)
        // PetscReal norm_local_before;
        // ierr = VecNorm(localVec, NORM_INFINITY, &norm_local_before); CHKERRQ(ierr);
        // LOG_ALLOW(GLOBAL, LOG_DEBUG,"Max norm '%s' (Local) BEFORE Ghost Update: %g\n", field_name, norm_local_before);
    }
 
    // --- 4. Perform the Global-to-Local Transfer (Ghost Update) ---
    LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Calling DMGlobalToLocalBegin/End for '%s'.\n", rank, field_name);
    ierr = DMGlobalToLocalBegin(dm, globalVec, INSERT_VALUES, localVec); CHKERRQ(ierr);
    ierr = DMGlobalToLocalEnd(dm, globalVec, INSERT_VALUES, localVec); CHKERRQ(ierr);
    ierr = RepairPeriodicNormalFaceGhosts(user, dm, localVec, dof, face_direction,
                                          component_staggered); CHKERRQ(ierr);
    LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Completed DMGlobalToLocalBegin/End for '%s'.\n", rank, field_name);
 
    // --- 5. Optional Debugging: Norm After Update ---
    // Use your logging convention check
    // if (get_log_level() >= LOG_LEVEL_DEBUG && is_function_allowed("UpdateLocalGhosts")) { // Example check
    if(get_log_level() == LOG_DEBUG && is_function_allowed(__func__)){ // Using your specific check
        PetscReal norm_local_after;
        ierr = VecNorm(localVec, NORM_INFINITY, &norm_local_after); CHKERRQ(ierr);
        LOG_ALLOW(GLOBAL, LOG_INFO,"Max norm '%s' (Local) AFTER Ghost Update: %g\n", field_name, norm_local_after);
 
        // --- 6. Optional Debugging: Specific Point Checks (Example for Ucat on Rank 0/1) ---
        //    (Keep this conditional if it's only for specific debug scenarios)
        if (field_id == FIELD_ID_UCAT) { // Only do detailed checks for Ucat for now
           PetscMPIInt rank_test;
           MPI_Comm_rank(PETSC_COMM_WORLD, &rank_test);
 
           // Get Local Info needed for indexing checks
           DMDALocalInfo info_check;
           ierr = DMDAGetLocalInfo(dm, &info_check); CHKERRQ(ierr); // Use the correct dm
 
           // Buffer for array pointer
           Cmpnts ***lUcat_arr_test = NULL;
           PetscErrorCode ierr_test = 0;
 
           LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Testing '%s' access immediately after ghost update...\n", rank_test, field_name);
           ierr_test = DMDAVecGetArrayDOFRead(dm, localVec, &lUcat_arr_test); // Use correct dm and localVec
 
           if (ierr_test) {
               LOG_ALLOW(LOCAL, LOG_ERROR, "Rank %d: ERROR %d getting '%s' array after ghost update!\n", rank_test, ierr_test, field_name);
           } else if (!lUcat_arr_test) {
                LOG_ALLOW(LOCAL, LOG_ERROR, "Rank %d: ERROR NULL pointer getting '%s' array after ghost update!\n", rank_test, field_name);
           }
           else {
               // Check owned interior point (e.g., first interior point)
               PetscInt k_int = info_check.zs + (info_check.zm > 1 ? 1 : 0); // Global k index (at least zs+1 if possible)
               PetscInt j_int = info_check.ys + (info_check.ym > 1 ? 1 : 0); // Global j index
               PetscInt i_int = info_check.xs + (info_check.xm > 1 ? 1 : 0); // Global i index
                // Ensure indices are within global bounds if domain is very small
               //if (k_int >= info_check.mz-1) k_int = info_check.mz-2; if (k_int < 1) k_int = 1;
               //if (j_int >= info_check.my-1) j_int = info_check.my-2; if (j_int < 1) j_int = 1;
           // if (i_int >= info_check.mx-1) i_int = info_check.mx-2; if (i_int < 1) i_int = 1;
           // clamp k_int to [1 .. mz-2] 
           if (k_int >= info_check.mz - 1) {
         k_int = info_check.mz - 2;
           }
           if (k_int < 1) {
         k_int = 1;
           }
 
           // clamp j_int to [1 .. my-2] 
           if (j_int >= info_check.my - 1) {
         j_int = info_check.my - 2;
           }
           if (j_int < 1) {
         j_int = 1;
           }
 
           // clamp i_int to [1 .. mx-2]
           if (i_int >= info_check.mx - 1) {
         i_int = info_check.mx - 2;
           }
           if (i_int < 1) {
         i_int = 1;
           }
 
               // Only attempt read if indices are actually owned (relevant for multi-rank)
               if (k_int >= info_check.zs && k_int < info_check.zs + info_check.zm &&
                   j_int >= info_check.ys && j_int < info_check.ys + info_check.ym &&
                   i_int >= info_check.xs && i_int < info_check.xs + info_check.xm)
               {
                  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);
                  Cmpnts test_val_owned_interior = lUcat_arr_test[k_int][j_int][i_int];
                  LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: SUCCESS reading owned interior: x=%g\n", rank_test, test_val_owned_interior.x);
               } else {
                  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);
               }
 
 
               // Check owned boundary point (e.g., first owned point)
               PetscInt k_bnd = info_check.zs; // Global k index
               PetscInt j_bnd = info_check.ys; // Global j index
               PetscInt i_bnd = info_check.xs; // Global i index
               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);
               Cmpnts test_val_owned_boundary = lUcat_arr_test[k_bnd][j_bnd][i_bnd];
               LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: SUCCESS reading owned boundary: x=%g\n", rank_test, test_val_owned_boundary.x);
 
 
               // Check ghost point (e.g., one layer below in k, if applicable)
               if (info_check.zs > 0) { // Only if there's a rank below
                   PetscInt k_ghost = info_check.zs - 1;
                   PetscInt j_ghost = info_check.ys; // Use start of owned y, simple example
                   PetscInt i_ghost = info_check.xs; // Use start of owned x, simple example
                   LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Attempting test read GHOST [%d][%d][%d] (Global)\n", rank_test, k_ghost, j_ghost, i_ghost);
                   Cmpnts test_val_ghost = lUcat_arr_test[k_ghost][j_ghost][i_ghost];
                   LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: SUCCESS reading ghost: x=%g\n", rank_test, test_val_ghost.x);
               } else {
                   LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Skipping ghost test read (zs=0).\n", rank_test);
               }
 
               // Restore the array
               ierr_test = DMDAVecRestoreArrayDOFRead(dm, localVec, &lUcat_arr_test);
               if(ierr_test){ LOG_ALLOW(LOCAL, LOG_ERROR, "Rank %d: ERROR %d restoring '%s' array after test read!\n", rank_test, ierr_test, field_name); }
               LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Finished testing '%s' access.\n", rank_test, field_name);
           }
        } // end if Ucat
    } // end debug logging check
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Rank %d: Completed ghost update for field '%s'.\n", rank, field_name);
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

SetupBoundaryConditions()

PetscErrorCode SetupBoundaryConditions

(

SimCtx *

simCtx

)

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

Parameters

simCtx

Simulation context controlling the operation.

Returns

PetscErrorCode 0 on success.

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

Local to this translation unit.

Definition at line 2678 of file setup.c.

{
    PetscErrorCode ierr;
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    LOG_ALLOW(GLOBAL,LOG_INFO, "--- Setting up Boundary Conditions ---\n");
    // --- Phase 1: Parse and initialize BC configuration for all blocks ---
    LOG_ALLOW(GLOBAL,LOG_INFO,"Parsing BC configuration files and initializing boundary condition data structures.\n");
    UserCtx *user_finest = simCtx->usermg.mgctx[simCtx->usermg.mglevels-1].user;
    for (PetscInt bi = 0; bi < simCtx->block_number; bi++) {
        LOG_ALLOW(GLOBAL,LOG_DEBUG, "  -> Processing Block %d:\n", bi);
 
    // --- Generate the filename for the current block ---
    const char *current_bc_filename = simCtx->bcs_files[bi];
    LOG_ALLOW(GLOBAL,LOG_DEBUG,"  -> Processing Block %d using config file '%s'\n", bi, current_bc_filename);
       // This will populate user_finest[bi].boundary_faces
 
        //ierr = ParseAllBoundaryConditions(&user_finest[bi],current_bc_filename); CHKERRQ(ierr);
 
        ierr = BoundarySystem_Initialize(&user_finest[bi], current_bc_filename); CHKERRQ(ierr);
    }
 
    // Propogate BC Configuration to coarser levels.
    ierr = PropagateBoundaryConfigToCoarserLevels(simCtx); CHKERRQ(ierr);
 
    // Validate the geometric contract before any metric consumes periodic geometry.
    for (PetscInt level = simCtx->usermg.mglevels - 1; level >= 0; level--) {
        UserCtx *level_users = simCtx->usermg.mgctx[level].user;
        for (PetscInt bi = 0; bi < simCtx->block_number; bi++) {
            ierr = ValidatePeriodicGeometry(&level_users[bi]); CHKERRQ(ierr);
        }
    }
 
    // --- Calculate Grid Metrics (requires BC configuration) ---
    // NOTE: This MUST be called here (after BC initialization but before inlet/outlet calculations) because:
    // 1. Periodic BC corrections in metric calculations need boundary_faces data to be populated
    // 2. Inlet/Outlet area calculations (below) require computed metrics (Csi, Eta, Zet) to be available
    // Previously this was in SetupGridAndSolvers, but that caused metrics to be computed without BC info.
    LOG_ALLOW(GLOBAL,LOG_INFO,"Computing grid metrics with boundary condition information.\n");
    ierr = CalculateAllGridMetrics(simCtx); CHKERRQ(ierr);
 
    // --- Phase 2: Calculate inlet/outlet properties (requires computed metrics) ---
    LOG_ALLOW(GLOBAL,LOG_INFO,"Calculating inlet and outlet face properties.\n");
    for (PetscInt bi = 0; bi < simCtx->block_number; bi++) {
        // Call the function to calculate the center of the inlet face & the inlet area, which may be used to calculate Boundary values.
        ierr = CalculateInletProperties(&user_finest[bi]); CHKERRQ(ierr);
 
        // Call the function to calculate the center of the outlet face & the outlet area, which may be used to calculate Boundary values.
        ierr = CalculateOutletProperties(&user_finest[bi]); CHKERRQ(ierr);
    }
 
    LOG_ALLOW(GLOBAL,LOG_INFO, "--- Boundary Conditions setup complete ---\n");   
 
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

InitializeSolutionConvergenceState()

PetscErrorCode InitializeSolutionConvergenceState

(

SimCtx *

simCtx

)

Allocates any runtime storage required by solution-convergence logging.

This helper is setup-owned because it allocates mode-specific buffers once the finest-level Eulerian vectors already exist. It does not compute or log any metrics; it only prepares storage for later use by LOG_SOLUTION_CONVERGENCE().

Allocation depends on the active solution-convergence mode:

• steady/transient: no extra storage beyond counters
• periodic_deterministic: one phase-aligned Ucat / P reference slot per configured period step and per block
• statistical_steady: rolling scalar histories sized for two adjacent windows of mean_speed and mean_ke

The helper also resets the recorded-sample counter so warmup behavior starts from a clean state each time solver setup is rebuilt.

Parameters

[in,out]

simCtx

Master simulation context owning the runtime storage.

Returns

PetscErrorCode 0 on success.

Allocates any runtime storage required by solution-convergence logging.

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

See also

InitializeSolutionConvergenceState()

Definition at line 96 of file setup.c.

{
    UserCtx *user = NULL;
    PetscInt history_capacity = 0;
 
    PetscFunctionBeginUser;
    if (!simCtx) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "SimCtx cannot be NULL.");
    if (simCtx->exec_mode != EXEC_MODE_SOLVER) PetscFunctionReturn(0);
    if (!simCtx->solutionConvergenceEnabled) PetscFunctionReturn(0);
    if (!simCtx->usermg.mgctx) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE,
                "Multigrid hierarchy must exist before initializing solution convergence storage.");
    }
 
    user = simCtx->usermg.mgctx[simCtx->usermg.mglevels - 1].user;
    if (!user) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE,
                "Finest-level UserCtx must exist before initializing solution convergence storage.");
    }
 
    simCtx->solutionConvergenceSamplesRecorded = 0;
 
    if (simCtx->solutionConvergenceMode == SOLUTION_CONVERGENCE_PERIODIC_DETERMINISTIC) {
        for (PetscInt bi = 0; bi < simCtx->block_number; ++bi) {
            PetscCall(PetscCalloc1((size_t)simCtx->solutionConvergencePeriodSteps,
                                   &user[bi].solutionConvergencePeriodicUcatRef));
            PetscCall(PetscCalloc1((size_t)simCtx->solutionConvergencePeriodSteps,
                                   &user[bi].solutionConvergencePeriodicPRef));
            for (PetscInt phase = 0; phase < simCtx->solutionConvergencePeriodSteps; ++phase) {
                PetscCall(VecDuplicate(user[bi].Ucat, &user[bi].solutionConvergencePeriodicUcatRef[phase]));
                PetscCall(VecSet(user[bi].solutionConvergencePeriodicUcatRef[phase], 0.0));
                PetscCall(VecDuplicate(user[bi].P, &user[bi].solutionConvergencePeriodicPRef[phase]));
                PetscCall(VecSet(user[bi].solutionConvergencePeriodicPRef[phase], 0.0));
            }
        }
    }
 
    if (simCtx->solutionConvergenceMode == SOLUTION_CONVERGENCE_STATISTICAL_STEADY) {
        history_capacity = 2 * simCtx->solutionConvergenceWindowSteps;
        PetscCall(PetscCalloc1((size_t)history_capacity, &simCtx->solutionConvergenceMeanSpeedHistory));
        PetscCall(PetscCalloc1((size_t)history_capacity, &simCtx->solutionConvergenceMeanKEHistory));
    }
 
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

DestroySolutionConvergenceState()

PetscErrorCode DestroySolutionConvergenceState

(

SimCtx *

simCtx

)

Frees any runtime storage allocated for solution-convergence logging.

This helper is the cleanup counterpart to InitializeSolutionConvergenceState(). It releases periodic reference buffers and statistical-history arrays, then resets the associated pointers and counters. It is safe to call during teardown even if a given mode never allocated optional storage.

Parameters

[in,out]

simCtx

Master simulation context owning the runtime storage.

Returns

PetscErrorCode 0 on success.

Frees any runtime storage allocated for solution-convergence logging.

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

See also

DestroySolutionConvergenceState()

Definition at line 148 of file setup.c.

{
    UserCtx *user = NULL;
 
    PetscFunctionBeginUser;
    if (!simCtx) PetscFunctionReturn(0);
    if (simCtx->fieldStatisticsWindows) {
        PetscCall(PetscFree(simCtx->fieldStatisticsWindows));
        simCtx->fieldStatisticsWindows = NULL;
        simCtx->fieldStatisticsWindowCount = 0;
    }
    if (!simCtx->usermg.mgctx) PetscFunctionReturn(0);
 
    user = simCtx->usermg.mgctx[simCtx->usermg.mglevels - 1].user;
    if (user) {
        for (PetscInt bi = 0; bi < simCtx->block_number; ++bi) {
            if (user[bi].solutionConvergencePeriodicUcatRef) {
                for (PetscInt phase = 0; phase < simCtx->solutionConvergencePeriodSteps; ++phase) {
                    if (user[bi].solutionConvergencePeriodicUcatRef[phase]) {
                        PetscCall(VecDestroy(&user[bi].solutionConvergencePeriodicUcatRef[phase]));
                    }
                }
                PetscCall(PetscFree(user[bi].solutionConvergencePeriodicUcatRef));
                user[bi].solutionConvergencePeriodicUcatRef = NULL;
            }
            if (user[bi].solutionConvergencePeriodicPRef) {
                for (PetscInt phase = 0; phase < simCtx->solutionConvergencePeriodSteps; ++phase) {
                    if (user[bi].solutionConvergencePeriodicPRef[phase]) {
                        PetscCall(VecDestroy(&user[bi].solutionConvergencePeriodicPRef[phase]));
                    }
                }
                PetscCall(PetscFree(user[bi].solutionConvergencePeriodicPRef));
                user[bi].solutionConvergencePeriodicPRef = NULL;
            }
        }
    }
 
    if (simCtx->solutionConvergenceMeanSpeedHistory) {
        PetscCall(PetscFree(simCtx->solutionConvergenceMeanSpeedHistory));
        simCtx->solutionConvergenceMeanSpeedHistory = NULL;
    }
    if (simCtx->solutionConvergenceMeanKEHistory) {
        PetscCall(PetscFree(simCtx->solutionConvergenceMeanKEHistory));
        simCtx->solutionConvergenceMeanKEHistory = NULL;
    }
    simCtx->solutionConvergenceSamplesRecorded = 0;
 
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

Allocate3DArrayScalar()

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

Allocates a 3D array of PetscReal values using PetscCalloc.

This function dynamically allocates memory for a 3D array of PetscReal values with dimensions nz (layers) x ny (rows) x nx (columns). It uses PetscCalloc1 to ensure the memory is zero-initialized.

The allocation is done in three steps:

1. Allocate an array of nz pointers (one for each layer).
2. Allocate a contiguous block for nz*ny row pointers and assign each layer’s row pointers.
3. Allocate a contiguous block for all nz*ny*nx PetscReal values.

This setup allows the array to be accessed as array[k][j][i], and the memory for the data is contiguous, which improves cache efficiency.

Parameters

[out]	array	Pointer to the 3D array to be allocated.
[in]	nz	Number of layers (z-direction).
[in]	ny	Number of rows (y-direction).
[in]	nx	Number of columns (x-direction).

Returns

PetscErrorCode 0 on success, nonzero on failure.

Allocates a 3D array of PetscReal values using PetscCalloc.

Local to this translation unit.

Definition at line 2742 of file setup.c.

{
  PetscErrorCode ierr;
  PetscReal      ***data;
  PetscReal      *dataContiguous;
  PetscInt       k, j;
 
  PetscFunctionBegin;
  /* Step 1: Allocate memory for an array of nz layer pointers (zero-initialized) */
  ierr = PetscCalloc1(nz, &data); CHKERRQ(ierr);
 
  /* Step 2: Allocate memory for all row pointers (nz * ny pointers) */
  ierr = PetscCalloc1(nz * ny, &data[0]); CHKERRQ(ierr);
  for (k = 1; k < nz; k++) {
    data[k] = data[0] + k * ny;
  }
 
  /* Step 3: Allocate one contiguous block for all data elements (nz*ny*nx) */
  ierr = PetscCalloc1(nz * ny * nx, &dataContiguous); CHKERRQ(ierr);
 
  /* Build the 3D pointer structure: each row pointer gets the correct segment of data */
  for (k = 0; k < nz; k++) {
    for (j = 0; j < ny; j++) {
      data[k][j] = dataContiguous + (k * ny + j) * nx;
      /* Memory is already zeroed by PetscCalloc1, so no manual initialization is needed */
    }
  }
  *array = data;
  PetscFunctionReturn(0);
}

Deallocate3DArrayScalar()

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

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

This function frees the memory allocated for a 3D array of PetscReal values. It assumes the memory was allocated using Allocate3DArrayScalar, which allocated three separate memory blocks: one for the contiguous data, one for the row pointers, and one for the layer pointers.

Parameters

[in]	array	Pointer to the 3D array to be deallocated.
[in]	nz	Number of layers (z-direction).
[in]	ny	Number of rows (y-direction).

Returns

PetscErrorCode 0 on success, nonzero on failure.

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

Local to this translation unit.

Definition at line 2777 of file setup.c.

{
  PetscErrorCode ierr;
  (void)nz;
  (void)ny;
 
  PetscFunctionBegin;
   if (!array || !array[0] || !array[0][0] ) { // Added more robust check
      LOG_ALLOW(GLOBAL, LOG_WARNING, "Deallocate3DArrayScalar called with potentially unallocated or NULL array.\n");
       if (array) {
           if (array[0]) { // Check if row pointers might exist
               // Cannot safely access array[0][0] if array[0] might be invalid/freed
               // Standard deallocation below assumes valid pointers.
                ierr = PetscFree(array[0]); CHKERRQ(ierr); // Free row pointers if they exist
           }
            ierr = PetscFree(array); CHKERRQ(ierr); // Free layer pointers if they exist
       }
       PetscFunctionReturn(0);
   }
 
  // --- Standard Deallocation (assuming valid allocation) ---
 
  /* 1. Free the contiguous block of PetscReal values.
     The starting address was stored in array[0][0]. */
  ierr = PetscFree(array[0][0]); CHKERRQ(ierr); // Free the ACTUAL DATA
 
  /* 2. Free the contiguous block of row pointers.
     The starting address was stored in array[0]. */
  ierr = PetscFree(array[0]); CHKERRQ(ierr); // Free the ROW POINTERS
 
  /* 3. Free the layer pointer array.
     The starting address is 'array' itself. */
  ierr = PetscFree(array); CHKERRQ(ierr); // Free the LAYER POINTERS
 
  PetscFunctionReturn(0);
}

Allocate3DArrayVector()

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

Allocates a contiguous 3D array of Cmpnts values.

The memory layout mirrors Allocate3DArrayScalar: layer pointers, row pointers, and contiguous payload are allocated such that indexing as array[k][j][i] is valid while keeping payload data contiguous.

Parameters

[out]	array	Pointer to the 3D vector array to allocate.
[in]	nz	Number of layers in the z-direction.
[in]	ny	Number of rows in the y-direction.
[in]	nx	Number of columns in the x-direction.

Returns

PetscErrorCode 0 on success, nonzero on failure.

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 2820 of file setup.c.

{
  PetscErrorCode ierr;
  Cmpnts         ***data;
  Cmpnts         *dataContiguous;
  PetscInt       k, j;
  PetscMPIInt    rank;
 
  PetscFunctionBegin;
 
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
  
  /* Step 1: Allocate memory for nz layer pointers (zeroed) */
  ierr = PetscCalloc1(nz, &data); CHKERRQ(ierr);
 
  LOG_ALLOW(LOCAL,LOG_DEBUG," [Rank %d] memory allocated for outermost layer (%d k-layer pointers).\n",rank,nz);
  
  /* Step 2: Allocate memory for all row pointers (nz * ny pointers) */
  ierr = PetscCalloc1(nz * ny, &data[0]); CHKERRQ(ierr);
  for (k = 1; k < nz; k++) {
    data[k] = data[0] + k * ny;
  }
 
  LOG_ALLOW(LOCAL,LOG_DEBUG,"[Rank %d] memory allocated for %dx%d row pointers.\n",rank,nz,ny);
  
  /* Step 3: Allocate one contiguous block for nz*ny*nx Cmpnts structures (zeroed) */
  ierr = PetscCalloc1(nz * ny * nx, &dataContiguous); CHKERRQ(ierr);
 
  LOG_ALLOW(GLOBAL,LOG_DEBUG,"[Rank %d] memory allocated for contigous block of %dx%dx%d Cmpnts structures).\n",rank,nz,ny,nx); 
 
  /* Build the 3D pointer structure for vector data */
  for (k = 0; k < nz; k++) {
    for (j = 0; j < ny; j++) {
      data[k][j] = dataContiguous + (k * ny + j) * nx;
      /* The PetscCalloc1 call has already initialized each Cmpnts to zero. */
    }
  }
 
  LOG_ALLOW(GLOBAL,LOG_DEBUG,"[Rank %d] 3D pointer structure for vector data created. \n",rank);
  
  *array = data;
  PetscFunctionReturn(0);
}

Deallocate3DArrayVector()

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

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

This function frees the memory allocated for a 3D array of Cmpnts structures. It assumes the memory was allocated using Allocate3DArrayVector, which created three separate memory blocks: one for the contiguous vector data, one for the row pointers, and one for the layer pointers.

Parameters

[in]	array	Pointer to the 3D array to be deallocated.
[in]	nz	Number of layers in the z-direction.
[in]	ny	Number of rows in the y-direction.

Returns

PetscErrorCode 0 on success, nonzero on failure.

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 2870 of file setup.c.

{
  PetscErrorCode ierr;
  (void)nz;
  (void)ny;
 
  PetscFunctionBegin;
  // If array is NULL or hasn't been allocated properly, just return.
  if (!array || !array[0] || !array[0][0] ) {
      LOG_ALLOW(GLOBAL, LOG_WARNING, "Deallocate3DArrayVector called with potentially unallocated or NULL array.\n");
      // Attempt to free what might exist, but be cautious
      if (array) {
          if (array[0]) { // Check if row pointers were allocated
             // We don't have a direct pointer to the contiguous data block
             // saved separately in this allocation scheme. The allocation relies
             // on array[0][0] pointing to it. If array[0] was freed first,
             // accessing array[0][0] is unsafe.
             // The allocation scheme where the contiguous data block is not
             // stored separately makes safe deallocation tricky if freeing
             // happens out of order or if parts are NULL.
 
             // A SAFER ALLOCATION/DEALLOCATION would store the data pointer separately.
             // Given the current allocation scheme, the order MUST be:
             // 1. Free the data block (pointed to by array[0][0])
             // 2. Free the row pointer block (pointed to by array[0])
             // 3. Free the layer pointer block (pointed to by array)
 
             // Let's assume the allocation was successful and pointers are valid.
             // Get pointer to the contiguous data block *before* freeing row pointers
             Cmpnts *dataContiguous = array[0][0];
             ierr = PetscFree(dataContiguous); CHKERRQ(ierr); // Free data block
 
             // Now free the row pointers block
             ierr = PetscFree(array[0]); CHKERRQ(ierr); // Free row pointers
 
          }
          // Finally, free the array of layer pointers
          ierr = PetscFree(array); CHKERRQ(ierr);
      }
      PetscFunctionReturn(0); // Return gracefully if input was NULL initially
  }
 
 
  // --- Standard Deallocation (assuming valid allocation) ---
 
  /* 1. Free the contiguous block of Cmpnts structures.
     The starting address was stored in array[0][0] by Allocate3DArrayVector. */
  ierr = PetscFree(array[0][0]); CHKERRQ(ierr); // Free the ACTUAL DATA
 
  /* 2. Free the contiguous block of row pointers.
     The starting address was stored in array[0]. */
  ierr = PetscFree(array[0]); CHKERRQ(ierr); // Free the ROW POINTERS
 
  /* 3. Free the layer pointer array.
     The starting address is 'array' itself. */
  ierr = PetscFree(array); CHKERRQ(ierr); // Free the LAYER POINTERS
 
  PetscFunctionReturn(0);
}

GetOwnedCellRange()

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

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

Ownership is defined by the rank owning the cell's origin node (min i,j,k corner).

Parameters

[in]	info_nodes	Pointer to the DMDALocalInfo struct obtained from the NODE-based DMDA (e.g., user->da or user->fda, assuming they have consistent nodal partitioning for defining cell origins).
[in]	dim	The dimension to compute the range for (0 for x/i, 1 for y/j, 2 for z/k).
[out]	xs_cell_global_out	Pointer to store the starting global cell index owned by this process.
[out]	xm_cell_local_out	Pointer to store the number of owned cells in this dimension.

Returns

PetscErrorCode 0 on success.

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 2936 of file setup.c.

{
    PetscErrorCode ierr = 0; // Standard PETSc error code, not explicitly set here but good practice.
    PetscInt xs_node_global_rank;   // Global index of the first node owned by this rank in the specified dimension.
    PetscInt num_nodes_owned_rank;  // Number of nodes owned by this rank in this dimension (local count, excluding ghosts).
    PetscInt GlobalNodesInDim_from_info; // Total number of DA points in this dimension, from DMDALocalInfo.
 
    PetscFunctionBeginUser;
 
    // --- 1. Input Validation ---
    if (!info_nodes || !xs_cell_global_out || !xm_cell_local_out) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Null pointer passed to GetOwnedCellRange.");
    }
 
    // --- 2. Extract Node Ownership and Global Dimension Information from DMDALocalInfo ---
    if (dim == 0) { // I-direction
        xs_node_global_rank      = info_nodes->xs;
        num_nodes_owned_rank     = info_nodes->xm;
        GlobalNodesInDim_from_info = info_nodes->mx;
    } else if (dim == 1) { // J-direction
        xs_node_global_rank      = info_nodes->ys;
        num_nodes_owned_rank     = info_nodes->ym;
        GlobalNodesInDim_from_info = info_nodes->my;
    } else if (dim == 2) { // K-direction
        xs_node_global_rank      = info_nodes->zs;
        num_nodes_owned_rank     = info_nodes->zm;
        GlobalNodesInDim_from_info = info_nodes->mz;
    } else {
      SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d in GetOwnedCellRange. Must be 0, 1, or 2.", dim);
    }
 
    // --- 3. Correct for User-Defined Ghost Node ---
    // Per the function's contract (@warning), the DA size includes an extra, non-physical
    // node. We subtract 1 to get the true number of physical nodes for cell calculations.
    const PetscInt physical_nodes_in_dim = GlobalNodesInDim_from_info - 1;
 
    // --- 4. Handle Edge Cases for Physical Domain Size ---
    // If the physical domain has 0 or 1 node, no cells can be formed.
    if (physical_nodes_in_dim <= 1) {
        *xs_cell_global_out = xs_node_global_rank; // Still report the rank's starting node
        *xm_cell_local_out  = 0;                    // But 0 cells
        PetscFunctionReturn(0);
    }
 
    // --- 5. Determine Cell Ownership Based on Corrected Node Ownership ---
    // The first cell this rank *could* define has its origin at the first node this rank owns.
    *xs_cell_global_out = xs_node_global_rank;
 
    // If the rank owns no nodes in this dimension, it can't form any cell origins.
    if (num_nodes_owned_rank == 0) {
        *xm_cell_local_out = 0;
    } else {
        // --- BUG FIX APPLIED HERE ---
        // The previous logic incorrectly assumed a cell's end node (N_{k+1}) must be on the
        // same rank as its origin node (N_k). The correct logic is to find the intersection
        // between the nodes this rank owns and the nodes that are valid origins globally.
 
        // The first node owned by the rank is its first potential origin.
        PetscInt first_owned_origin = xs_node_global_rank;
 
        // The absolute last node owned by this rank. Any node up to and including this one
        // is a potential cell origin from this rank's perspective.
        PetscInt last_node_owned_by_rank = xs_node_global_rank + num_nodes_owned_rank - 1;
 
        // The absolute last node in the entire PHYSICAL domain that can serve as a cell origin.
        // If there are `N` physical nodes (0 to N-1), this index is `N-2`.
        PetscInt last_possible_origin_global_idx = physical_nodes_in_dim - 2;
 
        // The actual last origin this rank can provide is the *minimum* of what it owns
        // and what is globally possible. This correctly handles both ranks in the middle of
        // the domain and the very last rank.
        PetscInt actual_last_origin_this_rank_can_form = PetscMin(last_node_owned_by_rank, last_possible_origin_global_idx);
 
        // If the first potential origin this rank owns is already beyond the actual last
        // origin it can form, then this rank forms no valid cell origins. This happens if
        // the rank only owns the very last physical node.
        if (first_owned_origin > actual_last_origin_this_rank_can_form) {
            *xm_cell_local_out = 0;
        } else {
            // The number of cells is the count of valid origins this rank owns.
            // (Count = Last Index - First Index + 1)
            *xm_cell_local_out = actual_last_origin_this_rank_can_form - first_owned_origin + 1;
        }
    }
 
    PetscFunctionReturn(ierr);
}

Here is the caller graph for this function:

ComputeAndStoreNeighborRanks()

PetscErrorCode ComputeAndStoreNeighborRanks

(

UserCtx *

user

)

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

This function retrieves the neighbor information from the primary DMDA (user->da) and stores the face neighbors (xm, xp, ym, yp, zm, zp) in the user->neighbors structure. It assumes a standard PETSc ordering for the neighbors array returned by DMDAGetNeighbors. Logs warnings if the assumed indices seem incorrect (e.g., center rank mismatch).

Parameters

[in,out]

user

Pointer to the UserCtx structure where neighbor info will be stored.

Returns

PetscErrorCode 0 on success, non-zero on failure.

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

Local to this translation unit.

Definition at line 3033 of file setup.c.

{
    PetscErrorCode    ierr;
    PetscMPIInt       rank;
    PetscMPIInt       size; // MPI communicator size
    const PetscMPIInt *neighbor_ranks_ptr; // Pointer to raw neighbor data from PETSc
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
    ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size); CHKERRQ(ierr); // Get MPI size for validation
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Rank %d: Computing DMDA neighbor ranks.\n", rank);
 
    if (!user || !user->da) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "UserCtx or user->da is NULL in ComputeAndStoreNeighborRanks.");
    }
 
    // Get the neighbor information from the DMDA
    // neighbor_ranks_ptr will point to an internal PETSc array of 27 ranks.
    ierr = DMDAGetNeighbors(user->da, &neighbor_ranks_ptr); CHKERRQ(ierr);
 
    // Log the raw values from DMDAGetNeighbors for boundary-relevant directions for debugging
    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",
                   rank,
                   neighbor_ranks_ptr[12], neighbor_ranks_ptr[14],
                   neighbor_ranks_ptr[10], neighbor_ranks_ptr[16],
                   neighbor_ranks_ptr[4],  neighbor_ranks_ptr[22],
                   (int)MPI_PROC_NULL);
 
    // PETSc standard indices for 3D face neighbors from the 27-point stencil:
    // Index = k_offset*9 + j_offset*3 + i_offset (where offsets -1,0,1 map to 0,1,2)
    // Center: (i_off=1, j_off=1, k_off=1) => 1*9 + 1*3 + 1 = 13
    // X-min:  (i_off=0, j_off=1, k_off=1) => 1*9 + 1*3 + 0 = 12
    // X-plus: (i_off=2, j_off=1, k_off=1) => 1*9 + 1*3 + 2 = 14
    // Y-min:  (i_off=1, j_off=0, k_off=1) => 1*9 + 0*3 + 1 = 10
    // Y-plus: (i_off=1, j_off=2, k_off=1) => 1*9 + 2*3 + 1 = 16
    // Z-min:  (i_off=1, j_off=1, k_off=0) => 0*9 + 1*3 + 1 = 4
    // Z-plus: (i_off=1, j_off=1, k_off=2) => 2*9 + 1*3 + 1 = 22
 
    if (neighbor_ranks_ptr[13] != rank) {
        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",
                  rank, neighbor_ranks_ptr[13], rank);
        // This warning is important. If the center isn't the current rank, the offsets are likely wrong.
        // However, PETSc should ensure this unless the DM is too small for a 3x3x3 stencil.
    }
 
    // Assign and sanitize each neighbor rank
    PetscMPIInt temp_neighbor;
 
    temp_neighbor = neighbor_ranks_ptr[12]; // xm
    if (temp_neighbor < 0  || temp_neighbor >= size) {
        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);
        user->neighbors.rank_xm = MPI_PROC_NULL;
    } else {
        user->neighbors.rank_xm = temp_neighbor;
    }
 
    temp_neighbor = neighbor_ranks_ptr[14]; // xp
    if (temp_neighbor < 0 || temp_neighbor >= size) {
        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);
        user->neighbors.rank_xp = MPI_PROC_NULL;
    } else {
        user->neighbors.rank_xp = temp_neighbor;
    }
 
    temp_neighbor = neighbor_ranks_ptr[10]; // ym
    if (temp_neighbor < 0 || temp_neighbor >= size) {
        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);
        user->neighbors.rank_ym = MPI_PROC_NULL;
    } else {
        user->neighbors.rank_ym = temp_neighbor;
    }
 
    temp_neighbor = neighbor_ranks_ptr[16]; // yp
    if (temp_neighbor < 0 || temp_neighbor >= size) {
        // The log for index 16 was "zm" in your output, should be yp
        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);
        user->neighbors.rank_yp = MPI_PROC_NULL;
    } else {
        user->neighbors.rank_yp = temp_neighbor;
    }
 
    temp_neighbor = neighbor_ranks_ptr[4]; // zm
    if (temp_neighbor < 0 || temp_neighbor >= size) {
        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);
        user->neighbors.rank_zm = MPI_PROC_NULL;
    } else {
        user->neighbors.rank_zm = temp_neighbor;
    }
 
    temp_neighbor = neighbor_ranks_ptr[22]; // zp
    if (temp_neighbor < 0 || temp_neighbor >= size) {
        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);
        user->neighbors.rank_zp = MPI_PROC_NULL;
    } else {
        user->neighbors.rank_zp = temp_neighbor;
    }
 
    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,
              user->neighbors.rank_xm, user->neighbors.rank_xp,
              user->neighbors.rank_ym, user->neighbors.rank_yp,
              user->neighbors.rank_zm, user->neighbors.rank_zp);
    PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT); // Ensure logs are flushed
 
    // Note: neighbor_ranks_ptr memory is managed by PETSc, do not free it.
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

SetDMDAProcLayout()

PetscErrorCode SetDMDAProcLayout	(	DM	dm,
		UserCtx *	user
	)

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

Reads the desired number of processors in x, y, and z directions using PETSc options (e.g., -dm_processors_x, -dm_processors_y, -dm_processors_z). If an option is not provided for a direction, PETSC_DECIDE is used for that direction. Applies the layout using DMDASetNumProcs.

Also stores the retrieved/decided values in user->procs_x/y/z if user context is provided.

Parameters

dm	The DMDA object to configure the layout for.
user	Pointer to the UserCtx structure (optional, used to store layout values).

Returns

PetscErrorCode 0 on success, non-zero on failure.

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

Local to this translation unit.

Definition at line 3149 of file setup.c.

{
    PetscErrorCode ierr;
    PetscMPIInt    size, rank;
    PetscInt       px = PETSC_DECIDE, py = PETSC_DECIDE, pz = PETSC_DECIDE;
    PetscBool      px_set = PETSC_FALSE, py_set = PETSC_FALSE, pz_set = PETSC_FALSE;
    SimCtx         *simCtx = user->simCtx;
 
    // Set no.of processors in direction 1
    if(simCtx->da_procs_x) {
      px_set = PETSC_TRUE;
      px = simCtx->da_procs_x;
    }
    // Set no.of processors in direction 2
    if(simCtx->da_procs_y) {
      py_set = PETSC_TRUE;
      py = simCtx->da_procs_y;
    }
    // Set no.of processors in direction 1
    if(simCtx->da_procs_z) {
      pz_set = PETSC_TRUE;
      pz = simCtx->da_procs_z;
    }
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size); CHKERRQ(ierr);
    ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank); CHKERRQ(ierr);
    LOG_ALLOW(GLOBAL, LOG_INFO, "Rank %d: Configuring DMDA processor layout for %d total processes.\n", rank, size);
 
    // --- Validate User Input (Optional but Recommended) ---
    // Check if specified processor counts multiply to the total MPI size
    if (px_set && py_set && pz_set) {
        if (px * py * pz != size) {
             SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_INCOMP,
                     "Specified processor layout %d x %d x %d = %d does not match MPI size %d",
                     px, py, pz, px * py * pz, size);
        }
         LOG_ALLOW(GLOBAL, LOG_INFO, "Using specified processor layout: %d x %d x %d\n", px, py, pz);
    } else if (px_set || py_set || pz_set) {
         // If only some are set, PETSC_DECIDE will be used for others
         LOG_ALLOW(GLOBAL, LOG_INFO, "Using partially specified processor layout: %d x %d x %d (PETSC_DECIDE for unspecified)\n", px, py, pz);
    } else {
        LOG_ALLOW(GLOBAL, LOG_INFO, "Using fully automatic processor layout (PETSC_DECIDE x PETSC_DECIDE x PETSC_DECIDE)\n");
    }
     // Additional checks: Ensure px, py, pz are positive if set
     if ((px_set && px <= 0) || (py_set && py <= 0) || (pz_set && pz <= 0)) {
         SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Specified processor counts must be positive.");
     }
 
 
    // --- Apply the layout to the DMDA ---
    ierr = DMDASetNumProcs(dm, px, py, pz); CHKERRQ(ierr);
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Rank %d: DMDASetNumProcs called with px=%d, py=%d, pz=%d.\n", rank, px, py, pz);
 
    // --- Store the values in UserCtx (Optional) ---
    // Note: If PETSC_DECIDE was used, PETSc calculates the actual values during DMSetUp.
    // We store the *requested* values here. To get the *actual* values used,
    // you would need to call DMDAGetInfo after DMSetUp.
    /*
    if (user) {
        user->procs_x = px;
        user->procs_y = py;
        user->procs_z = pz;
    }
    */
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

SetupDomainRankInfo()

PetscErrorCode SetupDomainRankInfo

(

SimCtx *

simCtx

)

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

This function orchestrates the following steps:

1. Compute and store the neighbor ranks in the user context.
2. Gather all local bounding boxes to rank 0.
3. Broadcast the complete bounding box list to all ranks.

The final result is that each rank has access to its immediate neighbors and the bounding box information of all ranks.

Parameters

[in,out]

simCtx

Pointer to initialized simulation context that owns all block UserCtx objects.

Returns

PetscErrorCode Returns 0 on success or non-zero PETSc error code.

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 3227 of file setup.c.

{
    PetscErrorCode ierr;
    PetscInt       nblk = simCtx->block_number;
    PetscInt       size = simCtx->size;
    BoundingBox    *final_bboxlist = NULL;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Starting full rank communication setup for %d block(s).\n", nblk);
 
    UserCtx *user_finest = simCtx->usermg.mgctx[simCtx->usermg.mglevels - 1].user;
 
    // --- Step 1: Compute neighbor ranks (unchanged) ---
    for (int bi = 0; bi < nblk; bi++) {
        ierr = ComputeAndStoreNeighborRanks(&user_finest[bi]); CHKERRQ(ierr);
    }
    LOG_ALLOW(GLOBAL, LOG_INFO, "Neighbor ranks computed and stored for all blocks.\n");
 
    // --- Step 2: Allocate the final, unified list on ALL ranks ---
    // Every rank will build this list in parallel.
    ierr = PetscMalloc1(size * nblk, &final_bboxlist); CHKERRQ(ierr);
 
    // --- Step 3: Loop through each block, gather then broadcast its bbox list ---
    for (int bi = 0; bi < nblk; bi++) {
        // This is a temporary pointer for the current block's list.
        BoundingBox *block_bboxlist = NULL;
 
        LOG_ALLOW(GLOBAL, LOG_INFO, "Processing bounding boxes for block %d...\n", bi);
 
        // A) GATHER: On rank 0, block_bboxlist is allocated and filled. On others, it's NULL.
        ierr = GatherAllBoundingBoxes(&user_finest[bi], &block_bboxlist); CHKERRQ(ierr);
        LOG_ALLOW(GLOBAL, LOG_DEBUG, "  -> Gather complete for block %d.\n", bi);
 
        // B) BROADCAST: On non-root ranks, block_bboxlist is allocated. Then, the data
        //    from rank 0 is broadcast to all ranks. After this call, ALL ranks have
        //    an identical, complete copy of the bounding boxes for the current block.
        ierr = BroadcastAllBoundingBoxes(&user_finest[bi], &block_bboxlist); CHKERRQ(ierr);
        LOG_ALLOW(GLOBAL, LOG_DEBUG, "  -> Broadcast complete for block %d.\n", bi);
 
        // C) ASSEMBLE: Every rank now copies the data for this block into the
        //    correct segment of its final, unified list.
        for (int r = 0; r < size; r++) {
            // The layout is [r0b0, r1b0, ..., r(size-1)b0, r0b1, r1b1, ...]
            final_bboxlist[bi * size + r] = block_bboxlist[r];
        }
        LOG_ALLOW(GLOBAL, LOG_DEBUG, "  -> Assembly into final list complete for block %d.\n", bi);
 
        // D) CLEANUP: Free the temporary list for this block on ALL ranks before the next iteration.
        //    Your helper functions use malloc, so we must use free.
        free(block_bboxlist);
    }
 
    // --- Step 4: Assign the final pointer and run the last setup step ---
    simCtx->bboxlist = final_bboxlist;
    LOG_ALLOW(GLOBAL, LOG_INFO, "Final unified bboxlist created on all ranks and stored in SimCtx.\n");
 
    ierr = SetupDomainCellDecompositionMap(&user_finest[0]); CHKERRQ(ierr);
    LOG_ALLOW(GLOBAL, LOG_INFO, "Domain Cell Composition set and broadcasted.\n");
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "SetupDomainRankInfo: Completed successfully.\n");
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

Contra2Cart()

PetscErrorCode Contra2Cart

(

UserCtx *

user

)

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

This function performs the transformation from a contravariant velocity representation (which is natural on a curvilinear grid) to a Cartesian (x,y,z) representation. For each interior computational cell owned by the rank, it performs the following:

1. It averages the contravariant velocity components (U¹, U², U³) from the surrounding faces to get an estimate of the contravariant velocity at the cell center.
2. It averages the metric vectors (Csi, Eta, Zet) from the surrounding faces to get an estimate of the metric tensor at the cell center. This tensor forms the transformation matrix.
3. It solves the linear system [MetricTensor] * [ucat] = [ucont] for the Cartesian velocity vector ucat = (u,v,w) using Cramer's rule.
4. The computed Cartesian velocity is stored only at independent owned cells in the global user->Ucat vector.

The function operates on local, ghosted versions of the input vectors (user->lUcont, user->lCsi, etc.) to ensure stencils are valid across processor boundaries.

Parameters

[in,out]

user

Pointer to the UserCtx structure. The function reads from user->lUcont, user->lCsi, user->lEta, user->lZet, user->lNvert and writes to the global user->Ucat vector.

Returns

PetscErrorCode 0 on success.

Note

• This function should be called AFTER user->lUcont and all local metric vectors (user->lCsi, etc.) have been populated with up-to-date ghost values via UpdateLocalGhosts.
• It only computes Ucat for interior cells (not on physical boundaries) and for cells not marked as solid/blanked by user->lNvert.
• The function performs no communication, does not synchronize periodic duplicate planes, and does not refresh user->lUcat.
• Before stencil consumption, the caller must synchronize periodic cell values with SynchronizePeriodicCellFields(user, 1, {FIELD_ID_UCAT}) and then call UpdateLocalGhosts(user, FIELD_ID_UCAT). The latter remains necessary for the fully nonperiodic case because periodic synchronization is then a no-op.

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

Local to this translation unit.

Definition at line 3300 of file setup.c.

{
    PetscErrorCode ierr;
    DMDALocalInfo  info;
    Cmpnts       ***lcsi_arr, ***leta_arr, ***lzet_arr; // Local metric arrays
    Cmpnts       ***lucont_arr;                       // Local contravariant velocity array
    Cmpnts       ***gucat_arr;                        // Global Cartesian velocity array
    PetscReal    ***lnvert_arr;                       // Local Nvert array
    PetscReal    ***laj_arr;                          // Local Jacobian Determinant inverse array
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Starting Contravariant-to-Cartesian velocity transformation.\n");
 
    // --- 1. Get DMDA Info and Check for Valid Inputs ---
    // All inputs (lUcont, lCsi, etc.) and outputs (Ucat) are on DMs from the UserCtx.
    // We get local info from fda, which governs the layout of most arrays here.
    ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
    if (!user->lUcont || !user->lCsi || !user->lEta || !user->lZet || !user->lNvert || !user->Ucat) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Contra2Cart requires lUcont, lCsi/Eta/Zet, lNvert, and Ucat to be non-NULL.");
    }
 
 
    // --- 2. Get Read-Only Array Access to Local Input Vectors (with ghosts) ---
    ierr = DMDAVecGetArrayRead(user->fda, user->lUcont, &lucont_arr); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lCsi,   &lcsi_arr);   CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lEta,   &leta_arr);   CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lZet,   &lzet_arr);   CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->da,  user->lNvert, &lnvert_arr); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->da,  user->lAj, &laj_arr); CHKERRQ(ierr);
 
    // --- 3. Get Write-Only Array Access to the Global Output Vector ---
    // We compute for local owned cells and write into the global vector.
    // PETSc handles mapping the global indices to the correct local memory locations.
    ierr = DMDAVecGetArray(user->fda, user->Ucat, &gucat_arr); CHKERRQ(ierr);
 
 
    // --- 4. Define Loop Bounds for INTERIOR Cells ---
    // We use adjusted bounds to avoid calculating Ucat on the physical domain boundaries,
    // as these are typically set explicitly by boundary condition functions.
    // The stencils use indices like i-1, j-1, k-1, so we must start loops at least at index 1.
    PetscInt i_start = (info.xs == 0) ? info.xs + 1 : info.xs;
    PetscInt i_end   = (info.xs + info.xm == info.mx) ? info.xs + info.xm - 1 : info.xs + info.xm;
 
    PetscInt j_start = (info.ys == 0) ? info.ys + 1 : info.ys;
    PetscInt j_end   = (info.ys + info.ym == info.my) ? info.ys + info.ym - 1 : info.ys + info.ym;
 
    PetscInt k_start = (info.zs == 0) ? info.zs + 1 : info.zs;
    PetscInt k_end   = (info.zs + info.zm == info.mz) ? info.zs + info.zm - 1 : info.zs + info.zm;
 
    // --- 5. Main Computation Loop ---
    // Loops over the GLOBAL indices of interior cells owned by this rank.
    for (PetscInt k_cell = k_start; k_cell < k_end; ++k_cell) {
        for (PetscInt j_cell = j_start; j_cell < j_end; ++j_cell) {
            for (PetscInt i_cell = i_start; i_cell < i_end; ++i_cell) {
 
                // Check if the cell is a fluid cell (not solid/blanked)
          //    if (lnvert_arr[k_cell][j_cell][i_cell] > 0.1) continue; // Skip solid/blanked cells
 
                // Transformation matrix [mat] is the metric tensor at the cell center,
                // estimated by averaging metrics from adjacent faces.
                PetscReal mat[3][3];
 
        //  PetscReal aj_center = laj_arr[k_cell+1][j_cell+1][i_cell+1];
        
                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;
                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;
                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;
 
                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;
                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;
                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;
 
                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;
                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;
                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;
    
                // Contravariant velocity vector `q` at the cell center,
                // estimated by averaging face-based contravariant velocities.
                PetscReal q[3];
                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
                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
                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
 
                // Solve the 3x3 system `mat * ucat = q` using Cramer's rule.
                PetscReal det = mat[0][0] * (mat[1][1] * mat[2][2] - mat[1][2] * mat[2][1]) -
                                mat[0][1] * (mat[1][0] * mat[2][2] - mat[1][2] * mat[2][0]) +
                                mat[0][2] * (mat[1][0] * mat[2][1] - mat[1][1] * mat[2][0]);
 
          if (PetscAbsReal(det) < 1.0e-18) {
              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);
                        }
 
                PetscReal det_inv = 1.0 / det;
 
                PetscReal det0 = q[0] * (mat[1][1] * mat[2][2] - mat[1][2] * mat[2][1]) -
                                 q[1] * (mat[0][1] * mat[2][2] - mat[0][2] * mat[2][1]) +
                                 q[2] * (mat[0][1] * mat[1][2] - mat[0][2] * mat[1][1]);
 
                PetscReal det1 = -q[0] * (mat[1][0] * mat[2][2] - mat[1][2] * mat[2][0]) +
                                  q[1] * (mat[0][0] * mat[2][2] - mat[0][2] * mat[2][0]) -
                                  q[2] * (mat[0][0] * mat[1][2] - mat[0][2] * mat[1][0]);
 
                PetscReal det2 = q[0] * (mat[1][0] * mat[2][1] - mat[1][1] * mat[2][0]) -
                                 q[1] * (mat[0][0] * mat[2][1] - mat[0][1] * mat[2][0]) +
                                 q[2] * (mat[0][0] * mat[1][1] - mat[0][1] * mat[1][0]);
 
                // Store computed Cartesian velocity in the GLOBAL Ucat array at the
                // array index corresponding to the cell's origin node.
                gucat_arr[k_cell][j_cell][i_cell].x = det0 * det_inv;
                gucat_arr[k_cell][j_cell][i_cell].y = det1 * det_inv;
                gucat_arr[k_cell][j_cell][i_cell].z = det2 * det_inv;
            }
        }
    }
 
    // --- 6. Restore Array Access ---
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lUcont, &lucont_arr); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi,   &lcsi_arr);   CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta,   &leta_arr);   CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet,   &lzet_arr);   CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->da,  user->lNvert, &lnvert_arr); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->da,  user->lAj, &laj_arr); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(user->fda, user->Ucat, &gucat_arr); CHKERRQ(ierr);
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Completed Contravariant-to-Cartesian velocity transformation. \n");
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

Cart2Contra()

PetscErrorCode Cart2Contra

(

UserCtx *

user

)

Convert the ghosted Cartesian velocity field to contravariant face fluxes.

Interpolates user->lUcat to each logical face and dots it with the corresponding face-area metric vector. Writes independent entries of user->Ucont. Periodic and physical endpoint values in user->lUcat must be finalized and its MPI ghosts current before entry. The resulting Ucont is not synchronized; callers must use the staggered-field workflow before stencil use.

Parameters

[in,out]

user

UserCtx for the block; Ucont is written in place.

Returns

PetscErrorCode 0 on success.

Convert the ghosted Cartesian velocity field to contravariant face fluxes.

Definition at line 3435 of file setup.c.

{
    PetscErrorCode ierr;
    DMDALocalInfo info;
    const Cmpnts ***ucat_arr, ***csi_arr, ***eta_arr, ***zet_arr;
    Cmpnts ***ucont_arr;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lUcat, &ucat_arr); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, &csi_arr); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lEta, &eta_arr); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lZet, &zet_arr); CHKERRQ(ierr);
    ierr = DMDAVecGetArray(user->fda, user->Ucont, &ucont_arr); CHKERRQ(ierr);
 
    const PetscInt i_start = PetscMax(info.xs, 1);
    const PetscInt j_start = PetscMax(info.ys, 1);
    const PetscInt k_start = PetscMax(info.zs, 1);
    const PetscInt i_end = PetscMin(info.xs + info.xm, info.mx - 1);
    const PetscInt j_end = PetscMin(info.ys + info.ym, info.my - 1);
    const PetscInt k_end = PetscMin(info.zs + info.zm, info.mz - 1);
 
    for (PetscInt k = k_start; k < k_end; k++) {
        for (PetscInt j = j_start; j < j_end; j++) {
            for (PetscInt i = i_start; i < i_end; i++) {
                const Cmpnts u_xi = {
                    0.5 * (ucat_arr[k][j][i].x + ucat_arr[k][j][i + 1].x),
                    0.5 * (ucat_arr[k][j][i].y + ucat_arr[k][j][i + 1].y),
                    0.5 * (ucat_arr[k][j][i].z + ucat_arr[k][j][i + 1].z)
                };
                const Cmpnts u_eta = {
                    0.5 * (ucat_arr[k][j][i].x + ucat_arr[k][j + 1][i].x),
                    0.5 * (ucat_arr[k][j][i].y + ucat_arr[k][j + 1][i].y),
                    0.5 * (ucat_arr[k][j][i].z + ucat_arr[k][j + 1][i].z)
                };
                const Cmpnts u_zeta = {
                    0.5 * (ucat_arr[k][j][i].x + ucat_arr[k + 1][j][i].x),
                    0.5 * (ucat_arr[k][j][i].y + ucat_arr[k + 1][j][i].y),
                    0.5 * (ucat_arr[k][j][i].z + ucat_arr[k + 1][j][i].z)
                };
                ucont_arr[k][j][i].x = csi_arr[k][j][i].x * u_xi.x + csi_arr[k][j][i].y * u_xi.y + csi_arr[k][j][i].z * u_xi.z;
                ucont_arr[k][j][i].y = eta_arr[k][j][i].x * u_eta.x + eta_arr[k][j][i].y * u_eta.y + eta_arr[k][j][i].z * u_eta.z;
                ucont_arr[k][j][i].z = zet_arr[k][j][i].x * u_zeta.x + zet_arr[k][j][i].y * u_zeta.y + zet_arr[k][j][i].z * u_zeta.z;
            }
        }
    }
 
    ierr = DMDAVecRestoreArray(user->fda, user->Ucont, &ucont_arr); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, &zet_arr); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, &eta_arr); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, &csi_arr); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lUcat, &ucat_arr); CHKERRQ(ierr);
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

UniformCart2Contra()

PetscErrorCode UniformCart2Contra	(	UserCtx *	user,
		PetscReal	u,
		PetscReal	v,
		PetscReal	w
	)

Populate contravariant fluxes from one uniform Cartesian velocity.

Parameters

[in,out]	user	UserCtx for the block; Ucont is written in place.
[in]	u	Physical x-velocity component.
[in]	v	Physical y-velocity component.
[in]	w	Physical z-velocity component.

Returns

PetscErrorCode 0 on success.

Definition at line 3504 of file setup.c.

{
    PetscErrorCode ierr;
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    DMDALocalInfo    info;
    Cmpnts         ***ucont_arr;
    const Cmpnts   ***csi_arr, ***eta_arr, ***zet_arr;
 
    ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
    ierr = DMDAVecGetArray(user->fda, user->Ucont, &ucont_arr);        CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, &csi_arr);       CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lEta, &eta_arr);       CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lZet, &zet_arr);       CHKERRQ(ierr);
 
    const PetscInt xs = info.xs, xe = info.xs + info.xm;
    const PetscInt ys = info.ys, ye = info.ys + info.ym;
    const PetscInt zs = info.zs, ze = info.zs + info.zm;
 
    for (PetscInt k = zs; k < ze; k++) {
        for (PetscInt j = ys; j < ye; j++) {
            for (PetscInt i = xs; i < xe; i++) {
                ucont_arr[k][j][i].x = csi_arr[k][j][i].x * u + csi_arr[k][j][i].y * v + csi_arr[k][j][i].z * w;
                ucont_arr[k][j][i].y = eta_arr[k][j][i].x * u + eta_arr[k][j][i].y * v + eta_arr[k][j][i].z * w;
                ucont_arr[k][j][i].z = zet_arr[k][j][i].x * u + zet_arr[k][j][i].y * v + zet_arr[k][j][i].z * w;
            }
        }
    }
 
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, &zet_arr);   CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, &eta_arr);   CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, &csi_arr);   CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(user->fda, user->Ucont, &ucont_arr);    CHKERRQ(ierr);
 
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Cart2Contra: set Ucont from Cartesian (%.3f, %.3f, %.3f).\n",
              (double)u, (double)v, (double)w);
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

SetupDomainCellDecompositionMap()

PetscErrorCode SetupDomainCellDecompositionMap

(

UserCtx *

user

)

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

This function is a critical part of the simulation setup. It determines the global cell ownership for each MPI rank and makes this information available to all other ranks. This "decomposition map" is essential for the robust "Walk and Handoff" particle migration strategy, allowing any rank to quickly identify the owner of a target cell.

The process involves:

1. Each rank gets its own node ownership information from the DMDA.
2. It converts this node information into cell ownership ranges using the GetOwnedCellRange helper function.
3. It participates in an MPI_Allgather collective operation to build a complete array (user->RankCellInfoMap) containing the ownership information for every rank.

This function should be called once during initialization after the primary DMDA (user->da) has been set up.

Parameters

[in,out]

user

Pointer to the UserCtx structure. The function will allocate and populate user->RankCellInfoMap and set user->num_ranks.

Returns

PetscErrorCode 0 on success, or a non-zero PETSc error code on failure. Errors can occur if input pointers are NULL or if MPI communication fails.

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

Local to this translation unit.

Definition at line 3551 of file setup.c.

{
    PetscErrorCode ierr;
    DMDALocalInfo  local_node_info;
    RankCellInfo   my_cell_info;
    PetscMPIInt    rank, size;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    // --- 1. Input Validation and MPI Info ---
    if (!user) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "UserCtx pointer is NULL in SetupDomainCellDecompositionMap.");
    }
    if (!user->da) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "user->da is not initialized in SetupDomainCellDecompositionMap.");
    }
 
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
    ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size); CHKERRQ(ierr);
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Setting up domain cell decomposition map for %d ranks.\n", size);
 
    // --- 2. Determine Local Cell Ownership ---
    // Get the local node ownership information from the primary DMDA.
    ierr = DMDAGetLocalInfo(user->da, &local_node_info); CHKERRQ(ierr);
 
    // Use the robust helper function to convert node ownership to cell ownership.
    // A cell's index is defined by its origin node.
     
    ierr = GetOwnedCellRange(&local_node_info, 0, &my_cell_info.xs_cell, &my_cell_info.xm_cell); CHKERRQ(ierr);
    ierr = GetOwnedCellRange(&local_node_info, 1, &my_cell_info.ys_cell, &my_cell_info.ym_cell); CHKERRQ(ierr);
    ierr = GetOwnedCellRange(&local_node_info, 2, &my_cell_info.zs_cell, &my_cell_info.zm_cell); CHKERRQ(ierr);
    
    // Log the calculated local ownership for debugging purposes.
    LOG_ALLOW(LOCAL, LOG_DEBUG, "[Rank %d] Owns cells: i[%d, %d), j[%d, %d), k[%d, %d)\n",
              rank, my_cell_info.xs_cell, my_cell_info.xs_cell + my_cell_info.xm_cell,
              my_cell_info.ys_cell, my_cell_info.ys_cell + my_cell_info.ym_cell,
              my_cell_info.zs_cell, my_cell_info.zs_cell + my_cell_info.zm_cell);
 
    // --- 3. Allocate and Distribute the Global Map ---
    // Allocate memory for the global map that will hold information from all ranks.
    ierr = PetscMalloc1(size, &user->RankCellInfoMap); CHKERRQ(ierr);
 
    // Perform the collective communication to gather the `RankCellInfo` struct from every rank.
    // Each rank sends its `my_cell_info` and receives the complete array in `user->RankCellInfoMap`.
    // We use MPI_BYTE to ensure portability across different systems and struct padding.
    ierr = MPI_Allgather(&my_cell_info, sizeof(RankCellInfo), MPI_BYTE,
                         user->RankCellInfoMap, sizeof(RankCellInfo), MPI_BYTE,
                         PETSC_COMM_WORLD); CHKERRQ(ierr);
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Domain cell decomposition map created and distributed successfully.\n");
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

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
	)

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

This is a standard binary search algorithm implemented as a PETSc-style helper function. It efficiently determines if a given key exists within a sorted array.

Parameters

[in]	n	The number of elements in the array.
[in]	arr	A pointer to the sorted array of PetscInt64 values to be searched.
[in]	key	The PetscInt64 value to search for.
[out]	found	A pointer to a PetscBool that will be set to PETSC_TRUE if the key is found, and PETSC_FALSE otherwise.

Returns

PetscErrorCode 0 on success, or a non-zero PETSc error code on failure.

Note

The input array arr must be sorted in ascending order for the algorithm to work correctly.

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 3616 of file setup.c.

{
    PetscInt low = 0, high = n - 1;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    // --- 1. Input Validation ---
    if (!found) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Output pointer 'found' is NULL in PetscBinarySearchInt64.");
    }
    if (n > 0 && !arr) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Input array 'arr' is NULL for n > 0.");
    }
    
    // Initialize output
    *found = PETSC_FALSE;
 
    // --- 2. Binary Search Algorithm ---
    while (low <= high) {
        // Use this form to prevent potential integer overflow on very large arrays
        PetscInt mid = low + (high - low) / 2;
 
        if (arr[mid] == key) {
            *found = PETSC_TRUE; // Key found!
            break;               // Exit the loop
        }
        
        if (arr[mid] < key) {
            low = mid + 1; // Search in the right half
        } else {
            high = mid - 1; // Search in the left half
        }
    }
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

ComputeDivergence()

PetscErrorCode ComputeDivergence

(

UserCtx *

user

)

Computes the discrete divergence of the contravariant velocity field.

This diagnostic/kernel routine evaluates continuity residuals on the local block and writes the resulting divergence field into the configured output vector(s).

Parameters

[in,out]

user

Block-level context containing velocity and metric fields.

Returns

PetscErrorCode 0 on success.

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 3685 of file setup.c.

{
  DM        da = user->da, fda = user->fda;
  DMDALocalInfo info = user->info;
 
  PetscInt ti = user->simCtx->step;
  
  PetscInt  xs = info.xs, xe = info.xs + info.xm;
  PetscInt      ys = info.ys, ye = info.ys + info.ym;
  PetscInt  zs = info.zs, ze = info.zs + info.zm;
  PetscInt  mx = info.mx, my = info.my, mz = info.mz;
 
  PetscInt  lxs, lys, lzs, lxe, lye, lze;
  PetscInt  i, j, k;
 
  Vec       Div;
  PetscReal ***div, ***aj, ***nvert,***p;
  Cmpnts    ***ucont;
  PetscReal maxdiv;
 
  lxs = xs; lxe = xe;
  lys = ys; lye = ye;
  lzs = zs; lze = ze;
 
  if (xs==0) lxs = xs+1;
  if (ys==0) lys = ys+1;
  if (zs==0) lzs = zs+1;
 
  if (xe==mx) lxe = xe-1;
  if (ye==my) lye = ye-1;
  if (ze==mz) lze = ze-1;
 
  PetscFunctionBeginUser;
  PROFILE_FUNCTION_BEGIN;
 
  DMDAVecGetArray(fda,user->lUcont, &ucont);
  DMDAVecGetArray(da, user->lAj, &aj);
  VecDuplicate(user->P, &Div);
  DMDAVecGetArray(da, Div, &div);
  DMDAVecGetArray(da, user->lNvert, &nvert);
   DMDAVecGetArray(da, user->P, &p);
  for (k=lzs; k<lze; k++) {
         for (j=lys; j<lye; j++){
      for (i=lxs; i<lxe; i++) {
        if (k==10 && j==10 && i==1){
          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]);
              }
 
        if (k==10 && j==10 && i==mx-3)
              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]);
        }
    }
    }
   DMDAVecRestoreArray(da, user->P, &p);
 
 
  for (k=lzs; k<lze; k++) {
    for (j=lys; j<lye; j++) {
      for (i=lxs; i<lxe; i++) {
    maxdiv = fabs((ucont[k][j][i].x - ucont[k][j][i-1].x +
               ucont[k][j][i].y - ucont[k][j-1][i].y +
               ucont[k][j][i].z - ucont[k-1][j][i].z)*aj[k][j][i]);
    if (nvert[k][j][i] + nvert[k+1][j][i] + nvert[k-1][j][i] +
        nvert[k][j+1][i] + nvert[k][j-1][i] +
        nvert[k][j][i+1] + nvert[k][j][i-1] > 0.1) maxdiv = 0.;
    div[k][j][i] = maxdiv;
 
      }
    }
  }
 
  if (zs==0) {
    k=0;
    for (j=ys; j<ye; j++) {
      for (i=xs; i<xe; i++) {
    div[k][j][i] = 0.;
      }
    }
  }
 
  if (ze == mz) {
    k=mz-1;
    for (j=ys; j<ye; j++) {
      for (i=xs; i<xe; i++) {
    div[k][j][i] = 0.;
      }
    }
  }
 
  if (xs==0) {
    i=0;
    for (k=zs; k<ze; k++) {
      for (j=ys; j<ye; j++) {
    div[k][j][i] = 0.;
      }
    }
  }
 
  if (xe==mx) {
    i=mx-1;
    for (k=zs; k<ze; k++) {
      for (j=ys; j<ye; j++) {
    div[k][j][i] = 0;
      }
    }
  }
 
  if (ys==0) {
    j=0;
    for (k=zs; k<ze; k++) {
      for (i=xs; i<xe; i++) {
    div[k][j][i] = 0.;
      }
    }
  }
 
  if (ye==my) {
    j=my-1;
    for (k=zs; k<ze; k++) {
      for (i=xs; i<xe; i++) {
    div[k][j][i] = 0.;
      }
    }
  }
  DMDAVecRestoreArray(da, Div, &div);
  PetscInt MaxFlatIndex;
  
  VecMax(Div, &MaxFlatIndex, &maxdiv);
 
  LOG_ALLOW(GLOBAL,LOG_INFO,"[Step %d]] The Maximum Divergence is %e at flat index %d.\n",ti,maxdiv,MaxFlatIndex); 
 
  user->simCtx->MaxDivFlatArg = MaxFlatIndex;
  user->simCtx->MaxDiv = maxdiv;
  
  for (k=zs; k<ze; k++) {
    for (j=ys; j<ye; j++) {
      for (i=xs; i<xe; i++) {
    if (Gidx(i,j,k,user) == MaxFlatIndex) {
      LOG_ALLOW(GLOBAL,LOG_INFO,"[Step %d] The Maximum Divergence(%e) is at location [%d][%d][%d]. \n", ti, maxdiv,k,j,i);
      user->simCtx->MaxDivz = k;
      user->simCtx->MaxDivy = j;
      user->simCtx->MaxDivx = i;
    }
      }
    }
  }
 
  
 DMDAVecRestoreArray(da, user->lNvert, &nvert);
 DMDAVecRestoreArray(fda, user->lUcont, &ucont);
 DMDAVecRestoreArray(da, user->lAj, &aj);
 VecDestroy(&Div);
 
 PROFILE_FUNCTION_END;
 PetscFunctionReturn(0);
}

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
	)

Initializes random number generators for assigning particle properties.

This function creates and configures separate PETSc random number generators for the x, y, and z coordinates.

Parameters

[in,out]	user	Pointer to the UserCtx structure containing simulation context.
[out]	randx	Pointer to store the RNG for the x-coordinate.
[out]	randy	Pointer to store the RNG for the y-coordinate.
[out]	randz	Pointer to store the RNG for the z-coordinate.

Returns

PetscErrorCode Returns 0 on success, non-zero on failure.

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 3851 of file setup.c.

                                                                                                                     {
    PetscErrorCode ierr;  // Error code for PETSc functions
    PetscMPIInt rank;
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
 
    // Initialize RNG for x-coordinate
    ierr = PetscRandomCreate(PETSC_COMM_SELF, randx); CHKERRQ(ierr);
    ierr = PetscRandomSetType((*randx), PETSCRAND48); CHKERRQ(ierr);
    ierr = PetscRandomSetInterval(*randx, user->bbox.min_coords.x, user->bbox.max_coords.x); CHKERRQ(ierr);
    /* Offsets reproduce the historical fixed seeds (12345, 67890, 54321 plus rank) at the
       default base seed, so existing placements are unchanged. */
    const unsigned long base_seed = (unsigned long)user->simCtx->particleRandomSeed;
    ierr = PetscRandomSetSeed(*randx, base_seed + (unsigned long)rank); CHKERRQ(ierr);
    ierr = PetscRandomSeed(*randx); CHKERRQ(ierr);
    LOG_ALLOW_SYNC(LOCAL,LOG_VERBOSE, "[Rank %d]Initialized RNG for X-axis.\n",rank);
 
    // Initialize RNG for y-coordinate
    ierr = PetscRandomCreate(PETSC_COMM_SELF, randy); CHKERRQ(ierr);
    ierr = PetscRandomSetType((*randy), PETSCRAND48); CHKERRQ(ierr);
    ierr = PetscRandomSetInterval(*randy, user->bbox.min_coords.y, user->bbox.max_coords.y); CHKERRQ(ierr);
    ierr = PetscRandomSetSeed(*randy, base_seed + 55545UL + (unsigned long)rank); CHKERRQ(ierr);
    ierr = PetscRandomSeed(*randy); CHKERRQ(ierr);
    LOG_ALLOW_SYNC(LOCAL,LOG_VERBOSE, "[Rank %d]Initialized RNG for Y-axis.\n",rank);
 
    // Initialize RNG for z-coordinate
    ierr = PetscRandomCreate(PETSC_COMM_SELF, randz); CHKERRQ(ierr);
    ierr = PetscRandomSetType((*randz), PETSCRAND48); CHKERRQ(ierr);
    ierr = PetscRandomSetInterval(*randz, user->bbox.min_coords.z, user->bbox.max_coords.z); CHKERRQ(ierr);
    ierr = PetscRandomSetSeed(*randz, base_seed + 41976UL + (unsigned long)rank); CHKERRQ(ierr);
    ierr = PetscRandomSeed(*randz); CHKERRQ(ierr);
    LOG_ALLOW_SYNC(LOCAL,LOG_VERBOSE, "[Rank %d]Initialized RNG for Z-axis.\n",rank);
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

InitializeLogicalSpaceRNGs()

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

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

This function creates and configures three separate PETSc random number generators, one for each logical dimension (i, j, k or xi, eta, zeta equivalent). Each RNG is configured to produce uniformly distributed real numbers in the interval [0.0, 1.0). These are typically used for selecting owned cells or generating intra-cell logical coordinates.

Parameters

[in]	base_seed	Base particle seed (SimCtx::particleRandomSeed); each stream adds a fixed offset and the rank.
[out]	rand_logic_i	Pointer to store the RNG for the i-logical dimension.
[out]	rand_logic_j	Pointer to store the RNG for the j-logical dimension.
[out]	rand_logic_k	Pointer to store the RNG for the k-logical dimension.

Returns

PetscErrorCode Returns 0 on success, non-zero on failure.

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

Local to this translation unit.

Definition at line 3895 of file setup.c.

                                                                                                                                               {
    PetscErrorCode ierr;
    PetscMPIInt rank;
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
 
    // --- RNG for i-logical dimension ---
    ierr = PetscRandomCreate(PETSC_COMM_SELF, rand_logic_i); CHKERRQ(ierr);
    ierr = PetscRandomSetType((*rand_logic_i), PETSCRAND48); CHKERRQ(ierr);
    ierr = PetscRandomSetInterval(*rand_logic_i, 0.0, 1.0); CHKERRQ(ierr); // Key change: [0,1)
    /* Offsets reproduce the historical seeds (202401..202403 plus rank) at the default base. */
    ierr = PetscRandomSetSeed(*rand_logic_i, (unsigned long)base_seed + 190056UL + (unsigned long)rank); CHKERRQ(ierr);
    ierr = PetscRandomSeed(*rand_logic_i); CHKERRQ(ierr);
    LOG_ALLOW(LOCAL,LOG_VERBOSE, "[Rank %d] Initialized RNG for i-logical dimension [0,1).\n",rank);
 
    // --- RNG for j-logical dimension ---
    ierr = PetscRandomCreate(PETSC_COMM_SELF, rand_logic_j); CHKERRQ(ierr);
    ierr = PetscRandomSetType((*rand_logic_j), PETSCRAND48); CHKERRQ(ierr);
    ierr = PetscRandomSetInterval(*rand_logic_j, 0.0, 1.0); CHKERRQ(ierr); // Key change: [0,1)
    ierr = PetscRandomSetSeed(*rand_logic_j, (unsigned long)base_seed + 190057UL + (unsigned long)rank); CHKERRQ(ierr);
    ierr = PetscRandomSeed(*rand_logic_j); CHKERRQ(ierr);
    LOG_ALLOW(LOCAL,LOG_VERBOSE, "[Rank %d] Initialized RNG for j-logical dimension [0,1).\n",rank);
 
    // --- RNG for k-logical dimension ---
    ierr = PetscRandomCreate(PETSC_COMM_SELF, rand_logic_k); CHKERRQ(ierr);
    ierr = PetscRandomSetType((*rand_logic_k), PETSCRAND48); CHKERRQ(ierr);
    ierr = PetscRandomSetInterval(*rand_logic_k, 0.0, 1.0); CHKERRQ(ierr); // Key change: [0,1)
    ierr = PetscRandomSetSeed(*rand_logic_k, (unsigned long)base_seed + 190058UL + (unsigned long)rank); CHKERRQ(ierr);
    ierr = PetscRandomSeed(*rand_logic_k); CHKERRQ(ierr);
    LOG_ALLOW(LOCAL,LOG_VERBOSE, "[Rank %d] Initialized RNG for k-logical dimension [0,1).\n",rank);
 
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

InitializeBrownianRNG()

PetscErrorCode InitializeBrownianRNG

(

SimCtx *

simCtx

)

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

Configures it for Uniform [0, 1) which is required for Box-Muller transformation. Seeded from particleRandomSeed, the rank, and StartStep, so identical inputs on the same rank count reproduce, and a restart draws a fresh sequence.

Parameters

[in,out]

simCtx

Pointer to the Simulation Context.

Returns

PetscErrorCode

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

Local to this translation unit.

Definition at line 3940 of file setup.c.

                                                     {
    PetscErrorCode ierr;
    PetscMPIInt rank;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
 
    // 1. Create the generator (stored in SimCtx, not UserCtx, as it is global physics)
    ierr = PetscRandomCreate(PETSC_COMM_WORLD, &simCtx->BrownianMotionRNG); CHKERRQ(ierr);
    ierr = PetscRandomSetType(simCtx->BrownianMotionRNG, PETSCRAND48); CHKERRQ(ierr);
 
    // 2. CRITICAL: Set interval to [0, 1). 
    // This is required for the Gaussian math to work.
    ierr = PetscRandomSetInterval(simCtx->BrownianMotionRNG, 0.0, 1.0); CHKERRQ(ierr);
 
    // 3. Seed from the configured base, the rank, and the start step. It was seeded from
    //    the wall clock, so no two runs of the same inputs agreed. The start step keeps a
    //    restart from replaying the increments the first segment already drew, which would
    //    correlate them; the same restart is still reproducible.
    unsigned long seed = (unsigned long)simCtx->particleRandomSeed
                       + (unsigned long)rank * 987654321UL
                       + (unsigned long)simCtx->StartStep * 1000003UL;
    ierr = PetscRandomSetSeed(simCtx->BrownianMotionRNG, seed); CHKERRQ(ierr);
    ierr = PetscRandomSeed(simCtx->BrownianMotionRNG); CHKERRQ(ierr);
 
    LOG_ALLOW(LOCAL, LOG_VERBOSE, "[Rank %d] Initialized Brownian Physics RNG.\n", rank);
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

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
	)

Transforms scalar derivatives from computational space to physical space.

using the chain rule. Formula: dPhi/dx = J * ( dPhi/dCsi * dCsi/dx + dPhi/dEta * dEta/dx + ... )

Parameters

	jacobian	Local mapping Jacobian determinant.
	csi_metrics	Physical gradient of the xi coordinate.
	eta_metrics	Physical gradient of the eta coordinate.
	zet_metrics	Physical gradient of the zeta coordinate.
	dPhi_dcsi	Scalar derivative in the xi direction.
	dPhi_deta	Scalar derivative in the eta direction.
	dPhi_dzet	Scalar derivative in the zeta direction.
[out]	gradPhi	Physical Cartesian scalar gradient.

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 3983 of file setup.c.

{
    // Gradient X component
    gradPhi->x = jacobian * (dPhi_dcsi * csi_metrics.x + dPhi_deta * eta_metrics.x + dPhi_dzet * zet_metrics.x);
    
    // Gradient Y component
    gradPhi->y = jacobian * (dPhi_dcsi * csi_metrics.y + dPhi_deta * eta_metrics.y + dPhi_dzet * zet_metrics.y);
    
    // Gradient Z component
    gradPhi->z = jacobian * (dPhi_dcsi * csi_metrics.z + dPhi_deta * eta_metrics.z + dPhi_dzet * zet_metrics.z);
}

Here is the caller graph for this function:

ComputeScalarFieldDerivatives()

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

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

Parameters

	user	Block metric context.
	i	Cell index in the xi direction.
	j	Cell index in the eta direction.
	k	Cell index in the zeta direction.
	field_data	Ghosted scalar values indexed as [k][j][i].
[out]	grad	Physical Cartesian gradient at the requested cell.

Returns

PetscErrorCode 0 on success.

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

Local to this translation unit.

Definition at line 4031 of file setup.c.

{
    PetscErrorCode ierr;
    Cmpnts    ***csi, ***eta, ***zet;
    PetscReal ***jac;
    PetscReal d_csi, d_eta, d_zet;
 
    PetscFunctionBeginUser;
 
    // 1. Get read-only access to metrics
    ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, &csi); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lEta, &eta); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lZet, &zet); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->da,  user->lAj,  &jac); CHKERRQ(ierr);
 
    // 2. Compute derivatives in computational space (Central Difference)
    //    Assumes ghosts are available at i+/-1
    d_csi = 0.5 * (field_data[k][j][i+1] - field_data[k][j][i-1]);
    d_eta = 0.5 * (field_data[k][j+1][i] - field_data[k][j-1][i]);
    d_zet = 0.5 * (field_data[k+1][j][i] - field_data[k-1][j][i]);
 
    // 3. Transform to physical space
    TransformScalarDerivativesToPhysical(jac[k][j][i], 
                                         csi[k][j][i], eta[k][j][i], zet[k][j][i],
                                         d_csi, d_eta, d_zet,
                                         grad);
 
    // 4. Restore arrays
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, &csi); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, &eta); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, &zet); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->da,  user->lAj,  &jac); CHKERRQ(ierr);
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

ComputeVectorFieldDerivatives()

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

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

This function orchestrates the calculation of spatial derivatives. It first computes the derivatives in computational space (d/dcsi, d/deta, d/dzet) using a central difference scheme and then transforms them into physical space (d/dx, d/dy, d/dz).

Parameters

	user	Block metric context.
	i	Cell index in the xi direction.
	j	Cell index in the eta direction.
	k	Cell index in the zeta direction.
	field_data	Ghosted Cartesian vector values indexed as [k][j][i].
[out]	dudx	Gradient of the x velocity component.
[out]	dvdx	Gradient of the y velocity component.
[out]	dwdx	Gradient of the z velocity component.

Returns

PetscErrorCode 0 on success.

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

Local to this translation unit.

Definition at line 4074 of file setup.c.

{
    PetscErrorCode ierr;
    Cmpnts ***csi, ***eta, ***zet;
    PetscReal ***jac;
    PetscFunctionBeginUser;
 
    // 1. Get read-only access to the necessary metric data arrays
    ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, &csi); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lEta, &eta); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lZet, &zet); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->da,  user->lAj,  &jac); CHKERRQ(ierr);
 
    // 2. Calculate derivatives in computational space using central differencing
    Cmpnts deriv_csi, deriv_eta, deriv_zet;
    deriv_csi.x = (field_data[k][j][i+1].x - field_data[k][j][i-1].x) * 0.5;
    deriv_csi.y = (field_data[k][j][i+1].y - field_data[k][j][i-1].y) * 0.5;
    deriv_csi.z = (field_data[k][j][i+1].z - field_data[k][j][i-1].z) * 0.5;
 
    deriv_eta.x = (field_data[k][j+1][i].x - field_data[k][j-1][i].x) * 0.5;
    deriv_eta.y = (field_data[k][j+1][i].y - field_data[k][j-1][i].y) * 0.5;
    deriv_eta.z = (field_data[k][j+1][i].z - field_data[k][j-1][i].z) * 0.5;
 
    deriv_zet.x = (field_data[k+1][j][i].x - field_data[k-1][j][i].x) * 0.5;
    deriv_zet.y = (field_data[k+1][j][i].y - field_data[k-1][j][i].y) * 0.5;
    deriv_zet.z = (field_data[k+1][j][i].z - field_data[k-1][j][i].z) * 0.5;
 
    // 3. Transform derivatives to physical space
    TransformDerivativesToPhysical(jac[k][j][i], csi[k][j][i], eta[k][j][i], zet[k][j][i],
                                   deriv_csi, deriv_eta, deriv_zet,
                                   dudx, dvdx, dwdx);
 
    // 4. Restore access to the PETSc data arrays
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, &csi); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, &eta); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, &zet); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->da,  user->lAj,  &jac); CHKERRQ(ierr);
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

LESConfigSetDefaults()

PetscErrorCode LESConfigSetDefaults

(

LESConfig *

config

)

Installs the shipped defaults for every LES closure parameter.

One definition of the defaults, shared by the control-file ingress and by test fixtures, so the two cannot drift apart. It lives with the rest of the context construction because defaults are ingress, not closure behaviour, and because the postprocessor builds a simulation context without linking the LES module.

Local averaging and clamped limiting are the defaults because they are the only settings valid on arbitrary geometry; cases with homogeneous directions or an interest in backscatter opt into the alternatives.

Parameters

[out]

config

Configuration block to populate.

Returns

PetscErrorCode 0 on success, or PETSC_ERR_ARG_NULL for a null argument.

Installs the shipped defaults for every LES closure parameter.

Full API contract is documented with the header declaration in include/setup.h.

See also

LESConfigSetDefaults()

Definition at line 206 of file setup.c.

{
    PetscFunctionBeginUser;
    PetscCheck(config != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
               "LES configuration destination cannot be NULL.");
 
    config->dynamic_frequency       = 1;
    config->constant_cs             = 0.03;
    // Vreman (2004) recommends c = 2.5 Cs^2, 0.07 for Cs = 0.17; Nicoud & Ducros (1999)
    // calibrate C_w near 0.5-0.6 against the same isotropic decay. Both are calibrations
    // from isotropic turbulence, not universal values.
    config->vreman_coefficient      = 0.07;
    config->wale_coefficient        = 0.5;
    config->filter_width_model      = LES_FILTER_WIDTH_CUBE_ROOT_VOLUME;
    config->test_filter_kernel      = LES_TEST_FILTER_VOLUME_WEIGHTED_BOX;
    config->test_filter_width_ratio = 2.0;
    // Local averaging makes no homogeneity claim, so it is the only default that is
    // correct on every geometry. Homogeneous and global averaging are opted into.
    config->averaging_mode          = LES_AVERAGING_LOCAL;
    config->averaging_direction[0]  = PETSC_FALSE;
    config->averaging_direction[1]  = PETSC_FALSE;
    config->averaging_direction[2]  = PETSC_FALSE;
    config->clip_mode               = LES_CLIP_CLAMP;
    // Roughly twice the physical Smagorinsky constant: high enough to catch a
    // diverging coefficient, low enough that it does not shape the distribution.
    config->max_cs                  = 0.3;
    config->min_viscosity_ratio     = 0.0;
    config->yoshizawa_ci            = 0.09;
    config->diagnostics_enabled     = PETSC_FALSE;
    config->diagnostics_cadence     = 1;
 
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

DestroyUserVectors()

PetscErrorCode DestroyUserVectors

(

UserCtx *

user

)

Destroys all PETSc Vec objects within a single UserCtx structure.

This helper function systematically destroys the Vec objects stored in a UserCtx. The vectors are organized into functional groups for clarity:

• Primary flow fields (Ucont, Ucat, P, Nvert)
• Solver work vectors (Phi)
• Time-stepping vectors (Ucont_o, Ucont_rm1, etc.)
• Grid metrics (cell-centered, face-centered, coordinates)
• Turbulence vectors (Nu_t, CS, lFriction_Velocity)
• Particle vectors (ParticleCount, Psi)
• Boundary condition vectors (Ubcs, Uch)
• Post-processing vectors (P_nodal, Qcrit)
• And more...

All destroys are protected with NULL checks to handle conditional allocations safely.

Parameters

[in,out]

user

Pointer to the UserCtx containing the vectors to destroy.

Returns

PetscErrorCode 0 on success.

Destroys all PETSc Vec objects within a single UserCtx structure.

Local to this translation unit.

Definition at line 4128 of file setup.c.

{
    PetscErrorCode ierr;
    PetscFunctionBeginUser;
 
    // --- Group A: Primary Flow Fields (Always allocated at all levels) ---
    if (user->Ucont) { ierr = VecDestroy(&user->Ucont); CHKERRQ(ierr); }
    if (user->lUcont) { ierr = VecDestroy(&user->lUcont); CHKERRQ(ierr); }
    if (user->Ucat) { ierr = VecDestroy(&user->Ucat); CHKERRQ(ierr); }
    if (user->lUcat) { ierr = VecDestroy(&user->lUcat); CHKERRQ(ierr); }
    if (user->P) { ierr = VecDestroy(&user->P); CHKERRQ(ierr); }
    if (user->lP) { ierr = VecDestroy(&user->lP); CHKERRQ(ierr); }
    if (user->Nvert) { ierr = VecDestroy(&user->Nvert); CHKERRQ(ierr); }
    if (user->lNvert) { ierr = VecDestroy(&user->lNvert); CHKERRQ(ierr); }
 
    // --- Group A2: Derived Flow Fields (Conditional) ---
    if(user->Diffusivity) {ierr = VecDestroy(&user->Diffusivity); CHKERRQ(ierr);}
    if(user->lDiffusivity){ierr = VecDestroy(&user->lDiffusivity); CHKERRQ(ierr);}
    if(user->DiffusivityGradient){ierr = VecDestroy(&user->DiffusivityGradient); CHKERRQ(ierr);}
    if(user->lDiffusivityGradient){ierr = VecDestroy(&user->lDiffusivityGradient); CHKERRQ(ierr);}
    
    // --- Group B: Solver Work Vectors (All levels) ---
    if (user->Phi) { ierr = VecDestroy(&user->Phi); CHKERRQ(ierr); }
    if (user->lPhi) { ierr = VecDestroy(&user->lPhi); CHKERRQ(ierr); }
 
    // --- Group C: Time-Stepping Vectors (Finest level only) ---
    if (user->Ucont_o) { ierr = VecDestroy(&user->Ucont_o); CHKERRQ(ierr); }
    if (user->Ucont_rm1) { ierr = VecDestroy(&user->Ucont_rm1); CHKERRQ(ierr); }
    if (user->Ucat_o) { ierr = VecDestroy(&user->Ucat_o); CHKERRQ(ierr); }
    if (user->P_o) { ierr = VecDestroy(&user->P_o); CHKERRQ(ierr); }
    if (user->Nvert_o) { ierr = VecDestroy(&user->Nvert_o); CHKERRQ(ierr); }
    if (user->lUcont_o) { ierr = VecDestroy(&user->lUcont_o); CHKERRQ(ierr); }
    if (user->lUcont_rm1) { ierr = VecDestroy(&user->lUcont_rm1); CHKERRQ(ierr); }
    if (user->lNvert_o) { ierr = VecDestroy(&user->lNvert_o); CHKERRQ(ierr); }
 
    // --- Group D: Grid Metrics - Face Centered (All levels) ---
    if (user->Csi) { ierr = VecDestroy(&user->Csi); CHKERRQ(ierr); }
    if (user->Eta) { ierr = VecDestroy(&user->Eta); CHKERRQ(ierr); }
    if (user->Zet) { ierr = VecDestroy(&user->Zet); CHKERRQ(ierr); }
    if (user->Aj) { ierr = VecDestroy(&user->Aj); CHKERRQ(ierr); }
    if (user->lCsi) { ierr = VecDestroy(&user->lCsi); CHKERRQ(ierr); }
    if (user->lEta) { ierr = VecDestroy(&user->lEta); CHKERRQ(ierr); }
    if (user->lZet) { ierr = VecDestroy(&user->lZet); CHKERRQ(ierr); }
    if (user->lAj) { ierr = VecDestroy(&user->lAj); CHKERRQ(ierr); }
 
    // --- Group E: Grid Metrics - Face Centered (All levels) ---
    if (user->ICsi) { ierr = VecDestroy(&user->ICsi); CHKERRQ(ierr); }
    if (user->IEta) { ierr = VecDestroy(&user->IEta); CHKERRQ(ierr); }
    if (user->IZet) { ierr = VecDestroy(&user->IZet); CHKERRQ(ierr); }
    if (user->JCsi) { ierr = VecDestroy(&user->JCsi); CHKERRQ(ierr); }
    if (user->JEta) { ierr = VecDestroy(&user->JEta); CHKERRQ(ierr); }
    if (user->JZet) { ierr = VecDestroy(&user->JZet); CHKERRQ(ierr); }
    if (user->KCsi) { ierr = VecDestroy(&user->KCsi); CHKERRQ(ierr); }
    if (user->KEta) { ierr = VecDestroy(&user->KEta); CHKERRQ(ierr); }
    if (user->KZet) { ierr = VecDestroy(&user->KZet); CHKERRQ(ierr); }
    if (user->IAj) { ierr = VecDestroy(&user->IAj); CHKERRQ(ierr); }
    if (user->JAj) { ierr = VecDestroy(&user->JAj); CHKERRQ(ierr); }
    if (user->KAj) { ierr = VecDestroy(&user->KAj); CHKERRQ(ierr); }
    if (user->lICsi) { ierr = VecDestroy(&user->lICsi); CHKERRQ(ierr); }
    if (user->lIEta) { ierr = VecDestroy(&user->lIEta); CHKERRQ(ierr); }
    if (user->lIZet) { ierr = VecDestroy(&user->lIZet); CHKERRQ(ierr); }
    if (user->lJCsi) { ierr = VecDestroy(&user->lJCsi); CHKERRQ(ierr); }
    if (user->lJEta) { ierr = VecDestroy(&user->lJEta); CHKERRQ(ierr); }
    if (user->lJZet) { ierr = VecDestroy(&user->lJZet); CHKERRQ(ierr); }
    if (user->lKCsi) { ierr = VecDestroy(&user->lKCsi); CHKERRQ(ierr); }
    if (user->lKEta) { ierr = VecDestroy(&user->lKEta); CHKERRQ(ierr); }
    if (user->lKZet) { ierr = VecDestroy(&user->lKZet); CHKERRQ(ierr); }
    if (user->lIAj) { ierr = VecDestroy(&user->lIAj); CHKERRQ(ierr); }
    if (user->lJAj) { ierr = VecDestroy(&user->lJAj); CHKERRQ(ierr); }
    if (user->lKAj) { ierr = VecDestroy(&user->lKAj); CHKERRQ(ierr); }
 
    // --- Group F: Cell/Face Coordinates and Grid Spacing (All levels) ---
    if (user->Cent) { ierr = VecDestroy(&user->Cent); CHKERRQ(ierr); }
    if (user->lCent) { ierr = VecDestroy(&user->lCent); CHKERRQ(ierr); }
    if (user->GridSpace) { ierr = VecDestroy(&user->GridSpace); CHKERRQ(ierr); }
    if (user->lGridSpace) { ierr = VecDestroy(&user->lGridSpace); CHKERRQ(ierr); }
    if (user->Centx) { ierr = VecDestroy(&user->Centx); CHKERRQ(ierr); }
    if (user->Centy) { ierr = VecDestroy(&user->Centy); CHKERRQ(ierr); }
    if (user->Centz) { ierr = VecDestroy(&user->Centz); CHKERRQ(ierr); }
    if (user->lCentx) { ierr = VecDestroy(&user->lCentx); CHKERRQ(ierr); }
    if (user->lCenty) { ierr = VecDestroy(&user->lCenty); CHKERRQ(ierr); }
    if (user->lCentz) { ierr = VecDestroy(&user->lCentz); CHKERRQ(ierr); }
 
    // --- Group G: Turbulence Model Vectors (Finest level, conditional on les) ---
    if (user->Nu_t) { ierr = VecDestroy(&user->Nu_t); CHKERRQ(ierr); }
    if (user->lNu_t) { ierr = VecDestroy(&user->lNu_t); CHKERRQ(ierr); }
    if (user->CS) { ierr = VecDestroy(&user->CS); CHKERRQ(ierr); }
    if (user->lCs) { ierr = VecDestroy(&user->lCs); CHKERRQ(ierr); }
    if (user->Nu_Wall) { ierr = VecDestroy(&user->Nu_Wall); CHKERRQ(ierr); }
    if (user->lNu_Wall) { ierr = VecDestroy(&user->lNu_Wall); CHKERRQ(ierr); }
    if (user->Friction_Velocity) { ierr = VecDestroy(&user->Friction_Velocity); CHKERRQ(ierr); }
    if (user->lFriction_Velocity) { ierr = VecDestroy(&user->lFriction_Velocity); CHKERRQ(ierr); }
 
    // --- Group H: Particle Vectors (Finest level, conditional on np > 0) ---
    if (user->ParticleCount) { ierr = VecDestroy(&user->ParticleCount); CHKERRQ(ierr); }
    if (user->lParticleCount) { ierr = VecDestroy(&user->lParticleCount); CHKERRQ(ierr); }
    if (user->Psi) { ierr = VecDestroy(&user->Psi); CHKERRQ(ierr); }
    if (user->lPsi) { ierr = VecDestroy(&user->lPsi); CHKERRQ(ierr); }
 
    // --- Group I: Boundary Condition Vectors (All levels) ---
    if (user->Bcs.Ubcs) { ierr = VecDestroy(&user->Bcs.Ubcs); CHKERRQ(ierr); }
    if (user->Bcs.Uch) { ierr = VecDestroy(&user->Bcs.Uch); CHKERRQ(ierr); }
 
    // --- Group J: Post-Processing Vectors (Finest level, postprocessor mode) ---
    if (user->P_nodal) { ierr = VecDestroy(&user->P_nodal); CHKERRQ(ierr); }
    if (user->Ucat_nodal) { ierr = VecDestroy(&user->Ucat_nodal); CHKERRQ(ierr); }
    if (user->Qcrit) { ierr = VecDestroy(&user->Qcrit); CHKERRQ(ierr); }
    if (user->lQcrit) { ierr = VecDestroy(&user->lQcrit); CHKERRQ(ierr); }
    if (user->Qcrit_nodal) { ierr = VecDestroy(&user->Qcrit_nodal); CHKERRQ(ierr); }
    if (user->Psi_nodal) { ierr = VecDestroy(&user->Psi_nodal); CHKERRQ(ierr); }
 
    // --- Group K: Interpolation Vectors (Lazy allocation) ---
    if (user->fieldStatisticsStorage) {
        SimCtx *stats_ctx = user->simCtx;
        const PetscInt window_count = stats_ctx ? stats_ctx->fieldStatisticsWindowCount : 0;
        for (PetscInt w = 0; w < window_count; ++w) {
            ierr = PicurvWindowStorageDestroy(&user->fieldStatisticsStorage[w]); CHKERRQ(ierr);
        }
        ierr = PetscFree(user->fieldStatisticsStorage); CHKERRQ(ierr);
        user->fieldStatisticsStorage = NULL;
    }
    if (user->CellScalarAtCorner) { ierr = VecDestroy(&user->CellScalarAtCorner); CHKERRQ(ierr); }
    if (user->lCellScalarAtCorner) { ierr = VecDestroy(&user->lCellScalarAtCorner); CHKERRQ(ierr); }
    if (user->CellVectorAtCorner) { ierr = VecDestroy(&user->CellVectorAtCorner); CHKERRQ(ierr); }
    if (user->PostScalar) { ierr = VecDestroy(&user->PostScalar); CHKERRQ(ierr); }
    if (user->lPostScalar) { ierr = VecDestroy(&user->lPostScalar); CHKERRQ(ierr); }
    if (user->PostScalarNodal) { ierr = VecDestroy(&user->PostScalarNodal); CHKERRQ(ierr); }
    if (user->PostVector) { ierr = VecDestroy(&user->PostVector); CHKERRQ(ierr); }
    if (user->lPostVector) { ierr = VecDestroy(&user->lPostVector); CHKERRQ(ierr); }
    if (user->PostVectorNodal) { ierr = VecDestroy(&user->PostVectorNodal); CHKERRQ(ierr); }
    if (user->lCellVectorAtCorner) { ierr = VecDestroy(&user->lCellVectorAtCorner); CHKERRQ(ierr); }
 
    // --- Group L: Implicit Solver Temporary Vectors (Destroyed after use, but check anyway) ---
    if (user->Rhs) { ierr = VecDestroy(&user->Rhs); CHKERRQ(ierr); }
    if (user->dUcont) { ierr = VecDestroy(&user->dUcont); CHKERRQ(ierr); }
    if (user->pUcont) { ierr = VecDestroy(&user->pUcont); CHKERRQ(ierr); }
 
    // --- Group M: Poisson Solver Vectors (kept for the run by PoissonSolver_Multigrid) ---
    if (user->B) { ierr = VecDestroy(&user->B); CHKERRQ(ierr); }
    if (user->R) { ierr = VecDestroy(&user->R); CHKERRQ(ierr); }
 
    LOG_ALLOW(LOCAL, LOG_DEBUG, "All vectors destroyed for UserCtx.\n");
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

DestroyUserContext()

PetscErrorCode DestroyUserContext

(

UserCtx *

user

)

Destroys all resources allocated within a single UserCtx structure.

This function cleans up all memory and PETSc objects associated with a single UserCtx (grid level). It calls the helper functions and destroys remaining objects in the proper dependency order:

1. Boundary conditions (handlers and their data)
2. All PETSc vectors (via DestroyUserVectors)
3. Matrix and solver objects (A, MR, MP, ksp, nullsp)
4. Application ordering (AO)
5. Distributed mesh objects (DMs) - most derived first
6. Raw PetscMalloc'd arrays (RankCellInfoMap)

This function should be called for each UserCtx in the multigrid hierarchy.

Parameters

[in,out]

user

Pointer to the UserCtx to be destroyed.

Returns

PetscErrorCode 0 on success.

Destroys all resources allocated within a single UserCtx structure.

Local to this translation unit.

Definition at line 4278 of file setup.c.

{
    PetscErrorCode ierr;
    PetscFunctionBeginUser;
 
    if (!user) {
        LOG_ALLOW(LOCAL, LOG_WARNING, "DestroyUserContext called with NULL user pointer.\n");
        PetscFunctionReturn(0);
    }
 
    LOG_ALLOW(LOCAL, LOG_INFO, "Destroying UserCtx at level %d...\n", user->thislevel);
 
    // --- Step 1: Destroy Boundary Condition System ---
    // This handles all BC handlers and their private data.
    ierr = BoundarySystem_Destroy(user); CHKERRQ(ierr);
    LOG_ALLOW(LOCAL, LOG_DEBUG, "  Boundary system destroyed.\n");
 
    // --- Step 2: Destroy All Vectors ---
    // Handles ~74 Vec objects with proper NULL checking.
    ierr = DestroyUserVectors(user); CHKERRQ(ierr);
    LOG_ALLOW(LOCAL, LOG_DEBUG, "  All vectors destroyed.\n");
 
    // --- Step 3: Destroy Matrix and Solver Objects ---
    // Destroy pressure-Poisson matrices and solver.
    if (user->A) {
        ierr = MatDestroy(&user->A); CHKERRQ(ierr);
        LOG_ALLOW(LOCAL, LOG_DEBUG, "  Matrix A destroyed.\n");
    }
    if (user->MR) {
        ierr = MatDestroy(&user->MR); CHKERRQ(ierr);
        LOG_ALLOW(LOCAL, LOG_DEBUG, "  Matrix MR destroyed.\n");
    }
    if (user->MP) {
        ierr = MatDestroy(&user->MP); CHKERRQ(ierr);
        LOG_ALLOW(LOCAL, LOG_DEBUG, "  Matrix MP destroyed.\n");
    }
    if (user->ksp) {
        ierr = KSPDestroy(&user->ksp); CHKERRQ(ierr);
        LOG_ALLOW(LOCAL, LOG_DEBUG, "  KSP solver destroyed.\n");
    }
    if (user->nullsp) {
        ierr = MatNullSpaceDestroy(&user->nullsp); CHKERRQ(ierr);
        LOG_ALLOW(LOCAL, LOG_DEBUG, "  MatNullSpace destroyed.\n");
    }
 
    // --- Step 4: Destroy Application Ordering ---
    if (user->ao) {
        ierr = AODestroy(&user->ao); CHKERRQ(ierr);
        LOG_ALLOW(LOCAL, LOG_DEBUG, "  AO destroyed.\n");
    }
 
    // --- Step 5: Destroy DM Objects ---
    // Destroy in reverse order of dependency: post_swarm, swarm, fda6, fda, da
    if (user->post_swarm) {
        ierr = DMDestroy(&user->post_swarm); CHKERRQ(ierr);
        LOG_ALLOW(LOCAL, LOG_DEBUG, "  post_swarm DM destroyed.\n");
    }
    if (user->swarm) {
        ierr = DMDestroy(&user->swarm); CHKERRQ(ierr);
        LOG_ALLOW(LOCAL, LOG_DEBUG, "  swarm DM destroyed.\n");
    }
    if (user->fda6) {
        ierr = DMDestroy(&user->fda6); CHKERRQ(ierr);
        LOG_ALLOW(LOCAL, LOG_DEBUG, "  fda6 DM destroyed.\n");
    }
    if (user->da) {
        ierr = DMDestroy(&user->da); CHKERRQ(ierr);
        LOG_ALLOW(LOCAL, LOG_DEBUG, "  da DM destroyed.\n");
    }
 
    // --- Step 6: Free PetscMalloc'd Arrays ---
    // Free arrays allocated with PetscMalloc1
    if (user->RankCellInfoMap) {
        ierr = PetscFree(user->RankCellInfoMap); CHKERRQ(ierr);
        user->RankCellInfoMap = NULL;
        LOG_ALLOW(LOCAL, LOG_DEBUG, "  RankCellInfoMap freed.\n");
    }
 
    LOG_ALLOW(LOCAL, LOG_INFO, "UserCtx at level %d fully destroyed.\n", user->thislevel);
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

FinalizeSimulation()

PetscErrorCode FinalizeSimulation

(

SimCtx *

simCtx

)

Main cleanup function for the entire simulation context.

This function is responsible for destroying ALL memory and PETSc objects allocated during the simulation, including:

• All UserCtx structures in the multigrid hierarchy (via DestroyUserContext)
• The multigrid management structures (UserMG, MGCtx array)
• All SimCtx-level objects (logviewer, dm_swarm, bboxlist, string arrays, etc.)
• The SimCtx allocation itself.

This function should be called ONCE at the end of the simulation, after all computation is complete, but BEFORE PetscFinalize().

Call order in main:

1. [Simulation runs]
2. ProfilingFinalize(simCtx);
3. FinalizeSimulation(simCtx); <- This function
4. PetscFinalize();

Parameters

[in,out]

simCtx

Pointer to the master SimulationContext to be destroyed. The pointer is invalid after this call returns.

Returns

PetscErrorCode 0 on success.

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 4368 of file setup.c.

{
    PetscErrorCode ierr;
    PetscFunctionBeginUser;
 
    if (!simCtx) {
        LOG_ALLOW(GLOBAL, LOG_WARNING, "FinalizeSimulation called with NULL SimCtx pointer.\n");
        PetscFunctionReturn(0);
    }
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "========================================\n");
    LOG_ALLOW(GLOBAL, LOG_INFO, "Beginning simulation memory cleanup...\n");
    LOG_ALLOW(GLOBAL, LOG_INFO, "========================================\n");
 
    // ============================================================================
    // PHASE 1: DESTROY MULTIGRID HIERARCHY (All UserCtx structures)
    // ============================================================================
 
    ierr = DestroySolutionConvergenceState(simCtx); CHKERRQ(ierr);
    ierr = DestroyFieldStatisticsConfig(simCtx); CHKERRQ(ierr);
    ierr = ParticleFieldPlanDestroy(&simCtx->particleFieldPlan); CHKERRQ(ierr);
 
    if (simCtx->usermg.mgctx) {
        LOG_ALLOW(GLOBAL, LOG_INFO, "Destroying multigrid hierarchy (%d levels)...\n",
                  simCtx->usermg.mglevels);
 
        // Destroy each UserCtx from finest to coarsest (reverse order is safer)
        for (PetscInt level = simCtx->usermg.mglevels - 1; level >= 0; level--) {
            UserCtx *user = simCtx->usermg.mgctx[level].user;
            if (user) {
                LOG_ALLOW(LOCAL, LOG_INFO, "  Destroying level %d of %d...\n",
                          level, simCtx->usermg.mglevels - 1);
                ierr = DestroyUserContext(user); CHKERRQ(ierr);
 
                // Free the UserCtx structure itself
                ierr = PetscFree(user); CHKERRQ(ierr);
                simCtx->usermg.mgctx[level].user = NULL;
            }
 
            // Destroy the MGCtx-level packer DM
            if (simCtx->usermg.mgctx[level].packer) {
                ierr = DMDestroy(&simCtx->usermg.mgctx[level].packer); CHKERRQ(ierr);
                LOG_ALLOW(LOCAL, LOG_DEBUG, "  MGCtx[%d].packer destroyed.\n", level);
            }
        }
 
        // Free the MGCtx array itself
        ierr = PetscFree(simCtx->usermg.mgctx); CHKERRQ(ierr);
        simCtx->usermg.mgctx = NULL;
        LOG_ALLOW(GLOBAL, LOG_INFO, "All multigrid levels destroyed.\n");
    }
 
    // ============================================================================
    // PHASE 2: DESTROY USERMG-LEVEL OBJECTS
    // ============================================================================
 
    if (simCtx->usermg.packer) {
        ierr = DMDestroy(&simCtx->usermg.packer); CHKERRQ(ierr);
        LOG_ALLOW(LOCAL, LOG_DEBUG, "UserMG.packer DM destroyed.\n");
    }
 
    if (simCtx->usermg.snespacker) {
        ierr = SNESDestroy(&simCtx->usermg.snespacker); CHKERRQ(ierr);
        LOG_ALLOW(LOCAL, LOG_DEBUG, "UserMG.snespacker SNES destroyed.\n");
    }
 
    // ============================================================================
    // PHASE 3: DESTROY SIMCTX-LEVEL OBJECTS
    // ============================================================================
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Destroying SimCtx-level objects...\n");
 
    // --- PetscViewer for logging ---
    if (simCtx->logviewer) {
        ierr = PetscViewerDestroy(&simCtx->logviewer); CHKERRQ(ierr);
        LOG_ALLOW(LOCAL, LOG_DEBUG, "  logviewer destroyed.\n");
    }
 
    // --- Particle System DM ---
    if (simCtx->dm_swarm) {
        ierr = DMDestroy(&simCtx->dm_swarm); CHKERRQ(ierr);
        LOG_ALLOW(LOCAL, LOG_DEBUG, "  dm_swarm destroyed.\n");
    }
 
    // --- BoundingBox List (Array of BoundingBox structs) ---
    if (simCtx->bboxlist) {
        ierr = PetscFree(simCtx->bboxlist); CHKERRQ(ierr);
        simCtx->bboxlist = NULL;
        LOG_ALLOW(LOCAL, LOG_DEBUG, "  bboxlist freed.\n");
    }
 
    // --- Boundary Condition Files (Array of strings) ---
    if (simCtx->bcs_files) {
        for (PetscInt i = 0; i < simCtx->num_bcs_files; i++) {
            if (simCtx->bcs_files[i]) {
                ierr = PetscFree(simCtx->bcs_files[i]); CHKERRQ(ierr);
            }
        }
        ierr = PetscFree(simCtx->bcs_files); CHKERRQ(ierr);
        simCtx->bcs_files = NULL;
        LOG_ALLOW(LOCAL, LOG_DEBUG, "  bcs_files array freed (%d files).\n", simCtx->num_bcs_files);
    }
 
    // --- Brownian Motion RNG ---
    if (simCtx->BrownianMotionRNG) {
        ierr = PetscRandomDestroy(&simCtx->BrownianMotionRNG); CHKERRQ(ierr);
        LOG_ALLOW(LOCAL, LOG_DEBUG, "  BrownianMotionRNG destroyed.\n");
    }
    // --- Post-Processing Parameters ---
    // pps is allocated with PetscNew and contains only static char arrays and basic types.
    // No internal dynamic allocations need to be freed.
    if (simCtx->pps) {
        ierr = PetscFree(simCtx->pps); CHKERRQ(ierr);
        simCtx->pps = NULL;
        LOG_ALLOW(LOCAL, LOG_DEBUG, "  PostProcessParams freed.\n");
    }
 
    // --- IBM/FSI Objects ---
    // Note: These are initialized to NULL and currently have no dedicated destroy functions.
    // If these modules are extended with cleanup routines, call them here.
    if (simCtx->ibm != NULL) {
        LOG_ALLOW(GLOBAL, LOG_WARNING, "  WARNING: simCtx->ibm is non-NULL but no destroy function exists. Potential memory leak.\n");
    }
    if (simCtx->ibmv != NULL) {
        LOG_ALLOW(GLOBAL, LOG_WARNING, "  WARNING: simCtx->ibmv is non-NULL but no destroy function exists. Potential memory leak.\n");
    }
    if (simCtx->fsi != NULL) {
        LOG_ALLOW(GLOBAL, LOG_WARNING, "  WARNING: simCtx->fsi is non-NULL but no destroy function exists. Potential memory leak.\n");
    }
 
    // --- Logging Allowed Functions (Array of strings) ---
    // Note: The logging system maintains its own copy via set_allowed_functions(),
    // so freeing simCtx->allowedFuncs will NOT affect LOG_ALLOW functionality.
    if (simCtx->allowedFuncs) {
        for (PetscInt i = 0; i < simCtx->nAllowed; i++) {
            if (simCtx->allowedFuncs[i]) {
                ierr = PetscFree(simCtx->allowedFuncs[i]); CHKERRQ(ierr);
            }
        }
        ierr = PetscFree(simCtx->allowedFuncs); CHKERRQ(ierr);
        simCtx->allowedFuncs = NULL;
        LOG_ALLOW(LOCAL, LOG_DEBUG, "  allowedFuncs array freed (%d functions).\n", simCtx->nAllowed);
    }
 
    // --- Profiling Critical Functions (Array of strings) ---
    if (simCtx->profilingSelectedFuncs) {
        for (PetscInt i = 0; i < simCtx->nProfilingSelectedFuncs; i++) {
            if (simCtx->profilingSelectedFuncs[i]) {
                ierr = PetscFree(simCtx->profilingSelectedFuncs[i]); CHKERRQ(ierr);
            }
        }
        ierr = PetscFree(simCtx->profilingSelectedFuncs); CHKERRQ(ierr);
        simCtx->profilingSelectedFuncs = NULL;
        LOG_ALLOW(LOCAL, LOG_DEBUG, "  profilingSelectedFuncs array freed (%d functions).\n", simCtx->nProfilingSelectedFuncs);
    }
 
    // ============================================================================
    // PHASE 4: FINAL SUMMARY
    // ============================================================================
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "========================================\n");
    LOG_ALLOW(GLOBAL, LOG_INFO, "Simulation cleanup completed successfully.\n");
    LOG_ALLOW(GLOBAL, LOG_INFO, "All PETSc objects have been destroyed.\n");
    LOG_ALLOW(GLOBAL, LOG_INFO, "========================================\n");
 
    ierr = PetscFree(simCtx); CHKERRQ(ierr);
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function: