2021-10-15 Sparse Linear Algebra¶
Last time¶
Note on p-Laplacian activity
Weak forms and jump conditions
2D wave equations via gas dynamics
Appropriate time integrators
Hamiltonian structure and energy drift
Today¶
Recap and questions
Sparse direct solvers
matrix orderings
impact on formulation
cost scaling
Why iterative solvers
Discussion with Ioana Fleming
using Plots
using LinearAlgebra
using SparseArrays
default(linewidth=4)
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
function laplacian_matrix(n)
"Laplacian in n×n periodic domain"
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
laplacian_matrix (generic function with 1 method)
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?
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/1RWWg/src/animation.jl:114
Accuracy, conservation of mass with RK4¶
thist, solhist = ode_rk_explicit((t, sol) -> A * sol, sol0, h=.05,
tfinal=1)
tfinal = thist[end]
M = exp(Matrix(A*tfinal))
sol_exact = M * sol0
sol_final = solhist[:, end]
norm(sol_final - sol_exact)
0.02015111748435016
mass = vec(sum(solhist[1:n^2, :], dims=1))
plot(thist[2:end], mass[2:end] .- mass[1])
Conservation of energy with RK4¶
Hamiltonian structure¶
We can express the total energy for our system as a sum of kinetic and potential energy.
where we identify \(\rho\) as a generalized position and \(\dot\rho\) as generalized momentum. Hamilton’s equations state that the equations of motion are
where we have used the weak form to associate \(\int \nabla v \cdot \nabla u = v^T L u\).
function energy(sol, n)
L = laplacian_matrix(n)
rho = sol[1:end÷2]
rhodot = sol[end÷2+1:end]
kinetic = .5 * norm(rhodot)^2
potential = .5 * rho' * L * rho
kinetic + potential
end
ehist = [energy(solhist[:,i], n) for i in 1:length(thist)]
plot(thist, ehist)
Velocity Verlet integrator¶
function wave_verlet(n, u0; tfinal=1., h=0.1)
L = laplacian_matrix(n)
u = copy(u0)
t = 0.
thist = [t]
uhist = [u0]
irho = 1:n^2
irhodot = n^2+1:2*n^2
accel = -L * u[irho]
while t < tfinal
tnext = min(t+h, tfinal)
h = tnext - t
u[irho] += h * u[irhodot] + h^2/2 * accel
accel_next = -L * u[irho]
u[irhodot] += h/2 * (accel + accel_next)
accel = accel_next
t = tnext
push!(thist, t)
push!(uhist, copy(u))
end
thist, hcat(uhist...)
end
wave_verlet (generic function with 1 method)
thist, solhist = wave_verlet(n, sol0, h=.04)
@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/1RWWg/src/animation.jl:114