Advection equation
Start from
\[
\partial_t \psi + u \partial_x \psi = 0.
\]
Weak form, \(x \in [0, 1]\), weight \(w(x)\):
\[
\partial_t \int w \psi + [u w \psi]_0^1 - \int u \psi \partial_x w = 0 \, .
\]
Basis functions \(w = \sum w_A N_A(x)\) etc:
\[
\underbrace{M_{AB}}_{\text{mass}} \partial_t \psi_A + \underbrace{K_{AB}}_{\text{stiffness}} \psi_A = \underbrace{F_B}_{\text{force}} \, .
\]
Key matrices
\[
\begin{aligned}
\partial_t \int w \psi + [u w \psi]_0^1 - \int u \psi \partial_x w &= 0 \, , \\
\underbrace{M_{AB}}_{\text{mass}} \partial_t \psi_A + \underbrace{K_{AB}}_{\text{stiffness}} \psi_A &= \underbrace{F_B}_{\text{force}} \, ,
\end{aligned}
\]
where
\[
\begin{aligned}
M_{AB} &= \int N_A N_B \, , &
K_{AB} &= -u \int N_A \partial_x N_B \, , \\
F_B &= \text{boundary terms} \, .
\end{aligned}
\]
Example
nodes, dx = np.linspace(0, 1, N_elements+1, retstep=True)
LM = np.zeros((2, N_elements), dtype=np.int64)
for e in range(N_elements):
...
# Global mass and stiffness matrix and force vector
M = np.zeros((N_elements, N_elements))
K = np.zeros((N_elements, N_elements))
F = np.zeros((N_elements,))
# Loop over elements
for e in range(N_elements):
m_e = mass(nodes[e:e+2])
k_e = stiffness_time(nodes[e:e+2], u)
for a in range(2):
A = LM[a, e]
for b in range(2):
B = LM[b, e]
if (A >= 0) and (B >= 0):
M[A, B] += m_e[a, b]
K[A, B] += k_e[a, b]
# Iterate
dt = 0.9*dx
t = 0
while t < t_end:
if t + dt > t_end:
dt = t_end - t + 1e-10
t += dt
# RK2 step 1
dpsidt = np.linalg.solve(M, F - K @ Psi[1:])
Psi1 = Psi.copy()
Psi1[1:] += dt * dpsidt
# RK2 step 2
dpsidt = np.linalg.solve(M, F - K @ Psi1[1:])
Psi1 = Psi.copy()
Psi[1:] = (Psi[1:] + Psi1[1:] + dt * dpsidt) / 2
Matrix structure
- Matrices are sparse: \(\sim N\) entries.
- Use sparse structures for efficiency.
- Changes setting/use patterns.
from scipy import sparse
K = sparse.lil_array((N_elements,
N_elements))
for e in range(N_elements):
...
K[A, B] += k_e[a, b]
Psi_interior = sparse.linalg.solve(
K.to_csr(), F)
