4.2 Matsuno Model: Deep Water Equations

"""
    Matsuno(params)

Matsuno formulation of the deep water wave equations.

The Matsuno formulation reformulates the problem to avoid computing
spatial derivatives of implicit quantities, leading to better numerical
stability and efficiency in spectral implementations.
"""
mutable struct Matsuno <: AbstractModel

    label    :: String
    datasize :: Int              # Number of components (2: η and φ)
    Γ        :: Array{Float64,1}  # |k| (dispersion relation modifier)
    Dx       :: Array{ComplexF64,1}  # i·k (spatial derivative in Fourier space)
    H        :: Array{ComplexF64,1}  # Hilbert transform (-i·sign(k))
    Π⅔       :: BitArray{1}       # Dealiasing filter (2/3 rule)
    ϵ        :: Float64           # Nonlinearity parameter (amplitude)
    
    # Work arrays for intermediate computations
    hnew     :: Vector{ComplexF64}
    unew     :: Vector{ComplexF64}
    I₀       :: Vector{ComplexF64}
    I₁       :: Vector{ComplexF64}
    I₂       :: Vector{ComplexF64}
    I₃       :: Vector{ComplexF64}

    Px       :: FFTW.FFTWPlan     # Pre-planned FFT for efficiency

    function Matsuno(param::NamedTuple)
        """
        Initialize Matsuno model solver
        
        Parameters:
        - param: named tuple with ϵ (nonlinearity), N (grid size)
        """
        
        label = "Matsuno"
        datasize = 2
        ϵ = param.ϵ
        
        mesh = Mesh(param)
        
        # Dispersion relation: |k|
        Γ = abs.(mesh.k)
        
        # Spatial derivative: ∂/∂x → i·k in Fourier space
        Dx = 1im * mesh.k
        
        # Hilbert transform: -i·sign(k)
        # Used in nonlinear terms
        H = -1im * sign.(mesh.k)
        
        # 2/3 dealiasing rule: keep only modes |k| < 2/3·k_max
        # This removes high-frequency aliasing artifacts from nonlinear products
        Π⅔ = Γ .< mesh.kmax * 2/3
        
        # Allocate work arrays
        hnew = zeros(ComplexF64, mesh.N)
        unew = zeros(ComplexF64, mesh.N)
        I₀ = zeros(ComplexF64, mesh.N)
        I₁ = zeros(ComplexF64, mesh.N)
        I₂ = zeros(ComplexF64, mesh.N)
        I₃ = zeros(ComplexF64, mesh.N)
        
        # Pre-plan FFT for efficiency
        # FFTW.MEASURE tries different algorithms to find the fastest
        Px = plan_fft(hnew; flags=FFTW.MEASURE)

        new(label, datasize, Γ, Dx, H, Π⅔, ϵ,
            hnew, unew, I₀, I₁, I₂, I₃, Px)
    end
end

4.2.1 Matsuno Right-Hand Side

"""
Compute right-hand side of the Matsuno equations

This is a complex spectral computation involving:
1. Fourier transforms
2. Multiplication by spectral operators
3. Nonlinear products
4. Inverse transforms

The computation is carefully organized to:
- Minimize redundant FFTs
- Use in-place operations
- Apply dealiasing to prevent aliasing errors

Parameters:
- U: [η̂, φ̂] (Fourier coefficients, modified in-place)
"""
function (m::Matsuno)(U::Array{ComplexF64,2})
    
    # === Compute η̂_t = -∂φ/∂x ===
    
    # unew = |k| * η̂
    Threads.@threads for i in eachindex(m.hnew)
        @inbounds m.unew[i] = m.Γ[i] * U[i,1]
    end
    
    # Inverse FFT to get |k| * ∂_x φ in physical space
    ldiv!(m.unew, m.Px, m.unew)
    
    # hnew = ∂_x η̂
    Threads.@threads for i in eachindex(m.hnew)
        @inbounds m.hnew[i] = m.Dx[i] * U[i,1]
    end
    
    # I₁ = ∂_x η (physical space)
    ldiv!(m.I₁, m.Px, m.hnew)
    
    # I₁ = (∂_x η) * (|k| * ∂_x φ) in physical space
    Threads.@threads for i in eachindex(m.unew)
        @inbounds m.unew[i] *= m.I₁[i]
    end
    
    # Forward transform the product
    mul!(m.I₁, m.Px, m.unew)
    
    # Apply dealiasing and finalize η̂_t contribution
    Threads.@threads for i in eachindex(m.hnew)
        @inbounds m.I₁[i] = m.I₁[i] * m.ϵ * m.Π⅔[i] - m.hnew[i]
    end
    
    # === Compute φ̂_t terms ===
    
    # Transform η and φ to physical space
    ldiv!(m.hnew, m.Px, view(U, :, 1))  # η in physical space
    ldiv!(m.unew, m.Px, view(U, :, 2))  # φ in physical space
    
    # First nonlinear term: η * (∂_x η) in physical space
    Threads.@threads for i in eachindex(m.hnew)
        @inbounds m.I₂[i] = m.hnew[i] * m.unew[i]
    end
    
    # Transform to Fourier space
    mul!(m.I₃, m.Px, m.I₂)
    
    # === Second part of φ̂_t ===
    
    # Apply Hilbert transform to φ̂
    Threads.@threads for i in eachindex(m.H)
        @inbounds U[i,1] = m.H[i] * U[i,2]
        @inbounds m.I₀[i] = m.Γ[i] * U[i,2]
    end
    
    # Transform |k| * φ̂ to physical space
    ldiv!(m.I₂, m.Px, m.I₀)
    
    # Multiply by η in physical space
    Threads.@threads for i in eachindex(m.hnew)
        @inbounds m.I₂[i] *= m.hnew[i]
    end
    
    # Transform back to Fourier space
    mul!(m.hnew, m.Px, m.I₂)
    
    # Combine all terms in φ̂_t
    Threads.@threads for i in eachindex(m.unew)
        @inbounds U[i,1] -= (m.I₃[i] * m.Dx[i] + m.hnew[i] * m.H[i]) * m.ϵ * m.Π⅔[i]
        @inbounds m.I₃[i] = m.unew[i]^2  # For kinetic energy term
    end
    
    # Transform kinetic energy term
    mul!(m.unew, m.Px, m.I₃)
    
    # Final φ̂_t equation
    Threads.@threads for i in eachindex(m.unew)
        @inbounds U[i,2] = m.I₁[i] - m.unew[i] * m.Dx[i] * m.ϵ / 2 * m.Π⅔[i]
    end

end

Key Computational Features:

• Spectral operators: Multiplication by \(\hat{\Gamma}\), \(D_x\), \(H\) in Fourier space
• Nonlinear terms: Products in physical space to avoid circular convolutions
• Dealiasing: The \(\Pi_{2/3}\) filter removes aliasing errors from the nonlinear products
• Work array reuse: Multiple intermediate computations use the same arrays to save memory