Skip to content

Solver Divergence & Stability

Physics-informed symbolic solvers can experience numerical instability if the initial conditions are unphysical, the mesh quality is poor, or the time-step is too large. This page provides guidelines for diagnosing and resolving solver divergence when using Shardian solvers.

1. Monitoring Residuals

A healthy CFD simulation should show monotonic residual convergence. Residuals for velocity components (\(U_x, U_y, U_z\)), pressure (\(p\)), and turbulence parameters (\(k, \omega\)) should fall below \(10^{-5}\) or \(10^{-6}\).

Time = 150
DILUPBiCGStab:  Solving for Ux, Initial residual = 4.21e-05, Final residual = 1.02e-06, No Iterations 2
DILUPBiCGStab:  Solving for Uy, Initial residual = 8.11e-05, Final residual = 1.45e-06, No Iterations 2
GAMG:  Solving for p, Initial residual = 3.82e-04, Final residual = 6.12e-06, No Iterations 5

If residuals begin to rise or fluctuate wildly: 1. Aero: Check the local maximum value of nut. If it exceeds your inflow viscosity by more than 5 orders of magnitude, the turbulence corrections are overproducing eddy viscosity. 2. Atmos: Check the WRF log files (rsl.error.0000). Look for extreme vertical velocity values: 

Max w-velocity =  35.2 m/s at grid cell (45, 12, 10)
If vertical velocity exceeds \(50\text{ m/s}\) in non-convective grids, the boundary layer coupling has diverged.

2. OpenFOAM (Shardian Aero) Stability Guidelines

If your OpenFOAM simulation aborts with floating-point errors (sigFpe), apply these adjustments:

A. Under-Relaxation Factors

Reduce relaxation factors in system/fvSolution during the initial phase of the run:

relaxationFactors
{
    fields
    {
        p               0.2; // Standard: 0.3
    }
    equations
    {
        U               0.5; // Standard: 0.7
        k               0.5; // Standard: 0.7
        omega           0.5; // Standard: 0.7
    }
}

Once the residuals stabilize (e.g., after 200 iterations), you can restore standard relaxation factors (0.3 for pressure, 0.7 for velocity and turbulence).

B. Bounding Turbulence Quantities

Ensure k and omega are bounded in system/fvSchemes using limited schemes:

divSchemes
{
    default         none;
    div(phi,U)      bounded Gauss linearUpwind grad(U);
    div(phi,k)      bounded Gauss upwind;
    div(phi,omega)  bounded Gauss upwind;
}

3. WRF (Shardian Atmos) Stability Guidelines

In mesoscale meteorology, grid-scale stability is dictated by the Courant-Friedrichs-Lewy (CFL) condition:

\[\text{CFL} = \frac{u \Delta t}{\Delta x} \le 1\]

If Shardian Atmos experiences CFL divergence (often labeled as cfl error or w-wind exceeding limits in WRF output):

A. Reduce the Time Step (\(\Delta t\))

A conservative rule of thumb for WRF is to set the time step in seconds to 4 times the grid spacing in kilometers:

\[\Delta t \approx 4 \cdot \Delta x_{\text{km}}\]

For example, if \(\Delta x = 3\text{ km}\), set \(\Delta t = 12\text{ s}\). If the boundary layer is highly unstable (high convective fluxes), reduce this to 3 times the grid spacing (\(\Delta t = 9\text{ s}\)).

B. Increase Vertical Grid Resolution

Avoid rapid changes in vertical grid spacing. Ensure a smooth stretching ratio (less than 1.2) when defining vertical \(\eta\) levels in namelist.input (eta_levels).