Periodic linear advection

New to Julia
1

I’m trying to numerically advect a scalar field β(x, y, t) only in the x-direction (no y-advection), mimicking a shifting pattern like β(x - ct). I’m using ModelingToolkit.jl to solve the PDE, but my code hangs after the discretize step without completing. I’m new to this and would appreciate help identifying what’s going wrong. My goal is to reproduce the expected x-translation behavior using the PDE formulation. Code is below — thanks in advance!


@variables x y t
@variables β(..)

∂t = Differential(t)
∂x = Differential(x)

# advection-diffusion constants
v = 100.0
R = 1.0
Lx = 160000.0
Ly = 160000.0
ω = 2π / Lx  # wave number for the sinusoidal field

# Initial condition function
β₀(x, y, t) = 1000.0 + 1000.0 * sin(ω * x) * sin(ω * y)

# advection equation (note: removed R coefficient from time derivative for pure advection)
eqn = ∂t(β(x, y, t)) ~ (-v) * ∂x(β(x, y, t))

# Periodic boundary conditions - this is key for getting the traveling wave behavior
bcs = [
    β(x, y, 0) ~ β₀(x, y, 0),           # Initial condition
    β(0, y, t) ~ β(Lx, y, t)
]

# space and time domains
dom = [x ∈ (0, Lx), y ∈ (0, Ly), t ∈ (0, 1800.0)]

# define PDE system
@named sys = PDESystem(eqn, bcs, dom, [x, y, t], [β(x, y, t)])
println("passed PDESystem")
# Discretization with periodic boundary conditions
dx = 2000.0  # Increased for better stability
dy = 2000.0


# Use periodic boundary conditions and upwind scheme for advection
disc = MOLFiniteDifference(
    [x => dx, y => dy], 
    t,
    advection_scheme = UpwindScheme(1),
    grid_align = CenterAlignedGrid()
)
println("passed MOLFiniteDifference")
# Apply periodic boundary conditions
prob = discretize(sys, disc)
println("passed discretize")
# Solve with appropriate time step
sol = solve(prob, Tsit5(), saveat=60.0)
println("passed solve")
## Extract grid and solution data
x_grid = sol[x]
y_grid = sol[y]
times = sol.t

# Get the solution data
all_sol_u_keys = collect(keys(sol.u))
beta_key_from_sol_u = all_sol_u_keys[1]
println("Key for beta data in sol.u: ", beta_key_from_sol_u)

beta_full_3D_array = sol.u[beta_key_from_sol_u]
initial_beta_values = beta_full_3D_array[:, :, 1]

nx = length(x_grid)
ny = length(y_grid)

println("Grid dimensions: nx = $nx, ny = $ny")
println("Time steps: ", length(times))

# Create animation
anim = @animate for i in 1:length(sol.t)
    current_beta_values = beta_full_3D_array[:, :, i]
    current_time = times[i]
    
    heatmap(
        x_grid ./ 1000, 
        y_grid ./ 1000, 
        current_beta_values',
        xlabel = "x (km)",
        ylabel = "y (km)",
        title = "β Field at t = $(round(current_time, digits=1)) seconds",
        color = :viridis,
        clims = (0, 2000),
        size = (800, 600),
        dpi = 150,
        aspect_ratio = :equal
    )
end

# Save the animation
gif(anim, "test.gif", fps = 8)
println("Animation saved as beta_advection_periodic.gif")```