setup.c File Reference

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

#include <ctype.h>
#include <errno.h>
#include "setup.h"

Include dependency graph for setup.c:

Go to the source code of this file.

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

Functions
PetscBool	RuntimeWalltimeGuardParsePositiveSeconds (const char *text, PetscReal *seconds_out)
	Implementation of RuntimeWalltimeGuardParsePositiveSeconds().
PetscErrorCode	InitializeSolutionConvergenceState (SimCtx *simCtx)
	Implementation of InitializeSolutionConvergenceState().
PetscErrorCode	DestroySolutionConvergenceState (SimCtx *simCtx)
	Implementation of DestroySolutionConvergenceState().
PetscErrorCode	CreateSimulationContext (int argc, char **argv, SimCtx **p_simCtx)
	Implementation of CreateSimulationContext().
static PetscErrorCode	PetscMkdirRecursive (const char *path)
	Internal helper implementation: PetscMkdirRecursive().
PetscErrorCode	SetupSimulationEnvironment (SimCtx *simCtx)
	Internal helper implementation: SetupSimulationEnvironment().
static PetscErrorCode	AllocateContextHierarchy (SimCtx *simCtx)
	Internal helper implementation: AllocateContextHierarchy().
static PetscErrorCode	SetupSolverParameters (SimCtx *simCtx)
	Internal helper implementation: SetupSolverParameters().
PetscErrorCode	SetupGridAndSolvers (SimCtx *simCtx)
	Implementation of SetupGridAndSolvers().
PetscErrorCode	CreateAndInitializeAllVectors (SimCtx *simCtx)
	Internal helper implementation: CreateAndInitializeAllVectors().
static PetscErrorCode	RepairPeriodicNormalFaceGhosts (UserCtx *user, DM dm, Vec local_vec, PetscInt dof, char face_direction, PetscBool component_staggered)
	Repairs the adjacent normal ghost layer for periodic face-staggered data.
PetscErrorCode	UpdateLocalGhosts (UserCtx *user, const char *fieldName)
	Internal helper implementation: UpdateLocalGhosts().
PetscErrorCode	SetupBoundaryConditions (SimCtx *simCtx)
	Internal helper implementation: SetupBoundaryConditions().
PetscErrorCode	Allocate3DArrayScalar (PetscReal ****array, PetscInt nz, PetscInt ny, PetscInt nx)
	Internal helper implementation: Allocate3DArrayScalar().
PetscErrorCode	Deallocate3DArrayScalar (PetscReal ***array, PetscInt nz, PetscInt ny)
	Internal helper implementation: Deallocate3DArrayScalar().
PetscErrorCode	Allocate3DArrayVector (Cmpnts ****array, PetscInt nz, PetscInt ny, PetscInt nx)
	Implementation of Allocate3DArrayVector().
PetscErrorCode	Deallocate3DArrayVector (Cmpnts ***array, PetscInt nz, PetscInt ny)
	Implementation of Deallocate3DArrayVector().
PetscErrorCode	GetOwnedCellRange (const DMDALocalInfo *info_nodes, PetscInt dim, PetscInt *xs_cell_global_out, PetscInt *xm_cell_local_out)
	Internal helper implementation: GetOwnedCellRange().
PetscErrorCode	ComputeAndStoreNeighborRanks (UserCtx *user)
	Internal helper implementation: ComputeAndStoreNeighborRanks().
PetscErrorCode	SetDMDAProcLayout (DM dm, UserCtx *user)
	Internal helper implementation: SetDMDAProcLayout().
PetscErrorCode	SetupDomainRankInfo (SimCtx *simCtx)
	Implementation of SetupDomainRankInfo().
PetscErrorCode	Contra2Cart (UserCtx *user)
	Internal helper implementation: Contra2Cart().
PetscErrorCode	Cart2Contra (UserCtx *user)
	Convert a spatially varying Cartesian velocity field to contravariant fluxes.
PetscErrorCode	UniformCart2Contra (UserCtx *user, PetscReal u, PetscReal v, PetscReal w)
	Convert a uniform Cartesian velocity (u,v,w) to contravariant fluxes in Ucont.
PetscErrorCode	SetupDomainCellDecompositionMap (UserCtx *user)
	Internal helper implementation: SetupDomainCellDecompositionMap().
PetscErrorCode	BinarySearchInt64 (PetscInt n, const PetscInt64 arr[], PetscInt64 key, PetscBool *found)
	Implementation of BinarySearchInt64().
static PetscInt	Gidx (PetscInt i, PetscInt j, PetscInt k, UserCtx *user)
	Internal helper implementation: Gidx().
PetscErrorCode	ComputeDivergence (UserCtx *user)
	Implementation of ComputeDivergence().
PetscErrorCode	InitializeRandomGenerators (UserCtx *user, PetscRandom *randx, PetscRandom *randy, PetscRandom *randz)
	Implementation of InitializeRandomGenerators().
PetscErrorCode	InitializeLogicalSpaceRNGs (PetscRandom *rand_logic_i, PetscRandom *rand_logic_j, PetscRandom *rand_logic_k)
	Internal helper implementation: InitializeLogicalSpaceRNGs().
PetscErrorCode	InitializeBrownianRNG (SimCtx *simCtx)
	Internal helper implementation: InitializeBrownianRNG().
void	TransformScalarDerivativesToPhysical (PetscReal jacobian, Cmpnts csi_metrics, Cmpnts eta_metrics, Cmpnts zet_metrics, PetscReal dPhi_dcsi, PetscReal dPhi_deta, PetscReal dPhi_dzet, Cmpnts *gradPhi)
	Implementation of TransformScalarDerivativesToPhysical().
static void	TransformDerivativesToPhysical (PetscReal jacobian, Cmpnts csi_metrics, Cmpnts eta_metrics, Cmpnts zet_metrics, Cmpnts deriv_csi, Cmpnts deriv_eta, Cmpnts deriv_zet, Cmpnts *dudx, Cmpnts *dvdx, Cmpnts *dwdx)
	Internal helper implementation: TransformDerivativesToPhysical().
PetscErrorCode	ComputeScalarFieldDerivatives (UserCtx *user, PetscInt i, PetscInt j, PetscInt k, PetscReal ***field_data, Cmpnts *grad)
	Internal helper implementation: ComputeScalarFieldDerivatives().
PetscErrorCode	ComputeVectorFieldDerivatives (UserCtx *user, PetscInt i, PetscInt j, PetscInt k, Cmpnts ***field_data, Cmpnts *dudx, Cmpnts *dvdx, Cmpnts *dwdx)
	Internal helper implementation: ComputeVectorFieldDerivatives().
PetscErrorCode	DestroyUserVectors (UserCtx *user)
	Internal helper implementation: DestroyUserVectors().
PetscErrorCode	DestroyUserContext (UserCtx *user)
	Internal helper implementation: DestroyUserContext().
PetscErrorCode	FinalizeSimulation (SimCtx *simCtx)
	Implementation of FinalizeSimulation().

Detailed Description

// Setup code for running any simulation

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

Definition in file setup.c.

Macro Definition Documentation

__FUNCT__ [1/30]

#define __FUNCT__   "CreateSimulationContext"

Definition at line 143 of file setup.c.

__FUNCT__ [2/30]

#define __FUNCT__   "PetscMkdirRecursive"

Definition at line 143 of file setup.c.

__FUNCT__ [3/30]

#define __FUNCT__   "SetupSimulationEnvironment"

Definition at line 143 of file setup.c.

__FUNCT__ [4/30]

#define __FUNCT__   "AllocateContextHeirarchy"

Definition at line 143 of file setup.c.

__FUNCT__ [5/30]

#define __FUNCT__   "SetupSolverParameters"

Definition at line 143 of file setup.c.

__FUNCT__ [6/30]

#define __FUNCT__   "SetupGridAndSolvers"

Definition at line 143 of file setup.c.

__FUNCT__ [7/30]

#define __FUNCT__   "CreateAndInitializeAllVectors"

Definition at line 143 of file setup.c.

__FUNCT__ [8/30]

#define __FUNCT__   "RepairPeriodicNormalFaceGhosts"

Definition at line 143 of file setup.c.

__FUNCT__ [9/30]

#define __FUNCT__   "UpdateLocalGhosts"

Definition at line 143 of file setup.c.

__FUNCT__ [10/30]

#define __FUNCT__   "SetupBoundaryConditions"

Definition at line 143 of file setup.c.

__FUNCT__ [11/30]

#define __FUNCT__   "GetOwnedCellRange"

Definition at line 143 of file setup.c.

__FUNCT__ [12/30]

#define __FUNCT__   "ComputeAndStoreNeighborRanks"

Definition at line 143 of file setup.c.

__FUNCT__ [13/30]

#define __FUNCT__   "SetDMDAProcLayout"

Definition at line 143 of file setup.c.

__FUNCT__ [14/30]

#define __FUNCT__   "SetupDomainRankInfo"

Definition at line 143 of file setup.c.

__FUNCT__ [15/30]

#define __FUNCT__   "Contra2Cart"

Definition at line 143 of file setup.c.

__FUNCT__ [16/30]

#define __FUNCT__   "Cart2Contra"

Definition at line 143 of file setup.c.

__FUNCT__ [17/30]

#define __FUNCT__   "UniformCart2Contra"

Definition at line 143 of file setup.c.

__FUNCT__ [18/30]

#define __FUNCT__   "SetupDomainCellDecompositionMap"

Definition at line 143 of file setup.c.

__FUNCT__ [19/30]

#define __FUNCT__   "BinarySearchInt64"

Definition at line 143 of file setup.c.

__FUNCT__ [20/30]

#define __FUNCT__   "ComputeDivergence"

Definition at line 143 of file setup.c.

__FUNCT__ [21/30]

#define __FUNCT__   "InitializeRandomGenerators"

Definition at line 143 of file setup.c.

__FUNCT__ [22/30]

#define __FUNCT__   "InitializeLogicalSpaceRNGs"

Definition at line 143 of file setup.c.

__FUNCT__ [23/30]

#define __FUNCT__   "InitializeBrownianRNG"

Definition at line 143 of file setup.c.

__FUNCT__ [24/30]

#define __FUNCT__   "TransformScalarDerivativesToPhysical"

Definition at line 143 of file setup.c.

__FUNCT__ [25/30]

#define __FUNCT__   "TransformDerivativesToPhysical"

Definition at line 143 of file setup.c.

__FUNCT__ [26/30]

#define __FUNCT__   "ComputeScalarFieldDerivatives"

Definition at line 143 of file setup.c.

__FUNCT__ [27/30]

#define __FUNCT__   "ComputeVectorFieldDerivatives"

Definition at line 143 of file setup.c.

__FUNCT__ [28/30]

#define __FUNCT__   "DestroyUserVectors"

Definition at line 143 of file setup.c.

__FUNCT__ [29/30]

#define __FUNCT__   "DestroyUserContext"

Definition at line 143 of file setup.c.

__FUNCT__ [30/30]

#define __FUNCT__   "FinalizeSimulation"

Definition at line 143 of file setup.c.

Function Documentation

RuntimeWalltimeGuardParsePositiveSeconds()

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

Implementation of RuntimeWalltimeGuardParsePositiveSeconds().

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

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

See also

RuntimeWalltimeGuardParsePositiveSeconds()

Definition at line 18 of file setup.c.

{
    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:

InitializeSolutionConvergenceState()

PetscErrorCode InitializeSolutionConvergenceState

(

SimCtx *

simCtx

)

Implementation of InitializeSolutionConvergenceState().

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 47 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->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

)

Implementation of DestroySolutionConvergenceState().

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

{
    UserCtx *user = NULL;
 
    PetscFunctionBeginUser;
    if (!simCtx || !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:

CreateSimulationContext()

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

Implementation of CreateSimulationContext().

Allocates and populates the master SimulationContext object.

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

See also

CreateSimulationContext()

Definition at line 151 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->euler_subdir,"eulerian");
    strcpy(simCtx->particle_subdir,"particles");
    simCtx->_io_context_buffer[0] = '\0';
    simCtx->current_io_directory = NULL;
 
    // --- Group 3: High-Level Physics & Model Selection Flags ---
    simCtx->immersed = 0; simCtx->movefsi = 0; simCtx->rotatefsi = 0;
    simCtx->sediment = 0; simCtx->rheology = 0; simCtx->invicid = 0;
    simCtx->TwoD = 0; simCtx->thin = 0; simCtx->moveframe = 0;
    simCtx->rotateframe = 0; simCtx->blank = 0;
    simCtx->dgf_x = 0; simCtx->dgf_y = 1; simCtx->dgf_z = 0;
    simCtx->dgf_ax = 1; simCtx->dgf_ay = 0; simCtx->dgf_az = 0;
    strcpy(simCtx->AnalyticalSolutionType,"TGV3D");
 
    // --- Group 4: Specific Simulation Case Flags --- (DEPRICATED)
    simCtx->cop=0; simCtx->fish=0; simCtx->fish_c=0; simCtx->fishcyl=0;
    simCtx->eel=0; simCtx->pizza=0; simCtx->turbine=0; simCtx->Pipe=0;
    simCtx->wing=0; simCtx->hydro=0; simCtx->MHV=0; simCtx->LV=0;
    simCtx->channelz = 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;
    simCtx->mom_resid_atol = 0.0; simCtx->mom_resid_rtol = 0.0;
    simCtx->imp_stol = 1.e-8;
    simCtx->mglevels = 3; simCtx->mg_MAX_IT = 30; simCtx->mg_idx = 1;
    simCtx->mg_preItr = 1; simCtx->mg_poItr = 1;
    simCtx->poisson = 0; simCtx->poisson_tol = 5.e-9;
    simCtx->STRONG_COUPLING = 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->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->NumberOfBodies = 1; 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;
 
    // --- Group 7: Grid, Domain, and Boundary Condition Settings ---
    simCtx->block_number = 1; simCtx->inletprofile = 1;
    simCtx->grid1d = 0; simCtx->Ogrid = 0;
    simCtx->i_periodic = 0; simCtx->j_periodic = 0; simCtx->k_periodic = 0;
    simCtx->blkpbc = 10; 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->grid_rotation_angle  = 0.0;
    simCtx->Croty = 0.0; simCtx->Crotz = 0.0;
    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->targetVolumetricFlux  = 0.0;
    simCtx->bulkVelocityCorrection = 0.0;
    simCtx->boundaryVelocityCorrection = 0.0;
    simCtx->AreaInSum = 0.0; simCtx->AreaOutSum = 0.0;
    simCtx->U_bc = 0.0; simCtx->ccc = 0;
    simCtx->ratio = 0.0;
    
    
    // --- Group 8: Turbulence Modeling (LES/RANS) ---
    simCtx->averaging = PETSC_FALSE; simCtx->les = NO_LES_MODEL; simCtx->rans = 0;
    simCtx->wallfunction = 0; simCtx->mixed = 0; simCtx->clark = 0;
    simCtx->dynamic_freq = 1; simCtx->max_cs = 0.5;
    simCtx->Const_CS = 0.03;
    simCtx->testfilter_ik = 0; simCtx->testfilter_1d = 0;
    simCtx->i_homo_filter = 0; simCtx->j_homo_filter = 0; simCtx->k_homo_filter = 0;
 
    // --- Group 9: Particle / DMSwarm Data & Settings ---
    simCtx->np = 0; simCtx->readFields = PETSC_FALSE;
    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->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;
    simCtx->rstart_fsi = PETSC_FALSE; simCtx->duplicate = 0;
 
    // --- 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->summationRHS = 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);
 
    // === 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);
    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,"-euler_subdir",simCtx->euler_subdir,sizeof(simCtx->euler_subdir),NULL);CHKERRQ(ierr);
    ierr = PetscOptionsGetString(NULL,NULL,"-particle_subdir",simCtx->particle_subdir,sizeof(simCtx->particle_subdir),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, "-sediment", &simCtx->sediment, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-rheology", &simCtx->rheology, 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, "-thin", &simCtx->thin, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-mframe", &simCtx->moveframe, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-rframe", &simCtx->rotateframe, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-blk", &simCtx->blank, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-dgf_z", &simCtx->dgf_z, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-dgf_y", &simCtx->dgf_y, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-dgf_x", &simCtx->dgf_x, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-dgf_az", &simCtx->dgf_az, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-dgf_ay", &simCtx->dgf_ay, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-dgf_ax", &simCtx->dgf_ax, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetString(NULL,NULL,"-analytical_type",simCtx->AnalyticalSolutionType,sizeof(simCtx->AnalyticalSolutionType),NULL);CHKERRQ(ierr);
 
    //  --- Group 4
    LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 4: Specific Simulation Case Flags \n");
    ierr = PetscOptionsGetInt(NULL, NULL, "-cop", &simCtx->cop, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-fish", &simCtx->fish, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-pizza", &simCtx->pizza, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-turbine", &simCtx->turbine, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-fishcyl", &simCtx->fishcyl, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-eel", &simCtx->eel, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-cstart", &simCtx->fish_c, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-wing", &simCtx->wing, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-mhv", &simCtx->MHV, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-hydro", &simCtx->hydro, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-lv", &simCtx->LV, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-Pipe", &simCtx->Pipe, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-Turbulent_Channel_z", &simCtx->channelz, NULL); CHKERRQ(ierr);
    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 = 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);
 
    // 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_max_it", &simCtx->mg_MAX_IT, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-mg_idx", &simCtx->mg_idx, 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, "-poisson_tol", &simCtx->poisson_tol, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-str", &simCtx->STRONG_COUPLING, 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 = PetscOptionsGetInt(NULL, NULL, "-no_of_bodies", &simCtx->NumberOfBodies, NULL); CHKERRQ(ierr);
    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
    ierr = PetscOptionsGetInt(NULL, NULL, "-inlet", &simCtx->inletprofile, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-Ogrid", &simCtx->Ogrid, NULL); CHKERRQ(ierr);
    // NOTE: channelz was not parsed, likely set programmatically. We will omit its parsing call.
    ierr = PetscOptionsGetInt(NULL, NULL, "-grid1d", &simCtx->grid1d, NULL); CHKERRQ(ierr);
    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);
    ierr = PetscOptionsGetInt(NULL, NULL, "-pbc_domain", &simCtx->blkpbc, NULL); CHKERRQ(ierr);
    // NOTE: pseudo_periodic was not parsed. We will omit its parsing call. 
    ierr = PetscOptionsGetReal(NULL, NULL, "-grid_rotation_angle", &simCtx->grid_rotation_angle, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-Croty", &simCtx->Croty, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-Crotz", &simCtx->Crotz, NULL); CHKERRQ(ierr);
    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);
     ierr = PetscOptionsGetReal(NULL, NULL, "-U_bc", &simCtx->U_bc, NULL); 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/RANS) \n");
    PetscInt temp_les_model;
    ierr = PetscOptionsGetInt(NULL, NULL, "-les", &temp_les_model, NULL); CHKERRQ(ierr);
    simCtx->les = (LESModelType)temp_les_model;
    ierr = PetscOptionsGetInt(NULL, NULL, "-rans", &simCtx->rans, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-wallfunction", &simCtx->wallfunction, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-mixed", &simCtx->mixed, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-clark", &simCtx->clark, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-dynamic_freq", &simCtx->dynamic_freq, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-max_cs", &simCtx->max_cs, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-const_cs", &simCtx->Const_CS, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-testfilter_ik", &simCtx->testfilter_ik, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-testfilter_1d", &simCtx->testfilter_1d, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-i_homo_filter", &simCtx->i_homo_filter, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-j_homo_filter", &simCtx->j_homo_filter, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-k_homo_filter", &simCtx->k_homo_filter, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetBool(NULL, NULL, "-averaging", &simCtx->averaging, NULL); 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);
    ierr = PetscOptionsGetBool(NULL, NULL, "-read_fields", &simCtx->readFields, 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 = InitializeBrownianRNG(simCtx); CHKERRQ(ierr);
    // --- Group 10
    LOG_ALLOW(GLOBAL,LOG_DEBUG, "Parsing Group 10: Immersed Boundary & FSI Data Object Pointers \n");
    ierr = PetscOptionsGetBool(NULL, NULL, "-rs_fsi", &simCtx->rstart_fsi, NULL); CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL, NULL, "-duplicate", &simCtx->duplicate, NULL); CHKERRQ(ierr);
 
    // --- 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);
 
    // === 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:

PetscMkdirRecursive()

static PetscErrorCode PetscMkdirRecursive ( const char *  path )

static

Internal helper implementation: PetscMkdirRecursive().

Local to this translation unit.

Definition at line 974 of file setup.c.

{
    PetscErrorCode ierr;
    char           tmp_path[PETSC_MAX_PATH_LEN];
    char           *p = NULL;
    size_t         len;
    PetscBool      exists;
 
    PetscFunctionBeginUser;
 
    // Create a mutable copy of the path
    len = strlen(path);
    if (len >= sizeof(tmp_path)) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Path is too long to process: %s", path);
    }
    strcpy(tmp_path, path);
 
    // If the path ends with a separator, remove it
    if (tmp_path[len - 1] == '/') {
        tmp_path[len - 1] = 0;
    }
 
    // Iterate through the path, creating each directory level
    for (p = tmp_path + 1; *p; p++) {
        if (*p == '/') {
            *p = 0; // Temporarily terminate the string
 
            // Check if this directory level exists
            ierr = PetscTestDirectory(tmp_path, 'r', &exists); CHKERRQ(ierr);
            if (!exists) {
                ierr = PetscMkdir(tmp_path); CHKERRQ(ierr);
            }
            
            *p = '/'; // Restore the separator
        }
    }
    
    // Create the final, full directory path
    ierr = PetscTestDirectory(tmp_path, 'r', &exists); CHKERRQ(ierr);
    if (!exists) {
        ierr = PetscMkdir(tmp_path); CHKERRQ(ierr);
    }
 
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

SetupSimulationEnvironment()

PetscErrorCode SetupSimulationEnvironment

(

SimCtx *

simCtx

)

Internal helper implementation: SetupSimulationEnvironment().

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

Local to this translation unit.

Definition at line 1026 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.
            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 Directories ---
        char path_buffer[PETSC_MAX_PATH_LEN];
 
        // 1. Check/Create the main output directory
        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);
        }
 
        // 2. Check/Create the Eulerian subdirectory
        ierr = PetscSNPrintf(path_buffer, sizeof(path_buffer), "%s/%s", simCtx->output_dir, simCtx->euler_subdir); CHKERRQ(ierr);
        LOG_ALLOW(GLOBAL, LOG_DEBUG, "Verifying Eulerian subdirectory: %s\n", path_buffer);
        ierr = PetscTestDirectory(path_buffer, 'r', &exists); CHKERRQ(ierr);
        if (!exists) {
            LOG_ALLOW(GLOBAL, LOG_INFO, "Eulerian subdirectory not found. Creating: %s\n", path_buffer);
            ierr = PetscMkdir(path_buffer); CHKERRQ(ierr);
        }
 
        // 3. Check/Create the Particle subdirectory if needed
        if (simCtx->np > 0) {
          ierr = PetscSNPrintf(path_buffer, sizeof(path_buffer), "%s/%s", simCtx->output_dir, simCtx->particle_subdir); CHKERRQ(ierr);
          LOG_ALLOW(GLOBAL, LOG_DEBUG, "Verifying Particle subdirectory: %s\n", path_buffer);
          ierr = PetscTestDirectory(path_buffer, 'r', &exists); CHKERRQ(ierr);
          if (!exists) {
              LOG_ALLOW(GLOBAL, LOG_INFO, "Particle subdirectory not found. Creating: %s\n", path_buffer);
              ierr = PetscMkdir(path_buffer); 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);
 
    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:

AllocateContextHierarchy()

static PetscErrorCode AllocateContextHierarchy ( SimCtx *  simCtx )

static

Internal helper implementation: AllocateContextHierarchy().

Local to this translation unit.

Definition at line 1215 of file setup.c.

{
    PetscErrorCode ierr;
    UserMG         *usermg = &simCtx->usermg;
    MGCtx          *mgctx;
    PetscInt       nblk = simCtx->block_number;
    PetscBool      found;
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Allocating context hierarchy for %d levels and %d blocks...\n", simCtx->mglevels, nblk);
 
    // Store the number of levels in the UserMG struct itself
    usermg->mglevels = simCtx->mglevels;
 
    // --- 1. Allocate the array of MGCtx structs ---
    ierr = PetscMalloc(usermg->mglevels * sizeof(MGCtx), &usermg->mgctx); CHKERRQ(ierr);
    // Zero-initialize to ensure all pointers (especially packer) are NULL
    ierr = PetscMemzero(usermg->mgctx, usermg->mglevels * sizeof(MGCtx)); CHKERRQ(ierr);
    mgctx = usermg->mgctx;
    LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Allocated MGCtx array of size %d.\n", simCtx->rank, usermg->mglevels);
    
    // --- 2. Parse semi-coarsening options (logic from MG_Initial) ---
    // These flags determine if a dimension is coarsened in the multigrid hierarchy.
    PetscInt *isc, *jsc, *ksc;
    ierr = PetscMalloc3(nblk, &isc, nblk, &jsc, nblk, &ksc); CHKERRQ(ierr);
    // Set defaults to FALSE (full coarsening)
    for (PetscInt i = 0; i < nblk; ++i) {
        isc[i] = 0; jsc[i] = 0; ksc[i] = 0;
    }
 
// Use a temporary variable for the 'count' argument to the parsing function.
    // This protects the original 'nblk' which is needed for the loop bounds.
    PetscInt n_opts_found = nblk;
    ierr = PetscOptionsGetIntArray(NULL, NULL, "-mg_i_semi", isc, &n_opts_found, &found); CHKERRQ(ierr);
 
    n_opts_found = nblk; // Reset the temp variable before the next call
    ierr = PetscOptionsGetIntArray(NULL, NULL, "-mg_j_semi", jsc, &n_opts_found, &found); CHKERRQ(ierr);
 
    n_opts_found = nblk; // Reset the temp variable before the next call
    ierr = PetscOptionsGetIntArray(NULL, NULL, "-mg_k_semi", ksc, &n_opts_found, &found); CHKERRQ(ierr);
 
    // --- 3. Loop over levels and blocks to allocate UserCtx arrays ---
    for (PetscInt level = 0; level < simCtx->mglevels; level++) {
 
        LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Setting up MG Level %d...\n", simCtx->rank, level);        
        // Allocate the array of UserCtx structs for this level
        ierr = PetscMalloc(nblk * sizeof(UserCtx), &mgctx[level].user); CHKERRQ(ierr);
        // It's good practice to zero out the memory to avoid uninitialized values
        ierr = PetscMemzero(mgctx[level].user, nblk * sizeof(UserCtx)); CHKERRQ(ierr);
        mgctx[level].thislevel = level;
 
        for (PetscInt bi = 0; bi < nblk; bi++) {
            UserCtx *currentUser = &mgctx[level].user[bi];
            LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d:   Initializing UserCtx for Level %d, Block %d.\n", simCtx->rank, level, bi);
        
            // --- CRITICAL STEP: Set the back-pointer to the master context ---
            currentUser->simCtx = simCtx;
 
            // Initialize other per-context values
            currentUser->thislevel = level;
            currentUser->_this = bi; //
            currentUser->mglevels = usermg->mglevels;
 
            // Assign semi-coarsening flags
            currentUser->isc = isc[bi];
            currentUser->jsc = jsc[bi];
            currentUser->ksc = ksc[bi];
 
            // Link to finer/coarser contexts for multigrid operations
            if (level > 0) {
                currentUser->user_c = &mgctx[level-1].user[bi];
        LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG, "Rank %d:     -> Linked to coarser context (user_c).\n", simCtx->rank);
            }
            if (level < usermg->mglevels - 1) {
                currentUser->user_f = &mgctx[level+1].user[bi];
        LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG, "Rank %d:     -> Linked to finer context (user_f).\n", simCtx->rank);
            }
        }
    }
 
    // Log a summary of the parsed flags on each rank.
    if (get_log_level() >= LOG_DEBUG && nblk > 0) {
        LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Final semi-coarsening configuration view:\n", simCtx->rank);
        for (PetscInt bi = 0; bi < nblk; ++bi) {
            LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d:   Block %d: i-semi=%d, j-semi=%d, k-semi=%d\n", simCtx->rank, bi, isc[bi], jsc[bi], ksc[bi]);
        }
    }
 
    // Clean up temporary arrays
    ierr = PetscFree3(isc, jsc, ksc); CHKERRQ(ierr);
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Context hierarchy allocation complete.\n");
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

SetupSolverParameters()

static PetscErrorCode SetupSolverParameters ( SimCtx *  simCtx )

static

Internal helper implementation: SetupSolverParameters().

Local to this translation unit.

Definition at line 1318 of file setup.c.

                                                           {
  
  PetscFunctionBeginUser;
  PROFILE_FUNCTION_BEGIN;
 
  LOG_ALLOW(GLOBAL,LOG_INFO, " -- Setting up solver parameters -- .\n");
 
  UserMG         *usermg = &simCtx->usermg;
  MGCtx          *mgctx = usermg->mgctx;
  PetscInt       nblk = simCtx->block_number;
 
  for (PetscInt level  = usermg->mglevels-1; level >=0; level--) {
    for (PetscInt bi = 0; bi < nblk; bi++) {
      UserCtx *user = &mgctx[level].user[bi];
      LOG_ALLOW_SYNC(LOCAL, LOG_DEBUG, "Rank %d: Setting up parameters for level %d, block %d\n", simCtx->rank, level, bi);
 
      user->assignedA = PETSC_FALSE;
      user->multinullspace = PETSC_FALSE;
    }
  }
  PROFILE_FUNCTION_END;
  PetscFunctionReturn(0);
}

Here is the caller graph for this function:

SetupGridAndSolvers()

PetscErrorCode SetupGridAndSolvers

(

SimCtx *

simCtx

)

Implementation of SetupGridAndSolvers().

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

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

See also

SetupGridAndSolvers()

Definition at line 1350 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);
    ierr = SetupSolverParameters(simCtx); CHKERRQ(ierr);
    ierr = InitializeSolutionConvergenceState(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

)

Internal helper implementation: CreateAndInitializeAllVectors().

Creates and initializes all PETSc Vec objects for all fields.

Local to this translation unit.

Definition at line 1387 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 || simCtx->rans) {
                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 LES/RANS model.\n");
                if(simCtx->les){ 
                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, "Smagorinsky constant (CS) vectors created for LES model.\n");
                }
 
                if(simCtx->wallfunction){
                  ierr = DMCreateLocalVector(user->fda,&user->lFriction_Velocity); CHKERRQ(ierr); ierr = VecSet(user->lFriction_Velocity,0.0);
                }
                  // 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)...
        
            }
        // --- 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);
        
      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);
 
                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->Psi_nodal = NULL;
        }
      }
    
  }
}
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "All simulation vectors created and initialized.\n");
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

RepairPeriodicNormalFaceGhosts()

static PetscErrorCode RepairPeriodicNormalFaceGhosts ( UserCtx *  user, DM  dm, Vec  local_vec, PetscInt  dof, char  face_direction, PetscBool  component_staggered  )

static

Repairs the adjacent normal ghost layer for periodic face-staggered data.

PETSc wraps every component with cell-style indexing. A face family instead needs its adjacent normal ghost shifted by one additional physical face. With the width-three periodic DMDA, the required value is available in the deeper PETSc ghost at -3 or n+2. Tangential ghosts retain PETSc's native wraparound.

Definition at line 1586 of file setup.c.

{
    DMDALocalInfo info;
    PetscInt xs, xe, ys, ye, zs, ze;
    PetscInt gxs, gxe, gys, gye, gzs, gze;
    PetscInt mx, my, mz;
 
    PetscFunctionBeginUser;
    if (!face_direction && !component_staggered) PetscFunctionReturn(0);
 
    PetscCall(DMDAGetLocalInfo(dm, &info));
    xs = info.xs; xe = info.xs + info.xm;
    ys = info.ys; ye = info.ys + info.ym;
    zs = info.zs; ze = info.zs + info.zm;
    gxs = info.gxs; gxe = info.gxs + info.gxm;
    gys = info.gys; gye = info.gys + info.gym;
    gzs = info.gzs; gze = info.gzs + info.gzm;
    mx = info.mx; my = info.my; mz = info.mz;
 
    if (component_staggered) {
        Cmpnts ***array;
        PetscCall(DMDAVecGetArray(dm, local_vec, &array));
 
        if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && xs == 0) {
            PetscCheck(gxs <= -3, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
                       "Periodic Ucont.x ghost repair requires DMDA stencil width at least 3.");
            for (PetscInt k = gzs; k < gze; k++) for (PetscInt j = gys; j < gye; j++)
                array[k][j][-1].x = array[k][j][-3].x;
        }
        if (user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC && xe == mx) {
            PetscCheck(gxe > mx + 2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
                       "Periodic Ucont.x ghost repair requires DMDA stencil width at least 3.");
            for (PetscInt k = gzs; k < gze; k++) for (PetscInt j = gys; j < gye; j++)
                array[k][j][mx].x = array[k][j][mx + 2].x;
        }
        if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && ys == 0) {
            PetscCheck(gys <= -3, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
                       "Periodic Ucont.y ghost repair requires DMDA stencil width at least 3.");
            for (PetscInt k = gzs; k < gze; k++) for (PetscInt i = gxs; i < gxe; i++)
                array[k][-1][i].y = array[k][-3][i].y;
        }
        if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC && ye == my) {
            PetscCheck(gye > my + 2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
                       "Periodic Ucont.y ghost repair requires DMDA stencil width at least 3.");
            for (PetscInt k = gzs; k < gze; k++) for (PetscInt i = gxs; i < gxe; i++)
                array[k][my][i].y = array[k][my + 2][i].y;
        }
        if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && zs == 0) {
            PetscCheck(gzs <= -3, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
                       "Periodic Ucont.z ghost repair requires DMDA stencil width at least 3.");
            for (PetscInt j = gys; j < gye; j++) for (PetscInt i = gxs; i < gxe; i++)
                array[-1][j][i].z = array[-3][j][i].z;
        }
        if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC && ze == mz) {
            PetscCheck(gze > mz + 2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
                       "Periodic Ucont.z ghost repair requires DMDA stencil width at least 3.");
            for (PetscInt j = gys; j < gye; j++) for (PetscInt i = gxs; i < gxe; i++)
                array[mz][j][i].z = array[mz + 2][j][i].z;
        }
 
        PetscCall(DMDAVecRestoreArray(dm, local_vec, &array));
        PetscFunctionReturn(0);
    }
 
    if (dof == 1) {
        PetscReal ***array;
        PetscCall(DMDAVecGetArray(dm, local_vec, &array));
 
        if (face_direction == 'i') {
            if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && xs == 0) {
                PetscCheck(gxs <= -3, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
                           "Periodic I-face ghost repair requires DMDA stencil width at least 3.");
                for (PetscInt k = gzs; k < gze; k++) for (PetscInt j = gys; j < gye; j++)
                    array[k][j][-1] = array[k][j][-3];
            }
            if (user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC && xe == mx) {
                PetscCheck(gxe > mx + 2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
                           "Periodic I-face ghost repair requires DMDA stencil width at least 3.");
                for (PetscInt k = gzs; k < gze; k++) for (PetscInt j = gys; j < gye; j++)
                    array[k][j][mx] = array[k][j][mx + 2];
            }
        } else if (face_direction == 'j') {
            if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && ys == 0) {
                PetscCheck(gys <= -3, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
                           "Periodic J-face ghost repair requires DMDA stencil width at least 3.");
                for (PetscInt k = gzs; k < gze; k++) for (PetscInt i = gxs; i < gxe; i++)
                    array[k][-1][i] = array[k][-3][i];
            }
            if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC && ye == my) {
                PetscCheck(gye > my + 2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
                           "Periodic J-face ghost repair requires DMDA stencil width at least 3.");
                for (PetscInt k = gzs; k < gze; k++) for (PetscInt i = gxs; i < gxe; i++)
                    array[k][my][i] = array[k][my + 2][i];
            }
        } else if (face_direction == 'k') {
            if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && zs == 0) {
                PetscCheck(gzs <= -3, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
                           "Periodic K-face ghost repair requires DMDA stencil width at least 3.");
                for (PetscInt j = gys; j < gye; j++) for (PetscInt i = gxs; i < gxe; i++)
                    array[-1][j][i] = array[-3][j][i];
            }
            if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC && ze == mz) {
                PetscCheck(gze > mz + 2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
                           "Periodic K-face ghost repair requires DMDA stencil width at least 3.");
                for (PetscInt j = gys; j < gye; j++) for (PetscInt i = gxs; i < gxe; i++)
                    array[mz][j][i] = array[mz + 2][j][i];
            }
        }
 
        PetscCall(DMDAVecRestoreArray(dm, local_vec, &array));
    } else {
        Cmpnts ***array;
        PetscCall(DMDAVecGetArray(dm, local_vec, &array));
 
        if (face_direction == 'i') {
            if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && xs == 0) {
                PetscCheck(gxs <= -3, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
                           "Periodic I-face ghost repair requires DMDA stencil width at least 3.");
                for (PetscInt k = gzs; k < gze; k++) for (PetscInt j = gys; j < gye; j++)
                    array[k][j][-1] = array[k][j][-3];
            }
            if (user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC && xe == mx) {
                PetscCheck(gxe > mx + 2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
                           "Periodic I-face ghost repair requires DMDA stencil width at least 3.");
                for (PetscInt k = gzs; k < gze; k++) for (PetscInt j = gys; j < gye; j++)
                    array[k][j][mx] = array[k][j][mx + 2];
            }
        } else if (face_direction == 'j') {
            if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && ys == 0) {
                PetscCheck(gys <= -3, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
                           "Periodic J-face ghost repair requires DMDA stencil width at least 3.");
                for (PetscInt k = gzs; k < gze; k++) for (PetscInt i = gxs; i < gxe; i++)
                    array[k][-1][i] = array[k][-3][i];
            }
            if (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC && ye == my) {
                PetscCheck(gye > my + 2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
                           "Periodic J-face ghost repair requires DMDA stencil width at least 3.");
                for (PetscInt k = gzs; k < gze; k++) for (PetscInt i = gxs; i < gxe; i++)
                    array[k][my][i] = array[k][my + 2][i];
            }
        } else if (face_direction == 'k') {
            if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && zs == 0) {
                PetscCheck(gzs <= -3, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
                           "Periodic K-face ghost repair requires DMDA stencil width at least 3.");
                for (PetscInt j = gys; j < gye; j++) for (PetscInt i = gxs; i < gxe; i++)
                    array[-1][j][i] = array[-3][j][i];
            }
            if (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC && ze == mz) {
                PetscCheck(gze > mz + 2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
                           "Periodic K-face ghost repair requires DMDA stencil width at least 3.");
                for (PetscInt j = gys; j < gye; j++) for (PetscInt i = gxs; i < gxe; i++)
                    array[mz][j][i] = array[mz + 2][j][i];
            }
        }
 
        PetscCall(DMDAVecRestoreArray(dm, local_vec, &array));
    }
 
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

UpdateLocalGhosts()

PetscErrorCode UpdateLocalGhosts	(	UserCtx *	user,
		const char *	fieldName
	)

Internal helper implementation: UpdateLocalGhosts().

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

Local to this translation unit.

Definition at line 1755 of file setup.c.

{
    PetscErrorCode ierr;
    PetscMPIInt    rank;
    Vec            globalVec = NULL;
    Vec            localVec = NULL;
    DM             dm = NULL; // The DM associated with this field pair
    PetscInt       dof = 0;
    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);
    LOG_ALLOW(GLOBAL, LOG_INFO, "Rank %d: Starting ghost update for field '%s'.\n", rank, fieldName);
 
    // --- 1. Identify the correct Vectors and DM ---
    if (strcmp(fieldName, "Coordinates") == 0) {
        ierr = DMGetCoordinates(user->da, &globalVec); CHKERRQ(ierr);
        ierr = DMGetCoordinatesLocal(user->da, &localVec); CHKERRQ(ierr);
        dm = user->fda;
        dof = 3;
    } else if (strcmp(fieldName, "Ucat") == 0) {
        globalVec = user->Ucat;
        localVec  = user->lUcat;
        dm        = user->fda;
        dof       = 3;
    } else if (strcmp(fieldName, "Ucont") == 0) {
        globalVec = user->Ucont;
        localVec  = user->lUcont;
        dm        = user->fda;
        dof       = 3;
        component_staggered = PETSC_TRUE;
    } else if (strcmp(fieldName, "Ucont_o") == 0) {
        globalVec = user->Ucont_o;
        localVec  = user->lUcont_o;
        dm        = user->fda;
        dof       = 3;
        component_staggered = PETSC_TRUE;
    } else if (strcmp(fieldName, "Ucont_rm1") == 0) {
        globalVec = user->Ucont_rm1;
        localVec  = user->lUcont_rm1;
        dm        = user->fda;
        dof       = 3;
        component_staggered = PETSC_TRUE;
    } else if (strcmp(fieldName, "P") == 0) {
        globalVec = user->P;
        localVec  = user->lP;
        dm        = user->da;
    } else if (strcmp(fieldName, "Nu_t") == 0 || strcmp(fieldName, "Eddy Viscosity") == 0) {
        globalVec = user->Nu_t;
        localVec  = user->lNu_t;
        dm        = user->da;
    } else if (strcmp(fieldName, "CS") == 0 || strcmp(fieldName, "Cs") == 0) {
        globalVec = user->CS;
        localVec  = user->lCs;
        dm        = user->da;
    } else if (strcmp(fieldName, "Diffusivity") == 0) {
        globalVec = user->Diffusivity;
        localVec  = user->lDiffusivity;
        dm        = user->da;
    } else if (strcmp(fieldName, "DiffusivityGradient") == 0) {
        globalVec = user->DiffusivityGradient;
        localVec  = user->lDiffusivityGradient;
        dm        = user->fda;
    } else if (strcmp(fieldName, "Csi") == 0) {
        globalVec = user->Csi;
        localVec  = user->lCsi;
        dm        = user->fda;
        dof       = 3;
        face_direction = 'i';
    } else if (strcmp(fieldName, "Eta") == 0) {
        globalVec = user->Eta;
        localVec  = user->lEta;
        dm        = user->fda;
        dof       = 3;
        face_direction = 'j';
    }  else if (strcmp(fieldName, "Zet") == 0) {
        globalVec = user->Zet;
        localVec  = user->lZet;
        dm        = user->fda;
        dof       = 3;
        face_direction = 'k';
    }else if (strcmp(fieldName, "Nvert") == 0) {
        globalVec = user->Nvert;
        localVec  = user->lNvert;
        dm        = user->da;
     // Add other fields as needed
    } else if (strcmp(fieldName, "Aj") == 0) {
        globalVec = user->Aj;
        localVec  = user->lAj;
        dm        = user->da;
    } else if (strcmp(fieldName, "Cent") == 0) {
        globalVec = user->Cent;
        localVec  = user->lCent;
        dm        = user->fda;
    }else if (strcmp(fieldName, "GridSpace") == 0) {
        globalVec = user->GridSpace;
        localVec  = user->lGridSpace;
        dm        = user->fda;
    }else if (strcmp(fieldName, "Centx") == 0) {
        globalVec = user->Centx;
        localVec  = user->lCentx;
        dm        = user->fda;
        dof       = 3;
        face_direction = 'i';
    }else if (strcmp(fieldName, "Centy") == 0) {
        globalVec = user->Centy;
        localVec  = user->lCenty;
        dm        = user->fda;
        dof       = 3;
        face_direction = 'j';
    }else if (strcmp(fieldName, "Centz") == 0) {
        globalVec = user->Centz;
        localVec  = user->lCentz;
        dm        = user->fda;
        dof       = 3;
        face_direction = 'k';
    }else if (strcmp(fieldName,"ICsi") == 0){
      globalVec = user->ICsi;
      localVec  = user->lICsi;
      dm        = user->fda;
      dof       = 3;
      face_direction = 'i';
    }else if (strcmp(fieldName,"IEta") == 0){
      globalVec = user->IEta;
      localVec  = user->lIEta;
      dm        = user->fda;
      dof       = 3;
      face_direction = 'i';
    }else if (strcmp(fieldName,"IZet") == 0){
      globalVec = user->IZet;
      localVec  = user->lIZet;
      dm        = user->fda;
      dof       = 3;
      face_direction = 'i';
    }else if (strcmp(fieldName,"JCsi") == 0){
      globalVec = user->JCsi;
      localVec  = user->lJCsi;
      dm        = user->fda;
      dof       = 3;
      face_direction = 'j';
    }else if (strcmp(fieldName,"JEta") == 0){
      globalVec = user->JEta;
      localVec  = user->lJEta;
      dm        = user->fda;
      dof       = 3;
      face_direction = 'j';
    }else if (strcmp(fieldName,"JZet") == 0){
      globalVec = user->JZet;
      localVec  = user->lJZet;
      dm        = user->fda;
      dof       = 3;
      face_direction = 'j';
    }else if (strcmp(fieldName,"KCsi") == 0){
      globalVec = user->KCsi;
      localVec  = user->lKCsi;
      dm        = user->fda;
      dof       = 3;
      face_direction = 'k';
    }else if (strcmp(fieldName,"KEta") == 0){
      globalVec = user->KEta;
      localVec  = user->lKEta;
      dm        = user->fda;
      dof       = 3;
      face_direction = 'k';
    }else if (strcmp(fieldName,"KZet") == 0){
      globalVec = user->KZet;
      localVec  = user->lKZet;
      dm        = user->fda;
      dof       = 3;
      face_direction = 'k';
    }else if (strcmp(fieldName,"IAj") == 0){
      globalVec = user->IAj;
      localVec  = user->lIAj;
      dm        = user->da;
      dof       = 1;
      face_direction = 'i';
    }else if (strcmp(fieldName,"JAj") == 0){
      globalVec = user->JAj;
      localVec  = user->lJAj;
      dm        = user->da;
      dof       = 1;
      face_direction = 'j';
    }else if (strcmp(fieldName,"KAj") == 0){
      globalVec = user->KAj;
      localVec  = user->lKAj;
      dm        = user->da;
      dof       = 1;
      face_direction = 'k';
    }else if (strcmp(fieldName,"Phi") == 0){ // Pressure correction term.
      globalVec = user->Phi;
      localVec  = user->lPhi;
      dm        = user->da;
    }else if (strcmp(fieldName,"Psi") == 0){ // Particle scalar property.
      globalVec = user->Psi;
      localVec  = user->lPsi;
      dm        = user->da;
    }else if (strcmp(fieldName,"Nvert_o") == 0){
      globalVec = user->Nvert_o;
      localVec  = user->lNvert_o;
      dm        = user->da;
    }else if (strcmp(fieldName,"ParticleCount") == 0){
      globalVec = user->ParticleCount;
      localVec  = user->lParticleCount;
      dm        = user->da;
    }else if (strcmp(fieldName,"K_Omega") == 0){
      globalVec = user->K_Omega;
      localVec  = user->lK_Omega;
      dm        = user->fda2;
    }else if (strcmp(fieldName,"K_Omega_o") == 0){
      globalVec = user->K_Omega_o;
      localVec  = user->lK_Omega_o;
      dm        = user->fda2;
    }else {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Field '%s' not recognized for ghost update.", fieldName);
    }
 
    // --- 2. Check if components were found ---
    if (!globalVec) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Global vector for field '%s' is NULL.", fieldName);
    }
    if (!localVec) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Local vector for field '%s' is NULL.", fieldName);
    }
    if (!dm) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "DM for field '%s' is NULL.", fieldName);
    }
 
    LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Identified components for '%s': DM=%p, GlobalVec=%p, LocalVec=%p.\n",
              rank, fieldName, (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", fieldName, 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", fieldName, 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, fieldName);
    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, fieldName);
 
    // --- 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", fieldName, 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 (strcmp(fieldName, "Ucat") == 0) { // 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, fieldName);
           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, fieldName);
           } else if (!lUcat_arr_test) {
                LOG_ALLOW(LOCAL, LOG_ERROR, "Rank %d: ERROR NULL pointer getting '%s' array after ghost update!\n", rank_test, fieldName);
           }
           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, fieldName); }
               LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Finished testing '%s' access.\n", rank_test, fieldName);
           }
        } // end if Ucat
    } // end debug logging check
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Rank %d: Completed ghost update for field '%s'.\n", rank, fieldName);
    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

)

Internal helper implementation: SetupBoundaryConditions().

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

Local to this translation unit.

Definition at line 2124 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:

Allocate3DArrayScalar()

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

Internal helper implementation: Allocate3DArrayScalar().

Allocates a 3D array of PetscReal values using PetscCalloc.

Local to this translation unit.

Definition at line 2188 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
	)

Internal helper implementation: Deallocate3DArrayScalar().

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

Local to this translation unit.

Definition at line 2223 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
	)

Implementation of Allocate3DArrayVector().

Allocates a contiguous 3D array of Cmpnts values.

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

See also

Allocate3DArrayVector()

Definition at line 2266 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
	)

Implementation of Deallocate3DArrayVector().

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

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

See also

Deallocate3DArrayVector()

Definition at line 2316 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
	)

Internal helper implementation: GetOwnedCellRange().

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

Local to this translation unit.

Definition at line 2382 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

)

Internal helper implementation: ComputeAndStoreNeighborRanks().

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

Local to this translation unit.

Definition at line 2479 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
	)

Internal helper implementation: SetDMDAProcLayout().

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

Local to this translation unit.

Definition at line 2595 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

)

Implementation of SetupDomainRankInfo().

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

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

See also

SetupDomainRankInfo()

Definition at line 2673 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

)

Internal helper implementation: Contra2Cart().

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

Local to this translation unit.

Definition at line 2746 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 a spatially varying Cartesian velocity field to contravariant fluxes.

Convert the ghosted Cartesian velocity field to contravariant face fluxes.

Definition at line 2881 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
	)

Convert a uniform Cartesian velocity (u,v,w) to contravariant fluxes in Ucont.

Populate contravariant fluxes from one uniform Cartesian velocity.

Computes the dot product of the physical velocity with each face-area vector: U^xi = csi · (u,v,w), U^eta = eta · (u,v,w), U^zeta = zet · (u,v,w). Writes to all owned nodes (xs..xe, ys..ye, zs..ze); boundary ghosts are overwritten later by ApplyBoundaryConditions.

Parameters

user	UserCtx with fda, Ucont, lCsi, lEta, lZet populated.
u	Physical Cartesian x-velocity.
v	Physical Cartesian y-velocity.
w	Physical Cartesian z-velocity.

Definition at line 2953 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

)

Internal helper implementation: SetupDomainCellDecompositionMap().

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

Local to this translation unit.

Definition at line 3000 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
	)

Implementation of BinarySearchInt64().

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

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

See also

BinarySearchInt64()

Definition at line 3065 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:

Gidx()

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

static

Internal helper implementation: Gidx().

Local to this translation unit.

Definition at line 3109 of file setup.c.

{
  PetscInt nidx;
  DMDALocalInfo info = user->info;
 
  PetscInt  mx = info.mx, my = info.my;
  
  AO ao;
  DMDAGetAO(user->da, &ao);
  nidx=i+j*mx+k*mx*my;
  
  AOApplicationToPetsc(ao,1,&nidx);
  
  return (nidx);
}

Here is the caller graph for this function:

ComputeDivergence()

PetscErrorCode ComputeDivergence

(

UserCtx *

user

)

Implementation of ComputeDivergence().

Computes the discrete divergence of the contravariant velocity field.

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

See also

ComputeDivergence()

Definition at line 3135 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
	)

Implementation of InitializeRandomGenerators().

Initializes random number generators for assigning particle properties.

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

See also

InitializeRandomGenerators()

Definition at line 3301 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);
    ierr = PetscRandomSetSeed(*randx, rank + 12345); CHKERRQ(ierr);  // Unique seed per rank
    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, rank + 67890); CHKERRQ(ierr);  // Unique seed per rank
    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, rank + 54321); CHKERRQ(ierr);  // Unique seed per rank
    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	(	PetscRandom *	rand_logic_i,
		PetscRandom *	rand_logic_j,
		PetscRandom *	rand_logic_k
	)

Internal helper implementation: InitializeLogicalSpaceRNGs().

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

Local to this translation unit.

Definition at line 3342 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)
    ierr = PetscRandomSetSeed(*rand_logic_i, rank + 202401); CHKERRQ(ierr); // Unique seed
    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, rank + 202402); 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, rank + 202403); 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

)

Internal helper implementation: InitializeBrownianRNG().

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

Local to this translation unit.

Definition at line 3386 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 based on Rank to ensure spatial randomness
    // Multiplying by a large prime helps separate the streams significantly
    unsigned long seed = (unsigned long)rank * 987654321 + (unsigned long)time(NULL);
    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
	)

Implementation of TransformScalarDerivativesToPhysical().

Transforms scalar derivatives from computational space to physical space.

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

See also

TransformScalarDerivativesToPhysical()

Definition at line 3425 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:

TransformDerivativesToPhysical()

static void TransformDerivativesToPhysical ( PetscReal  jacobian, Cmpnts  csi_metrics, Cmpnts  eta_metrics, Cmpnts  zet_metrics, Cmpnts  deriv_csi, Cmpnts  deriv_eta, Cmpnts  deriv_zet, Cmpnts *  dudx, Cmpnts *  dvdx, Cmpnts *  dwdx  )

static

Internal helper implementation: TransformDerivativesToPhysical().

Local to this translation unit.

Definition at line 3450 of file setup.c.

{
    // Derivatives of the first component (u)
    dudx->x = jacobian * (deriv_csi.x * csi_metrics.x + deriv_eta.x * eta_metrics.x + deriv_zet.x * zet_metrics.x);
    dudx->y = jacobian * (deriv_csi.x * csi_metrics.y + deriv_eta.x * eta_metrics.y + deriv_zet.x * zet_metrics.y);
    dudx->z = jacobian * (deriv_csi.x * csi_metrics.z + deriv_eta.x * eta_metrics.z + deriv_zet.x * zet_metrics.z);
    // Derivatives of the second component (v)
    dvdx->x = jacobian * (deriv_csi.y * csi_metrics.x + deriv_eta.y * eta_metrics.x + deriv_zet.y * zet_metrics.x);
    dvdx->y = jacobian * (deriv_csi.y * csi_metrics.y + deriv_eta.y * eta_metrics.y + deriv_zet.y * zet_metrics.y);
    dvdx->z = jacobian * (deriv_csi.y * csi_metrics.z + deriv_eta.y * eta_metrics.z + deriv_zet.y * zet_metrics.z);
    // Derivatives of the third component (w)
    dwdx->x = jacobian * (deriv_csi.z * csi_metrics.x + deriv_eta.z * eta_metrics.x + deriv_zet.z * zet_metrics.x);
    dwdx->y = jacobian * (deriv_csi.z * csi_metrics.y + deriv_eta.z * eta_metrics.y + deriv_zet.z * zet_metrics.y);
    dwdx->z = jacobian * (deriv_csi.z * csi_metrics.z + deriv_eta.z * eta_metrics.z + deriv_zet.z * 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
	)

Internal helper implementation: ComputeScalarFieldDerivatives().

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

Local to this translation unit.

Definition at line 3474 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:

ComputeVectorFieldDerivatives()

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

Internal helper implementation: ComputeVectorFieldDerivatives().

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

Local to this translation unit.

Definition at line 3517 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:

DestroyUserVectors()

PetscErrorCode DestroyUserVectors

(

UserCtx *

user

)

Internal helper implementation: DestroyUserVectors().

Destroys all PETSc Vec objects within a single UserCtx structure.

Local to this translation unit.

Definition at line 3571 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/rans) ---
    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->lFriction_Velocity) { ierr = VecDestroy(&user->lFriction_Velocity); CHKERRQ(ierr); }
    if (user->K_Omega) { ierr = VecDestroy(&user->K_Omega); CHKERRQ(ierr); }
    if (user->lK_Omega) { ierr = VecDestroy(&user->lK_Omega); CHKERRQ(ierr); }
    if (user->K_Omega_o) { ierr = VecDestroy(&user->K_Omega_o); CHKERRQ(ierr); }
    if (user->lK_Omega_o) { ierr = VecDestroy(&user->lK_Omega_o); 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->Psi_nodal) { ierr = VecDestroy(&user->Psi_nodal); CHKERRQ(ierr); }
 
    // --- Group K: Interpolation Vectors (Lazy allocation) ---
    if (user->CellFieldAtCorner) { ierr = VecDestroy(&user->CellFieldAtCorner); CHKERRQ(ierr); }
    if (user->lCellFieldAtCorner) { ierr = VecDestroy(&user->lCellFieldAtCorner); CHKERRQ(ierr); }
 
    // --- Group L: Statistical Averaging Vectors (If allocated) ---
    if (user->Ucat_sum) { ierr = VecDestroy(&user->Ucat_sum); CHKERRQ(ierr); }
    if (user->Ucat_cross_sum) { ierr = VecDestroy(&user->Ucat_cross_sum); CHKERRQ(ierr); }
    if (user->Ucat_square_sum) { ierr = VecDestroy(&user->Ucat_square_sum); CHKERRQ(ierr); }
    if (user->P_sum) { ierr = VecDestroy(&user->P_sum); CHKERRQ(ierr); }
 
    // --- Group M: 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 N: Poisson Solver Vectors (Destroyed after solve, but check anyway) ---
    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 caller graph for this function:

DestroyUserContext()

PetscErrorCode DestroyUserContext

(

UserCtx *

user

)

Internal helper implementation: DestroyUserContext().

Destroys all resources allocated within a single UserCtx structure.

Local to this translation unit.

Definition at line 3709 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->C) {
        ierr = MatDestroy(&user->C); CHKERRQ(ierr);
        LOG_ALLOW(LOCAL, LOG_DEBUG, "  Matrix C 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, fda2, 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->fda2) {
        ierr = DMDestroy(&user->fda2); CHKERRQ(ierr);
        LOG_ALLOW(LOCAL, LOG_DEBUG, "  fda2 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");
    }
    if (user->KSKE) {
        ierr = PetscFree(user->KSKE); CHKERRQ(ierr);
        user->KSKE = NULL;
        LOG_ALLOW(LOCAL, LOG_DEBUG, "  KSKE array 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

)

Implementation of FinalizeSimulation().

Main cleanup function for the entire simulation context.

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

See also

FinalizeSimulation()

Definition at line 3808 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);
 
    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: