Numerical Integration

Numerical Integration

Goal

Provide a reliable workflow to select integrators, set tolerances, and manage adaptive time stepping for time-dependent simulations.

Requirements

• Python 3.8+
• NumPy (for some scripts)
• No heavy dependencies for core functionality

Inputs to Gather

Decision Guidance

Choosing an Integrator

Is the problem stiff?
├── YES → Is Jacobian available?
│   ├── YES → Use Rosenbrock or BDF
│   └── NO → Use BDF with numerical Jacobian
└── NO → Is high accuracy needed?
    ├── YES → Use RK45 or DOP853
    └── NO → Use RK4 or Adams-Bashforth

Stiff vs Non-Stiff Detection

Script Outputs (JSON Fields)

Workflow

1. Classify stiffness - Check eigenvalue spread or use stiffness_detector 2. Choose tolerances - See references/tolerance_guidelines.md 3. Select integrator - Run scripts/integrator_selector.py 4. Compute error norms - Use scripts/error_norm.py for step acceptance 5. Adapt step size - Use scripts/adaptive_step_controller.py 6. Plan IMEX/splitting - If mixed stiff/nonstiff, use scripts/imex_split_planner.py 7. Validate convergence - Repeat with tighter tolerances

Conversational Workflow Example

User: I'm solving the Allen-Cahn equation with a stiff double-well potential. What integrator should I use?

Agent workflow: 1. Check integrator options:

   python3 scripts/integrator_selector.py --stiff --jacobian-available --accuracy high --json

2. Plan the IMEX splitting (diffusion implicit, reaction explicit):

   python3 scripts/imex_split_planner.py --stiff-terms diffusion --nonstiff-terms reaction --coupling weak --json

3. Recommend: Use IMEX-BDF2 with diffusion term implicit, double-well reaction explicit.

Pre-Integration Checklist

• [ ] Identify stiffness and dominant time scales
• [ ] Set rtol/atol consistent with physics and units
• [ ] Confirm integrator compatibility with stiffness
• [ ] Use error norm to accept/reject steps
• [ ] Verify convergence with tighter tolerance run

CLI Examples

# Select integrator for stiff problem with Jacobian
python3 scripts/integrator_selector.py --stiff --jacobian-available --accuracy high --json

# Compute scaled error norm
python3 scripts/error_norm.py --error 0.01,0.02 --solution 1.0,2.0 --rtol 1e-3 --atol 1e-6 --json

# Adaptive step control with PI controller
python3 scripts/adaptive_step_controller.py --dt 1e-2 --error-norm 0.8 --order 4 --controller pi --json

# Plan IMEX splitting
python3 scripts/imex_split_planner.py --stiff-terms diffusion,elastic --nonstiff-terms reaction --coupling strong --json

# Estimate splitting error
python3 scripts/splitting_error_estimator.py --dt 1e-4 --scheme strang --commutator-norm 50 --target-error 1e-6 --json

Error Handling

Interpretation Guidance

Error Norm Values

Controller Selection

IMEX Strategy

Limitations

• No automatic stiffness detection: Use stiffness_detector from numerical-stability
• Splitting assumes separability: Terms must be cleanly separable
• Jacobian requirement: Some methods need analytical or numerical Jacobian

References

• references/method_catalog.md - Integrator options and properties
• references/tolerance_guidelines.md - Choosing rtol/atol
• references/error_control.md - Error norm and adaptation formulas
• references/imex_guidelines.md - Stiff/non-stiff splitting
• references/splitting_catalog.md - Operator splitting patterns
• references/multiphase_field_patterns.md - Phase-field specific splits

Version History

• v1.1.0 (2024-12-24): Enhanced documentation, decision guidance, examples
• v1.0.0: Initial release with 5 integration scripts

Error Control

Comprehensive guide for adaptive step size control in time integration.

Error Norm Computation

Scaled Error Norm

The standard approach scales the error by tolerance:

scale_i = atol_i + rtol × max(|y_n,i|, |y_{n+1},i|)
err_scaled_i = error_i / scale_i

Norm Types

RMS Norm (Recommended Default):

error_norm = sqrt( (1/n) × Σᵢ (error_i / scale_i)² )

Max Norm (Conservative):

error_norm = max_i |error_i / scale_i|

When to Use Each Norm

Step Acceptance Criterion

IF error_norm ≤ 1.0:
    ACCEPT step
    Possibly increase dt
ELSE:
    REJECT step
    Decrease dt
    Retry with smaller step

Near-Acceptance Zone

Step Size Controllers

I-Controller (Integral/Proportional)

Simplest controller based on current error only:

dt_new = dt × safety × error_norm^(-1/(p+1))

Where:

• safety = 0.8 - 0.9 (prevents oscillation)
• p = method order

Characteristics:

• Simple implementation
• May oscillate if error varies rapidly
• Suitable for smooth problems

PI-Controller (Proportional-Integral)

Uses current and previous error for smoother adaptation:

dt_new = dt × safety × error_norm^(-α) × error_prev^(β)

Standard coefficients:

α = 0.7 / (p + 1)
β = 0.4 / (p + 1)

Alternative (Gustafsson):

α = 1 / (p + 1)
β = 1 / (p + 1) × (dt_prev / dt)

Characteristics:

• Smoother step size evolution
• Reduces oscillations
• Industry standard choice

PID-Controller

Adds derivative term for aggressive adaptation:

dt_new = dt × safety × error_norm^(-α) × error_prev^(β) × error_prev2^(γ)

Typical coefficients:

α = 0.49 / (p + 1)
β = 0.34 / (p + 1)
γ = 0.10 / (p + 1)

Characteristics:

• Most aggressive adaptation
• Best for rapidly changing dynamics
• Risk of oscillation if poorly tuned

Controller Comparison

Safety Factors and Limits

Safety Factor

safety = 0.8 - 0.9 (typical)

Purpose:

• Prevents accepting steps at exact tolerance boundary
• Accounts for error estimate uncertainty
• Reduces rejected steps due to fluctuations

Tuning:

• Conservative (stiff/sensitive): safety = 0.8
• Aggressive (smooth): safety = 0.95

Step Size Limits

dt_min = ε_machine × |t| × 100  (prevent stall)
dt_max = T_end / 10  (or physics-based limit)

Max step increase per step:

dt_new ≤ factor_max × dt  (typically factor_max = 2-5)

Max step decrease per step:

dt_new ≥ factor_min × dt  (typically factor_min = 0.1-0.2)

Limiting Formulas

# Apply all limits
dt_new = dt × safety × error_norm^(-1/(p+1))
dt_new = max(dt_new, dt × factor_min)
dt_new = min(dt_new, dt × factor_max)
dt_new = max(dt_new, dt_min)
dt_new = min(dt_new, dt_max)

Error Estimation Methods

Embedded Methods

Use two methods of different orders in one step:

y_high = O(h^{p+1}) method result
y_low  = O(h^p) embedded result
error  = y_high - y_low

Advantage: Error estimate is "free" from the integration.

Richardson Extrapolation

Compare full step with two half-steps:

y_full = one step of size h
y_half = two steps of size h/2
error = (y_half - y_full) / (2^p - 1)

Advantage: Works with any method. Disadvantage: Doubles function evaluations.

Milne Device (Multi-step)

Compare predictor and corrector:

y_predict = Adams-Bashforth (explicit)
y_correct = Adams-Moulton (implicit)
error = (y_correct - y_predict) × C

Where C depends on method orders.

Handling Rejected Steps

Standard Approach

1. Compute y_new and error estimate
2. If error_norm > 1:
   a. Reduce dt using controller
   b. Discard y_new
   c. Retry from y_current with smaller dt
3. If many consecutive rejections:
   a. Check for stiffness
   b. Consider method change

Rejection Statistics

After Rejection Strategy

Conservative:

dt_new = dt × 0.5  (halve step)

Controller-based:

dt_new = dt × safety × error_norm^(-1/(p+1))

Minimum reduction:

dt_new = max(dt × 0.1, dt_min)  (never reduce by more than 10×)

Special Situations

Output Points

When output is required at specific times:

1. Integrate to t_out with adaptive stepping
2. Last step may need adjustment to hit t_out exactly
3. Options:
   a. Interpolate from nearby steps (dense output)
   b. Force exact step to t_out
   c. Accept small overstep and interpolate back

Discontinuities

When a discontinuity is detected (event, switch):

1. Locate discontinuity time t_disc (bisection/Newton)
2. Integrate exactly to t_disc
3. Apply discontinuity (jump conditions)
4. Restart integrator with small initial step
5. Rebuild multi-step history if needed

Stiffness Detection

Monitor for implicit stiffness transition:

IF step rejections increase rapidly:
   Estimate stiffness ratio
   IF ratio > threshold (e.g., 1000):
       Switch to implicit method

Implementation Example

def adaptive_step(y, t, dt, f, tol, order, controller='pi'):
    """One adaptive step with PI control."""
    safety = 0.9
    factor_min, factor_max = 0.2, 5.0

    # Error history for PI controller
    global error_prev

    while True:
        # Take step and estimate error
        y_new, error = step_with_error(y, t, dt, f)

        # Compute scaled error norm
        scale = tol['atol'] + tol['rtol'] * np.maximum(np.abs(y), np.abs(y_new))
        error_norm = np.sqrt(np.mean((error / scale)**2))

        if error_norm <= 1.0:
            # Accept step
            if controller == 'pi':
                factor = safety * error_norm**(-0.7/(order+1)) * error_prev**(0.4/(order+1))
            else:  # 'i' controller
                factor = safety * error_norm**(-1/(order+1))

            factor = np.clip(factor, factor_min, factor_max)
            dt_new = dt * factor
            error_prev = error_norm

            return y_new, t + dt, dt_new, True
        else:
            # Reject step
            factor = safety * error_norm**(-1/(order+1))
            factor = max(factor, factor_min)
            dt = dt * factor
            # Loop continues with smaller dt

Tuning Guidelines

Start Conservative

1. Begin with safety = 0.8, PI controller 2. Use default α, β coefficients 3. Set factor_max = 2.0

Tune for Efficiency

1. If few rejections, increase safety to 0.9 2. If smooth evolution, increase factor_max to 5.0 3. If oscillating dt, reduce α or switch to PI

Monitor Health

Track and report:

• Total step count
• Rejected step count
• Minimum dt reached
• Maximum error norm encountered

IMEX Guidelines

Comprehensive guide for Implicit-Explicit time integration methods.

When to Use IMEX

Problem Structure

IMEX methods are designed for problems of the form:

du/dt = f_stiff(u) + f_nonstiff(u)

Where:

• f_stiff: Requires implicit treatment for stability
• f_nonstiff: Can be handled explicitly (cheaper)

Canonical Examples

Decision Criteria

Use IMEX when: 1. Clear separation exists between stiff and non-stiff terms 2. Stiff term is linear or easy to solve implicitly 3. Non-stiff term dominates cost if treated implicitly 4. Explicit dt limit would be too restrictive for stiff part

Avoid IMEX when: 1. All terms have similar stiffness 2. Terms are strongly coupled (nonlinear interaction) 3. Stiff term is highly nonlinear 4. Splitting error is unacceptable

IMEX Scheme Classes

IMEX-Runge-Kutta

Combines explicit and implicit RK methods:

Stage i:
  U_i = u_n + dt × Σⱼ aᵢⱼᴱ f_explicit(Uⱼ) + dt × Σⱼ aᵢⱼᴵ f_implicit(Uⱼ)

Update:
  u_{n+1} = u_n + dt × Σᵢ bᵢᴱ f_explicit(Uᵢ) + dt × Σᵢ bᵢᴵ f_implicit(Uᵢ)

Common Schemes:

IMEX-Multistep (SBDF)

Semi-implicit BDF combines BDF for stiff with Adams-Bashforth for non-stiff:

SBDF1:

(u_{n+1} - u_n) / dt = f_implicit(u_{n+1}) + f_explicit(u_n)

SBDF2:

(3u_{n+1} - 4u_n + u_{n-1}) / (2dt) = f_implicit(u_{n+1}) + 2f_explicit(u_n) - f_explicit(u_{n-1})

IMEX-Peer Methods

Multi-stage multi-step methods with good parallel properties:

• Each stage can be computed in parallel
• Combine strengths of RK and multistep
• Emerging class with ongoing research

Splitting Strategies

Term Assignment

Semi-implicit Linearization

When nonlinear terms must be implicit, linearize:

Original: du/dt = -∇·(D(u)∇u)

Linearized: du/dt = -∇·(D(u_n)∇u_{n+1})
             ↑            ↑         ↑
          implicit    lagged    new value

Lagging strategies:

• Lag coefficient by one step (first-order in time for coefficient)
• Extrapolate coefficient (higher order)
• Iterate until converged (fully implicit)

Stability Analysis

Combined Stability Region

The effective stability region is the intersection of: 1. Implicit method's region for stiff term 2. Explicit method's region for non-stiff term

For stability: both λ_stiff × dt and λ_nonstiff × dt must be in respective regions

Stability Plots

For IMEX methods, plot stability in (λ_E × dt, λ_I × dt) plane:

• Horizontal axis: explicit eigenvalue × dt
• Vertical axis: implicit eigenvalue × dt
• Shaded region: stable combinations

CFL for Explicit Part

Even with IMEX, the explicit part imposes a CFL constraint:

dt ≤ C × dx / |v_max|  (for advection)

The implicit treatment of diffusion removes:

dt ≤ dx² / (2D)  (this constraint is gone!)

Accuracy Considerations

Order Conditions

For IMEX-RK, both tableaux must satisfy order conditions:

• Individual ERK and DIRK must be consistent
• Coupling conditions between tableaux
• Stiff accuracy (optional but helpful)

Splitting Error

IMEX introduces splitting error beyond method truncation error:

Total error = O(dt^p) + splitting error

Splitting error depends on:

• Commutator [f_stiff, f_nonstiff]
• Method coupling (how stages interact)

For additive IMEX-RK: splitting error = O(dt^p) when properly designed.

Order Reduction

Near stability boundaries:

• Classical order reduction for stiff ODEs applies
• Some IMEX schemes are "stiff accurate" to mitigate this

Implementation Guidelines

Linear Solver Requirements

The implicit part requires solving:

(I - dt × γ × J_implicit) × δu = RHS

Where J_implicit is the Jacobian of f_implicit.

Requirements: 1. Jacobian or Jacobian-vector product 2. Linear solver (direct or iterative) 3. Preconditioner for large systems

Jacobian Reuse

For IMEX-RK with SDIRK (same diagonal):

• Same Jacobian matrix at all implicit stages
• Factor once per step, reuse for all stages
• Significant cost savings

For varying Jacobian:

• Recompute each stage OR
• Lag Jacobian across stages (lose some order)

Iteration Strategy

Fixed number of iterations:

• Predictable cost
• May not converge fully
• OK if splitting dominates error

Converge to tolerance:

• Fully implicit behavior
• Variable cost per step
• Use when high accuracy needed

Memory Management

Practical Examples

Advection-Diffusion

∂u/∂t + v · ∇u = D ∇²u

f_explicit = -v · ∇u  (advection)
f_implicit = D ∇²u    (diffusion)

dt constraint: CFL for advection only
dt ≤ dx / |v_max|  (instead of min with dx²/D)

Allen-Cahn Equation

∂φ/∂t = M(ε²∇²φ - φ(φ²-1))

Option 1 (linear implicit):
  f_implicit = M ε² ∇²φ
  f_explicit = -M φ(φ²-1)

Option 2 (stabilized explicit):
  f_implicit = M ε² ∇²φ - M s φ  (add stabilizing term)
  f_explicit = -M φ(φ²-1) + M s φ  (subtract same term)

  Choose s to make explicit part bounded.

Cahn-Hilliard Equation

∂φ/∂t = M ∇²(f'(φ) - ε² ∇²φ)

Convex-concave splitting:
  f_implicit = M ∇²(c₁ φ - ε² ∇²φ)  (convex part)
  f_explicit = M ∇²(f'(φ) - c₁ φ)   (concave part)

  c₁ chosen so that f''(φ) - c₁ ≤ 0

Coupling Strength

Weak Coupling

Terms interact slowly:

Example: Heat conduction + slow reaction
Solution: Simple splitting works well

Moderate Coupling

Terms have some interaction:

Example: Fast diffusion + moderate reaction
Solution: IMEX-RK with good coupling coefficients

Strong Coupling

Terms are tightly coupled:

Example: Phase-field with strong anisotropy
Solution: Fully coupled IMEX-RK or iteration

Verification

Order Verification

1. Run with dt, dt/2, dt/4 2. Compute error vs exact (if known) or Richardson 3. Confirm error ∝ dt^p

Splitting Error Check

1. Run pure implicit (full problem) 2. Run IMEX with same dt 3. Difference reveals splitting error

Stability Boundary

1. Gradually increase dt 2. Find maximum stable dt 3. Compare with theoretical prediction

Method Catalog

Comprehensive reference for time integration methods in ODE/PDE simulations.

Explicit Methods (Non-Stiff)

Runge-Kutta Family

RK45 (Dormand-Prince)

• Default choice for non-stiff problems
• Embedded 4th-order method for error estimation
• FSAL (First Same As Last) optimization
• Recommended tolerances: rtol=1e-3, atol=1e-6

DOP853

• 8th-order with 5th and 3rd-order error estimates
• Excellent for high-precision requirements
• More expensive per step but fewer steps needed
• Use when rtol < 1e-6 is required

Adams-Bashforth (Multi-step)

Advantages:

• Only one function evaluation per step
• Efficient for smooth, non-stiff problems

Disadvantages:

• Requires startup procedure
• Less robust for discontinuous forcing
• Needs variable-step modification for adaptivity

Implicit Methods (Stiff)

BDF (Backward Differentiation Formulas)

When to use:

• Large eigenvalue spread (> 100)
• Chemistry with fast/slow reactions
• Diffusion-dominated problems

Considerations:

• Requires nonlinear solver (Newton)
• Jacobian needed (analytical or numerical)
• Order reduction near stability boundary

Radau IIA

Properties:

• L-stable (strong damping of stiff modes)
• Excellent for very stiff problems
• Superconvergent at endpoints

Use when:

• BDF order reduction is problematic
• DAE (differential-algebraic) systems
• Very stiff chemistry

Rosenbrock Methods

Characteristics:

• Linearly implicit (one Jacobian factorization per step)
• No nonlinear iteration needed
• Excellent for moderate stiffness

Advantages over BDF:

• No iteration convergence issues
• Fixed number of Jacobian evaluations
• Better for time-varying Jacobians

Structure-Preserving Methods

Symplectic Integrators

For Hamiltonian systems: dp/dt = -∂H/∂q, dq/dt = ∂H/∂p

Use for:

• Long-time molecular dynamics
• Orbital mechanics
• Oscillatory systems with energy conservation

Properties:

• Exactly conserve symplectic structure
• Near-conservation of Hamiltonian for exponentially long times
• Time-reversible (symmetric methods)

Geometric Integrators

IMEX Methods

For problems with mixed stiff/non-stiff terms: du/dt = f_stiff(u) + f_nonstiff(u)

Common IMEX-RK Schemes

IMEX-BDF

Selection Guide

Quick Decision Table

Cost Comparison

Stability Regions

Explicit Methods

• RK4: Extends to Re(λh) ≈ -2.8 on real axis
• RK45: Similar to RK4
• Adams-Bashforth: Smaller regions, decreasing with order

Implicit Methods

• BDF1-2: A-stable (entire left half-plane)
• BDF3-5: A(α)-stable (wedge-shaped regions)
• Radau: L-stable (A-stable + stiff decay)

Implementation Notes

Jacobian Handling

1. Analytical: Most accurate, requires code derivation 2. Automatic differentiation: Accurate, some overhead 3. Numerical (finite difference): Simple, may be inaccurate for stiff 4. Jacobian-free (GMRES): For very large systems

Step Size Limits

• Minimum dt: Floating-point precision, typically 1e-15 * t_current
• Maximum dt: Physical time scale or output frequency
• Safety factor: Typically 0.8-0.9 for adaptive methods

Multiphase-Field Splitting Patterns

Specialized integration patterns for phase-field and multi-order-parameter models.

Phase-Field Model Classes

Allen-Cahn (Non-Conserved)

Model A kinetics for non-conserved order parameters:

∂φ/∂t = -M × δF/δφ = -M × (f'(φ) - ε²∇²φ)

Where:

• φ: order parameter (e.g., phase indicator)
• M: mobility
• f(φ): bulk free energy density
• ε: interface width parameter

Cahn-Hilliard (Conserved)

Model B kinetics for conserved order parameters:

∂c/∂t = ∇·(M∇μ)
μ = δF/δc = f'(c) - ε²∇²c

Where:

• c: conserved field (e.g., concentration)
• μ: chemical potential
• M: mobility (may be M(c))

Multi-Order-Parameter

Multiple coupled order parameters:

∂φᵢ/∂t = -Mᵢ × δF/δφᵢ + Σⱼ coupling_ij(φⱼ)

Examples:

• Grain growth: φᵢ for each grain
• Eutectic: φ_α, φ_β, c
• Polycrystal: φᵢ + θᵢ (orientation)

Allen-Cahn Splitting

Standard Double-Well

f(φ) = (1/4)(φ² - 1)² = (1/4)φ⁴ - (1/2)φ² + 1/4

IMEX Split (Linear Implicit):

Implicit: ε²∇²φ           (stiffest term)
Explicit: -φ(φ² - 1)      (bulk driving force)

Stability: Explicit part bounded for |φ| ≤ 1 + O(ε)

Stabilized Splitting

Add/subtract linear stabilizing term:

∂φ/∂t = M[ε²∇²φ - sφ] + M[-φ³ + φ + sφ]
         └─implicit─┘   └───explicit───┘

Choose s > max|f''(φ)| = 3 for double-well.

Effect: Explicit term becomes contractive, unconditionally stable.

Convex-Concave Splitting (Eyre)

Split f(φ) into convex and concave parts:

f(φ) = f_convex(φ) + f_concave(φ)

f_convex = (c/2)φ²  where c ≥ max f''
f_concave = f(φ) - (c/2)φ²

Scheme:

(φ_{n+1} - φ_n)/dt = M[ε²∇²φ_{n+1} - f'_convex(φ_{n+1}) - f'_concave(φ_n)]

Properties:

• Unconditionally energy stable
• Unique solvability
• First-order in time

Second-Order Convex Splitting

Using Crank-Nicolson for spatial:

(φ_{n+1} - φ_n)/dt = M[ε²∇²(φ_{n+1}+φ_n)/2
                      - f'_convex(φ_{n+1})
                      - f'_concave(φ_n)]

Or BDF2:

(3φ_{n+1} - 4φ_n + φ_{n-1})/(2dt) = M[ε²∇²φ_{n+1}
                                     - f'_convex(φ_{n+1})
                                     - f'_concave(2φ_n - φ_{n-1})]

Cahn-Hilliard Splitting

Fourth-Order Challenge

Cahn-Hilliard involves ∇⁴:

∂c/∂t = M∇²μ = M∇²(f'(c) - ε²∇²c)

Direct discretization:

• Explicit: dt ~ dx⁴ (very restrictive!)
• Implicit: fourth-order operator

Mixed Formulation

Introduce chemical potential as separate variable:

∂c/∂t = M∇²μ
μ = f'(c) - ε²∇²c

Two second-order equations instead of one fourth-order.

IMEX for Mixed Form

Implicit: -ε²∇²c (in μ equation)
          ∇²μ (in c equation)
Explicit: f'(c)

Linear system per step:

[I      ε²L  ] [c_{n+1}]   [c_n + dt×M×L×μ_{n+1}]
[-1   I     ] [μ_{n+1}] = [f'(c_n)              ]

Where L = discrete Laplacian.

Convex Splitting for Cahn-Hilliard

f(c) = f_convex(c) + f_concave(c)

(c_{n+1} - c_n)/dt = M∇²[f'_convex(c_{n+1}) + f'_concave(c_n) - ε²∇²c_{n+1}]

Properties:

• Energy decreasing: F[c_{n+1}] ≤ F[c_n]
• Mass conserving: ∫c dx = constant
• Unconditionally stable

Scalar Auxiliary Variable (SAV)

Introduce r(t) = √(F[c] + C) where C ensures positivity:

dc/dt = M∇²μ
μ = r(t)/√(F+C) × δF/δc - ε²∇²c
dr/dt = (1/2√(F+C)) ∫ (δF/δc) × (dc/dt) dx

Advantages:

• Linear implicit solve
• Energy stable
• Second-order possible

Coupled Multi-Phase Systems

Two-Phase with Concentration

φ: phase field, c: concentration

∂φ/∂t = -M_φ × δF/δφ
∂c/∂t = ∇·(M_c ∇(δF/δc))

Splitting strategy:

Step 1: Evolve φ with c fixed (Allen-Cahn)
Step 2: Evolve c with φ fixed (diffusion with phase-dependent mobility)

Coupling terms:

• Free energy: f(φ,c) couples fields
• Mobility: M_c(φ) may depend on phase

Multi-Grain Systems

N order parameters φ₁, ..., φ_N with constraint:

Σᵢ φᵢ² = 1  (at interfaces)

Sequential update:

For i = 1 to N:
    Update φᵢ with φⱼ (j≠i) fixed
    Enforce constraint

Projection method:

Update all φᵢ without constraint
Project: φᵢ → φᵢ / √(Σⱼ φⱼ²)

Phase-Field Crystal

∂ψ/∂t = ∇²[(r + (1+∇²)²)ψ + ψ³]

Splitting:

Implicit: (1+∇²)²ψ (sixth-order!)
Explicit: rψ + ψ³

Requires careful treatment due to high-order derivatives.

Anisotropy Handling

Weakly Anisotropic

Interface energy: γ(n) = γ₀(1 + ε₄cos(4θ))

Gradient energy:

W(∇φ) = (1/2)|∇φ|² × a(n)²

IMEX approach:

Implicit: isotropic part (1/2)|∇φ|²
Explicit: anisotropic correction

Strongly Anisotropic

When a(n)² can become negative (missing orientations):

Regularization:

W_reg = (1/2)|∇φ|² × a(n)² + (β/2)|∇∇φ|²

Corner regularization term requires implicit treatment.

Faceted Interfaces

For crystalline anisotropy (Wulff shapes with corners):

Use Willmore regularization or
Level-set with crystalline curvature

Specialized schemes needed (not standard IMEX).

Adaptive Mesh Refinement (AMR)

Where to Refine

Refine near interfaces where gradients are large:

Criterion: |∇φ| > threshold
        or |φ - 0.5| < δ

Time Stepping with AMR

Subcycling:

Fine levels: dt_fine = dt_coarse / ratio
Synchronize at coarse step

Splitting with AMR:

1. Coarse step for bulk
2. Fine steps for interface region
3. Synchronize solutions at boundaries

Practical Guidelines

Time Step Selection

With IMEX: can exceed explicit limits by 10-100×.

Convergence Checking

1. Energy decay: dF/dt ≤ 0 (thermodynamically consistent) 2. Mass conservation: ∫c dx = const (for Cahn-Hilliard) 3. Interface width: Check ε² matches expected profile 4. Steady state: φ converges to equilibrium

Common Pitfalls

Interface Resolution

Rule of thumb:

dx ≤ ε / 3  (at least 3 points across interface)
Better: dx ≤ ε / 5  (5-6 points)

For IMEX stability:

dt ≤ C × dx² / (M × ε²)  (explicit part)

With stabilization, C can be O(1) instead of O(0.1).

Verification Benchmarks

Allen-Cahn

1. Shrinking circle: R(t) = √(R₀² - 2Mt) 2. Traveling wave: Compare to analytical profile 3. Equilibrium: Flat interface, f'(φ) = 0

Cahn-Hilliard

1. Spinodal decomposition: Compare coarsening rate 2. Ostwald ripening: R ~ t^{1/3} for droplets 3. Equilibrium: μ = constant, phase fractions correct

Multi-Order-Parameter

1. Grain growth: Normal grain growth law 2. Triple junction: 120° angles for isotropic 3. Wetting: Contact angle vs surface energies

Operator Splitting Catalog

Comprehensive reference for operator splitting methods in PDEs and multiphysics.

Fundamental Concepts

Problem Setup

Given:

du/dt = A(u) + B(u)

Where A and B are operators (may be differential, algebraic, or mixed).

Goal: Solve A and B separately, combine solutions.

Why Split?

First-Order Methods

Lie-Trotter Splitting

Step 1: Solve du/dt = A(u) for time dt → u*
Step 2: Solve du/dt = B(u) for time dt starting from u* → u_{n+1}

Properties:

• Order: 1 (error = O(dt))
• Simplest implementation
• Non-symmetric (order matters)

Error:

Error ≈ (dt/2) × [A, B] × u

Where [A, B] = AB - BA is the commutator.

Choosing Order (A then B vs B then A)

Second-Order Methods

Strang Splitting

Step 1: Solve du/dt = A(u) for time dt/2 → u*
Step 2: Solve du/dt = B(u) for time dt → u**
Step 3: Solve du/dt = A(u) for time dt/2 → u_{n+1}

Properties:

• Order: 2 (error = O(dt²))
• Symmetric (time-reversible)
• Standard choice for moderate accuracy

Error:

Error ≈ (dt²/24) × [[A, B], A + B] × u

Involves nested commutators.

Marchuk-Strang Alternating

Step n (even): A/2 → B → A/2
Step n (odd):  B/2 → A → B/2

Alternating between Strang orderings:

• Averages directional bias
• Can improve symmetry for anisotropic problems

Higher-Order Methods

Triple-Jump (Order 4)

γ = 1 / (2 - 2^{1/3}) ≈ 1.3512

Sub-step 1: Strang with dt × γ
Sub-step 2: Strang with dt × (1 - 2γ)  [negative!]
Sub-step 3: Strang with dt × γ

Note: Requires a negative time step!

• Problematic for irreversible operators (diffusion)
• Works for Hamiltonian systems

Yoshida (Order 4, Symmetric)

For symmetric operators only:

w₀ = -2^{1/3} / (2 - 2^{1/3})
w₁ = 1 / (2 - 2^{1/3})

Sequence: S(w₁ dt) → S(w₀ dt) → S(w₁ dt)

Where S is Strang splitting.

Forest-Ruth (Order 4)

For Hamiltonian splitting H = T(p) + V(q):

θ = 1 / (2 - 2^{1/3})
Sequence: V(θ/2) → T(θ) → V((1-θ)/2) → T(1-2θ) → V((1-θ)/2) → T(θ) → V(θ/2)

Directional Splitting (ADI)

2D Diffusion

∂u/∂t = D(∂²u/∂x² + ∂²u/∂y²)

Split by direction:

A: ∂u/∂t = D ∂²u/∂x²  (x-diffusion)
B: ∂u/∂t = D ∂²u/∂y²  (y-diffusion)

Douglas-Gunn (Unconditionally Stable)

Step 1: (I - dt/2 × D_xx) u* = (I + dt/2 × D_xx + dt × D_yy) u_n
Step 2: (I - dt/2 × D_yy) u_{n+1} = u* - dt/2 × D_yy u_n

Properties:

• Second-order accurate
• Unconditionally stable for diffusion
• Tridiagonal solves only

Peaceman-Rachford

Step 1: (I - dt/2 × D_xx) u* = (I + dt/2 × D_yy) u_n
Step 2: (I - dt/2 × D_yy) u_{n+1} = (I + dt/2 × D_xx) u*

Properties:

• Second-order
• Symmetric form
• Classic ADI method

3D Extension

For 3D: ∂u/∂t = D(∂²u/∂x² + ∂²u/∂y² + ∂²u/∂z²)

Douglas-Rachford-Gunn:

(I - dt/2 × D_xx) u* = (I + dt/2 × D_xx + dt × D_yy + dt × D_zz) u_n
(I - dt/2 × D_yy) u** = u* - dt/2 × D_yy u_n
(I - dt/2 × D_zz) u_{n+1} = u** - dt/2 × D_zz u_n

Physics-Based Splitting

Advection-Reaction

∂u/∂t + v·∇u = R(u)

A: ∂u/∂t + v·∇u = 0  (pure advection)
B: du/dt = R(u)       (pure reaction, ODE at each point)

Advantages:

• Use specialized advection scheme (upwind, WENO)
• Use stiff ODE solver for reaction
• Each solver optimized for its physics

Diffusion-Reaction

∂u/∂t = D∇²u + R(u)

A: ∂u/∂t = D∇²u    (diffusion)
B: du/dt = R(u)    (reaction)

Common Approach:

• Implicit for diffusion (larger dt)
• CVODE/LSODA for stiff reaction

Full Advection-Diffusion-Reaction

∂u/∂t + v·∇u = D∇²u + R(u)

Three-way split:
A: Advection
B: Diffusion
C: Reaction

Strang-like: A/2 → B/2 → C → B/2 → A/2

Splitting Error Analysis

Commutator Size

Splitting error depends on how much operators "fail to commute":

[A, B] = AB - BA

Error Estimation

A posteriori estimate:

Compare: Strang (dt) vs 2 × Strang (dt/2)
Error ≈ (result_dt - result_dt/2) / 3

Adaptive splitting:

• If error small: increase dt
• If error large: decrease dt or refine splitting

Error Accumulation

Over long times:

Total error ≈ (T_final / dt) × local_splitting_error
           = O(dt^{p-1}) for order-p method

Special Techniques

Balanced Splitting

For problems with conservation laws:

Ensure: ∫ u dx conserved by each sub-step

Modify operators to preserve conservation:

• Conservative discretization for each piece
• Flux-form splitting

Stabilized Splitting

Add and subtract stabilizing terms:

Original: du/dt = A(u)
Modified: du/dt = [A(u) + Su] - Su
          Split:   implicit      explicit

Choose S to improve stability of explicit part.

Iterative Splitting

Iterate between operators until converged:

While not converged:
    Solve A with current B(u)
    Solve B with updated A(u)

Converges to coupled solution (removes splitting error).

Implementation Patterns

Modular Code Structure

def splitting_step(u, dt, A_solver, B_solver, method='strang'):
    if method == 'lie':
        u = A_solver(u, dt)
        u = B_solver(u, dt)
    elif method == 'strang':
        u = A_solver(u, dt/2)
        u = B_solver(u, dt)
        u = A_solver(u, dt/2)
    return u

State Management

For multi-step splittings:

class SplitSolver:
    def __init__(self):
        self.u_history = []  # For multi-step methods

    def step(self, u, dt):
        # Store for multi-step
        self.u_history.append(u.copy())
        if len(self.u_history) > 3:
            self.u_history.pop(0)
        # Splitting step...

Parallel Splitting

When operators act on independent domains:

A: operates on spatial points 0:N/2
B: operates on spatial points N/2:N
→ Solve A and B in parallel!

Choosing a Method

Decision Table

Red Flags

Verification Tests

Manufactured Solutions

1. Choose exact u(x,t) 2. Compute source S = ∂u/∂t - A(u) - B(u) 3. Solve with splitting + source S 4. Compare to exact

Order Verification

1. Run with dt, dt/2, dt/4 2. Measure error vs reference 3. Compute order = log(e₁/e₂) / log(2) 4. Should match theoretical order

Conservation Test

1. Compute ∫ u dx at t=0 2. Run splitting simulation 3. Check ∫ u dx drift over time 4. Should be O(dt^p) or machine precision

Tolerance Guidelines

Comprehensive guide for setting integration tolerances in time-dependent simulations.

Fundamental Concepts

Relative vs Absolute Tolerance

Local error estimate ≤ atol + rtol × |y|

Error Interpretation

• rtol = 1e-3: ~3 significant digits accuracy
• rtol = 1e-6: ~6 significant digits accuracy
• atol = 1e-10: Values below 1e-10 treated as "zero"

Default Starting Points

By Application Type

By Problem Type

Multicomponent Systems

The Scaling Problem

When variables have vastly different magnitudes:

Temperature: O(1000 K)
Concentration: O(1e-6 mol/L)
Velocity: O(10 m/s)

A single atol cannot work for all!

Solution: Component-wise atol

# Python/SciPy example
atol = [1e-3, 1e-12, 1e-5]  # T, C, V
rtol = 1e-4  # Same for all

Determining Component atol

Nondimensionalization Alternative

Scale all variables to O(1):
T* = T / T_ref
C* = C / C_ref
Then use atol = 1e-6 for all

Advantages:

• Simpler tolerance specification
• Better numerical conditioning
• Easier debugging

Dimensional Analysis Approach

Physical Scale Identification

1. Identify characteristic scales

• Length: L (domain size or feature size)
• Time: τ (diffusion time, reaction time)
• Velocity: U (imposed or derived)
• Temperature: ΔT (driving temperature difference)

2. Compute derived scales

• Concentration change: ΔC = C_initial - C_eq
• Energy: ρ × c_p × ΔT × L³

3. Set atol

• atol_T = 1e-3 × ΔT
• atol_C = 1e-3 × ΔC
• atol_U = 1e-3 × U

Diagnosing Tolerance Issues

Symptoms and Solutions

Diagnostic Workflow

1. Run with default (rtol=1e-3, atol=1e-6)
2. Check:
   - Step rejection rate < 10%? If not, loosen
   - Conservation errors acceptable? If not, tighten
   - Physical variables positive? If not, reduce atol
3. Tighten by 10× and re-run
4. If results change significantly, tolerances were too loose

Conservation and Tolerance

Mass Conservation

For conserved quantities (mass, energy):

Relative conservation error ≈ rtol × (integration time / characteristic time)

Example:

• rtol = 1e-4
• Integration: 1000 time units
• Characteristic time: 1 unit
• Expected drift: ~1e-4 × 1000 = 0.1 (10% error!)

Solution: Use tighter rtol for long integrations or specialized conservative schemes.

Energy Drift in Hamiltonian Systems

For long-time simulations, prefer structure-preserving methods.

Stiff Problems

Special Considerations

1. Stiff components need tighter atol

• Fast modes decay to small values quickly
• Large atol masks stiff dynamics

2. Watch for order reduction

• Near stability boundary, effective order drops
• May need tighter rtol to compensate

3. Jacobian accuracy matters

• Numerical Jacobian errors ∝ sqrt(eps)
• May need tighter tolerances with numerical Jacobian

Recommended Settings for Stiff

Practical Examples

Phase-Field Simulation

# Two-phase Allen-Cahn
# phi: order parameter [-1, 1]
# T: temperature [300-2000 K]

rtol = 1e-4
atol_phi = 1e-6  # Order parameter precision
atol_T = 1e-1    # 0.1 K absolute accuracy

atol = [atol_phi, atol_T]

Reactive Flow

# Combustion with trace species
# Major species: O(1)
# Minor species: O(1e-6)
# Temperature: O(1000)

rtol = 1e-5  # Tighter for stiff chemistry
atol_major = 1e-8
atol_minor = 1e-12  # Detect small species
atol_T = 1e-2

atol = [..., atol_minor, ..., atol_T]

Solidification

# Phase-field solidification
# phi: solid fraction [0, 1]
# C: concentration [0.01 - 0.05]
# T: temperature [1700-1800 K]

rtol = 1e-4
atol = [1e-6,   # phi (interface precision)
        1e-6,   # C (microsegregation)
        0.01]   # T (0.01 K accuracy)

Tolerance Tightening Protocol

Convergence Study

1. Run with baseline tolerances (tol₀) 2. Run with tol₀/10 3. Compare key metrics (peak values, integrals, times) 4. If difference < acceptable threshold, tol₀ is adequate 5. If not, repeat with tighter tolerances

What to Compare

Common Mistakes

1. Same atol for all components

Problem: Variables with small physical values are ignored. Solution: Use component-wise atol scaled to each variable.

2. atol = 0

Problem: Division by zero when solution passes through zero. Solution: Always use small positive atol.

3. rtol too tight for noisy data

Problem: Excessive step rejection, slow progress. Solution: Match rtol to data precision.

4. Ignoring units

Problem: atol = 1e-6 means nothing without units. Solution: Express atol in simulation units: "1e-6 mol/L" or "0.1 K".