Optimize Elsasser RHS Laplacian computations

[anjor]

Context

Currently, z_plus_rhs() and z_minus_rhs() compute the Laplacian twice:

1. Perpendicular Laplacian: laplacian(z, kx, ky, kz=None) for Poisson bracket
2. Full 3D Laplacian: laplacian(z, kx, ky, kz) for dissipation

Potential Optimization

The perpendicular Laplacian is a subset of the full Laplacian:

∇²⊥z = -(kx² + ky²)|z|²
∇²z = -(kx² + ky² + kz²)|z|²

We could compute the full Laplacian once and extract the perpendicular part.

Current Performance Impact

• Low priority: Laplacian computation is cheap (just array multiplication)
• Current approach prioritizes code clarity
• No bottleneck identified in profiling yet

When to Optimize

Only if profiling shows Laplacian computation is a significant fraction (>5%) of runtime.

Implementation Notes

If optimizing:

1. Add separate_laplacians: bool = True parameter
2. When False, compute full Laplacian and derive perpendicular part
3. Benchmark to verify actual speedup
4. Keep current implementation as fallback for clarity

Related

• PR feat: Implement Alfvén dynamics with Elsasser variables (Issue #6) #40 (original implementation)
• Code review comment Step 2: Core Spectral Infrastructure #2 about performance optimization

Labels

enhancement, performance, optimization

[anjor]

Status Check - Status unchanged (low priority optimization)

Current State

• No optimization implemented
• z_plus_rhs() and z_minus_rhs() still compute two Laplacians separately
• Original issue correctly identified this as "low priority"

Context Update

The implementation is working well in production:

• ✅ examples/decaying_turbulence.py runs efficiently
• ✅ examples/orszag_tang.py performs well
• ✅ No performance bottlenecks reported

Profiling Needed

Per original issue:

Only if profiling shows Laplacian computation is a significant fraction (>5%) of runtime.

We haven't profiled yet (Issue #35 would provide this data).

Current Usage Pattern

# In physics.py z_plus_rhs():
laplacian_perp = laplacian(z_plus, kx, ky, kz=None)  # For Poisson bracket
laplacian_full = laplacian(z_plus, kx, ky, kz)       # For dissipation

# Cost:
# laplacian_perp: -(kx² + ky²) * z_plus  → 3 array ops
# laplacian_full: -(kx² + ky² + kz²) * z_plus  → 4 array ops

Savings from optimization: ~1 array operation per RHS call. This is tiny compared to FFT cost (6 FFTs per Poisson bracket).

Recommendation

Original issue assessment remains correct:

Current approach prioritizes code clarity

Benefits of optimization: ~1-2% speedup (estimate)
Cost of optimization: Slightly more complex code
Verdict: Not worth it unless profiling shows it matters

Decision Tree

1. Implement Issue Add performance benchmark for Poisson bracket #35 (Performance benchmarks) ✓ First priority
2. Profile typical simulation:
   • If Laplacian is <5% of runtime → close this issue
   • If Laplacian is >5% of runtime → implement optimization
3. Currently expecting <1% (FFTs dominate)

Related Issues

• Issue Add performance benchmark for Poisson bracket #35 - Need profiling data first
• Issue Consider caching strategy for repeated Poisson bracket calls #37 - Similar "investigate before optimizing" issue
• Issue Optimize diagnostics spectrum calculation performance #52 - Found actual bottleneck (Python loops, not Laplacian)

Priority: Very low - Wait for profiling. Likely will close as "not worth optimizing".

Keeping open pending profiling, but expect to close. 🔬