4.1 Julia implementation

using FFTW, LinearAlgebra, Plots

4.1.1 Abstract Types and Interfaces

abstract type AbstractModel end
abstract type InitialData end

Design Pattern: Abstract types define interfaces that different models and initial conditions can implement. This allows:

4.1.2 Mesh Structure

"""
1D Fourier spectral mesh for periodic domain

Stores both physical space grid (x) and frequency space grid (k).
Uses periodic boundary conditions on [-L, L].

Create a periodic mesh on [-L, L] with N grid points
    
Parameters (from named tuple):

- L: half-width of domain (domain is [-L, L])
- N: number of grid points
"""
struct Mesh
    
    N    :: Int64              # Number of grid points
    x    :: Vector{Float64}    # Physical space grid points
    kmin :: Float64            # Minimum wavenumber
    kmax :: Float64            # Maximum wavenumber
    k    :: Vector{Float64}    # Fourier wavenumber grid

    function Mesh(param::NamedTuple)
        
        xmin = -Float64(param.L)
        xmax = Float64(param.L)
        N = param.N
        
        # Physical grid spacing
        dx = (xmax - xmin) / N
        
        # Create physical space grid (exclude right endpoint for periodicity)
        x = LinRange(xmin, xmax, N+1)[1:end-1]
        
        # Fourier wavenumber spacing
        dk = 2π / (N * dx)
        
        # Wavenumber bounds
        kmin = -N / 2 * dk
        kmax = (N / 2 - 1) * dk
        
        # Fourier wavenumber grid
        # Standard FFT ordering: [0, 1, 2, ..., -n/2, ..., -2, -1]
        k = dk .* fftfreq(N, N)

        new(N, x, kmin, kmax, k)
    end

end

Key Design Features:

4.1.3 Time Stepping Structure

"""
Time stepping parameters and metadata
"""
struct Times
    Nt   :: Int              # Number of time steps
    dt   :: Float64          # Time step size

    function Times(dt, tfin)
        """
        Create time stepping array
        
        Parameters:
        - dt: time step size
        - tfin: final simulation time
        """
        # Calculate total number of steps
        Nt = length(0:dt:tfin)
        new(Nt, dt)
    end

end

4.1.4 RK4 Time Integrator

"""
    RK4(params, model)

Runge-Kutta fourth-order time integrator for ODEs.

The RK4 method has global error O(dt⁴), making it very accurate
while remaining reasonably efficient.

The method for dy/dt = f(t,y) is:
  k₁ = f(tⁿ, yⁿ)
  k₂ = f(tⁿ + dt/2, yⁿ + dt/2·k₁)
  k₃ = f(tⁿ + dt/2, yⁿ + dt/2·k₂)
  k₄ = f(tⁿ + dt, yⁿ + dt·k₃)
  yⁿ⁺¹ = yⁿ + dt/6·(k₁ + 2k₂ + 2k₃ + k₄)

"""
mutable struct RK4
    
    Uhat :: Array{ComplexF64,2}    # Stage values for RK4
    dU   :: Array{ComplexF64,2}    # Accumulated derivative

    function RK4(param::NamedTuple, model::AbstractModel)
        """
        Initialize RK4 solver
        
        Parameters:
        - param: named tuple with parameters (including N)
        - model: the model defining the dynamics
        """
        
        n = param.N
        
        # Allocate work arrays
        Uhat = zeros(ComplexF64, (n, 2))  # For stages
        dU = zeros(ComplexF64, (n, 2))    # For accumulation

        new(Uhat, dU)
    end

end

4.1.4.1 RK4 Step Implementation

function step!(s::RK4,
               f!::AbstractModel,
               U::Array{ComplexF64,2},
               dt::Float64)
    """
    Perform one RK4 time step
    
    Parameters:
    - s: RK4 solver state
    - f!: right-hand side function (model dynamics)
    - U: current solution state [η, φ]
    - dt: time step size
    
    The function is modified in-place, updating U to the next time level.
    """
    
    # === Stage 1: k₁ = f(tⁿ, yⁿ) ===
    # Copy current state to work array
    Threads.@threads for i in eachindex(U)
        @inbounds s.Uhat[i] = U[i]
    end
    
    # Evaluate right-hand side
    f!(s.Uhat)
    
    # === Stage 2: k₂ = f(tⁿ + dt/2, yⁿ + dt/2·k₁) ===
    # Copy and accumulate k₁
    Threads.@threads for i in eachindex(U)
        @inbounds s.dU[i] = s.Uhat[i]              # dU = k₁
        @inbounds s.Uhat[i] = U[i] + dt/2 * s.Uhat[i]  # New y value
    end
    
    # Evaluate right-hand side
    f!(s.Uhat)
    
    # === Stage 3: k₃ = f(tⁿ + dt/2, yⁿ + dt/2·k₂) ===
    # Add 2k₂ and prepare for stage 3
    Threads.@threads for i in eachindex(U)
        @inbounds s.dU[i] += 2 * s.Uhat[i]        # dU = k₁ + 2k₂
        @inbounds s.Uhat[i] = U[i] + dt/2 * s.Uhat[i]
    end
    
    # Evaluate right-hand side
    f!(s.Uhat)
    
    # === Stage 4: k₄ = f(tⁿ + dt, yⁿ + dt·k₃) ===
    # Add 2k₃ and prepare for stage 4
    Threads.@threads for i in eachindex(U)
        @inbounds s.dU[i] += 2 * s.Uhat[i]         # dU = k₁ + 2k₂ + 2k₃
        @inbounds s.Uhat[i] = U[i] + dt * s.Uhat[i]
    end
    
    # Evaluate right-hand side
    f!(s.Uhat)
    
    # === Final update: yⁿ⁺¹ = yⁿ + dt/6·(k₁ + 2k₂ + 2k₃ + k₄) ===
    Threads.@threads for i in eachindex(U)
        @inbounds s.dU[i] += s.Uhat[i]              # dU = k₁ + 2k₂ + 2k₃ + k₄
        @inbounds U[i] += dt/6 * s.dU[i]            # Final update
    end

end

Algorithmic Notes:

4.1.5 Initial Conditions: Bell Curve

"""
Bell curve (Gaussian-like) initial condition for surface elevation.

This represents a smooth, localized wave packet.

Create bell curve initial condition

Parameters:
- p: named tuple with mesh parameters
- theta: shape parameter (higher = sharper peak)

Formula: η(x) = exp(-(|x|^θ) * ln(2))

This creates a smooth hump that is initially at rest (u = 0).
"""
struct BellCurve <: InitialData
    
    h :: Vector{ComplexF64}    # Surface elevation η
    u :: Vector{ComplexF64}    # Velocity potential φ

    function BellCurve(p::NamedTuple, theta::Real)
        
        mesh = Mesh(p)
        
        # Create smooth elevation profile
        h = zeros(ComplexF64, mesh.N)
        h .= exp.(-(abs.(mesh.x).^theta) .* log(2))
        
        # Initially at rest
        u = zeros(ComplexF64, mesh.N)

        new(h, u)
    end

end

Physical Interpretation:



CC BY-NC-SA 4.0 Pierre Navaro