using Pkg
Pkg.activate(".")
Pkg.instantiate()
Activating project at `c:\Users\610mo\Desktop\NumericalNavier`
using IJulia
using LaTeXStrings
using FileIO
using ImageIO
using Images
using Plots
using Base64
using LinearAlgebra
using ImageMagick
using Statistics
function display_image(path)
image = load(path)
display(image)
end
function display_gif(path)
display(MIME("text/html"), """<img src=$path>""")
end
function meshgrid(x, y)
X = [i for i in x, j in 1:length(y)]
Y = [j for i in 1:length(x), j in y]
return X, Y
end
meshgrid (generic function with 1 method)
Navier Stokes for Numerical Blokes¶
The Navier-Stokes equations are partial differential equations which describe the motion of viscous fluid substances, developed from 1822 (Navier) to 1842-1850 (Stokes).¶
These equations have shown to be notoriously hard to solve analytically. Thus, the goal of the project is the following:¶
- reformulate the problem into one that can be solved numerically
- analyze the convergence of the solution
- analyze its well-possendess
In this notebook we still focus on the Navier-Stokes equations for incompressible 2D fluid, which are formulated as such:¶
display_image("images/navier-stokes-2d.png")
Let us begin by addressing objective 1:The reformulation of the problem into one that can be solved numerically¶
We choose to solve this problem by using finite-difference-methods (FDM) which are a class of numerical techniquesforsolving differenial equations by approximating derivatvies with finite differences. Both the spatial domain and time domain (if applicable) are discretized, and the values of the solution at the end poins of the intervals are approximated by solving algebraic equations.
To give you a gist of what this means, lets consider the following simple example: $$ u'(x) = 3u(x) + 2\\ \Rightarrow u(x+h) \approx u(x) + h(3u(x)+2) $$
The discretization above is known as "forward difference." There are many different discretization schemes and we will analyze the effictiveness of the following on solving our (PDE):
- forward time, backward space
- forward time, center space (Lax-Friedrichs)
- forward time, forward space
To get a feel for how discretization works, we will first focus on Diffusion in 1D, then Burger's Equation in 2D, and then return to our initial (PDE). Burger's equation is a convection-diffusion equation whose analytical solutions are well known (and thus it is easy to calculate the error of our solution).
In 1D, the diffusion equation reads¶
$$\frac{\partial u}{\partial t} = \nu \frac{\partial^2 u}{\partial x^2}$$
Discretizing the first partial is trivial, but how can we do so for the second? Consider the Taylor expansion of $u_{i+1}$ and $u_{i-1}$ around $u_i$.
display_image("images/taylor-expansion.png")
Adding the two, we see that
display_image("images/2nd-partial-discritization.png")
From here, we obtain
display_image("images/diffusion-1d.png")
And thus our equation is simply
display_image("images/diffusion-solved.png")
Given an initial condition, this equation can be solved easily. Let's do so now. For demonstrative purposes (which will be explained later on), we define our initial conditions to be a hat function.
function animate_wave_propagation_diffusion_1d(nx=81, max_nt=100, sigma=.3, nu=0.2)
dx = 2.0 / (nx-1) # Spatial step size
dt = sigma*dx^2/nu
u = ones(nx) # Initial condition, wave at rest
# Set initial wave pulse
u[Int(floor(0.5 / dx)) + 1 : Int(floor(1 / dx)) + 1] .= 2
u_values = zeros(max_nt+1, nx) # Matrix to store u values at each time step
u_values[1, :] = u # Store initial condition
anim = @animate for nt in 1:max_nt
un = copy(u) # Copy the current state of u into un
for i in 2:nx-1 # Iterate over spatial grid points, leaving out the first point for boundary condition
# Finite difference method for convection equation
u[i] = un[i] + nu * dt / dx^2 * (un[i+1] - 2 * un[i] + un[i-1])
end
u_values[nt+1, :] = u # Store u values at each time step
x = range(0, stop=1, length=nx) # Spatial domain
plot(x, u, legend=false, ylims=(0, 3)) # Plot current state with fixed y-axis limits
end
return u_values, anim
end
animate_wave_propagation_diffusion_1d (generic function with 5 methods)
u_values, anim = animate_wave_propagation_diffusion_1d();
gif(anim, "gifs/wave_propagation_diffusion_1d_81.gif", fps = 10)
┌ Info: Saved animation to c:\Users\610mo\Desktop\NumericalNavier\gifs\wave_propagation_diffusion_1d_81.gif └ @ Plots C:\Users\610mo\.julia\packages\Plots\ju9dp\src\animation.jl:156
Let's look at it again but with a slightly smaller step size.
u_values, anim = animate_wave_propagation_diffusion_1d(31);
gif(anim, "gifs/wave_propagation_diffusion_1d_31.gif", fps = 10)
┌ Info: Saved animation to c:\Users\610mo\Desktop\NumericalNavier\gifs\wave_propagation_diffusion_1d_31.gif └ @ Plots C:\Users\610mo\.julia\packages\Plots\ju9dp\src\animation.jl:156
The result is as expected (almost), but I bet you couldn't spot the problem with your naked eye. We will address this issue later on when it becomes more apparent. (Hint: how does step size relate to numerical stability).
The 1D Diffusion Equation is a relatively simple BVP whose solution is somewhat (not really at all tbh but compared to others) straightforward (ask me after class if you really want to know how it's obtained). It's closed form solution can be written as $$\begin{equation} u(x, t) = 1 + \sum_{n=1}^{\infty} \left[ \frac{2}{n \pi} \left(\cos\left(\frac{n \pi}{4}\right) - \cos\left(\frac{n \pi}{2}\right) \right) \sin\left(\frac{n \pi x}{2}\right)\right] e^{-\left(\frac{n \pi}{2}\right)^2 \nu t} \end{equation}$$ Let's code it up and see how well we did.
function diffusion_1d_u(x, t, D, n_max)
u_xt = 1.0
for n in 1:n_max
B_n = 2 * (cos(pi * n / 4) - cos(pi * n / 2)) / (pi * n)
u_xt += B_n * sin(n * pi * x / 2) * exp(-D * (n * pi / 2)^2 * t)
end
return u_xt
end
diffusion_1d_u (generic function with 1 method)
# Simulation parameters
nx = 81
L = 2.0
max_nt = 100
nu = 0.2
sigma = 0.3
n = 50 # Number of Fourier terms
dx = L / (nx - 1)
dt = sigma * dx^2 / nu
# Initial condition
u = ones(nx)
u[Int(floor(0.5 / dx)) + 1 : Int(floor(1 / dx)) + 1] .= 2
# Setup for storing results
u_values = zeros(max_nt+1, nx)
u_values[1, :] = u
# Time stepping
for nt in 1:max_nt
un = copy(u)
for i in 2:nx-1
u[i] = un[i] + nu * dt / dx^2 * (un[i+1] - 2 * un[i] + un[i-1])
end
u_values[nt+1, :] = u
end
# Compute the closed-form solution at each time step
x = range(0, stop=L, length=nx)
closed_form_solutions = [diffusion_1d_u(xi, t, nu, n) for t in 0:dt:(max_nt*dt), xi in x]
# Calculate errors
errors = abs.(u_values .- closed_form_solutions)
max_error = maximum(errors)
avg_error = mean(errors)
rmse = sqrt(mean(errors.^2))
mape = mean(abs.(u_values .- closed_form_solutions) ./ closed_form_solutions)*100
# Output error metrics
println("Maximum error: $max_error")
println("Average error: $avg_error")
println("RMSE: $rmse")
println("MAPE: $mape")
# Plotting the error evolution
plot(0:max_nt, maximum(errors, dims=2), xlabel="Time step", ylabel="Maximum error", title="Error evolution over time", legend=false)
Maximum error: 0.500179281114093 Average error: 0.012530598188992187 RMSE: 0.022172727441353953 MAPE: 0.8820982088348815
Not bad! With only 81 discretization points, we got a MAPE of less than 1%. With more discretization points, it can only go down from here.
There's more! CFL is a condition which must be met for the numerical stability of our simulation. It rquires that our step size remains proportional to our discretization such that the <property> which is diffusing does not overstep the bounds of spatially adjacent cells. That is, the timestep is upper bounded. In explicit cases, this bound is usually 1. That is, the "Courant Number", defined as
$$C = u\frac{dt}{dx} \le1$$
where $u$ represents velocity, is lte to 1. Our discretization is static and challenging to recompute at each step in the simulation, as it would require the lengthening or shortening of various data structures. Our velocity is, however, dynamic in various cases, leading to a Courant number which may vary over the course of the simulation. We can therefore dynamically recompute the timestep size to remain within this bound.
Below is a demonstration on 1D convection:
nx = 81 #position step count
@show dx = 1.0 / 40 #space each position step covers
total_time = 1 #1 second of total time
initial_dt = .01
elapsed_time = 0
u = ones(nx)*1
u[11:21].= 2;
anim = @animate while elapsed_time < total_time
max_u = maximum(u) #find fastest moving <property>
dt = 0.9 * dx / max_u #max_u/dx <= 1*dt, compute dynamic timestep
global elapsed_time = dt + elapsed_time
Courant = max_u*dt/dx
@show elapsed_time
@show Courant
u2 = copy(u);
plot(u, ylimit=[0,3])
for space in 2:nx
u[space] = u2[space]-u2[space]*dt/dx*(u2[space]-u2[space-1]);
end
end
gif(anim, "./gifs/waveOverTime.gif", fps = 10)
dx = 1.0 / 40 = 0.025 elapsed_time = 0.011250000000000001 Courant = 0.9 elapsed_time = 0.022500000000000003 Courant = 0.9 elapsed_time = 0.03375 Courant = 0.9 elapsed_time = 0.045000000000000005 Courant = 0.9 elapsed_time = 0.05625000000000001 Courant = 0.9 elapsed_time = 0.0675 Courant = 0.9 elapsed_time = 0.07875 Courant = 0.9 elapsed_time = 0.09 Courant = 0.9 elapsed_time = 0.10124999999999999 Courant = 0.9 elapsed_time = 0.11249999999999999 Courant = 0.9 elapsed_time = 0.12374999999999999 Courant = 0.9 elapsed_time = 0.13500000000236761 Courant = 0.9 elapsed_time = 0.14625000003947036 Courant = 0.9 elapsed_time = 0.15750000036405654 Courant = 0.9 elapsed_time = 0.16875000247291805 Courant = 0.9 elapsed_time = 0.1800000138952206 Courant = 0.9 elapsed_time = 0.1912500687750817 Courant = 0.9 elapsed_time = 0.20250031169103203 Courant = 0.9 elapsed_time = 0.2137513266415551 Courant = 0.9 elapsed_time = 0.22500539049050494 Courant = 0.9 elapsed_time = 0.2362710673614378 Courant = 0.9 elapsed_time = 0.24757862362206068 Courant = 0.9 elapsed_time = 0.2590206367153264 Courant = 0.9 elapsed_time = 0.2708094802587342 Courant = 0.9 elapsed_time = 0.2826082224428979 Courant = 0.9 elapsed_time = 0.2945001766307793 Courant = 0.9 elapsed_time = 0.3067097842350265 Courant = 0.9 elapsed_time = 0.31894508727433074 Courant = 0.9 elapsed_time = 0.3312658523823315 Courant = 0.9 elapsed_time = 0.3438797233165884 Courant = 0.9 elapsed_time = 0.35650603011814963 Courant = 0.9 elapsed_time = 0.36921080460871836 Courant = 0.9 elapsed_time = 0.38217707493558495 Courant = 0.9 elapsed_time = 0.3951585683178615 Courant = 0.9 elapsed_time = 0.4082070901301655 Courant = 0.9 elapsed_time = 0.42147772684728746 Courant = 0.9 elapsed_time = 0.43478514419870573 Courant = 0.9 elapsed_time = 0.44814475513226243 Courant = 0.9 elapsed_time = 0.4616829073789103 Courant = 0.9 elapsed_time = 0.4752920631721925 Courant = 0.9 elapsed_time = 0.4889374364398899 Courant = 0.9 elapsed_time = 0.5027176900626145 Courant = 0.9 elapsed_time = 0.5166085934001614 Courant = 0.9 elapsed_time = 0.5305205616152779 Courant = 0.9 elapsed_time = 0.5445270957949846 Courant = 0.9 elapsed_time = 0.5586832611409489 Courant = 0.9 elapsed_time = 0.5728474306456454 Courant = 0.9 elapsed_time = 0.5870715898438413 Courant = 0.9 elapsed_time = 0.6014551484073926 Courant = 0.9 elapsed_time = 0.6158605459384817 Courant = 0.9 elapsed_time = 0.6302986223674919 Courant = 0.9 elapsed_time = 0.644840005133007 Courant = 0.9 elapsed_time = 0.6594771044181736 Courant = 0.9 elapsed_time = 0.6741272243104618 Courant = 0.9 elapsed_time = 0.6888376309393197 Courant = 0.9 elapsed_time = 0.703692104428226 Courant = 0.9 elapsed_time = 0.7185515880928749 Courant = 0.9 elapsed_time = 0.733440371181698 Courant = 0.9 elapsed_time = 0.7484152334438278 Courant = 0.9 elapsed_time = 0.763480499140763 Courant = 0.9 elapsed_time = 0.7785547567467803 Courant = 0.9 elapsed_time = 0.7936738985141376 Courant = 0.9 elapsed_time = 0.8089006166347136 Courant = 0.9 elapsed_time = 0.8241628864503017 Courant = 0.9 elapsed_time = 0.8394429018080939 Courant = 0.9 elapsed_time = 0.8547806027093788 Courant = 0.9 elapsed_time = 0.8702300770323148 Courant = 0.9 elapsed_time = 0.8856806586187549 Courant = 0.9 elapsed_time = 0.9011560207603788 Courant = 0.9 elapsed_time = 0.9166980158093698 Courant = 0.9 elapsed_time = 0.9323228949019486 Courant = 0.9 elapsed_time = 0.9479527415260856 Courant = 0.9 elapsed_time = 0.9636123332295983 Courant = 0.9 elapsed_time = 0.979343627716045 Courant = 0.8999999999999999 elapsed_time = 0.9951358013560668 Courant = 0.9000000000000002 elapsed_time = 1.0109355731354688 Courant = 0.9000000000000002
┌ Info: Saved animation to c:\Users\610mo\Desktop\NumericalNavier\gifs\waveOverTime.gif └ @ Plots C:\Users\610mo\.julia\packages\Plots\ju9dp\src\animation.jl:156
Lets repeat the diffusion experiment above but see if we can obtain better error values or the same error in fewer timesteps using dynamic timestepping, over the same window of real undescretized time!
The courant bound and courant number for diffusion is known to be: $$C = v\frac{dt}{dx^{2}}\le\frac{1}{2}$$
Since v, or "nu" is constant, this does not need to be dynamically recomputed, maping for a less interesting situation overall. Still, let's see what happens...
# Simulation parameters
nx = 81
L = 2.0
max_nt = 100
nu = 0.2
nt = 100
sigma = 0.3
n = 50 # Number of Fourier terms
dx = L / (nx - 1)
dt = (dx^2)/(2*nu) #set dt to exactly the max value within CFL
Courant = nu*dt/dx^2
elapsed_time = 0
total_time = dt*nt
@show total_time
@show Courant
# Initial condition
u = ones(nx)
u[Int(floor(0.5 / dx)) + 1 : Int(floor(1 / dx)) + 1] .= 2
# Setup for storing results
u_values = zeros(max_nt+1, nx)
u_values[1, :] = u
# Compute the closed-form solution at each time step
x = range(0, stop=L, length=nx)
closed_form_solutions = [diffusion_1d_u(xi, t, nu, n) for t in 0:dt:(total_time), xi in x]
# Time stepping
anim = @animate while elapsed_time < total_time
global elapsed_time = dt + elapsed_time
un = copy(u)
plot(x, [u, diffusion_1d_u.(x, elapsed_time, nu, n)])
for i in 2:nx-1
u[i] = un[i] + nu * dt / dx^2 * (un[i+1] - 2 * un[i] + un[i-1])
end
u_values[nt+1, :] = u
end
# Calculate errors
errors = abs.(u_values .- closed_form_solutions)
max_error = maximum(errors)
avg_error = mean(errors)
rmse = sqrt(mean(errors.^2))
mape = mean(abs.(u_values .- closed_form_solutions) ./ closed_form_solutions)*100
# Output error metrics
println("Maximum error: $max_error")
println("Average error: $avg_error")
println("RMSE: $rmse")
println("MAPE: $mape")
# Plotting the error evolution
plot(0:max_nt, maximum(errors, dims=2), xlabel="Time step", ylabel="Maximum error", title="Error evolution over time", legend=false)
total_time = 0.15625000000000003 Courant = 0.5 Maximum error: 2.004569246031523 Average error: 1.2221701979675559 RMSE: 1.271895368161698 MAPE: 98.04843393414171
gif(anim, "./gifs/courantDiffusion.gif", fps=20)
┌ Info: Saved animation to c:\Users\610mo\Desktop\NumericalNavier\gifs\courantDiffusion.gif └ @ Plots C:\Users\610mo\.julia\packages\Plots\ju9dp\src\animation.jl:156
# Simulation parameters
nx = 81
L = 2.0
max_nt = 100
nu = 0.2
nt = 100
sigma = 0.3
n = 50 # Number of Fourier terms
dx = L / (nx - 1)
dt = 0.95 * (dx^2)/(2*nu) #set dt to exactly the max value within CFL
Courant = nu*dt/dx^2
elapsed_time = 0
total_time = dt*nt
@show total_time
@show Courant
# Initial condition
u = ones(nx)
u[Int(floor(0.5 / dx)) + 1 : Int(floor(1 / dx)) + 1] .= 2
# Setup for storing results
u_values = zeros(max_nt+1, nx)
u_values[1, :] = u
# Compute the closed-form solution at each time step
x = range(0, stop=L, length=nx)
closed_form_solutions = [diffusion_1d_u(xi, t, nu, n) for t in 0:dt:(total_time), xi in x]
# Time stepping
anim = @animate while elapsed_time < total_time
global elapsed_time = dt + elapsed_time
un = copy(u)
plot(x, [u, diffusion_1d_u.(x, elapsed_time, nu, n)])
for i in 2:nx-1
u[i] = un[i] + nu * dt / dx^2 * (un[i+1] - 2 * un[i] + un[i-1])
end
u_values[nt+1, :] = u
end
# Calculate errors
errors = abs.(u_values .- closed_form_solutions)
max_error = maximum(errors)
avg_error = mean(errors)
rmse = sqrt(mean(errors.^2))
mape = mean(abs.(u_values .- closed_form_solutions) ./ closed_form_solutions)*100
# Output error metrics
println("Maximum error: $max_error")
println("Average error: $avg_error")
println("RMSE: $rmse")
println("MAPE: $mape")
# Plotting the error evolution
plot(0:max_nt, maximum(errors, dims=2), xlabel="Time step", ylabel="Maximum error", title="Error evolution over time", legend=false)
total_time = 0.14843750000000003 Courant = 0.475 Maximum error: 2.0050663402841864 Average error: 1.2222434342386443 RMSE: 1.2726990615149791 MAPE: 98.04813134051798
gif(anim, "./gifs/courantDiffusion90CFL.gif", fps=20)
#much better. Still bad, but much better...
┌ Info: Saved animation to c:\Users\610mo\Desktop\NumericalNavier\gifs\courantDiffusion90CFL.gif └ @ Plots C:\Users\610mo\.julia\packages\Plots\ju9dp\src\animation.jl:156
Finally, let's see what happens when the CFL gets bad, say, 1.05% of CMax.
# Simulation parameters
nx = 81
L = 2.0
max_nt = 100
nu = 0.2
nt = 100
sigma = 0.3
n = 50 # Number of Fourier terms
dx = L / (nx - 1)
dt = 1.05 * (dx^2)/(2*nu) #set dt to exactly the max value within CFL
Courant = nu*dt/dx^2
elapsed_time = 0
total_time = dt*nt
@show total_time
@show Courant
# Initial condition
u = ones(nx)
u[Int(floor(0.5 / dx)) + 1 : Int(floor(1 / dx)) + 1] .= 2
# Setup for storing results
u_values = zeros(max_nt+1, nx)
u_values[1, :] = u
# Compute the closed-form solution at each time step
x = range(0, stop=L, length=nx)
closed_form_solutions = [diffusion_1d_u(xi, t, nu, n) for t in 0:dt:(total_time), xi in x]
# Time stepping
anim = @animate while elapsed_time < total_time
global elapsed_time = dt + elapsed_time
un = copy(u)
plot(x, [u, diffusion_1d_u.(x, elapsed_time, nu, n)])
for i in 2:nx-1
u[i] = un[i] + nu * dt / dx^2 * (un[i+1] - 2 * un[i] + un[i-1])
end
u_values[nt+1, :] = u
end
# Calculate errors
errors = abs.(u_values .- closed_form_solutions)
max_error = maximum(errors)
avg_error = mean(errors)
rmse = sqrt(mean(errors.^2))
mape = mean(abs.(u_values .- closed_form_solutions) ./ closed_form_solutions)*100
# Output error metrics
println("Maximum error: $max_error")
println("Average error: $avg_error")
println("RMSE: $rmse")
println("MAPE: $mape")
# Plotting the error evolution
plot(0:max_nt, maximum(errors, dims=2), xlabel="Time step", ylabel="Maximum error", title="Error evolution over time", legend=false)
total_time = 0.16406250000000003 Courant = 0.525 Maximum error: 320.71735522166585 Average error: 2.8678999462281425 RMSE: 20.95971580995941 MAPE: 216.45079983899404
gif(anim, "./gifs/courantDiffusion1.5CFL.gif", fps=20)
#Terrible
┌ Info: Saved animation to c:\Users\610mo\Desktop\NumericalNavier\gifs\courantDiffusion1.5CFL.gif └ @ Plots C:\Users\610mo\.julia\packages\Plots\ju9dp\src\animation.jl:156
Discussing CFL further.
In our computation,
u[i] = un[i] + nu * dt / dx^2 * (un[i+1] - 2 * un[i] + un[i-1])
Let's look at a single step, and say nu*dt/dx^2 is 0.6, not 0.5, the known upper bound. This means that
u[i] = un[i] + 0.6(un[i+1] - 2 * un[i] + un[i-1])
That is, the amount of <property> at a point x will change by 0.6(right - 2 current + left) Let's say u = [0,0,0, 1,2,3, 2,1,0, 0,0] u at i+1 = [0,0,.6,1,0,1.8,2,1,.6,0,0] Note the 1.8 created at the former peak of our diffusing <property>. This, to a reasonable observer, should be a 2. That is, physically speaking, the center of a droplet of honey does not collapse before its walls do. Its walls must diffuse proportionately into it, but here, it loses more material than it should. This in turn generates smaller peaks out of its adjacent points, which will also lose more property/material than they should in ensuing steps, causing this chain of noisy exploding values.
This oscillatory behavior can be more intriguingly explained through Von Neumann analysis. By assuming error changes as a function of time, whether growing or shrinking, one may perform a fourier decomposition on error as defined by analytical - numerical solution at each timestep t.
The fourier decomposition leads us to a solution of the form: $$u_n[i] = \xi^n e^{ikx_i}$$ where $\xi$ is the growth factor per time step, $k$ is the wave number, and $x_i$ is the spatial location. Substituting this assumed solution into the discretized equation gives: $$\xi e^{ikx_i} = e^{ikx_i} + \nu \frac{\Delta t}{\Delta x^2} (e^{ikx_{i+1}} - 2 e^{ikx_i} + e^{ikx_{i-1}})$$
Using $$e^{ikx_{i\pm1}} = e^{ik(x_i \pm \Delta x)} = e^{ikx_i} e^{\pm ik\Delta x}$$ simplifies to:
$$\xi = 1 + \nu \frac{\Delta t}{\Delta x^2} (e^{ik\Delta x} - 2 + e^{-ik\Delta x})$$
$$\xi = 1 + \nu \frac{\Delta t}{\Delta x^2} (2 \cos(k \Delta x) - 2)$$
This cosine step is a bit quirky and was transcribed from Wikipedia with an admittedly limited understanding without a further careful review of fourier analysis, which we are unprepared to do. This simplifies further to: $$\xi = 1 - 4 \nu \frac{\Delta t}{\Delta x^2} \sin^2\left(\frac{k \Delta x}{2}\right)$$ For stability, the magnitude of $\xi$ must be less than or equal to 1: $$|1 - 4 \nu \frac{\Delta t}{\Delta x^2} \sin^2\left(\frac{k \Delta x}{2}\right)| \leq 1$$
Considering the worst-case scenario, where $$\sin^2\left(\frac{k \Delta x}{2}\right) = 1$$, we get: $$|1 - 4 \nu \frac{\Delta t}{\Delta x^2}| \leq 1$$ which leads to the CFL condition: $$0 \leq \nu \frac{\Delta t}{\Delta x^2} \leq \frac{1}{2} $$
That's pretty bad. Lots of flickering, likely caused by numerical instability of the CFL being precisely met. Let's instead set dt to 95% of the CFL.
Let's move on to Burger's Equation¶
In 2D, Burger's Equation reads $$ \frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} = \nu \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right) $$ $$ \frac{\partial v}{\partial t} + u\frac{\partial v}{\partial x} + v\frac{\partial v}{\partial y} = \nu \left( \frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial y^2} \right) $$
As before, we discretize each term obtaining
display_image("images/burgers-discretized.png")
Solving for our unknowns,
display_image("images/burgers-solved.png")
Burgers¶
###variable declarations
nx = 41
ny = 41
nt = 120
dx = 2 / (nx - 1)
dy = 2 / (ny - 1)
sigma = .0009
nu = .001
dt = sigma * dx * dy / nu
x = 0:dx:2
y = 0:dy:2
u = ones(ny, nx) # create a 1xn vector of 1's
v = ones(ny, nx)
un = ones(ny, nx)
vn = ones(ny, nx)
###Assign initial conditions
##set hat function I.C. : u(.5<=x<=1 && .5<=y<=1 ) is 2
u[Int(.5 / dy):Int(1 / dy + 1),Int(.5 / dx):Int(1 / dx + 1)] .= 2;
##set hat function I.C. : u(.5<=x<=1 && .5<=y<=1 ) is 2
v[Int(.5 / dy):Int(1 / dy + 1),Int(.5 / dx):Int(1 / dx + 1)] .= 2;
heatmap(x, y, u, color=cgrad([:lightblue,:white,:red]), clim=(1,2))
heatmap(x, y, v, color=cgrad([:lightblue,:white,:red]), clim=(1,2))
anim = @animate for n in 1:nt+1
un = copy(u)
vn = copy(v)
u[2:end-1, 2:end-1] .= (un[2:end-1, 2:end-1] -
dt / dx .* un[2:end-1, 2:end-1] .*
(un[2:end-1, 2:end-1] .- un[2:end-1, 1:end-2]) -
dt / dy .* vn[2:end-1, 2:end-1] .*
(un[2:end-1, 2:end-1] .- un[1:end-2, 2:end-1]) +
nu * dt / dx^2 .*
(un[2:end-1,3:end] - 2 .* un[2:end-1, 2:end-1] .+ un[2:end-1, 1:end-2]) +
nu * dt / dy^2 .*
(un[3:end, 2:end-1] - 2 .* un[2:end-1, 2:end-1] .+ un[1:end-2, 2:end-1]))
v[2:end-1, 2:end-1] .= (vn[2:end-1, 2:end-1] -
dt / dx .* un[2:end-1, 2:end-1] .*
(vn[2:end-1, 2:end-1] .- vn[2:end-1, 1:end-2]) -
dt / dy .* vn[2:end-1, 2:end-1] .*
(vn[2:end-1, 2:end-1] .- vn[1:end-2, 2:end-1]) +
nu * dt / dx^2 .*
(vn[2:end-1, 3:end] - 2 .* vn[2:end-1, 2:end-1] .+ vn[2:end-1, 1:end-2]) +
nu * dt / dy^2 .*
(vn[3:end, 2:end-1] - 2 .* vn[2:end-1, 2:end-1] .+ vn[1:end-2, 2:end-1]))
u[1, :] .= 1
u[end, :] .= 1
u[:, 1] .= 1
u[:, end] .= 1
v[1, :] .= 1
v[end, :] .= 1
v[:, 1] .= 1
v[:, end] .= 1
heatmap(x, y, u, color=cgrad([:lightblue,:white,:red]), clim=(1,2))
#plot(x, y, u, st=:surface, color=cgrad([:black,:white]), zlimits=[0.99, 2.01])
end
gif(anim, "gifs/burgers_2d.gif", fps = 60)
┌ Info: Saved animation to c:\Users\610mo\Desktop\NumericalNavier\gifs\burgers_2d.gif └ @ Plots C:\Users\610mo\.julia\packages\Plots\ju9dp\src\animation.jl:156
As before, lets compare with the known closed form solution. $$\alpha = \frac{1}{4 \nu t}, \quad \beta = -\frac{1}{4 \nu}$$ $$u(x, y, t, \nu) = v(x, y, t, \nu) = \frac{1 + \exp\left(\beta (x + y - t)\right)}{1 + \exp\left(\beta (x + y - t)\right) + \exp\left(\alpha (x^2 + y^2)\right)} $$
# True solution
function burgers_closed_form(x, y, t, ν)
α = 1 / (4 * ν * t)
β = -1 / (4 * ν)
u_true = @. (1 + exp(β * (x + y - t))) / (1 + exp(β * (x + y - t)) + exp(α * (x^2 + y^2)))
v_true = @. (1 + exp(β * (x + y - t))) / (1 + exp(β * (x + y - t)) + exp(α * (x^2 + y^2)))
return u_true, v_true
end
# Error analysis
err_u = zeros(ny, nx)
err_v = zeros(ny, nx)
# Create gif
@gif for n in 1:nt+1
t = (n - 1) * dt
for i in 1:ny
for j in 1:nx
u_true, v_true = burgers_closed_form(x[j], y[i], t, nu)
err_u[i, j] = abs(u[i, j] - u_true)
err_v[i, j] = abs(v[i, j] - v_true)
end
end
p1 = heatmap(x, y, err_u', title="Error in u at t = $(round(t, digits=2))", xlabel="x", ylabel="y", color=:viridis, clim=(0, maximum(err_u)))
p2 = heatmap(x, y, err_v', title="Error in v at t = $(round(t, digits=2))", xlabel="x", ylabel="y", color=:viridis, clim=(0, maximum(err_v)))
plot(p1, p2, layout=(1, 2), size=(800, 400))
end
┌ Info: Saved animation to c:\Users\610mo\Desktop\NumericalNavier\tmp.gif └ @ Plots C:\Users\610mo\.julia\packages\Plots\ju9dp\src\animation.jl:156
Finally, we are equipped to take on our original (PDE)¶
Recall what our equations are:
display_image("images/navier-stokes.png")
We discretize as usual:
display_image("images/navier-stokes-discrete.png")
Rearranging to get the explicit form we want,
display_image("images/navier-stokes-solved.png")
Navier Stokes Channel Flow and Test Case Experiments¶
These have been designed in accordance with the test cases in our midterm report!
The below test cases and implementations of channel flow attempt to demonstrate known physical circumstances which can be evaluated against our numerical methods.
Function to build the boundary conditions¶
function build_up_b(rho, dt, dx, dy, u, v)
b = zeros(size(u))
b[2:end-1, 2:end-1] .= (rho * (1 / dt * ((u[2:end-1, 3:end] - u[2:end-1, 1:end-2]) / (2 * dx) +
(v[3:end, 2:end-1] - v[1:end-2, 2:end-1]) / (2 * dy)) -
((u[2:end-1, 3:end] - u[2:end-1, 1:end-2]) / (2 * dx)).^2 -
2 * ((u[3:end, 2:end-1] - u[1:end-2, 2:end-1]) / (2 * dy) .*
(v[2:end-1, 3:end] - v[2:end-1, 1:end-2]) / (2 * dx)) -
((v[3:end, 2:end-1] - v[1:end-2, 2:end-1]) / (2 * dy)).^2))
# Periodic BC Pressure @ x = 2
b[2:end-1, end] .= (rho * (1 / dt * ((u[2:end-1, 1] - u[2:end-1, end-1]) / (2 * dx) +
(v[3:end, end] - v[1:end-2, end]) / (2 * dy)) -
((u[2:end-1, 1] - u[2:end-1, end-1]) / (2 * dx)).^2 -
2 * ((u[3:end, end] - u[1:end-2, end]) / (2 * dy) .*
(v[2:end-1, 1] - v[2:end-1, end-1]) / (2 * dx)) -
((v[3:end, end] - v[1:end-2, end]) / (2 * dy)).^2))
# Periodic BC Pressure @ x = 0
b[2:end-1, 1] .= (rho * (1 / dt * ((u[2:end-1, 2] - u[2:end-1, end]) / (2 * dx) +
(v[3:end, 1] - v[1:end-2, 1]) / (2 * dy)) -
((u[2:end-1, 2] - u[2:end-1, end]) / (2 * dx)).^2 -
2 * ((u[3:end, 1] - u[1:end-2, 1]) / (2 * dy) .*
(v[2:end-1, 2] - v[2:end-1, end]) / (2 * dx)) -
((v[3:end, 1] - v[1:end-2, 1]) / (2 * dy)).^2))
return b
end
build_up_b (generic function with 1 method)
Solve for pressure field¶
function pressure_poisson_periodic(p, dx, dy, b, nit)
for q in 1:nit
pn = copy(p)
p[2:end-1, 2:end-1] .= ((pn[2:end-1, 3:end] + pn[2:end-1, 1:end-2]) * dy^2 +
(pn[3:end, 2:end-1] + pn[1:end-2, 2:end-1]) * dx^2) /
(2 * (dx^2 + dy^2)) -
dx^2 * dy^2 / (2 * (dx^2 + dy^2)) * b[2:end-1, 2:end-1]
# Periodic BC Pressure @ x = 2
p[2:end-1, end] .= ((pn[2:end-1, 1] + pn[2:end-1, end-1]) * dy^2 +
(pn[3:end, end] + pn[1:end-2, end]) * dx^2) /
(2 * (dx^2 + dy^2)) -
dx^2 * dy^2 / (2 * (dx^2 + dy^2)) * b[2:end-1, end]
# Periodic BC Pressure @ x = 0
p[2:end-1, 1] .= ((pn[2:end-1, 2] + pn[2:end-1, end]) * dy^2 +
(pn[3:end, 1] + pn[1:end-2, 1]) * dx^2) /
(2 * (dx^2 + dy^2)) -
dx^2 * dy^2 / (2 * (dx^2 + dy^2)) * b[2:end-1, 1]
# Wall boundary conditions, pressure
p[end, :] .= p[end-1, :] # dp/dy = 0 at y = 2
p[1, :] .= p[2, :] # dp/dy = 0 at y = 0
end
return p
end
pressure_poisson_periodic (generic function with 1 method)
nx = 41
ny = 41
nt = 10
nit = 50
c = 1
dx = 2 / (nx - 1)
dy = 2 / (ny - 1)
x = LinRange(0, 2, nx)
y = LinRange(0, 2, ny)
X, Y = meshgrid(x,y)
##physical variables
rho = 1
nu = .1
F = 1
dt = .01
#initial conditions
u = zeros(ny, nx)
un = zeros(ny, nx)
v = zeros(ny, nx)
vn = zeros(ny, nx)
p = ones(ny, nx)
pn = ones(ny, nx)
b = zeros(ny, nx);
function calculate_cfl(u, v, dt, dx, dy)
u_max = maximum(abs.(u))
v_max = maximum(abs.(v))
cfl = (u_max * dt / dx) + (v_max * dt / dy)
return cfl
end
calculate_cfl (generic function with 1 method)
udiff = 1.0
stepcount = 0
max_cfl = 0.5
while udiff > 0.001
un = copy(u)
vn = copy(v)
b = build_up_b(rho, dt, dx, dy, u, v)
p = pressure_poisson_periodic(p, dx, dy, b, nit)
@. u[2:end-1, 2:end-1] = (un[2:end-1, 2:end-1]-
un[2:end-1, 2:end-1] * dt / dx *
(un[2:end-1, 2:end-1] - un[2:end-1, 1:end-2]) -
vn[2:end-1, 2:end-1] * dt / dy *
(un[2:end-1, 2:end-1] - un[1:end-2, 2:end-1]) -
dt / (2 * rho * dx) * (p[2:end-1, 3:end] - p[2:end-1, 1:end-2]) +
nu * (dt / dx^2 *
(un[2:end-1, 3:end] - 2 * un[2:end-1, 2:end-1] + un[2:end-1, 1:end-2]) +
dt / dy^2 *
(un[3:end, 2:end-1] - 2 * un[2:end-1, 2:end-1] + un[1:end-2, 2:end-1])))
@. v[2:end-1,2:end-1] = (vn[2:end-1, 2:end-1] -
un[2:end-1, 2:end-1] * dt / dx *
(vn[2:end-1, 2:end-1] - vn[2:end-1, 1:end-2]) -
vn[2:end-1, 2:end-1] * dt / dy *
(vn[2:end-1, 2:end-1] - vn[1:end-2, 2:end-1]) -
dt / (2 * rho * dy) * (p[3:end, 2:end-1] - p[1:end-2, 2:end-1]) +
nu * (dt / dx^2 *
(vn[2:end-1, 3:end] - 2 * vn[2:end-1, 2:end-1] + vn[2:end-1, 1:end-2]) +
dt / dy^2 *
(vn[3:end, 2:end-1] - 2 * vn[2:end-1, 2:end-1] + vn[1:end-2, 2:end-1])))
u[1, :] .= 0
u[end, :] .= 0
v[1, :] .= 0
v[end, :] .= 0
udiff = abs.(sum(u) - sum(un)) / sum(u)
stepcount += 1
cfl = calculate_cfl(u, v, dt, dx, dy)
if cfl > max_cfl
println("CFL number is too large: $cfl.")
break
end
end
quiver(X[1:3:end, 1:3:end], Y[1:3:end, 1:3:end], quiver=(u[1:3:end, 1:3:end], v[1:3:end, 1:3:end]))
Convergence and Well-possedness¶
Now that we have obtained a numerical approximation to our solution, and are satisfied with our "error", lets discuss project objectives 2 & 3. Note that we put "error" in quotation marks because we are yet to prove that our solution converges as we further refine our discretization. The "closed-form" solutions we compared to are numerical approximations obtained from well known julia packages. However, if we cannot prove convergence, who's to say their numerical approximation is any better than ours. Additionally, let's say that we do converge to something, how do we know that we are converging to an approximation of the closed form solution?
A numerical method is said to be consistent is the discretization errors induced by the method decrease to zero as the discretization parameters approach 0. In the context of the local truncation error, we require $$\lim_{\Delta t, \Delta x \to 0}LTE = 0$$ Checking that this condition holds rigourosly can sometimes be challenging, Lagrange Remainder Theorem can help.
Conceptually, stability means that the numerical errors do not grow uncontrollably as the computation advances in time (CFL helps but is not enough).
Why are both needed?
- Consistency without stability might mean accurately representing the differential equation, but small errors could grow exponentially, leading to useless results.
- Stability without consistency might control error growth well, but the errors might not necessarily diminish with finer discretizations, potentially leading to a stable but inaccurate solution
^Lax-Richtmyer equivalence theorem:^ For a consistent finite difference method for a well-posed linear initial value problem, the method is convergent if and only if it is stable.
Stuff becomes harder to prove when we get to BVPs and introduce nonlinearity.
TEST CASES¶
# initializing different boundary conditions
function loopingBoundary(u, v) # velocity on one border is transfered to the other border
u[1, :] .= u[end, :]
u[:, 1] .= u[:, end]
v[1, :] .= v[end, :]
v[:, 1] .= v[:, end]
end
function wallBoundary(u, v) # velocity component on each border is set to zero, simulate no-slip wall
u[1, :] .= 0
u[end, :] .= 0
u[:, 1] .= 0
u[:, end] .= 0
v[1, :] .= 0
v[end, :] .= 0
v[:, 1] .= 0
v[:, end] .= 0
end
function inflowLeftOutflowRight(u, v, flow=1) # right flowing vector field, x component of vector field has constant inflow, all others are zero
u[1, :] .= flow
u[end, :] .= flow
u[:, 1] .= 0
u[:, end] .= 0
v[1, :] .= 0
v[end, :] .= 0
v[:, 1] .= 0
v[:, end] .= 0
end
inflowLeftOutflowRight (generic function with 2 methods)
function process(boundary, u, v, nx, ny)
nt = 10
nit = 50
c = 1
dx = 2 / (nx - 1)
dy = 2 / (ny - 1)
x = LinRange(0, 2, nx)
y = LinRange(0, 2, ny)
X, Y = meshgrid(x,y)
##physical variables
rho = 1
nu = .1
F = 1
dt = .01
#initial conditions
un = zeros(ny, nx)
vn = zeros(ny, nx)
p = ones(ny, nx)
pn = ones(ny, nx)
b = zeros(ny, nx);
udiff = 1.0
stepcount = 0
max_cfl = 0.5
while udiff > 0.001
un = copy(u)
vn = copy(v)
b = build_up_b(rho, dt, dx, dy, u, v)
p = pressure_poisson_periodic(p, dx, dy, b, nit)
@. u[2:end-1, 2:end-1] = (un[2:end-1, 2:end-1]-
un[2:end-1, 2:end-1] * dt / dx *
(un[2:end-1, 2:end-1] - un[2:end-1, 1:end-2]) -
vn[2:end-1, 2:end-1] * dt / dy *
(un[2:end-1, 2:end-1] - un[1:end-2, 2:end-1]) -
dt / (2 * rho * dx) * (p[2:end-1, 3:end] - p[2:end-1, 1:end-2]) +
nu * (dt / dx^2 *
(un[2:end-1, 3:end] - 2 * un[2:end-1, 2:end-1] + un[2:end-1, 1:end-2]) +
dt / dy^2 *
(un[3:end, 2:end-1] - 2 * un[2:end-1, 2:end-1] + un[1:end-2, 2:end-1])))
@. v[2:end-1,2:end-1] = (vn[2:end-1, 2:end-1] -
un[2:end-1, 2:end-1] * dt / dx *
(vn[2:end-1, 2:end-1] - vn[2:end-1, 1:end-2]) -
vn[2:end-1, 2:end-1] * dt / dy *
(vn[2:end-1, 2:end-1] - vn[1:end-2, 2:end-1]) -
dt / (2 * rho * dy) * (p[3:end, 2:end-1] - p[1:end-2, 2:end-1]) +
nu * (dt / dx^2 *
(vn[2:end-1, 3:end] - 2 * vn[2:end-1, 2:end-1] + vn[2:end-1, 1:end-2]) +
dt / dy^2 *
(vn[3:end, 2:end-1] - 2 * vn[2:end-1, 2:end-1] + vn[1:end-2, 2:end-1])))
boundary(u, v)
udiff = abs.(sum(u) - sum(un)) / sum(u)
stepcount += 1
cfl = calculate_cfl(u, v, dt, dx, dy)
if cfl > max_cfl
println("CFL number is too large: $cfl.")
break
end
end
return u, v
end
process (generic function with 1 method)
# uniform density with looping boundary condition
nx = 41
ny = 41
u = zeros(ny, nx)
v = zeros(ny, nx)
u, v = process(loopingBoundary, u, v, nx, ny)
quiver(X[1:3:end, 1:3:end], Y[1:3:end, 1:3:end], quiver=(u[1:3:end, 1:3:end], v[1:3:end, 1:3:end]))
# uniform density with wall boundary condition
nx = 41
ny = 41
u = zeros(ny, nx)
v = zeros(ny, nx)
u, v = process(wallBoundary, u, v, nx, ny)
quiver(X[1:3:end, 1:3:end], Y[1:3:end, 1:3:end], quiver=(u[1:3:end, 1:3:end], v[1:3:end, 1:3:end]))
# centered peak with wall boundary condition
nx = 41
ny = 41
u = zeros(ny, nx)
v = zeros(ny, nx)
u[Int(.9 / dy):Int(1.1 / dy),Int(0.9 / dx):Int(1.1 / dx)] .= 2;
v[Int(.9 / dy):Int(1.1 / dy + 1),Int(.9 / dx):Int(1.1 / dx + 1)] .= 2;
u, v = process(wallBoundary, u, v, nx, ny)
quiver(X[1:3:end, 1:3:end], Y[1:3:end, 1:3:end], quiver=(u[1:3:end, 1:3:end], v[1:3:end, 1:3:end]))
CFL number is too large: 0.6639384061761726.
# uniform density with inflow left and outflow right, wall conditions on the top and bottom
nx = 41
ny = 41
u = zeros(ny, nx)
v = zeros(ny, nx)
u, v = process(inflowLeftOutflowRight, u, v, nx, ny)
quiver(X[1:3:end, 1:3:end], Y[1:3:end, 1:3:end], quiver=(u[1:3:end, 1:3:end], v[1:3:end, 1:3:end]))
CFL number is too large: 1.052899407053174.