so below is the code snippet i have for Julia:
using Printf
using Plots
gr()
# Parameters definition
Ny = 21 # Number of points in y-direction
Aspect_ratio = 10 # Ratio of channel length to width
Re = 100 # Reynolds number
dt = 0.008 # Time step
N_steps = 2000 # Total number of time steps
Plot_freq = 20 # Plot save frequency
Poisson_max_iter = 1000 # Total number of pressure-poisson iterations
tol = 1e-4 # Tolerance for pressure-poisson
# Function to compute an error in L2 norm
function l2_err(u1, u2)
return sqrt(sum((u1 .- u2).^2))
end
# Generating the 2D grid
function ndgrid(x_range, y_range)
Nx = length(x_range)
Ny = length(y_range)
coordinates_x = zeros(Nx, Ny)
coordinates_y = zeros(Nx, Ny)
for i in 1:Nx
for j in 1:Ny
coordinates_x[i, j] = x_range[i]
coordinates_y[i, j] = y_range[j]
end
end
return coordinates_x, coordinates_y
end
# Pressure-Poisson using the Jacobi method
function jacobi_method(p0, b, dx, dy, tolerance, max_it)
pnew = copy(p0)
it = 0
err = 1.0
while err > tolerance && it < max_it
p = copy(pnew)
# Calculate p_new for interior points
pnew[2:end-1, 2:end-1] = 0.25 * (
p[3:end, 2:end-1] + p[1:end-2, 2:end-1] + p[2:end-1, 3:end] + p[2:end-1, 1:end-2] -
dx^2 * b[2:end-1, 2:end-1]
)
# Apply boundary conditions
pnew[1, 2:end-1] .= pnew[2, 2:end-1]
pnew[end, 2:end-1] .= -pnew[end-1, 2:end-1]
pnew[:, 1] .= pnew[:, 2]
pnew[:, end] .= pnew[:, end-1]
err = l2_err(pnew, p)
it += 1
end
return pnew
end
# Computing diffusion, convection, and pressure gradient terms
function compute_diffusion(u, v, dx, dy, Re)
diffusion_x = (u[3:end, 2:end-1] - 2u[2:end-1, 2:end-1] + u[1:end-2, 2:end-1]) / dx^2 +
(u[2:end-1, 3:end] - 2u[2:end-1, 2:end-1] + u[2:end-1, 1:end-2]) / dy^2
diffusion_y = (v[3:end, 2:end-1] - 2v[2:end-1, 2:end-1] + v[1:end-2, 2:end-1]) / dx^2 +
(v[2:end-1, 3:end] - 2v[2:end-1, 2:end-1] + v[2:end-1, 1:end-2]) / dy^2
return diffusion_x / Re, diffusion_y / Re
end
function compute_convection(u, v, dx, dy)
convection_x = u[2:end-1, 2:end-1] .* (u[3:end, 2:end-1] - u[1:end-2, 2:end-1]) / (2dx) +
v[2:end-1, 2:end-1] .* (u[2:end-1, 3:end] - u[2:end-1, 1:end-2]) / (2dy)
convection_y = u[2:end-1, 2:end-1] .* (v[3:end, 2:end-1] - v[1:end-2, 2:end-1]) / (2dx) +
v[2:end-1, 2:end-1] .* (v[2:end-1, 3:end] - v[2:end-1, 1:end-2]) / (2dy)
return convection_x, convection_y
end
function compute_pressure_gradient(pressure, dx, dy)
pressure_gradient_x = (pressure[3:end, 2:end-1] - pressure[1:end-2, 2:end-1]) / (2dx)
pressure_gradient_y = (pressure[2:end-1, 3:end] - pressure[2:end-1, 1:end-2]) / (2dy)
return pressure_gradient_x, pressure_gradient_y
end
# Incompressible Navier Stokes Function
function main()
# Make an equidistant grid of dx=dy=0.1
cell_length = 2.0 / (Ny - 1)
Nx = (Ny - 1) * Aspect_ratio + 1
# Create grid
x_range = range(0.0, stop=2.0 * Aspect_ratio, length=Nx)
y_range = range(-1.0, stop=1.0, length=Ny)
# Generate a 2D grid
coordinates_x, coordinates_y = ndgrid(x_range, y_range)
# Initial conditions
u_prev = ones(Nx, Ny + 1) # u^*(t^* = 0) = 1
u_prev[:, 1] .= -u_prev[:, 2] # No-slip boundary at bottom wall
u_prev[:, end] .= -u_prev[:, end-1] # No-slip boundary at top wall
v_prev = zeros(Nx + 1, Ny) # v^*(t^* = 0) = 0
pressure_prev = zeros(Nx + 1, Ny + 1) # p^*(t^* = 0) = 0
# Pre-Allocate some arrays
u_int = zeros(size(u_prev))
u_next = zeros(size(u_prev))
v_int = zeros(size(v_prev))
v_next = zeros(size(v_prev))
# Initialize empty array to store frames
frames = []
for iter in 1:N_steps
# Compute the diffusion, convection and pressure gradient terms
diffusion_x, diffusion_y = compute_diffusion(u_prev, v_prev, cell_length, cell_length, Re)
convection_x, convection_y = compute_convection(u_prev, v_prev, cell_length, cell_length)
pressure_gradient_x, pressure_gradient_y = compute_pressure_gradient(pressure_prev, cell_length, cell_length)
# Update intermediate velocities
u_int[2:end-1, 2:end-1] = u_prev[2:end-1, 2:end-1] + dt * (
diffusion_x - convection_x - pressure_gradient_x
)
v_int[2:end-1, 2:end-1] = v_prev[2:end-1, 2:end-1] + dt * (
diffusion_y - convection_y - pressure_gradient_y
)
# Apply BC on u
u_int[:, 1] .= -u_int[:, 2]
u_int[:, end] .= -u_int[:, end-1]
u_int[1, :] .= 1.0
u_int[end, :] .= u_int[end-1, :]
# Apply BC on v
v_int[:, 1] .= -v_int[:, 2]
v_int[:, end] .= -v_int[:, end-1]
v_int[1, :] .= 0.0
v_int[end, :] .= 0.0
# Compute the divergence as it will be the rhs of the Pressure-Poisson
divergence = (u_int[2:end-1, 2:end-1] - u_int[1:end-2, 2:end-1]) / cell_length +
(v_int[2:end-1, 2:end-1] - v_int[2:end-1, 1:end-2]) / cell_length
pressure_poisson_rhs = divergence / dt
# Solve the pressure correction poisson problem
pressure_correction_prev = zeros(size(pressure_prev))
pressure_correction_next = jacobi_method(pressure_correction_prev, pressure_poisson_rhs,
cell_length, cell_length, tol, Poisson_max_iter)
# Update the pressure
pressure_next = pressure_prev + pressure_correction_next
# Correct the velocities to be incompressible
pressure_correction_gradient_x, pressure_correction_gradient_y = compute_pressure_gradient(pressure_correction_next, cell_length, cell_length)
u_next[2:end-1, 2:end-1] .= u_int[2:end-1, 2:end-1] - dt * pressure_correction_gradient_x
v_next[2:end-1, 2:end-1] .= v_int[2:end-1, 2:end-1] - dt * pressure_correction_gradient_y
# Enforce mass conservation
u_next[1, 2:end-1] .= 1.0
inflow_mass_rate_next = sum(u_next[1, 2:end-1])
outflow_mass_rate_next = sum(u_next[end-1, 2:end-1])
u_next[end, 2:end-1] .= u_next[end-1, 2:end-1] * inflow_mass_rate_next / outflow_mass_rate_next
# Again apply BC
u_next[:, 1] .= -u_next[:, 2]
u_next[:, end] .= -u_next[:, end-1]
v_next[1, 2:end-1] .= 0.0
v_next[end, 2:end-1] .= 0.0
v_next[:, 1] .= -v_next[:, 2]
v_next[:, end] .= -v_next[:, end-1]
# Compute error between u_prev and u_next
err_u = l2_err(u_prev, u_next)
println("Iteration $iter, Error in u: $err_u")
if err_u < 1e-8
println("Converged at iteration $iter")
break
end
# Advance in time
u_prev .= u_next
v_prev .= v_next
pressure_prev .= pressure_next
# Visualization
if iter % Plot_freq == 0
u_centered = (u_next[:, 2:end] + u_next[:, 1:end-1]) / 2
# Velocity plot
p1 = plot(coordinates_x, coordinates_y, u_centered, st=:surface, title="Velocity Field", xlabel="x", ylabel="y")
# Pressure plot
p2 = plot(coordinates_x, coordinates_y, pressure_next, st=:surface, title="Pressure Field", xlabel="x", ylabel="y")
# Combine both plots
p = plot(p1, p2, layout = @layout [a; b])
# Push the current frame to the frames array
push!(frames, deepcopy(p))
end
end
# Create the animation from the stored frames
anim = @animate for frame in frames
plot(frame)
end
return anim
end
# Initialize and run the solver
animation = main()
but it always returns the error
{
"name": "DimensionMismatch",
"message": "DimensionMismatch: arrays could not be broadcast to a common size; got a dimension with lengths 200 and 199",
"stack": "DimensionMismatch: arrays could not be broadcast to a common size; got a dimension with lengths 200 and 199
Stacktrace:
[1] _bcs1
@ ./broadcast.jl:555 [inlined]
[2] _bcs
@ ./broadcast.jl:549 [inlined]
[3] broadcast_shape
@ ./broadcast.jl:543 [inlined]
[4] combine_axes
@ ./broadcast.jl:524 [inlined]
[5] instantiate
@ ./broadcast.jl:306 [inlined]
[6] materialize
@ ./broadcast.jl:903 [inlined]
[7] compute_convection(u::Matrix{Float64}, v::Matrix{Float64}, dx::Float64, dy::Float64)
@ Main ~/Documents/filepath.ipynb:74
[8] main()
@ Main ~/Documents/filepath.ipynb:123
[9] top-level scope
@ ~/Documents/filepath.ipynb:221"
}
i’ve tried to ask my peers about this problem but either they don’t have it or are just as clueless as i am
i’m looking forward to suggestions 😀
i just expect the code to function without any warnings/errors, at least for now
New contributor
CFDn00b is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.
The error message says that either in your computation of convection_x or convection_y (whichever one is on the correct line), you are trying to do element-wise operations on arrays of different lengths. This is caused by the fact that uprev has size of Nx by Ny + 1, while vprev has size Nx + 1 by Ny. To fix this, you need to figure out what size everything actually needs to be for the math to line up.