Investigate 64³ resolution instability anomaly in Alfvénic cascade

[anjor]

Problem

The 64³ resolution exhibits anomalous instability requiring 10× stronger dissipation (η=20.0) compared to expected scaling from 32³ (η~1.5). This contrasts with 128³ which follows expected scaling (η=2.0).

Evidence from PR #81

Systematic parameter search showed 7 failed attempts before finding stable parameters:

Attempt	η	r	amp	Result
1	0.5	4	1.0	Blowup (~1.4×10⁶)
2	0.5	4	0.1	Blowup
3	5.0	4	0.1	Catastrophic (1.6×10²⁷)
4	1.5	2	0.8	NaN at t=2.1 τ_A
5	2.5	2	0.05	NaN at t=14.8 τ_A
6	8.0	2	0.05	NaN at t=20.4 τ_A
7	20.0	2	0.05	Stable to 40 τ_A
8	20.0	2	0.01	Stable to 100 τ_A ✅

Current workaround: η=20.0 with gentle forcing (amplitude=0.01) produces correct k⊥^(-5/3) spectra but is not a fundamental solution.

Possible Root Causes

1. Wavenumber Resonance

• Forcing at k=1,2 may resonate with dealiased modes at k~21 (2/3 cutoff of N=64)
• Specific mode interactions at this resolution?
• Check if 48³ or 80³ exhibit similar behavior

2. Critical Balance Violation

• τ_nl ~ τ_A condition may fail at specific k⊥ range for 64³
• Need to verify eddy turnover time vs Alfvén time across scales
• Compare (k⊥v⊥)^(-1) vs (k∥v_A)^(-1) at different resolutions

3. Dealiasing Issues

• 2/3 cutoff at k_max = (N-1)//3 may have resolution-dependent artifacts
• For N=64: k_max=21, for N=128: k_max=42
• Could be specific to how nonlinear terms alias at this boundary

4. CFL/Courant Number

• Max CFL = max(|u|)·dt/dx may exceed stability threshold
• Need to log max velocity and check if approaching CFL limit
• Timestep dt~0.005 may be inappropriate for 64³ specifically

5. Aliasing Error Accumulation

• Energy may leak from inertial range to Nyquist modes
• Accumulation over time causes instability
• Exponential growth near k_max that normalized dissipation can't control

Recommended Investigation

Diagnostic Logging (High Priority)

Add to alfvenic_cascade_benchmark.py:

# After each timestep
max_velocity = jnp.max(jnp.abs(irfftn(state.z_plus + state.z_minus)))
cfl = max_velocity * dt / dx
max_nonlinear = jnp.max(jnp.abs(poisson_bracket(...)))

print(f"Max CFL: {cfl:.4f}, Max nonlinear: {max_nonlinear:.2e}")

Resolution Sweep (Medium Priority)

Run at intermediate resolutions:

• 48³: Between 32³ and 64³
• 80³: Between 64³ and 128³
• Check if anomaly is specific to N=64 or gradual transition

Wavenumber Analysis (Medium Priority)

# After forcing
forcing_spectrum = jnp.abs(state.z_plus)**2
plt.semilogy(k_perp, forcing_spectrum, label='After forcing')
# Check for unexpected energy at high k

Critical Balance Check (Low Priority)

# Compute eddy turnover time
tau_nl = 1 / (k_perp * v_perp)
tau_A = 1 / (k_parallel * v_A)
print(f"τ_nl/τ_A ratio: {tau_nl/tau_A}")  # Should be ~1

Success Criteria

• Identify root cause of 64³ instability
• Either fix instability OR document as intrinsic limitation
• Determine if intermediate resolutions (48³, 80³) have similar issues
• If fixable, reduce η from 20.0 to expected ~1.5

References

• PR Implement normalized hyper-dissipation for Alfvénic cascade benchmark #81: Original discovery and parameter search
• examples/alfvenic_cascade_benchmark.py:116-125: Current workaround with inline comments
• Original GANDALF thesis Section 2.6.3: Alfvénic cascade parameters

Related Code

• examples/alfvenic_cascade_benchmark.py:115-125 (64³ parameters)
• src/krmhd/timestepping.py:444-456 (normalized hyper-dissipation)
• src/krmhd/physics.py:z_plus_rhs() (nonlinear terms)

[anjor]

New Insight: Forcing Amplitude Also Resolution-Dependent

Code review revealed another resolution-dependent anomaly that may be related:

Forcing Amplitude Scaling

if args.resolution == 32:
    force_amplitude = 0.05   # Baseline
elif args.resolution == 64:
    force_amplitude = 0.01   # 5× weaker than 32³ AND 128³
else:  # 128³
    force_amplitude = 0.05   # Same as 32³

Key observation: 64³ requires both:

1. 10× stronger dissipation (η=20 vs η~2)
2. 5× weaker forcing (0.01 vs 0.05)

This suggests the instability may be forcing-driven rather than dissipation-limited.

Hypothesis: Forcing-Dissipation Resonance

The anomaly could stem from specific wavenumber interactions:

• Forcing at k=1,2 (large scales)
• Cascade to k~10-15 (intermediate scales at 64³)
• Specific resonance with 2/3 dealiasing boundary at k=21

At 64³ specifically, this intermediate range may have insufficient damping relative to energy injection rate.

Recommended Tests

1. Decouple forcing and dissipation:

   # Test: 64³ with force_amplitude=0.05 (like 32³ and 128³)
   #       but keep η=2.0 (like 128³)
   # If unstable → dissipation issue
   # If stable → forcing issue

2. Test forcing-free decay:

   # Run decaying turbulence at 64³ with η=2.0, no forcing
   # If stable → proves forcing interaction is the culprit

3. Check energy injection at specific k:

   # After forcing, compute spectrum at k=10-15
   # Compare injection rate vs dissipation rate at these modes

4. Test different forcing modes:

   # Force only mode 1 (no mode 2)
   # Force modes 3-4 instead of 1-2
   # Check if instability persists

Updated Root Cause Candidates

1. Forcing-cascade resonance (NEW - HIGH PRIORITY)

   • Energy injection at k=1,2 cascades to k~10-15
   • At 64³, this range has specific spectral spacing
   • Insufficient dissipation relative to local cascade rate

2. Wavenumber-dependent CFL violation

   • Nonlinear CFL condition at k~10-15 exceeds linear estimate
   • Requires adaptive timestep or stricter safety factor

3. Critical balance breakdown

   • τ_nl ≈ τ_A assumption violated at 64³ intermediate scales
   • Need to verify eddy turnover time at k~10-15

4. Dealiasing ineffectiveness

   • Energy leaks past 2/3 cutoff at specific resolution
   • Accumulates in dealiased region

5. Forcing amplitude units issue

   • force_amplitude may have implicit resolution dependence
   • Check if 0.01 vs 0.05 is physically meaningful or empirical

The forcing amplitude discrepancy is a major clue—it suggests the problem isn't purely numerical but involves physics-numerical coupling.