momentumsolvers.h File Reference

#include "variables.h"
#include "logging.h"
#include "rhs.h"
#include "Boundaries.h"

Include dependency graph for momentumsolvers.h:

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

Go to the source code of this file.

Data Structures
struct	MomStabilityReport
	Diagnostic report produced by ComputeMomentumStabilityEstimate(). More...

Macros
#define	MomentumSolver_DualTime_Picard_RK4   MomentumSolver_DualTime_Picard_JamesonRK

Enumerations
enum	MomStabCandidate { MOM_STAB_CAND_B = 0 , MOM_STAB_CAND_C = 1 , MOM_STAB_CAND_D = 2 }
	Convective-estimate candidate selector. More...
enum	MomStabLimiter { MOM_STAB_LIMITER_TIME = 0 , MOM_STAB_LIMITER_CONVECTION = 1 , MOM_STAB_LIMITER_VISCOSITY = 2 }
	Dominant stiffness contributor at the controlling cell. More...

Functions
PetscErrorCode	MomentumSolver_Explicit_RungeKutta4 (UserCtx *user, IBMNodes *ibm, FSInfo *fsi)
	Advances the momentum equations using an explicit 4th-order Runge-Kutta scheme.
PetscErrorCode	MomentumSolver_NewtonKrylov (UserCtx *user, IBMNodes *ibm, FSInfo *fsi)
	Solves one physical momentum step with matrix-free Newton–Krylov.
PetscErrorCode	MomentumSolver_DualTime_Picard_JamesonRK (UserCtx *user, IBMNodes *ibm, FSInfo *fsi)
	Solves the momentum equations using dual-time Picard iteration with Jameson RK smoothing.
PetscBool	MomentumUsesBDF2 (SimCtx *simCtx)
	Returns whether the current physical step uses the BDF2 discretization.
PetscReal	MomentumBDFCoefficient (SimCtx *simCtx)
	Returns the BDF physical-time coefficient a0 for the current step.
PetscErrorCode	ComputeTotalResidual (UserCtx *user)
	Computes the shared spatial-plus-BDF momentum residual in user->Rhs.
PetscErrorCode	ComputeMomentumStabilityEstimate (UserCtx *user, PetscInt block_number, PetscReal dt, MomStabCandidate candidate, MomStabilityReport *rep)
	Compute the momentum pseudo-time stability estimate (shadow/diagnostic).
PetscInt	MomCellActiveRows (PetscReal ***nvert, PetscInt k, PetscInt j, PetscInt i, PetscInt mx, PetscInt my, PetscInt mz, PetscBool np_x1, PetscBool np_y1, PetscBool np_z1, PetscInt twoD)
	Active staggered-momentum row mask for a cell (exposed for unit testing).

---

Data Structure Documentation

MomStabilityReport

struct MomStabilityReport

Diagnostic report produced by ComputeMomentumStabilityEstimate().

This is a PRACTICAL CONSERVATIVE STABILITY ESTIMATE (operator-scaled pseudo-time estimate), not a proven spectral radius. All lambda_* are global maxima in [1/s].

Definition at line 165 of file momentumsolvers.h.

Data Fields
PetscReal	lambda
PetscReal	lambda_t
PetscReal	lambda_c
PetscReal	lambda_v
PetscReal	lambda_B
PetscReal	lambda_C
PetscReal	lambda_D
PetscInt	ci
PetscInt	cj
PetscInt	ck
PetscInt	cblock
PetscInt	cclass
PetscInt	one_sided
PetscInt	active_cells
PetscBool	estimate_incomplete
MomStabLimiter	limiter

Macro Definition Documentation

MomentumSolver_DualTime_Picard_RK4

#define MomentumSolver_DualTime_Picard_RK4   MomentumSolver_DualTime_Picard_JamesonRK

Deprecated:

Use MomentumSolver_DualTime_Picard_JamesonRK().

Definition at line 94 of file momentumsolvers.h.

Enumeration Type Documentation

MomStabCandidate

enum MomStabCandidate

Convective-estimate candidate selector.

See ComputeMomentumStabilityEstimate().

B: six-face transport scale (f_c * Aj * sum|U_f| / 2). C: B + frozen-advector discrete-divergence diagonal term. D: C + nonlinear velocity-gradient row-norm term (lambda_grad_u).

Enumerator
MOM_STAB_CAND_B
MOM_STAB_CAND_C
MOM_STAB_CAND_D

Definition at line 144 of file momentumsolvers.h.

             {
    MOM_STAB_CAND_B = 0,
    MOM_STAB_CAND_C = 1,
    MOM_STAB_CAND_D = 2
} MomStabCandidate;

MomStabLimiter

enum MomStabLimiter

Dominant stiffness contributor at the controlling cell.

Enumerator
MOM_STAB_LIMITER_TIME
MOM_STAB_LIMITER_CONVECTION
MOM_STAB_LIMITER_VISCOSITY

Definition at line 153 of file momentumsolvers.h.

             {
    MOM_STAB_LIMITER_TIME       = 0,
    MOM_STAB_LIMITER_CONVECTION = 1,
    MOM_STAB_LIMITER_VISCOSITY  = 2
} MomStabLimiter;

Function Documentation

MomentumSolver_Explicit_RungeKutta4()

PetscErrorCode MomentumSolver_Explicit_RungeKutta4 ( UserCtx *  user, IBMNodes *  ibm, FSInfo *  fsi  )

extern

Advances the momentum equations using an explicit 4th-order Runge-Kutta scheme.

Parameters

user	Array of UserCtx structs for all blocks.
ibm	(Optional) Pointer to IBM data. Pass NULL if disabled.
fsi	(Optional) Pointer to FSI data. Pass NULL if disabled.

Returns

PetscErrorCode 0 on success.

Note

Testing status: The explicit RK path remains on the near-term backlog for direct positive-path bespoke coverage; today it is weaker than the dual-time path in the test surface.

Advances the momentum equations using an explicit 4th-order Runge-Kutta scheme.

Local to this translation unit.

Definition at line 93 of file momentumsolvers.c.

{
    PetscErrorCode ierr;
    (void)fsi;
    // --- Context Acquisition ---
    // Get the master simulation context from the first block's UserCtx.
    // This is the bridge to access all former global variables.
    SimCtx *simCtx = user[0].simCtx;
    PetscReal dt = simCtx->dt;
    //PetscReal st = simCtx->st;
    PetscInt  istage;
    PetscReal alfa[] = {0.25, 1.0/3.0, 0.5, 1.0};
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Executing explicit momentum solver (Runge-Kutta) for %d block(s).\n",simCtx->block_number);
 
    // --- 1. Pre-Loop Initialization (Legacy Logic) ---
    // This block prepares boundary conditions and allocates the RHS vector for all blocks
    // before the main RK loop begins. This logic is preserved from the original code.
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Preparing all blocks for Runge-Kutta solve...\n");
    for (PetscInt bi = 0; bi < simCtx->block_number; bi++) {
 
        ierr = ApplyBoundaryConditions(&user[bi]); CHKERRQ(ierr);
 
        /*
        // Immersed boundary interpolation (if enabled)
        if (simCtx->immersed) {
            LOG_ALLOW(LOCAL, LOG_DEBUG, "  Performing pre-RK IBM interpolation for block %d.\n", bi);
            for (PetscInt ibi = 0; ibi < simCtx->NumberOfBodies; ibi++) {
                // The 'ibm' and 'fsi' pointers are passed directly from FlowSolver
                ierr = ibm_interpolation_advanced(&user[bi], &ibm[ibi], ibi, 1); CHKERRQ(ierr);
            }
        }
        */
        
        // Allocate the persistent RHS vector for this block's context
        ierr = VecDuplicate(user[bi].Ucont, &user[bi].Rhs); CHKERRQ(ierr);
    }
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Pre-loop initialization complete.\n");
 
    // --- 2. Main Runge-Kutta Loop ---
    // The legacy code had an outer `pseudot` loop that only ran once. We preserve it.
    for (PetscInt pseudot = 0; pseudot < 1; pseudot++) {
        // Loop over each block to perform the RK stages
        for (PetscInt bi = 0; bi < simCtx->block_number; bi++) {
            for (istage = 0; istage < 4; istage++) {
                LOG_ALLOW(LOCAL, LOG_DEBUG, "  Block %d, RK Stage %d (alpha=%.4f)...\n", bi, istage, alfa[istage]);
 
                // a. Calculate the Right-Hand Side (RHS) of the momentum equation.
                ierr = ComputeRHS(&user[bi], user[bi].Rhs); CHKERRQ(ierr);
 
                // b. Advance Ucont to the next intermediate stage using the RK coefficient.
                //    Ucont_new = Ucont_old + alpha * dt * RHS
                ierr = VecWAXPY(user[bi].Ucont, alfa[istage] * dt, user[bi].Rhs, user[bi].Ucont_o); CHKERRQ(ierr);
 
                // c. Synchronize periodic endpoints and local ghosts for the new intermediate Ucont.
                const char *staggered_fields[] = {"Ucont"};
                ierr = SynchronizePeriodicStaggeredFields(&user[bi], 1, staggered_fields); CHKERRQ(ierr);
 
                // d. Re-apply boundary conditions for the new intermediate velocity.
                //    This is crucial for the stability and accuracy of the multi-stage scheme.
                ierr = ApplyBoundaryConditions(&user[bi]); CHKERRQ(ierr);
                
            } // End of RK stages for one block
 
            /*
            // Final IBM Interpolation for the block (if enabled)
            if (simCtx->immersed) {
                LOG_ALLOW(LOCAL, LOG_DEBUG, "  Performing post-RK IBM interpolation for block %d.\n", bi);
                for (PetscInt ibi = 0; ibi < simCtx->NumberOfBodies; ibi++) {
                    ierr = ibm_interpolation_advanced(&user[bi], &ibm[ibi], ibi, 1); CHKERRQ(ierr);
                }
            }
            */
 
        } // End loop over blocks
 
        // --- 3. Inter-Block Communication (Legacy Logic) ---
        // This is called after all blocks have completed their RK stages.
        if (simCtx->block_number > 1) {
      //    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Updating multi-block interfaces after RK stages.\n");
        //   ierr = Block_Interface_U(user); CHKERRQ(ierr);
        }
 
    } // End of pseudo-time loop
 
    // --- 4. Cleanup ---
    // Destroy the RHS vectors that were created at the start of this function.
    for (PetscInt bi = 0; bi < simCtx->block_number; bi++) {
        ierr = VecDestroy(&user[bi].Rhs); CHKERRQ(ierr);
    }
    
    LOG_ALLOW(GLOBAL, LOG_INFO, "Runge-Kutta solve completed for all blocks.\n");
 
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

MomentumSolver_NewtonKrylov()

PetscErrorCode MomentumSolver_NewtonKrylov	(	UserCtx *	user,
		IBMNodes *	ibm,
		FSInfo *	fsi
	)

Solves one physical momentum step with matrix-free Newton–Krylov.

Version one uses a finite-difference matrix-free Jacobian, GMRES, and no preconditioner. All PETSc solver objects are local to this call. Rows removed by legacy boundary residual enforcement are made explicit: conditioned normal rows use X-Uconditioned, untouched dummy/tangential rows use X, and periodic duplicates use Xdup-Xrepresentative. Unsupported masked, interface, and component-disabled rows are rejected before setup.

Parameters

user	Single-block momentum context.
ibm	Must be NULL; immersed boundaries are not supported in version one.
fsi	Must be NULL; moving-body coupling is not supported in version one.

Returns

0 on convergence, PETSC_ERR_CONV_FAILED after rollback on nonconvergence.

Solves one physical momentum step with matrix-free Newton–Krylov.

Parameters

user	Single-block momentum context.
ibm	Must be NULL in version one.
fsi	Must be NULL in version one.

Returns

0 on convergence or an error after restoring the canonical entry state.

Definition at line 461 of file momentum_newton_krylov.c.

{
    PetscErrorCode ierr = PETSC_SUCCESS, cleanup_ierr;
    SimCtx        *simCtx;
    SNES           snes = NULL;
    Mat            J = NULL;
    Vec            solution = NULL, entry_backup = NULL;
    KSP            ksp = NULL;
    PC             pc = NULL;
    const char    *pc_type = NULL;
    PetscBool      pc_is_none = PETSC_FALSE;
    PetscBool      restore_entry = PETSC_FALSE;
    PetscBool      rhs_created = PETSC_FALSE;
    PetscBool      solve_started = PETSC_FALSE;
    PetscBool      committed = PETSC_FALSE;
    SNESConvergedReason reason = SNES_CONVERGED_ITERATING;
    PetscInt       nonlinear_its = 0, function_evals = 0, linear_its = 0;
    PetscReal      final_norm = PETSC_MAX_REAL;
    MomentumNewtonKrylovContext ctx = {0};
    const char *staggered_fields[] = {"Ucont"};
 
    PetscFunctionBeginUser;
    PetscCall(MomentumNewtonKrylov_Validate(user));
    PetscCheck(ibm == NULL && fsi == NULL, PETSC_COMM_WORLD, PETSC_ERR_SUP,
               "Newton Krylov version one does not accept IBM or FSI objects.");
    PetscCheck(user->Rhs == NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE,
               "Newton Krylov requires UserCtx::Rhs to be unallocated on entry.");
    simCtx = user->simCtx;
 
    ierr = VecDuplicate(user->Ucont, &solution); if (ierr) goto cleanup;
    ierr = VecDuplicate(user->Ucont, &entry_backup); if (ierr) goto cleanup;
    ierr = VecDuplicate(user->Ucont, &user->Rhs); if (ierr) goto cleanup;
    rhs_created = PETSC_TRUE;
    ierr = SNESCreate(PetscObjectComm((PetscObject)user->Ucont), &snes); if (ierr) goto cleanup;
 
    ierr = SynchronizePeriodicStaggeredFields(user, 1, staggered_fields); if (ierr) goto cleanup;
    ierr = ApplyBoundaryConditions(user); if (ierr) goto cleanup;
    ierr = VecCopy(user->Ucont, entry_backup); if (ierr) goto cleanup;
    restore_entry = PETSC_TRUE;
    ierr = VecCopy(user->Ucont, solution); if (ierr) goto cleanup;
 
    ctx.user = user;
    ierr = SNESSetOptionsPrefix(snes, "mom_nk_"); if (ierr) goto cleanup;
    ierr = SNESSetType(snes, SNESNEWTONLS); if (ierr) goto cleanup;
    ierr = SNESSetDM(snes, user->fda); if (ierr) goto cleanup;
    ierr = SNESSetFunction(snes, NULL, MomentumNewtonKrylov_FormResidual, &ctx); if (ierr) goto cleanup;
    ierr = MatCreateSNESMF(snes, &J); if (ierr) goto cleanup;
    ierr = SNESSetJacobian(snes, J, J, MatMFFDComputeJacobian, NULL); if (ierr) goto cleanup;
    ierr = SNESGetKSP(snes, &ksp); if (ierr) goto cleanup;
    ierr = KSPSetType(ksp, KSPGMRES); if (ierr) goto cleanup;
    ierr = KSPGetPC(ksp, &pc); if (ierr) goto cleanup;
    ierr = PCSetType(pc, PCNONE); if (ierr) goto cleanup;
    ierr = SNESSetFromOptions(snes); if (ierr) goto cleanup;
    ierr = SNESMonitorSet(snes, MomentumNewtonKrylov_Monitor, &ctx, NULL); if (ierr) goto cleanup;
    ierr = PCGetType(pc, &pc_type); if (ierr) goto cleanup;
    ierr = PetscStrcmp(pc_type, PCNONE, &pc_is_none); if (ierr) goto cleanup;
    if (!pc_is_none) {
        LOG(GLOBAL, LOG_ERROR,
            "Newton Krylov version one requires PCNONE; option processing selected PC type '%s'.\n",
            pc_type ? pc_type : "(unset)");
        ierr = PETSC_ERR_SUP;
        goto cleanup;
    }
 
    MomentumNewtonKrylov_OpenHistory(&ctx);
    solve_started = PETSC_TRUE;
    ierr = SNESSolve(snes, NULL, solution);
    if (ierr) goto cleanup;
    ierr = SNESGetConvergedReason(snes, &reason); if (ierr) goto cleanup;
    ierr = SNESGetIterationNumber(snes, &nonlinear_its); if (ierr) goto cleanup;
    ierr = SNESGetNumberFunctionEvals(snes, &function_evals); if (ierr) goto cleanup;
    ierr = SNESGetLinearSolveIterations(snes, &linear_its); if (ierr) goto cleanup;
    ierr = SNESGetFunctionNorm(snes, &final_norm); if (ierr) goto cleanup;
 
    if (reason > 0) {
        ierr = VecCopy(solution, user->Ucont); if (ierr) goto cleanup;
        simCtx->mom_last_converged = PETSC_TRUE;
    } else {
        ierr = VecCopy(entry_backup, user->Ucont); if (ierr) goto cleanup;
        simCtx->mom_last_converged = PETSC_FALSE;
    }
    ierr = SynchronizePeriodicStaggeredFields(user, 1, staggered_fields); if (ierr) goto cleanup;
    ierr = ApplyBoundaryConditions(user); if (ierr) goto cleanup;
    restore_entry = PETSC_FALSE;
    committed = (PetscBool)(reason > 0);
 
    LOG_ALLOW(GLOBAL, LOG_INFO,
              "Newton Krylov momentum solve: reason=%s (%d), Newton iterations=%d, residual evaluations=%d, Krylov iterations=%d, final norm=%.6e, state=%s.\n",
              SNESConvergedReasons[reason], (PetscInt)reason, nonlinear_its, function_evals,
              linear_its, (double)final_norm, reason > 0 ? "committed" : "rolled back");
    if (reason <= 0) ierr = PETSC_ERR_CONV_FAILED;
 
cleanup:
    /* A PETSc solve error can bypass the normal statistics path. Query whatever
       SNES retained without replacing the primary error so a failed attempt is
       still represented in the structured log. */
    if (solve_started && snes) {
        (void)SNESGetConvergedReason(snes, &reason);
        (void)SNESGetIterationNumber(snes, &nonlinear_its);
        (void)SNESGetNumberFunctionEvals(snes, &function_evals);
        (void)SNESGetLinearSolveIterations(snes, &linear_its);
        (void)SNESGetFunctionNorm(snes, &final_norm);
    }
    if (restore_entry && entry_backup) {
        cleanup_ierr = VecCopy(entry_backup, user->Ucont);
        if (!ierr) ierr = cleanup_ierr;
        cleanup_ierr = SynchronizePeriodicStaggeredFields(user, 1, staggered_fields);
        if (!ierr) ierr = cleanup_ierr;
        cleanup_ierr = ApplyBoundaryConditions(user);
        if (!ierr) ierr = cleanup_ierr;
        simCtx->mom_last_converged = PETSC_FALSE;
        committed = PETSC_FALSE;
    }
    if (ctx.history_file) {
        (void)fclose(ctx.history_file);
        ctx.history_file = NULL;
    }
    if (solve_started) {
        MomentumNewtonKrylov_WriteSummary(&ctx, reason, nonlinear_its, function_evals,
                                          linear_its, final_norm, committed);
    }
    if (rhs_created) {
        cleanup_ierr = VecDestroy(&user->Rhs);
        if (!ierr) ierr = cleanup_ierr;
    }
    cleanup_ierr = VecDestroy(&entry_backup); if (!ierr) ierr = cleanup_ierr;
    cleanup_ierr = VecDestroy(&solution); if (!ierr) ierr = cleanup_ierr;
    cleanup_ierr = MatDestroy(&J); if (!ierr) ierr = cleanup_ierr;
    cleanup_ierr = SNESDestroy(&snes); if (!ierr) ierr = cleanup_ierr;
    PetscFunctionReturn(ierr);
}

Here is the call graph for this function:

MomentumSolver_DualTime_Picard_JamesonRK()

PetscErrorCode MomentumSolver_DualTime_Picard_JamesonRK	(	UserCtx *	user,
		IBMNodes *	ibm,
		FSInfo *	fsi
	)

Solves the momentum equations using dual-time Picard iteration with Jameson RK smoothing.

---

GLOSSARY & THEORETICAL BASIS

1. METHODOLOGY: Dual-Time Stepping (Pseudo-Time Integration) We aim to solve the implicit BDF equation: R_spatial(U) + dU/dt_physical = 0. We do this by introducing a fictitious "Pseudo-Time" (tau) and iterating to steady state: dU/d(tau) = - [ R_spatial(U) + BDF_Terms(U) ] When dU/d(tau) -> 0, the physical time step is satisfied.
2. ALGORITHM: Fixed-Point Iteration with Explicit Runge-Kutta This is technically a Fixed-Point iteration on the operator: U_new = U_old + pseudo_dtau * alfa_stage * Total_Residual(U_old) where pseudo_dtau = pseudo_cfl / lambda_max is the spectral-radius-based pseudo-time step (lambda_max = global max convective spectral radius). This makes pseudo_cfl a true dimensionless Courant number, independent of the physical time step dt. We use a 4-Stage Explicit RK scheme (Jameson-Schmidt-Turkel coeffs) to smooth errors.
3. STABILITY: Adaptive Pseudo-CFL Trial Acceptance and Rollback If a pseudo-time trial causes excessive residual growth, the solver restores the

previous accepted state, reduces the global pseudo-CFL, and retries.

VARIABLE MAPPING

– Physics Variables (Legacy Names Kept) – ti : Physical Time Step Index. dt : Physical Time Step size (Delta t). st : Pseudo-Time Step size (Delta tau). alfa : Runge-Kutta stage coefficients {1/4, 1/3, 1/2, 1}. – Convergence & Solver Control (Renamed) – pseudo_iter : Counter for the inner dual-time loop. pseudo_dtau : Adaptive pseudo-time step [physical time], = pseudo_cfl / lambda_max. lambda_max : Global max convective spectral radius [1/s] from the current field. delta_sol_norm : The L_inf norm of the change in solution (dU).

resid_norm : The L_inf norm of the Total Residual (RHS).

Parameters

user	Primary UserCtx input for the operation.
ibm	Parameter ibm passed to MomentumSolver_DualTime_Picard_JamesonRK().
fsi	Parameter fsi passed to MomentumSolver_DualTime_Picard_JamesonRK().

Returns

PetscErrorCode 0 on success.

Note

Testing status: This solver is covered primarily through runtime smoke and orchestration tests. A smaller direct invariant-style positive-path harness remains part of the next-gap backlog.

Solves the momentum equations using dual-time Picard iteration with Jameson RK smoothing.

Local to this translation unit.

Definition at line 742 of file momentumsolvers.c.

{
    (void)ibm;
    (void)fsi;
    // --- CONTEXT ACQUISITION BLOCK ---
    SimCtx *simCtx = user[0].simCtx;
    const PetscInt block_number = simCtx->block_number;
    
    // Legacy Names (Physics/Time) - Kept for brevity in formulas
    const PetscInt ti = simCtx->step;
    const PetscReal dt = simCtx->dt;       // Physical Time Step
    /* cfl: initial pseudo-CFL Courant number (dimensionless, flow-independent since Phase 3).
     * Actual dtau = cfl / lambda_max, computed after the spectral radius estimate below. */
    const PetscReal cfl = simCtx->pseudo_cfl;
    const PetscReal alfa[] = {0.25, 1.0/3.0, 0.5, 1.0}; // Jameson RK smoothing coefficients
    Vec *pRhs; // Per-block backup of Rhs at last accepted pseudo-state.
 
    // State flags
    PetscBool force_restart = PETSC_FALSE;
 
    // Renamed Solver Parameters
    const PetscInt  max_pseudo_steps = simCtx->mom_max_pseudo_steps; // Max Dual-Time Iterations
    const PetscReal tol_abs_delta    = simCtx->mom_atol; // Stop if |dU| < tol
    const PetscReal tol_rtol_delta   = simCtx->mom_rtol; // Stop if |dU|/|dU0| < tol
    // --- END CONTEXT ACQUISITION BLOCK ---
 
    PetscErrorCode ierr;
    PetscMPIInt    rank;
    PetscInt       istage, pseudo_iter;
    PetscInt       accepted_iter = 0, rejected_iter = 0, recovery_streak = 0;
    PetscReal      ts, te, cput;
 
    // --- Global Convergence Metrics ---
    PetscReal global_norm_delta = 10.0;
    PetscReal global_rel_delta  = 1.0;
    PetscReal global_norm_resid = 1.0;
    PetscReal global_rel_resid  = 1.0;
    PetscReal smoothed_trial_ratio = 1.0; /* EMA-smoothed ratio; initialized neutral */
 
    // --- Spectral-radius-based pseudo-time step ---
    /* lambda_max: global maximum spectral radius [1/s], computed once per physical timestep.
     * pseudo_dtau[bi]: physical pseudo-time step dtau = pseudo_cfl / lambda_max [time].
     *   Unlike the old scheme (dtau = pseudo_cfl * dt), this is independent of dt and makes
     *   pseudo_cfl a true Courant number. (2.83 is only the scalar imaginary-axis RK limit,
     *   not a generally-stable value for the actual operator -- to be characterized in A4.)
     * dtau_min / dtau_max: per-timestep bounds derived from simCtx->min/max_pseudo_cfl. */
    PetscReal  lambda_max = 0.0;
    PetscReal  dtau_min, dtau_max;
 
    // --- Local Metric Arrays (Per Block) ---
    PetscReal *delta_sol_norm_init, *delta_sol_norm_prev, *delta_sol_norm_curr, *delta_sol_rel_curr;
    PetscReal *resid_norm_init,     *resid_norm_prev,     *resid_norm_curr,     *resid_rel_curr;
    PetscReal *pseudo_dtau;          /* per-block pseudo-time step dtau [physical time, NOT dtau/dt] */
    PetscReal *trial_ratio_log;        /* per-block step-to-step residual ratio (for diagnostics) */
    PetscReal  last_accepted_resid;    /* last accepted |R|, for final summary */
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
    LOG_ALLOW(GLOBAL, LOG_INFO, "Executing Dual-Time Momentum Solver for %d block(s).\n", block_number);
 
    // --- Allocate metric arrays ---
    ierr = PetscMalloc2(block_number, &delta_sol_norm_init, block_number, &delta_sol_norm_prev); CHKERRQ(ierr);
    ierr = PetscMalloc2(block_number, &delta_sol_norm_curr, block_number, &delta_sol_rel_curr); CHKERRQ(ierr);
    ierr = PetscMalloc2(block_number, &resid_norm_init,     block_number, &resid_norm_prev); CHKERRQ(ierr);
    ierr = PetscMalloc2(block_number, &resid_norm_curr,     block_number, &pseudo_dtau); CHKERRQ(ierr);
    ierr = PetscMalloc1(block_number, &resid_rel_curr); CHKERRQ(ierr);
    ierr = PetscMalloc1(block_number, &trial_ratio_log); CHKERRQ(ierr);
    ierr = PetscMalloc1(block_number, &pRhs); CHKERRQ(ierr);
    last_accepted_resid = 0.0;
 
    ierr = PetscTime(&ts); CHKERRQ(ierr);
 
    // --- 1. Pre-Loop Initialization ---
    
    if (block_number > 1) {
      // ierr = Block_Interface_U(user); CHKERRQ(ierr);
    }
 
    for (PetscInt bi = 0; bi < block_number; bi++) {
        ierr = ApplyBoundaryConditions(&user[bi]); CHKERRQ(ierr);
        
        // Immersed boundary interpolation (if enabled)
 
        // Allocate workspace
        ierr = VecDuplicate(user[bi].Ucont, &user[bi].Rhs); CHKERRQ(ierr);
        ierr = VecDuplicate(user[bi].Rhs, &pRhs[bi]); CHKERRQ(ierr);
        ierr = VecDuplicate(user[bi].Ucont, &user[bi].dUcont); CHKERRQ(ierr);
        ierr = VecDuplicate(user[bi].Ucont, &user[bi].pUcont); CHKERRQ(ierr);
        
        // Initialize Backup (pUcont) with current state
        ierr = VecCopy(user[bi].Ucont, user[bi].pUcont); CHKERRQ(ierr);
 
        // --- Calculate Initial Total Residual (Spatial + Temporal) ---
        ierr = ComputeTotalResidual(&user[bi]); CHKERRQ(ierr);
 
        // Backup Rhs with current state
        ierr = VecCopy(user[bi].Rhs, pRhs[bi]); CHKERRQ(ierr);
 
        // Compute Initial Norms
        ierr = VecNorm(user[bi].Rhs, NORM_INFINITY, &resid_norm_init[bi]); CHKERRQ(ierr);
        
        // Initialize history for backtracking logic
        resid_norm_prev[bi]     = resid_norm_init[bi]; /* meaningful ratio on first iteration */
        delta_sol_norm_prev[bi] = 1000.0;
        /* pseudo_dtau[bi] is set after ComputeGlobalSpectralRadiusEstimate below. */
 
        LOG_ALLOW(GLOBAL,LOG_INFO," Block %d | Max RHS = %.6f | initial pseudo-CFL = %.4f .\n", bi, resid_norm_init[bi], cfl);
    }
 
    /* --- Spectral Radius Estimate (computed once per physical timestep) ---
     * lambda_max [1/s] is the global maximum of (face_flux_sum * inverse_cell_volume)
     * over all owned cells across all MPI ranks. See ComputeGlobalSpectralRadiusEstimate().
     *
     * dtau = pseudo_cfl / lambda_max   (pseudo_cfl is now a true Courant number)
     * dtau_min/max are derived from the user-specified min/max_pseudo_cfl bounds.
     *
     * The BDF2 lower bound inside the helper ensures lambda_max >= 1.5/dt, so dtau
     * remains finite even when the velocity field is zero (startup / stagnant ICs). */
    ierr = ComputeGlobalSpectralRadiusEstimate(user, block_number, dt, &lambda_max); CHKERRQ(ierr);
    dtau_min = simCtx->min_pseudo_cfl / lambda_max;
    dtau_max = simCtx->max_pseudo_cfl / lambda_max;
 
    /* Initialise per-block pseudo-time step from the warm-start CFL (carried from last timestep). */
    for (PetscInt bi = 0; bi < block_number; bi++) {
        pseudo_dtau[bi] = cfl / lambda_max;   /* dtau [physical time], NOT a fraction of dt */
    }
 
    /* --- SHADOW MODE: new operator-scaled stability estimate (diagnostic only) ---
     * Computes the Workstream-A estimate alongside the legacy one and logs a compact
     * comparison. The legacy estimate above continues to drive production dtau; the
     * shadow estimate changes nothing. Enable with -mom_stability_shadow. */
    {
        PetscBool shadow = PETSC_FALSE;
        ierr = PetscOptionsGetBool(NULL, NULL, "-mom_stability_shadow", &shadow, NULL); CHKERRQ(ierr);
        if (shadow) {
            MomStabilityReport rep;
            ierr = ComputeMomentumStabilityEstimate(user, block_number, dt, MOM_STAB_CAND_C, &rep); CHKERRQ(ierr);
            const char *lim = (rep.limiter==MOM_STAB_LIMITER_TIME) ? "time"
                            : (rep.limiter==MOM_STAB_LIMITER_CONVECTION) ? "convection" : "viscosity";
            /* Detailed shadow comparison at DEBUG only (gated by level + function allowlist);
               never printed unconditionally. Enable in tests via the logging controls. */
            LOG_ALLOW(GLOBAL, LOG_DEBUG,
                "Momentum scale [shadow]: legacy_dtau=%.4e new_dtau=%.4e ratio=%.4f limiter=%s | "
                "lambda_legacy=%.4e lambda_new=%.4e (B=%.4e C=%.4e D=%.4e) | "
                "lt=%.4e lc=%.4e lv=%.4e | cell=(%d,%d,%d) blk=%d class=%d onesided=%d\n",
                cfl / lambda_max, cfl / rep.lambda, rep.lambda / lambda_max, lim,
                lambda_max, rep.lambda, rep.lambda_B, rep.lambda_C, rep.lambda_D,
                rep.lambda_t, rep.lambda_c, rep.lambda_v,
                (int)rep.ci, (int)rep.cj, (int)rep.ck, (int)rep.cblock, (int)rep.cclass, (int)rep.one_sided);
        }
    }
 
    LOG_ALLOW(GLOBAL, LOG_INFO,
        "Dual-time solver: lambda_max=%.4e [1/s]  dtau_init=%.4e  dtau range [%.4e, %.4e]  "
        "CFL range [%.4f, %.4f]  rejection_threshold=%.3f (EMA alpha=%.2f)  "
        "growth=%.3f  reduction=%.3f  max_accepted=%d.\n",
        lambda_max, cfl / lambda_max, dtau_min, dtau_max,
        simCtx->min_pseudo_cfl, simCtx->max_pseudo_cfl,
        simCtx->mom_dt_jameson_residual_norm_noise_allowance_factor,
        simCtx->mom_ratio_ema_alpha,
        simCtx->pseudo_cfl_growth_factor, simCtx->pseudo_cfl_reduction_factor,
        max_pseudo_steps);
 
    // --- 2. Main Pseudo-Time Iteration Loop ---
    pseudo_iter = 0;
    PetscBool residual_convergence_enabled =
        (PetscBool)(simCtx->mom_resid_atol > 0.0 || simCtx->mom_resid_rtol > 0.0);
    PetscBool converged = PETSC_FALSE;
    PetscBool last_trial_nonfinite = PETSC_FALSE;
    /* Rejected iterations no longer consume the accepted-iteration budget.
       A hard cap of 3x max_pseudo_steps guards against infinite rejection loops. */
    const PetscInt max_total_attempts = max_pseudo_steps * 3;
    while (!converged && accepted_iter < max_pseudo_steps && pseudo_iter < max_total_attempts)
    {
        pseudo_iter++;
        force_restart = PETSC_FALSE;
 
        // Every attempt starts from the last globally accepted state.
        for (PetscInt bi = 0; bi < block_number; bi++) {
            ierr = VecCopy(user[bi].pUcont, user[bi].Ucont); CHKERRQ(ierr);
            ierr = VecCopy(pRhs[bi], user[bi].Rhs); CHKERRQ(ierr);
            const char *staggered_fields[] = {"Ucont"};
            ierr = SynchronizePeriodicStaggeredFields(&user[bi], 1, staggered_fields); CHKERRQ(ierr);
            ierr = ApplyBoundaryConditions(&user[bi]); CHKERRQ(ierr);
        }
    
        for (PetscInt bi = 0; bi < block_number; bi++) {
            
            // === 4-Stage Jameson RK Smoothing Loop ===
            for (istage = 0; istage < 4; istage++) {
                
                LOG_ALLOW(GLOBAL,LOG_TRACE," Pseudo-Iter: %d | RK-Stage: %d\n", pseudo_iter, istage);
                
                /* RK Update: U_new = U_old + (dtau * alpha) * Residual
                 * pseudo_dtau[bi] IS dtau [physical time] — no further multiplication by dt.
                 * (In the old scheme dtau = pseudo_cfl * dt; here dtau = pseudo_cfl / lambda_max.) */
                ierr = VecWAXPY(user[bi].Ucont,
                                pseudo_dtau[bi] * alfa[istage],
                                user[bi].Rhs,
                                user[bi].pUcont); CHKERRQ(ierr);
 
                // Sync Ghosts & Re-apply BCs for intermediate stage
                const char *staggered_fields[] = {"Ucont"};
                ierr = SynchronizePeriodicStaggeredFields(&user[bi], 1, staggered_fields); CHKERRQ(ierr);
                
                ierr = ApplyBoundaryConditions(&user[bi]); CHKERRQ(ierr);
 
                // --- Re-calculate Total Residual for next stage ---
                ierr = ComputeTotalResidual(&user[bi]); CHKERRQ(ierr);
      
            } // End RK Stages
 
            // Immersed boundary interpolation (if enabled)
 
            // === Convergence Metrics Calculation ===
            
            // Calculate dU = U_current - U_backup
            ierr = VecWAXPY(user[bi].dUcont, -1.0, user[bi].pUcont, user[bi].Ucont); CHKERRQ(ierr);
 
            // Compute Infinity Norms
            ierr = VecNorm(user[bi].dUcont, NORM_INFINITY, &delta_sol_norm_curr[bi]); CHKERRQ(ierr);
            ierr = VecNorm(user[bi].Rhs, NORM_INFINITY, &resid_norm_curr[bi]); CHKERRQ(ierr);
        
            // Normalize relative metrics
            if (pseudo_iter == 1) {
                delta_sol_norm_init[bi] = delta_sol_norm_curr[bi];
                delta_sol_rel_curr[bi]  = 1.0;
                resid_rel_curr[bi]      = 1.0;
                // resid_norm_init[bi] set correctly before the loop — do not overwrite
            } else {
                if (delta_sol_norm_init[bi] > 1.0e-10) 
                    delta_sol_rel_curr[bi] = delta_sol_norm_curr[bi] / delta_sol_norm_init[bi];
                else 
                    delta_sol_rel_curr[bi] = 0.0;
 
                if(resid_norm_init[bi] > 1.0e-10) 
                    resid_rel_curr[bi] = resid_norm_curr[bi] / resid_norm_init[bi];
                else 
                    resid_rel_curr[bi] = 0.0;
            }
          
            /* Pre-compute per-block step-to-step trial ratio for later file logging. */
            {
                const PetscReal resid_floor_log = 1.0e-30;
                if (resid_norm_prev[bi] > resid_floor_log)
                    trial_ratio_log[bi] = resid_norm_curr[bi] / resid_norm_prev[bi];
                else if (resid_norm_curr[bi] <= resid_floor_log)
                    trial_ratio_log[bi] = 0.0;
                else
                    trial_ratio_log[bi] = PETSC_MAX_REAL;
            }
        } // End loop over blocks
 
        // --- Update Global Convergence Criteria ---
        global_norm_delta = -1.0e20; 
        global_rel_delta  = -1.0e20; 
        global_norm_resid = -1.0e20;
        global_rel_resid  = -1.0e20;
        
        for (PetscInt bi = 0; bi < block_number; bi++) {
            global_norm_delta = PetscMax(delta_sol_norm_curr[bi], global_norm_delta);
            global_rel_delta  = PetscMax(delta_sol_rel_curr[bi],  global_rel_delta);
            global_norm_resid = PetscMax(resid_norm_curr[bi],      global_norm_resid);
            global_rel_resid  = PetscMax(resid_rel_curr[bi],      global_rel_resid);
        }
        ierr = PetscTime(&te); CHKERRQ(ierr);
        cput = te - ts;
        LOG_ALLOW(GLOBAL, LOG_INFO, "  Pseudo-Iter(k) %d: |dUk|=%e, |dUk|/|dU0| = %e, |Rk|/|R0| = %e, CPU=%.2fs\n",
              pseudo_iter, global_norm_delta, global_rel_delta, global_rel_resid, cput);
 
        // === Adaptive Pseudo-CFL Trial Acceptance and Rollback ===
        const PetscReal resid_floor = 1.0e-30;
        PetscReal global_trial_ratio = 0.0;
        PetscBool global_nonfinite = PETSC_FALSE;
        for (PetscInt bi = 0; bi < block_number; bi++) {
            PetscReal ratio;
            if (resid_norm_prev[bi] > resid_floor) {
                ratio = resid_norm_curr[bi] / resid_norm_prev[bi];
            } else if (resid_norm_curr[bi] <= resid_floor) {
                ratio = 0.0;
            } else {
                ratio = PETSC_MAX_REAL;
            }
            global_trial_ratio = PetscMax(global_trial_ratio, ratio);
            global_nonfinite = (PetscBool)(global_nonfinite ||
                PetscIsInfOrNanReal(delta_sol_norm_curr[bi]) ||
                PetscIsInfOrNanReal(resid_norm_curr[bi]) ||
                PetscIsInfOrNanReal(ratio));
        }
 
        PetscReal reduced_trial_ratio;
        PetscMPIInt local_nonfinite = global_nonfinite ? 1 : 0;
        PetscMPIInt reduced_nonfinite = 0;
        ierr = MPI_Allreduce(&global_trial_ratio, &reduced_trial_ratio, 1, MPIU_REAL, MPI_MAX, PETSC_COMM_WORLD); CHKERRQ(ierr);
        ierr = MPI_Allreduce(&local_nonfinite, &reduced_nonfinite, 1, MPI_INT, MPI_LOR, PETSC_COMM_WORLD); CHKERRQ(ierr);
        global_trial_ratio = reduced_trial_ratio;
        global_nonfinite = reduced_nonfinite ? PETSC_TRUE : PETSC_FALSE;
        last_trial_nonfinite = global_nonfinite;
 
        /* EMA-smooth the trial ratio to tolerate occasional non-monotonic residual bumps.
           Non-finite trials bypass smoothing and always trigger rollback. */
        if (!global_nonfinite) {
            smoothed_trial_ratio = simCtx->mom_ratio_ema_alpha * global_trial_ratio
                                 + (1.0 - simCtx->mom_ratio_ema_alpha) * smoothed_trial_ratio;
        }
        LOG_ALLOW(GLOBAL, LOG_DEBUG,
                  "    [k=%d] raw_ratio=%.4e smoothed_ratio=%.4e (threshold=%.3f) | "
                  "|R_prev|=%.6e | |R_curr|=%.6e | CFL=%.6f\n",
                  pseudo_iter, global_trial_ratio, smoothed_trial_ratio,
                  simCtx->mom_dt_jameson_residual_norm_noise_allowance_factor,
                  resid_norm_prev[0], resid_norm_curr[0], pseudo_dtau[0]);
 
        if (simCtx->no_pseudo_cfl_backtrack) {
            /* Diagnostic mode: commit every finite trial; only non-finite triggers rollback. */
            force_restart = global_nonfinite;
        } else {
            force_restart = (PetscBool)(global_nonfinite ||
                smoothed_trial_ratio > simCtx->mom_dt_jameson_residual_norm_noise_allowance_factor);
        }
 
        if (force_restart) {
            /* old_dtau: dtau used for this (rejected) trial [physical time].
             * next_dtau: reduced dtau for the retry, clamped to dtau_min.
             * Both are physical-time values; the effective CFL = dtau * lambda_max. */
            PetscReal old_dtau  = pseudo_dtau[0];
            PetscReal next_dtau = PetscMax(dtau_min, old_dtau * simCtx->pseudo_cfl_reduction_factor);
            rejected_iter++;
            recovery_streak = 0;
            for (PetscInt bi = 0; bi < block_number; bi++) {
                ierr = VecCopy(user[bi].pUcont, user[bi].Ucont); CHKERRQ(ierr);
                ierr = VecCopy(pRhs[bi], user[bi].Rhs); CHKERRQ(ierr);
                const char *staggered_fields[] = {"Ucont"};
                ierr = SynchronizePeriodicStaggeredFields(&user[bi], 1, staggered_fields); CHKERRQ(ierr);
                ierr = ApplyBoundaryConditions(&user[bi]); CHKERRQ(ierr);
                pseudo_dtau[bi] = next_dtau;
            }
            LOG_ALLOW(GLOBAL, LOG_WARNING,
                      "  Trial %d REJECTED (raw_ratio=%.4e, smoothed=%.4e, nonfinite=%d); "
                      "dtau %.4e -> %.4e (cfl_eff %.4f -> %.4f)%s\n",
                      pseudo_iter, global_trial_ratio, smoothed_trial_ratio, (int)global_nonfinite,
                      old_dtau, next_dtau, old_dtau * lambda_max, next_dtau * lambda_max,
                      (old_dtau == next_dtau) ? " [AT FLOOR — no dtau reduction]" : "");
            /* Log rejected trial to file before rolling back. */
            if (!rank) {
                for (PetscInt bi = 0; bi < block_number; bi++) {
                    FILE *f;
                    char filen[PETSC_MAX_PATH_LEN + 128];
                    ierr = PetscSNPrintf(filen, sizeof(filen),
                        "%s/Momentum_Solver_DualTime_Picard_Jameson_RK_History_Block_%1d.log",
                        simCtx->log_dir, bi); CHKERRQ(ierr);
                    if (simCtx->step == simCtx->StartStep + 1 && pseudo_iter == 1 && !simCtx->continueMode)
                        f = fopen(filen, "w");
                    else
                        f = fopen(filen, "a");
                    if (simCtx->continueMode && simCtx->step == simCtx->StartStep + 1 && pseudo_iter == 1)
                        PetscFPrintf(PETSC_COMM_WORLD, f, "# Continuation from step %" PetscInt_FMT "\n", simCtx->StartStep);
                    /* Log format: dtau [physical time] + cfl_eff [dimensionless Courant number = dtau*lambda_max] */
                    PetscFPrintf(PETSC_COMM_WORLD, f,
                        "Step: %d | PseudoIter(k): %d | dtau: %.6e | cfl_eff: %.4f | |dUk|: %le | "
                        "|dUk|/|dU0|: %le | |Rk|: %le | |Rk|/|R0|: %le | "
                        "trial_ratio: %le | smoothed_ratio: %le | status: rejected | "
                        "dtau_after: %.6e | cfl_eff_after: %.4f\n",
                        (int)ti, (int)pseudo_iter, old_dtau, old_dtau * lambda_max,
                        delta_sol_norm_curr[bi], delta_sol_rel_curr[bi],
                        resid_norm_curr[bi], resid_rel_curr[bi],
                        trial_ratio_log[bi], smoothed_trial_ratio,
                        next_dtau, next_dtau * lambda_max);
                    fclose(f);
                }
            }
            if (old_dtau <= dtau_min) {
                if (global_nonfinite) {
                    SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_CONV_FAILED,
                            "Momentum solver produced a non-finite trial at minimum pseudo-CFL.");
                }
                /* Finite rejection at the dtau floor: reducing further is impossible.
                   Break instead of retrying bit-identically until max_pseudo_steps. */
                LOG_ALLOW(GLOBAL, LOG_WARNING,
                          "  Trial %d REJECTED (ratio=%.4e) at minimum dtau (%.4e, cfl_eff=%.4f) "
                          "with no further reduction possible; breaking retry loop.\n",
                          pseudo_iter, global_trial_ratio, old_dtau, old_dtau * lambda_max);
                break;
            }
            continue;
        }
 
        accepted_iter++;
        last_trial_nonfinite = PETSC_FALSE;
        last_accepted_resid  = resid_norm_curr[0]; /* track for final summary */
        for (PetscInt bi = 0; bi < block_number; bi++) {
            ierr = VecCopy(user[bi].Ucont, user[bi].pUcont); CHKERRQ(ierr);
            ierr = VecCopy(user[bi].Rhs, pRhs[bi]); CHKERRQ(ierr);
            resid_norm_prev[bi] = resid_norm_curr[bi];
            delta_sol_norm_prev[bi] = delta_sol_norm_curr[bi];
        }
 
        /* old_dtau: dtau used for this (accepted) trial.
         * next_dtau: dtau for the next trial, grown/reduced based on convergence rate.
         * Bounds (dtau_min, dtau_max) are derived from simCtx->min/max_pseudo_cfl / lambda_max. */
        PetscReal old_dtau  = pseudo_dtau[0]; /* save before any update */
        PetscReal next_dtau = old_dtau;
        if (global_trial_ratio < 0.90) {
            /* Residual decreasing fast: grow dtau to accelerate convergence. */
            next_dtau *= simCtx->pseudo_cfl_growth_factor;
            recovery_streak = 0;
        } else if (global_trial_ratio <= 1.0) {
            /* Residual barely decreasing: grow only after 3 clean consecutive trials. */
            recovery_streak++;
            if (recovery_streak >= 3) {
                next_dtau *= simCtx->pseudo_cfl_growth_factor;
                recovery_streak = 0;
            }
        } else {
            /* Accepted but residual slightly grew (within EMA noise allowance): reduce dtau. */
            next_dtau *= simCtx->pseudo_cfl_reduction_factor;
            recovery_streak = 0;
        }
        next_dtau = PetscMin(next_dtau, dtau_max);
        next_dtau = PetscMax(next_dtau, dtau_min);
        for (PetscInt bi = 0; bi < block_number; bi++) pseudo_dtau[bi] = next_dtau;
 
        LOG_ALLOW(GLOBAL, LOG_DEBUG,
                  "  Trial %d ACCEPTED (raw_ratio=%.4e, smoothed=%.4e); |dU|=%.6e | "
                  "dtau %.4e -> %.4e (cfl_eff %.4f -> %.4f)\n",
                  pseudo_iter, global_trial_ratio, smoothed_trial_ratio, global_norm_delta,
                  old_dtau, next_dtau, old_dtau * lambda_max, next_dtau * lambda_max);
 
        /* --- Post-decision file logging (one row per accepted trial) --- */
        if (!rank) {
            for (PetscInt bi = 0; bi < block_number; bi++) {
                FILE *f;
                char filen[PETSC_MAX_PATH_LEN + 128];
                ierr = PetscSNPrintf(filen, sizeof(filen),
                    "%s/Momentum_Solver_DualTime_Picard_Jameson_RK_History_Block_%1d.log",
                    simCtx->log_dir, bi); CHKERRQ(ierr);
                if (simCtx->step == simCtx->StartStep + 1 && pseudo_iter == 1 && !simCtx->continueMode)
                    f = fopen(filen, "w");
                else
                    f = fopen(filen, "a");
                if (simCtx->continueMode && simCtx->step == simCtx->StartStep + 1 && pseudo_iter == 1)
                    PetscFPrintf(PETSC_COMM_WORLD, f, "# Continuation from step %" PetscInt_FMT "\n", simCtx->StartStep);
                /* Log format: dtau [physical time] + cfl_eff [dimensionless Courant number = dtau*lambda_max] */
                PetscFPrintf(PETSC_COMM_WORLD, f,
                    "Step: %d | PseudoIter(k): %d | dtau: %.6e | cfl_eff: %.4f | |dUk|: %le | "
                    "|dUk|/|dU0|: %le | |Rk|: %le | |Rk|/|R0|: %le | "
                    "trial_ratio: %le | smoothed_ratio: %le | status: accepted | "
                    "dtau_after: %.6e | cfl_eff_after: %.4f\n",
                    (int)ti, (int)pseudo_iter, old_dtau, old_dtau * lambda_max,
                    delta_sol_norm_curr[bi], delta_sol_rel_curr[bi],
                    resid_norm_curr[bi], resid_rel_curr[bi],
                    trial_ratio_log[bi], smoothed_trial_ratio,
                    next_dtau, next_dtau * lambda_max);
                fclose(f);
            }
        }
 
        PetscBool update_pass;
        PetscBool residual_pass = PETSC_TRUE;
        if (residual_convergence_enabled) {
            update_pass = (PetscBool)(global_norm_delta <= tol_abs_delta && global_rel_delta <= tol_rtol_delta);
            residual_pass = (PetscBool)(
                (simCtx->mom_resid_atol > 0.0 && global_norm_resid <= simCtx->mom_resid_atol) ||
                (simCtx->mom_resid_rtol > 0.0 && global_rel_resid <= simCtx->mom_resid_rtol));
        } else {
            update_pass = (PetscBool)(global_norm_delta <= tol_abs_delta && global_rel_delta <= tol_rtol_delta);
        }
        converged = (PetscBool)(update_pass && residual_pass);
 
        if (block_number > 1) {
             // ierr = Block_Interface_U(user); CHKERRQ(ierr);
        }
    } // End while loop
 
    if (last_trial_nonfinite) {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_CONV_FAILED,
                "Momentum solver exhausted its attempt limit while recovering from a non-finite trial.");
    }
 
    /* Convert the final dtau back to a CFL number for the warm-start of the next physical timestep.
     * simCtx->pseudo_cfl stores a dimensionless Courant number (dtau * lambda_max); the next
     * timestep's lambda_max may differ as the flow evolves, so storing the CFL number (not dtau)
     * ensures the warm-start adapts correctly to the new spectral radius. */
    PetscReal next_dtau_start  = pseudo_dtau[0];
    PetscReal next_cfl_warmstart = next_dtau_start * lambda_max;  /* CFL = dtau * lambda_max */
 
    if (!rank) {
        LOG_ALLOW(GLOBAL, LOG_INFO,
                  "  Step %d finished: dtau=%.4e  cfl_eff=%.4f  lambda_max=%.4e [1/s]. "
                  "Next step warm-starts at cfl=%.4f.\n",
                  simCtx->step, next_dtau_start, next_cfl_warmstart, lambda_max, next_cfl_warmstart);
    }
 
    /* Store the CFL number (not dtau) so it remains meaningful across physical timesteps. */
    simCtx->pseudo_cfl        = next_cfl_warmstart;
    simCtx->mom_last_lambda_max = lambda_max;    
 
    /* Expose the last accepted finite state. On converged or floor-break exits pUcont==Ucont
       already (VecCopy would be an identity), but ghosts may have gone stale; always re-sync
       and re-apply BCs. Only restore the solution vector itself when nothing was ever accepted. */
    if (accepted_iter == 0) {
        for (PetscInt bi = 0; bi < block_number; bi++) {
            ierr = VecCopy(user[bi].pUcont, user[bi].Ucont); CHKERRQ(ierr);
            ierr = VecCopy(pRhs[bi], user[bi].Rhs); CHKERRQ(ierr);
        }
    }
    for (PetscInt bi = 0; bi < block_number; bi++) {
        const char *staggered_fields[] = {"Ucont"};
        ierr = SynchronizePeriodicStaggeredFields(&user[bi], 1, staggered_fields); CHKERRQ(ierr);
        ierr = ApplyBoundaryConditions(&user[bi]); CHKERRQ(ierr);
    }
 
    simCtx->mom_last_converged = converged;
    if (!converged) {
        PetscPrintf(PETSC_COMM_WORLD,
            "[WARNING] Momentum solver step %d: reached %d total attempts (%d accepted, %d rejected) "
            "without convergence; continuing from last accepted finite state.\n",
            (int)ti, pseudo_iter, accepted_iter, rejected_iter);
    }
    if (accepted_iter == 0) {
        PetscPrintf(PETSC_COMM_WORLD,
            "[WARNING] Momentum solver step %d: no pseudo-time trials were accepted; "
            "retaining physical-step entry state.\n", (int)ti);
    }
    LOG_ALLOW(GLOBAL, LOG_INFO,
              "Momentum solver finished: %d attempts (%d accepted, %d rejected) of %d max accepted / %d hard cap, "
              "converged=%s, last accepted |R|=%.6e, last accepted |dU|=%.6e, "
              "next_dtau=%.4e (cfl=%.4f).\n",
              pseudo_iter, accepted_iter, rejected_iter, max_pseudo_steps, max_total_attempts,
              converged ? "yes" : "no",
              (accepted_iter > 0) ? last_accepted_resid : resid_norm_init[0],
              (accepted_iter > 0) ? delta_sol_norm_prev[0] : 0.0,
              next_dtau_start, next_cfl_warmstart);
 
    // --- Final Cleanup ---
    for (PetscInt bi = 0; bi < block_number; bi++) {
        ierr = VecDestroy(&user[bi].Rhs); CHKERRQ(ierr);
        ierr = VecDestroy(&user[bi].dUcont); CHKERRQ(ierr);
        ierr = VecDestroy(&user[bi].pUcont); CHKERRQ(ierr);
        ierr = VecDestroy(&pRhs[bi]); CHKERRQ(ierr);
    }
    ierr = PetscFree(pRhs); CHKERRQ(ierr);
 
    ierr = PetscFree2(delta_sol_norm_init, delta_sol_norm_prev);CHKERRQ(ierr);
    ierr = PetscFree2(delta_sol_norm_curr, delta_sol_rel_curr);CHKERRQ(ierr);
    ierr = PetscFree2(resid_norm_init,     resid_norm_prev);CHKERRQ(ierr);
    ierr = PetscFree2(resid_norm_curr,     pseudo_dtau); CHKERRQ(ierr);
    ierr = PetscFree(resid_rel_curr); CHKERRQ(ierr);
    ierr = PetscFree(trial_ratio_log); CHKERRQ(ierr);
    
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

MomentumUsesBDF2()

PetscBool MomentumUsesBDF2

(

SimCtx *

simCtx

)

Returns whether the current physical step uses the BDF2 discretization.

Single source of truth for the BDF1/BDF2 selection. The predicate is identical to the one historically inlined in ComputeTotalResidual(): BDF2 when COEF_TIME_ACCURACY > 1.1 AND step != StartStep AND step != 1, otherwise BDF1 (startup step and the first step after a restart).

Parameters

simCtx

Master simulation context (reads step, StartStep).

Returns

PETSC_TRUE for BDF2, PETSC_FALSE for BDF1.

Returns whether the current physical step uses the BDF2 discretization.

See header.

This reproduces exactly the predicate historically inlined in ComputeTotalResidual(): BDF1 on the cold-start first step (ti == 1) and on the first step after a restart (ti == tistart); BDF2 thereafter (when the second- order time accuracy is enabled). Restart-history validity is a separate concern and is intentionally not addressed here.

Definition at line 19 of file momentumsolvers.c.

{
    const PetscInt ti      = simCtx->step;
    const PetscInt tistart = simCtx->StartStep;
    return (PetscBool)(COEF_TIME_ACCURACY > 1.1 && ti != tistart && ti != 1);
}

Here is the caller graph for this function:

MomentumBDFCoefficient()

PetscReal MomentumBDFCoefficient

(

SimCtx *

simCtx

)

Returns the BDF physical-time coefficient a0 for the current step.

a0 = 1.5 (== COEF_TIME_ACCURACY) for BDF2, a0 = 1.0 for BDF1. Used both as the leading coefficient of the physical-time term in the residual and as the additive temporal contribution lambda_t = a0/dt in the momentum stability estimate, keeping the two numerically consistent.

Parameters

simCtx

Master simulation context.

Returns

a0 in {1.0, 1.5}.

Returns the BDF physical-time coefficient a0 for the current step.

See header.

Definition at line 31 of file momentumsolvers.c.

{
    return MomentumUsesBDF2(simCtx) ? COEF_TIME_ACCURACY : 1.0;
}

Here is the call graph for this function:

Here is the caller graph for this function:

ComputeTotalResidual()

PetscErrorCode ComputeTotalResidual

(

UserCtx *

user

)

Computes the shared spatial-plus-BDF momentum residual in user->Rhs.

Parameters

user	Block context with an allocated Rhs vector.

Returns

PetscErrorCode 0 on success.

Computes the shared spatial-plus-BDF momentum residual in user->Rhs.

Definition at line 42 of file momentumsolvers.c.

{
    PetscErrorCode ierr;
    SimCtx *simCtx = user->simCtx;
    
    // Extract Time Parameters from Context.
    // BDF step selection now lives in MomentumUsesBDF2() (shared with the stability estimate).
    const PetscReal dt      = simCtx->dt;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    // 1. Calculate Spatial Terms (stored in user->Rhs)
    //    Rhs = -Div(Flux) + Viscous + Source
    ierr = ComputeRHS(user, user->Rhs); CHKERRQ(ierr);
 
    // 2. Add Physical Time Derivative Terms (BDF Discretization)
    //    The equation solved is: dU/dtau = RHS_Spatial + RHS_Temporal
    
    /* BDF order + leading coefficient from the shared helpers (single source of truth,
       shared with the momentum stability estimate). a0 = 1.5 (==COEF_TIME_ACCURACY) for
       BDF2, 1.0 for BDF1; using a0 for the current-state term keeps the residual
       numerically identical to the historical inlined coefficients. */
    const PetscBool use_bdf2 = MomentumUsesBDF2(simCtx);
    const PetscReal a0       = MomentumBDFCoefficient(simCtx);
    if (use_bdf2) {
        // --- BDF2 (Second Order Backward Difference) ---
        // (a0*U^{n} - 2.0*U^{n-1} + 0.5*U^{n-2}) / dt = RHS_Spatial(U^{n})
        ierr = VecAXPY(user->Rhs, -a0/dt,  user->Ucont);     CHKERRQ(ierr);
        ierr = VecAXPY(user->Rhs, +2.0/dt, user->Ucont_o);   CHKERRQ(ierr);
        ierr = VecAXPY(user->Rhs, -0.5/dt, user->Ucont_rm1); CHKERRQ(ierr);
    } else {
        // --- BDF1 (First Order / Euler Implicit) ---
        // (a0*U^{n} - U^{n-1}) / dt = RHS_Spatial(U^{n}), with a0 = 1.0
        ierr = VecAXPY(user->Rhs, -a0/dt,  user->Ucont);     CHKERRQ(ierr);
        ierr = VecAXPY(user->Rhs, +1.0/dt, user->Ucont_o);   CHKERRQ(ierr);
    }
 
    // 3. Enforce Boundary Conditions on the Residual
    ierr = EnforceRHSBoundaryConditions(user); CHKERRQ(ierr);
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

ComputeMomentumStabilityEstimate()

PetscErrorCode ComputeMomentumStabilityEstimate	(	UserCtx *	user,
		PetscInt	block_number,
		PetscReal	dt,
		MomStabCandidate	candidate,
		MomStabilityReport *	rep
	)

Compute the momentum pseudo-time stability estimate (shadow/diagnostic).

Conservative, operator-scaled estimate: lambda = max_cell (a0/dt + lambda_c + lambda_nu), over active, non-solid cells, blocks, and MPI ranks, where lambda_c already includes the per-direction QUICK scheme factors. Read-only; performs no halo exchange, but does perform global scalar collectives (see implementation). This is a PRACTICAL CONSERVATIVE estimate, not a proven spectral radius. See the Workstream-A design.

Call-site preconditions (NOT enforced internally): lUcont, lUcat, lNu_t, lNvert must be fresh; lAj, face Jacobians and face metrics are static after grid init.

Parameters

[in]	user	Array of UserCtx (one per block).
[in]	block_number	Number of blocks.
[in]	dt	Physical time step.
[in]	candidate	Convective candidate driving rep->lambda (B, C or D).
[out]	rep	Filled diagnostic report (global maxima + breakdown).

Returns

PetscErrorCode 0 on success.

Compute the momentum pseudo-time stability estimate (shadow/diagnostic).

See header.

Operator-scaled estimate lambda = max_cell (a0/dt + lambda_c + lambda_nu) over active, non-solid interior cells, blocks and MPI ranks (lambda_c already carries the per-direction QUICK scheme factors – it is NOT multiplied by f_c again here). Read-only; no halo exchange, 5 global scalar collectives. NOT a proven spectral radius. Drives nothing in shadow mode.

Definition at line 396 of file momentumsolvers.c.

{
    PetscErrorCode ierr;
    PetscMPIInt rank;
    PetscFunctionBeginUser;
 
    /* ---- input validation (no silent fallback) ---- */
    PetscCheck(user != NULL, PETSC_COMM_WORLD, PETSC_ERR_ARG_NULL, "user is NULL");
    PetscCheck(rep  != NULL, PETSC_COMM_WORLD, PETSC_ERR_ARG_NULL, "rep is NULL");
    PetscCheck(block_number > 0, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "block_number must be > 0");
    PetscCheck(candidate >= MOM_STAB_CAND_B && candidate <= MOM_STAB_CAND_D,
               PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "unsupported stability candidate enum");
    PetscCheck(PetscIsNormalReal(dt) && dt > 0.0, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE,
               "dt must be finite and positive");
 
    SimCtx         *simCtx   = user[0].simCtx;
    PetscCheck(simCtx != NULL, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "user[0].simCtx is NULL");
    const PetscReal a0       = MomentumBDFCoefficient(simCtx);
    const PetscBool centered = (PetscBool)(simCtx->les || simCtx->central);
    const PetscBool has_nut  = (PetscBool)(simCtx->les || simCtx->rans);
    const PetscBool inviscid = (PetscBool)simCtx->invicid;
    PetscCheck(inviscid || (PetscIsNormalReal(simCtx->ren) && simCtx->ren > 0.0),
               PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Reynolds number must be finite and positive when viscous");
    const PetscReal lambda_t = a0 / dt;
    const PetscReal nu_mol   = inviscid ? 0.0 : 1.0 / simCtx->ren;
    const PetscInt  twoD     = simCtx->TwoD;
 
    /* Shadow-mode completeness flag. The estimate does NOT cover the Clark nonlinear stress
       Jacobian, nor RANS eddy-viscosity sign behaviour (not verified sign-definite).
       Body forces in the supported configs read a per-timestep-frozen scalar
       (simCtx->bulkVelocityCorrection) inside ComputeRHS, so within a pseudo-solve they are a
       constant forcing with ZERO velocity Jacobian (consistent with the frozen-pressure
       treatment) -- they do not make the estimate incomplete. */
    rep->estimate_incomplete = (PetscBool)(simCtx->clark || simCtx->rans);
 
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
 
    /* Per-candidate local maxima; selected candidate also tracks the controlling cell. */
    PetscReal locmaxB = 0.0, locmaxC = 0.0, locmaxD = 0.0, sel_max = 0.0;
    PetscReal sel_lc = 0.0, sel_lv = 0.0;
    PetscInt  sel_ci = -1, sel_cj = -1, sel_ck = -1, sel_blk = -1, sel_class = 0, sel_os = 0;
    PetscInt  local_active = 0;
 
    for (PetscInt bi = 0; bi < block_number; bi++) {
        UserCtx *u = &user[bi];
        DMDALocalInfo info = u->info;
        const PetscInt mx = info.mx, my = info.my, mz = info.mz;
        const PetscInt xs = info.xs, xe = xs + info.xm;
        const PetscInt ys = info.ys, ye = ys + info.ym;
        const PetscInt zs = info.zs, ze = zs + info.zm;
        /* Interior range == where ComputeRHS evaluates the residual; boundary faces
           (0 and m-1) are RHS-masked, so excluding them implements active-row masking
           for the normal boundary rows. The last physical face (m-2) stays IN with the
           full cell estimate (its tangential rows remain active). */
        const PetscInt lxs = (xs==0)?xs+1:xs, lxe = (xe==mx)?xe-1:xe;
        const PetscInt lys = (ys==0)?ys+1:ys, lye = (ye==my)?ye-1:ye;
        const PetscInt lzs = (zs==0)?zs+1:zs, lze = (ze==mz)?ze-1:ze;
 
        /* Non-periodic boundary flags (QUICK stencil is degraded near these). */
        const PetscBool npx0 = (PetscBool)(u->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC);
        const PetscBool npx1 = (PetscBool)(u->boundary_faces[BC_FACE_POS_X].mathematical_type != PERIODIC);
        const PetscBool npy0 = (PetscBool)(u->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC);
        const PetscBool npy1 = (PetscBool)(u->boundary_faces[BC_FACE_POS_Y].mathematical_type != PERIODIC);
        const PetscBool npz0 = (PetscBool)(u->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC);
        const PetscBool npz1 = (PetscBool)(u->boundary_faces[BC_FACE_POS_Z].mathematical_type != PERIODIC);
        const PetscBool wxn = (PetscBool)(u->boundary_faces[BC_FACE_NEG_X].mathematical_type == WALL);
        const PetscBool wxp = (PetscBool)(u->boundary_faces[BC_FACE_POS_X].mathematical_type == WALL);
        const PetscBool wyn = (PetscBool)(u->boundary_faces[BC_FACE_NEG_Y].mathematical_type == WALL);
        const PetscBool wyp = (PetscBool)(u->boundary_faces[BC_FACE_POS_Y].mathematical_type == WALL);
        const PetscBool wzn = (PetscBool)(u->boundary_faces[BC_FACE_NEG_Z].mathematical_type == WALL);
        const PetscBool wzp = (PetscBool)(u->boundary_faces[BC_FACE_POS_Z].mathematical_type == WALL);
 
        Cmpnts ***ucont, ***ucat = NULL;
        Cmpnts ***icsi, ***ieta, ***izet, ***jcsi, ***jeta, ***jzet, ***kcsi, ***keta, ***kzet;
        Cmpnts ***csi = NULL, ***eta = NULL, ***zet = NULL;
        PetscReal ***aj, ***iaj, ***jaj, ***kaj, ***nvert, ***nut = NULL;
 
        ierr = DMDAVecGetArrayRead(u->fda, u->lUcont, &ucont); CHKERRQ(ierr);
        ierr = DMDAVecGetArrayRead(u->da,  u->lAj,   &aj);     CHKERRQ(ierr);
        ierr = DMDAVecGetArrayRead(u->da,  u->lIAj,  &iaj);    CHKERRQ(ierr);
        ierr = DMDAVecGetArrayRead(u->da,  u->lJAj,  &jaj);    CHKERRQ(ierr);
        ierr = DMDAVecGetArrayRead(u->da,  u->lKAj,  &kaj);    CHKERRQ(ierr);
        ierr = DMDAVecGetArrayRead(u->fda, u->lICsi, &icsi);   CHKERRQ(ierr);
        ierr = DMDAVecGetArrayRead(u->fda, u->lIEta, &ieta);   CHKERRQ(ierr);
        ierr = DMDAVecGetArrayRead(u->fda, u->lIZet, &izet);   CHKERRQ(ierr);
        ierr = DMDAVecGetArrayRead(u->fda, u->lJCsi, &jcsi);   CHKERRQ(ierr);
        ierr = DMDAVecGetArrayRead(u->fda, u->lJEta, &jeta);   CHKERRQ(ierr);
        ierr = DMDAVecGetArrayRead(u->fda, u->lJZet, &jzet);   CHKERRQ(ierr);
        ierr = DMDAVecGetArrayRead(u->fda, u->lKCsi, &kcsi);   CHKERRQ(ierr);
        ierr = DMDAVecGetArrayRead(u->fda, u->lKEta, &keta);   CHKERRQ(ierr);
        ierr = DMDAVecGetArrayRead(u->fda, u->lKZet, &kzet);   CHKERRQ(ierr);
        ierr = DMDAVecGetArrayRead(u->da,  u->lNvert, &nvert); CHKERRQ(ierr);
        if (has_nut) { ierr = DMDAVecGetArrayRead(u->da, u->lNu_t, &nut); CHKERRQ(ierr); }
        /* Cartesian velocity + cell metrics: always read so candidate D's gradient term
           (and hence rep->lambda_D) is meaningful regardless of the selected candidate.
           Precondition: lUcat fresh (Contra2Cart has run in the preceding residual eval). */
        ierr = DMDAVecGetArrayRead(u->fda, u->lUcat, &ucat); CHKERRQ(ierr);
        ierr = DMDAVecGetArrayRead(u->fda, u->lCsi,  &csi);  CHKERRQ(ierr);
        ierr = DMDAVecGetArrayRead(u->fda, u->lEta,  &eta);  CHKERRQ(ierr);
        ierr = DMDAVecGetArrayRead(u->fda, u->lZet,  &zet);  CHKERRQ(ierr);
 
        /* Error flag for the cell loop: on a bad metric/estimate we record the location and
           jump to block_cleanup so EVERY acquired array is restored before the error returns. */
        PetscErrorCode cell_err = 0;
        PetscInt       be_i = -1, be_j = -1, be_k = -1;
        const char    *be_what = NULL;
 
        for (PetscInt k = lzs; k < lze; k++) {
            for (PetscInt j = lys; j < lye; j++) {
                for (PetscInt i = lxs; i < lxe; i++) {
                    /* ---- active-row mask: skip only locations with no active unknown ---- */
                    const PetscInt rows = MomCellActiveRows(nvert, k, j, i, mx, my, mz,
                                                            npx1, npy1, npz1, twoD);
                    if (rows == 0) continue;
                    local_active++;
                    const PetscReal Ajc = aj[k][j][i];
                    if (!(PetscIsNormalReal(Ajc) && Ajc > 0.0)) {
                        cell_err = PETSC_ERR_FP; be_i = i; be_j = j; be_k = k;
                        be_what = "non-finite/non-positive cell inverse-Jacobian";
                        goto block_cleanup;
                    }
 
                    /* ---- directional QUICK modification (per direction; both faces; 2nd neighbour) ---- */
                    const PetscBool mod_x = MomQuickDirModified(nvert, k, j, i, i, mx, npx0, npx1,
                                                               info.gxs, info.gxs+info.gxm, 'i');
                    const PetscBool mod_y = MomQuickDirModified(nvert, k, j, i, j, my, npy0, npy1,
                                                               info.gys, info.gys+info.gym, 'j');
                    const PetscBool mod_z = MomQuickDirModified(nvert, k, j, i, k, mz, npz0, npz1,
                                                               info.gzs, info.gzs+info.gzm, 'k');
 
                    /* ---- classify cell: interior / physical-boundary / IB-adjacent ---- */
                    const PetscBool bnd = (PetscBool)((npx0 && i<=1) || (npx1 && i>=mx-2) ||
                                                      (npy0 && j<=1) || (npy1 && j>=my-2) ||
                                                      (npz0 && k<=1) || (npz1 && k>=mz-2));
                    const PetscBool ib  = MomCellHasSolidNeighbor(nvert, k, j, i);
                    const PetscInt cls = ib ? 2 : (bnd ? 1 : 0);
 
                    /* ---- convective: six contravariant face fluxes, DIRECTIONAL QUICK factors ---- */
                    const PetscReal Uxp = ucont[k][j][i].x,   Uxm = ucont[k][j][i-1].x;
                    const PetscReal Uyp = ucont[k][j][i].y,   Uym = ucont[k][j-1][i].y;
                    const PetscReal Uzp = ucont[k][j][i].z,   Uzm = ucont[k-1][j][i].z;
                    const PetscReal divU = PetscAbsReal((Uxp-Uxm)+(Uyp-Uym)+(Uzp-Uzm));
                    const PetscReal fx = centered ? 1.0 : (mod_x ? 2.5 : (4.0/3.0));
                    const PetscReal fy = centered ? 1.0 : (mod_y ? 2.5 : (4.0/3.0));
                    const PetscReal fz = centered ? 1.0 : (mod_z ? 2.5 : (4.0/3.0));
 
                    /* per-direction transport scale; C's divergence term stays unscaled by f_d. */
                    const PetscReal lcB = 0.5 * Ajc * (
                          fx * (PetscAbsReal(Uxp)+PetscAbsReal(Uxm))
                        + fy * (PetscAbsReal(Uyp)+PetscAbsReal(Uym))
                        + fz * (PetscAbsReal(Uzp)+PetscAbsReal(Uzm)) );
                    const PetscReal lcC = lcB + 0.5 * Ajc * divU;
                    /* lambda_grad_u = max_i sum_j |du_i/dx_j| from fresh Cartesian velocity
                       (the nonlinear zero-order Jacobian term delta_u . grad u).
                       Physical gradient: d/dx_dir = Ajc*(Csi.dir d/dxi + Eta.dir d/deta + Zet.dir d/dzeta). */
                    const Cmpnts duc = { 0.5*(ucat[k][j][i+1].x-ucat[k][j][i-1].x),
                                         0.5*(ucat[k][j][i+1].y-ucat[k][j][i-1].y),
                                         0.5*(ucat[k][j][i+1].z-ucat[k][j][i-1].z) };
                    const Cmpnts due = { 0.5*(ucat[k][j+1][i].x-ucat[k][j-1][i].x),
                                         0.5*(ucat[k][j+1][i].y-ucat[k][j-1][i].y),
                                         0.5*(ucat[k][j+1][i].z-ucat[k][j-1][i].z) };
                    const Cmpnts duz = { 0.5*(ucat[k+1][j][i].x-ucat[k-1][j][i].x),
                                         0.5*(ucat[k+1][j][i].y-ucat[k-1][j][i].y),
                                         0.5*(ucat[k+1][j][i].z-ucat[k-1][j][i].z) };
                    const Cmpnts C = csi[k][j][i], E = eta[k][j][i], Z = zet[k][j][i];
                    #define MOM_ROW(cmp) ( \
                        PetscAbsReal(Ajc*(C.x*duc.cmp + E.x*due.cmp + Z.x*duz.cmp)) + \
                        PetscAbsReal(Ajc*(C.y*duc.cmp + E.y*due.cmp + Z.y*duz.cmp)) + \
                        PetscAbsReal(Ajc*(C.z*duc.cmp + E.z*due.cmp + Z.z*duz.cmp)) )
                    const PetscReal lgrad = PetscMax(MOM_ROW(x), PetscMax(MOM_ROW(y), MOM_ROW(z)));
                    #undef MOM_ROW
                    const PetscReal lcD = lcC + lgrad;
 
                    /* ---- viscous: six faces, face Jacobians, full metric rows ---- */
                    PetscReal lv = 0.0;
                    if (!inviscid) {
                        /* xi+ (i), xi- (i-1) */
                        for (PetscInt s = 0; s < 2; s++) {
                            const PetscInt fi = (s==0) ? i : i-1;
                            if (!(PetscIsNormalReal(iaj[k][j][fi]) && iaj[k][j][fi] > 0.0)) {
                                cell_err = PETSC_ERR_FP; be_i=i; be_j=j; be_k=k;
                                be_what = "non-finite/non-positive xi-face inverse-Jacobian"; goto block_cleanup; }
                            PetscReal nuf = nu_mol;
                            if (has_nut) {
                                PetscReal nt = 0.5*(nut[k][j][fi] + nut[k][j][fi+1]);
                                if ((wxn && fi==0) || (wxp && fi==mx-2)) nt = 0.0;
                                nuf += nt;
                            }
                            lv += PetscAbsReal(nuf) * Ajc * iaj[k][j][fi]
                                * MomFaceGabs(icsi[k][j][fi], ieta[k][j][fi], izet[k][j][fi]);
                        }
                        /* eta+ (j), eta- (j-1) */
                        for (PetscInt s = 0; s < 2; s++) {
                            const PetscInt fj = (s==0) ? j : j-1;
                            if (!(PetscIsNormalReal(jaj[k][fj][i]) && jaj[k][fj][i] > 0.0)) {
                                cell_err = PETSC_ERR_FP; be_i=i; be_j=j; be_k=k;
                                be_what = "non-finite/non-positive eta-face inverse-Jacobian"; goto block_cleanup; }
                            PetscReal nuf = nu_mol;
                            if (has_nut) {
                                PetscReal nt = 0.5*(nut[k][fj][i] + nut[k][fj+1][i]);
                                if ((wyn && fj==0) || (wyp && fj==my-2)) nt = 0.0;
                                nuf += nt;
                            }
                            /* eta-face normal is the eta metric (jeta) -> pass as N. */
                            lv += PetscAbsReal(nuf) * Ajc * jaj[k][fj][i]
                                * MomFaceGabs(jeta[k][fj][i], jcsi[k][fj][i], jzet[k][fj][i]);
                        }
                        /* zeta+ (k), zeta- (k-1) */
                        for (PetscInt s = 0; s < 2; s++) {
                            const PetscInt fk = (s==0) ? k : k-1;
                            if (!(PetscIsNormalReal(kaj[fk][j][i]) && kaj[fk][j][i] > 0.0)) {
                                cell_err = PETSC_ERR_FP; be_i=i; be_j=j; be_k=k;
                                be_what = "non-finite/non-positive zeta-face inverse-Jacobian"; goto block_cleanup; }
                            PetscReal nuf = nu_mol;
                            if (has_nut) {
                                PetscReal nt = 0.5*(nut[fk][j][i] + nut[fk+1][j][i]);
                                if ((wzn && fk==0) || (wzp && fk==mz-2)) nt = 0.0;
                                nuf += nt;
                            }
                            /* zeta-face normal is the zeta metric (kzet) -> pass as N. */
                            lv += PetscAbsReal(nuf) * Ajc * kaj[fk][j][i]
                                * MomFaceGabs(kzet[fk][j][i], kcsi[fk][j][i], keta[fk][j][i]);
                        }
                        lv *= 4.0;   /* 2 (two-face row-sum) x 2 (longitudinal full-stress) */
                    }
 
                    /* one-sided viscous cross-derivative branch near a solid: conservative x2,
                       applied ONCE per cell, after a complete 3x3x3 solid-band check. */
                    PetscInt one_sided = 0;
                    if (!inviscid && MomCellUsesOneSidedViscousStencil(nvert, k, j, i)) {
                        one_sided = 1;
                        lv *= 2.0;
                    }
 
                    /* ---- candidate totals and running maxima ---- */
                    const PetscReal lB = lambda_t + lcB + lv;
                    const PetscReal lC = lambda_t + lcC + lv;
                    const PetscReal lD = lambda_t + lcD + lv;
                    if (PetscIsInfOrNanReal(lD)) {
                        cell_err = PETSC_ERR_FP; be_i=i; be_j=j; be_k=k;
                        be_what = "non-finite stability estimate"; goto block_cleanup;
                    }
                    locmaxB = PetscMax(locmaxB, lB);
                    locmaxC = PetscMax(locmaxC, lC);
                    locmaxD = PetscMax(locmaxD, lD);
 
                    const PetscReal lsel  = (candidate==MOM_STAB_CAND_B)?lB:(candidate==MOM_STAB_CAND_C)?lC:lD;
                    const PetscReal lcsel = (candidate==MOM_STAB_CAND_B)?lcB:(candidate==MOM_STAB_CAND_C)?lcC:lcD;
                    if (lsel > sel_max) {
                        sel_max = lsel; sel_lc = lcsel; sel_lv = lv;
                        sel_ci = i; sel_cj = j; sel_ck = k; sel_blk = bi;
                        sel_class = cls; sel_os = one_sided;
                    }
                }
            }
        }
 
        /* Single cleanup point: reached on normal completion AND on any cell-loop error.
           Every successful GetArrayRead above has exactly one matching restore here. */
        block_cleanup:
        ierr = DMDAVecRestoreArrayRead(u->fda, u->lUcont, &ucont); CHKERRQ(ierr);
        ierr = DMDAVecRestoreArrayRead(u->da,  u->lAj,   &aj);     CHKERRQ(ierr);
        ierr = DMDAVecRestoreArrayRead(u->da,  u->lIAj,  &iaj);    CHKERRQ(ierr);
        ierr = DMDAVecRestoreArrayRead(u->da,  u->lJAj,  &jaj);    CHKERRQ(ierr);
        ierr = DMDAVecRestoreArrayRead(u->da,  u->lKAj,  &kaj);    CHKERRQ(ierr);
        ierr = DMDAVecRestoreArrayRead(u->fda, u->lICsi, &icsi);   CHKERRQ(ierr);
        ierr = DMDAVecRestoreArrayRead(u->fda, u->lIEta, &ieta);   CHKERRQ(ierr);
        ierr = DMDAVecRestoreArrayRead(u->fda, u->lIZet, &izet);   CHKERRQ(ierr);
        ierr = DMDAVecRestoreArrayRead(u->fda, u->lJCsi, &jcsi);   CHKERRQ(ierr);
        ierr = DMDAVecRestoreArrayRead(u->fda, u->lJEta, &jeta);   CHKERRQ(ierr);
        ierr = DMDAVecRestoreArrayRead(u->fda, u->lJZet, &jzet);   CHKERRQ(ierr);
        ierr = DMDAVecRestoreArrayRead(u->fda, u->lKCsi, &kcsi);   CHKERRQ(ierr);
        ierr = DMDAVecRestoreArrayRead(u->fda, u->lKEta, &keta);   CHKERRQ(ierr);
        ierr = DMDAVecRestoreArrayRead(u->fda, u->lKZet, &kzet);   CHKERRQ(ierr);
        ierr = DMDAVecRestoreArrayRead(u->da,  u->lNvert, &nvert); CHKERRQ(ierr);
        if (has_nut) { ierr = DMDAVecRestoreArrayRead(u->da, u->lNu_t, &nut); CHKERRQ(ierr); }
        ierr = DMDAVecRestoreArrayRead(u->fda, u->lUcat, &ucat); CHKERRQ(ierr);
        ierr = DMDAVecRestoreArrayRead(u->fda, u->lCsi,  &csi);  CHKERRQ(ierr);
        ierr = DMDAVecRestoreArrayRead(u->fda, u->lEta,  &eta);  CHKERRQ(ierr);
        ierr = DMDAVecRestoreArrayRead(u->fda, u->lZet,  &zet);  CHKERRQ(ierr);
 
        /* Now that this block's arrays are restored, surface any cell-loop error. */
        PetscCheck(cell_err == 0, PETSC_COMM_SELF, cell_err, "%s at (%" PetscInt_FMT
                   ",%" PetscInt_FMT ",%" PetscInt_FMT ")", be_what ? be_what : "error", be_i, be_j, be_k);
    }
 
    /* ---- global reductions (no ghost/halo exchange; scalar collectives only) ----
       5 collectives total: one 3-element MAX (B/C/D), one SUM (active-cell count),
       one MIN (portable owner selection), and two broadcasts (the controlling-cell
       real breakdown and integer indices). */
    PetscReal loc3[3] = { locmaxB, locmaxC, locmaxD }, glo3[3];
    ierr = MPI_Allreduce(loc3, glo3, 3, MPIU_REAL, MPI_MAX, PETSC_COMM_WORLD); CHKERRQ(ierr);
 
    PetscInt global_active = 0;
    ierr = MPI_Allreduce(&local_active, &global_active, 1, MPIU_INT, MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
 
    /* The selected candidate's global max equals glo3[candidate] (no extra reduction).
       Portable owner pick: lowest rank whose local sel_max equals the global max. */
    const PetscReal sel_global = glo3[(int)candidate];
    PetscMPIInt nproc, claim, owner;
    ierr = MPI_Comm_size(PETSC_COMM_WORLD, &nproc); CHKERRQ(ierr);
    claim = (sel_max == sel_global) ? rank : nproc;
    ierr = MPI_Allreduce(&claim, &owner, 1, MPI_INT, MPI_MIN, PETSC_COMM_WORLD); CHKERRQ(ierr);
    if (owner == nproc) owner = 0;   /* no active cell anywhere: deterministic fallback owner */
 
    PetscReal dbuf[2] = { sel_lc, sel_lv };
    PetscInt  ibuf[6] = { sel_ci, sel_cj, sel_ck, sel_blk, sel_class, sel_os };
    ierr = MPI_Bcast(dbuf, 2, MPIU_REAL, owner, PETSC_COMM_WORLD); CHKERRQ(ierr);
    ierr = MPI_Bcast(ibuf, 6, MPIU_INT,  owner, PETSC_COMM_WORLD); CHKERRQ(ierr);
 
    rep->lambda_t = lambda_t;
    rep->lambda_c = dbuf[0];
    rep->lambda_v = dbuf[1];
    rep->lambda   = sel_global;
    rep->lambda_B = glo3[0]; rep->lambda_C = glo3[1]; rep->lambda_D = glo3[2];
    rep->ci = ibuf[0]; rep->cj = ibuf[1]; rep->ck = ibuf[2];
    rep->cblock = ibuf[3]; rep->cclass = ibuf[4]; rep->one_sided = ibuf[5];
    rep->active_cells = global_active;
 
    /* Only a genuinely empty active set (all cells masked) may fall back to the temporal
       term alone. A non-finite or non-positive estimate with active cells is a hard error. */
    if (global_active == 0) {
        rep->lambda = lambda_t;
    } else {
        PetscCheck(PetscIsNormalReal(rep->lambda) && rep->lambda > 0.0, PETSC_COMM_WORLD, PETSC_ERR_FP,
                   "momentum stability estimate is non-finite or non-positive with %" PetscInt_FMT
                   " active cells", global_active);
    }
 
    /* Dominant limiter at the controlling cell. */
    if (rep->lambda_t >= rep->lambda_c && rep->lambda_t >= rep->lambda_v) rep->limiter = MOM_STAB_LIMITER_TIME;
    else if (rep->lambda_c >= rep->lambda_v)                              rep->limiter = MOM_STAB_LIMITER_CONVECTION;
    else                                                                  rep->limiter = MOM_STAB_LIMITER_VISCOSITY;
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

MomCellActiveRows()

PetscInt MomCellActiveRows	(	PetscReal ***	nvert,
		PetscInt	k,
		PetscInt	j,
		PetscInt	i,
		PetscInt	mx,
		PetscInt	my,
		PetscInt	mz,
		PetscBool	np_x1,
		PetscBool	np_y1,
		PetscBool	np_z1,
		PetscInt	twoD
	)

Active staggered-momentum row mask for a cell (exposed for unit testing).

Returns a 3-bit mask (xi=1, eta=2, zeta=4). A solid cell yields 0 (all inactive); a positive solid neighbour or a positive non-periodic physical face clears the corresponding normal row; TwoD (1/2/3) clears the homogeneous direction's row.

Parameters

nvert	Local nvert array (ghosted).
k	Cell k index.
j	Cell j index.
i	Cell i index.
mx	Global x dimension.
my	Global y dimension.
mz	Global z dimension.
np_x1	True if the positive-x face is non-periodic.
np_y1	True if the positive-y face is non-periodic.
np_z1	True if the positive-z face is non-periodic.
twoD	TwoD homogeneous-direction selector (0 none, 1 xi, 2 eta, 3 zeta).

Returns

3-bit active-row mask (0 when the location carries no active unknown).

Definition at line 368 of file momentumsolvers.c.

{
    if (MOM_SKIP_SOLID(nvert[k][j][i])) return 0;     /* solid cell: all rows inactive */
    PetscInt bits = 0x7;
    if (MOM_SKIP_SOLID(nvert[k][j][i+1])) bits &= ~0x1;  /* positive xi neighbour solid   */
    if (MOM_SKIP_SOLID(nvert[k][j+1][i])) bits &= ~0x2;  /* positive eta neighbour solid  */
    if (MOM_SKIP_SOLID(nvert[k+1][j][i])) bits &= ~0x4;  /* positive zeta neighbour solid */
    if (np_x1 && i == mx-2) bits &= ~0x1;                /* positive non-periodic xi face  */
    if (np_y1 && j == my-2) bits &= ~0x2;
    if (np_z1 && k == mz-2) bits &= ~0x4;
    if      (twoD == 1) bits &= ~0x1;                    /* TwoD homogeneous direction: xi  */
    else if (twoD == 2) bits &= ~0x2;                    /*                             eta  */
    else if (twoD == 3) bits &= ~0x4;                    /*                             zeta */
    return bits;
}

Here is the caller graph for this function: