1.1 Julia implementation

1.1.1 Import Dependencies

using FFTW
using LinearAlgebra
using Plots
using BenchmarkTools

• FFTW: Fast Fourier Transform library for computing spatial derivatives efficiently
• LinearAlgebra: Vector and matrix operations (used for norm())
• Plots: Visualization and animation
• BenchmarkTools: Accurate timing of functions (handles JIT compilation)

1.1.2 Cartesian Mesh Structure

"""
    Mesh(xmin, xmax, nx, ymin, ymax, ny)

2D uniform Cartesian mesh for periodic domain with Fourier wavenumber arrays.

This structure stores:
- Physical space grid points (x, y)
- Fourier wavenumber grids (kx, ky)
- Grid dimensions (nx, ny)

The mesh uses periodic boundary conditions on [-π, π]².
        
Parameters:

  - `xmin, xmax`: x-domain boundaries
  - `nx`: number of grid points in x
  - `ymin, ymax`: y-domain boundaries
  - `ny`: number of grid points in y

"""
struct Mesh
    
    nx   :: Int64              # Number of grid points in x-direction
    ny   :: Int64              # Number of grid points in y-direction
    x    :: Vector{Float64}    # Physical space x-coordinates
    y    :: Vector{Float64}    # Physical space y-coordinates
    kx   :: Vector{Float64}    # Fourier wavenumbers in x-direction
    ky   :: Vector{Float64}    # Fourier wavenumbers in y-direction
    
    function Mesh(xmin, xmax, nx, ymin, ymax, ny)
        
        # Create physical space grid points
        # We exclude the right endpoint to enforce periodic boundary conditions
        # This ensures f(x_max) = f(x_min) for periodicity
        x = LinRange(xmin, xmax, nx+1)[1:end-1]
        y = LinRange(ymin, ymax, ny+1)[1:end-1]
        
        # Create Fourier wavenumber grids
        # Standard FFT ordering: [0, 1, 2, ..., -n/2, ..., -2, -1]
        # This matches FFTW's output convention
        kx = 2π / (xmax - xmin) * fftfreq(nx, nx)
        ky = 2π / (ymax - ymin) * fftfreq(ny, ny)

        new(nx, ny, x, y, kx, ky)
    end
end

Design Notes:

• Periodicity: Excluding the endpoint ensures the grid wraps around
• Wavenumber ordering: Follows FFTW convention for efficiency
• Memory layout: Arrays are stored in column-major order (Julia standard)

1.1.3 Exact Solution Function

"""
Compute the exact solution at a given time

The exact solution is computed by rotating the initial Gaussian:

f(t, x, y) = f₀(x cos(t) + y sin(t), -x sin(t) + y cos(t))

where f₀(x, y) = exp(-(x-1)²/0.1) × exp(-(y-1)²/0.1) is a Gaussian
centered at (1, 1) with width parameter 0.1.

Parameters:
- time: current time t
- mesh: 2D mesh structure
- shift: center location of the Gaussian (default: 1)

Returns:
- f: exact solution field (2D array)
"""
function exact(time, mesh::Mesh; shift=1)
    
    f = zeros(Float64, (mesh.nx, mesh.ny))
    
    # Iterate over all grid points
    for (i, x) in enumerate(mesh.x), (j, y) in enumerate(mesh.y)
        
        # Apply rotation matrix (inverse rotation to get the initial position)
        # If current position is (x, y), where was it at t=0?
        # Answer: rotate backward by -t
        xn = cos(time) * x - sin(time) * y
        yn = sin(time) * x + cos(time) * y
        
        # Evaluate the rotated Gaussian
        # This is the initial Gaussian evaluated at the rotated coordinates
        f[i, j] = exp(-(xn - shift)^2 / 0.1) * exp(-(yn - shift)^2 / 0.1)
    end

    f
end

Physical Interpretation:

• The function computes where the initial Gaussian would be at time \(t\) after rotating by angle \(t\)
• We store a Gaussian at \((1, 1)\) initially; at time \(t\) it’s rotated around the origin
• This allows us to compare our numerical solution with the true answer

1.1.4 Visualization of Initial Condition

mesh = Mesh(-π, π, 128, -π, π, 128)
f = exact(0.0, mesh)
contour(mesh.x, mesh.y, f; aspect_ratio=:equal, clims=(0., 1.))

This shows a smooth Gaussian bump that will rotate during the simulation.

1.1.5 Animation: Visualizing the Solution

function animation(tf, nt)
    """
    Create an animated visualization of the rotating Gaussian
    
    Parameters:
    - tf: final time
    - nt: number of frames
    """
    
    nx, ny = 64, 64
    xmin, xmax = -π, π
    ymin, ymax = -π, π
    mesh = Mesh(xmin, xmax, nx, ymin, ymax, ny)
    f = zeros(Float64, (nx, ny))
    dt = tf / nt
    t = 0
    
    @gif for n = 1:nt
        
        # Compute exact solution at current time
        f .= exact(t, mesh)
        t += dt
        
        # Plot the solution
        p = contour(mesh.x, mesh.y, f)
        plot!(p[1]; clims=(0., 1.), aspect_ratio=:equal)
        
        # Draw a circle of radius √2 to show the rotation center
        plot!(sqrt(2) .* cos.(-π:0.1:π+0.1), 
              sqrt(2) .* sin.(-π:0.1:π+0.1),
              label="Radius √2")
        
    end
end

animation(2π, 100)

The circle shows the rotation path of the Gaussian center.