2022-10-05 Factorization
Contents
2022-10-05 Factorization#
Last time#
FD methods in 2D
Cost profile
The need for fast algebraic solvers
Today#
Wave equation and Hamiltonians
Symplectic integrators
Sparse direct solvers
matrix orderings
impact on formulation
cost scaling
using Plots
default(linewidth=3)
using LinearAlgebra
using SparseArrays
function my_spy(A)
cmax = norm(vec(A), Inf)
s = max(1, ceil(120 / size(A, 1)))
spy(A, marker=(:square, s), c=:diverging_rainbow_bgymr_45_85_c67_n256, clims=(-cmax, cmax))
end
function plot_stability(Rz, title; xlims=(-2, 2), ylims=(-2, 2))
x = LinRange(xlims[1], xlims[2], 100)
y = LinRange(ylims[1], ylims[2], 100)
heatmap(x, y, (x, y) -> abs(Rz(x + 1im*y)), c=:bwr, clims=(0, 2), aspect_ratio=:equal, title=title)
end
struct RKTable
A::Matrix
b::Vector
c::Vector
function RKTable(A, b)
s = length(b)
A = reshape(A, s, s)
c = vec(sum(A, dims=2))
new(A, b, c)
end
end
function rk_stability(z, rk)
s = length(rk.b)
1 + z * rk.b' * ((I - z*rk.A) \ ones(s))
end
rk4 = RKTable([0 0 0 0; .5 0 0 0; 0 .5 0 0; 0 0 1 0], [1, 2, 2, 1] / 6)
function ode_rk_explicit(f, u0; tfinal=1., h=0.1, table=rk4)
u = copy(u0)
t = 0.
n, s = length(u), length(table.c)
fY = zeros(n, s)
thist = [t]
uhist = [u0]
while t < tfinal
tnext = min(t+h, tfinal)
h = tnext - t
for i in 1:s
ti = t + h * table.c[i]
Yi = u + h * sum(fY[:,1:i-1] * table.A[i,1:i-1], dims=2)
fY[:,i] = f(ti, Yi)
end
u += h * fY * table.b
t = tnext
push!(thist, t)
push!(uhist, u)
end
thist, hcat(uhist...)
end
ode_rk_explicit (generic function with 1 method)
Gas equations of state#
There are many ways to describe a gas
Name |
variable |
units |
---|---|---|
pressure |
\(p\) |
force/area |
density |
\(\rho\) |
mass/volume |
temperature |
\(T\) |
Kelvin |
(specific) internal energy |
\(e\) |
energy/mass |
entropy |
\(s\) |
energy/Kelvin |
Equation of state#
Ideal gas#
pressure(rho, T) = rho*T
contour(LinRange(0, 2, 30), LinRange(0, 2, 30), pressure, xlabel="\$\\rho\$", ylabel="\$T\$")
Conservation equations#
Mass#
Let \(\mathbf u\) be the fluid velocity. The mass flux (mass/time) moving through an area \(A\) is
If mass is conserved in a volume \(V\) with surface \(A\), then the total mass inside the volume must evolve as
where we have applied the divergence theorem. Dropping the integrals over arbitrary volumes, we have the evolution equation for conservation of mass.
Momentum#
The momentum \(\rho \mathbf u\) has a flux that includes
convection \(\rho \mathbf u \otimes \mathbf u\)
this is saying that each component of momentum is carried along in the vector velocity field
pressure \(p I\)
viscous \(-\boldsymbol\tau\)
A similar integral principle leads to the momentum equation
Simplifications#
Ignore viscous stress tensor \(\boldsymbol \tau\)
Ignore energy equation (not yet written) and assume constant temperature
\(p = a^2 \rho\) where \(a\) is the acoustic wave speed
Linearization#
Split each state variable into a mean state and a small fluctuation
\(\rho = \bar\rho + \tilde\rho\)
\(u = \bar u + \tilde u\)
Let \(\widetilde{\rho u} = (\bar\rho + \tilde\rho) (\bar u + \tilde u) - \bar\rho\bar u \approx \tilde \rho \bar u + \bar\rho \tilde u\), where we have dropped the second order term \(\tilde \rho\tilde u\) because both are assumed small.
We consider background state \(\bar u = 0\) and constant \(\bar\rho(x,y,t) = \bar\rho\). Then
Two forms of acoustic wave equation#
Divide the momentum equation through by background density and dropping the tildes yields the standard form.
Examine second equation
Let’s differentiate the first equation,
Note: we had to assume these derivatives exist!
We can reduce this to a first order system as
Question#
How is the problem size different?
What might we be concerned about in choosing the second formulation?
Laplacian in periodic domain#
function laplacian_matrix(n)
h = 2 / n
rows = Vector{Int64}()
cols = Vector{Int64}()
vals = Vector{Float64}()
wrap(i) = (i + n - 1) % n + 1
idx(i, j) = (wrap(i)-1)*n + wrap(j)
stencil_diffuse = [-1, -1, 4, -1, -1] / h^2
for i in 1:n
for j in 1:n
append!(rows, repeat([idx(i,j)], 5))
append!(cols, [idx(i-1,j), idx(i,j-1), idx(i,j), idx(i+1,j), idx(i,j+1)])
append!(vals, stencil_diffuse)
end
end
sparse(rows, cols, vals)
end
cond(Matrix(laplacian_matrix(5)))
2.9959163385932148e16
L = laplacian_matrix(10)
ev = eigvals(Matrix(L))
scatter(real(ev), imag(ev))
Wave operator#
function wave_matrix(n; a=1)
Z = spzeros(n^2, n^2)
L = laplacian_matrix(n)
[Z I; -a^2*L Z]
end
wave_matrix(2)
8Ă—8 SparseMatrixCSC{Float64, Int64} with 16 stored entries:
â‹… â‹… â‹… â‹… 1.0 â‹… â‹… â‹…
â‹… â‹… â‹… â‹… â‹… 1.0 â‹… â‹…
â‹… â‹… â‹… â‹… â‹… â‹… 1.0 â‹…
â‹… â‹… â‹… â‹… â‹… â‹… â‹… 1.0
-4.0 2.0 2.0 â‹… â‹… â‹… â‹… â‹…
2.0 -4.0 â‹… 2.0 â‹… â‹… â‹… â‹…
2.0 â‹… -4.0 2.0 â‹… â‹… â‹… â‹…
â‹… 2.0 2.0 -4.0 â‹… â‹… â‹… â‹…
A = wave_matrix(8; a=2) * .1
ev = eigvals(Matrix(A))
plot_stability(z -> rk_stability(z, rk4), "RK4", xlims=(-4, 4), ylims=(-4, 4))
scatter!(real(ev), imag(ev), color=:black)
Question: would forward Euler work?#
Example 2D wave solver with RK4#
n = 20
A = wave_matrix(n)
x = LinRange(-1, 1, n+1)[1:end-1]
y = x
rho0 = vec(exp.(-9*((x .+ 1e-4).^2 .+ y'.^2)))
sol0 = vcat(rho0, zero(rho0))
thist, solhist = ode_rk_explicit((t, sol) -> A * sol, sol0, h=.02)
size(solhist)
(800, 51)
@gif for tstep in 1:length(thist)
rho = solhist[1:n^2, tstep]
contour(x, y, reshape(rho, n, n), title="\$ t = $(thist[tstep])\$")
end
┌ Info: Saved animation to
│ fn = /home/jed/cu/numpde/slides/tmp.gif
â”” @ Plots /home/jed/.julia/packages/Plots/lW9ll/src/animation.jl:137