# 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#

$\rho, e \mapsto p, T$

## Ideal gas#

(23)#\begin{align} p &= \rho R T & e &= e(T) \end{align}
$p = (\gamma - 1) \rho e$
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

$\int_A \rho \mathbf u \cdot \hat{\mathbf n} .$

If mass is conserved in a volume $$V$$ with surface $$A$$, then the total mass inside the volume must evolve as

$\int_V \rho_t = \left( \int_V \rho \right)_t = - \underbrace{\int_A \rho\mathbf u \cdot \hat{\mathbf n}}_{\int_V \nabla\cdot (\rho\mathbf u)},$

where we have applied the divergence theorem. Dropping the integrals over arbitrary volumes, we have the evolution equation for conservation of mass.

$\rho_t + \nabla\cdot (\rho \mathbf u) = 0$

## 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

$(\rho \mathbf u)_t + \nabla\cdot\big[ \rho \mathbf u \otimes \mathbf u + p I - \boldsymbol \tau \big] = 0$

# 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

$\begin{split}\begin{pmatrix} \rho \\ \rho \mathbf u \end{pmatrix}_t + \nabla\cdot \begin{bmatrix} \rho \mathbf u \\ \rho \mathbf u \otimes \mathbf u + \rho a^2 I \end{bmatrix} = 0 \end{split}$

## 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

$\begin{split}\begin{pmatrix} \tilde \rho \\ \bar\rho \mathbf{\tilde u} \end{pmatrix}_t + \nabla\cdot \begin{bmatrix} \bar\rho \mathbf{\tilde u} \\ \tilde \rho a^2 I \end{bmatrix} = 0 \end{split}$

# Two forms of acoustic wave equation#

Divide the momentum equation through by background density and dropping the tildes yields the standard form.

$\begin{split}\begin{pmatrix} \rho \\ \mathbf u \end{pmatrix}_t + \nabla\cdot \begin{bmatrix} \bar\rho \mathbf u \\ \rho \frac{a^2}{\bar\rho} I \end{bmatrix} = 0 .\end{split}$

Examine second equation

$\frac{a^2}{\bar\rho} \nabla\cdot\big[ \rho I \big] = \frac{a^2}{\bar\rho} \nabla \rho$
and thus \begin{pmatrix} \rho \ \mathbf u \end{pmatrix}_t +

(24)#$\begin{bmatrix} & \bar\rho \nabla\cdot \\ \frac{a^2}{\bar\rho} \nabla & \\ \end{bmatrix}$
(25)#$\begin{pmatrix} \rho \\ \mathbf u \end{pmatrix}$

Letâ€™s differentiate the first equation,

$\rho_{tt} + \bar\rho\nabla\cdot(\mathbf u_t) = 0$
and substitute in the second equation
$\rho_{tt} = a^2 \nabla\cdot(\nabla \rho)$

• Note: we had to assume these derivatives exist!

We can reduce this to a first order system as

$\begin{split}\begin{pmatrix} \rho \\ \dot \rho \end{pmatrix}_t + \begin{bmatrix} & -I \\ -a^2 \nabla\cdot\nabla & \end{bmatrix} \begin{pmatrix} \rho \\ \dot\rho \end{pmatrix} = 0\end{split}$

## 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#

$\begin{split}\begin{pmatrix} \rho \\ \dot \rho \end{pmatrix}_t = \begin{bmatrix} & I \\ a^2 \nabla\cdot\nabla & \end{bmatrix} \begin{pmatrix} \rho \\ \dot\rho \end{pmatrix}\end{split}$
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)


# 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