Elliptical PDE Solvers ¶ Motivation ¶ In many PDEs, we have the famous elliptical system of the form:
∇ 2 u = f L u = F \begin{aligned}
\nabla^2u &= f \\
\mathbf{Lu} &=\mathbf{F}
\end{aligned} ∇ 2 u Lu = f = F where ∇ 2 \nabla^2 ∇ 2 f f f u u u L \mathbf{L} L F \mathbf{F} F u \mathbf{u} u ∇ 2 \nabla^2 ∇ 2
u = L − 1 F \mathbf{u} = \mathbf{L}^{-1}\mathbf{F} u = L − 1 F Typically we have to use some sort of linear solver in order to do the inverse.
There are traditional inversion schemes, , e.g. Jacobi, Successive OverRelaxation (SOR), Optimized SOR, etc, but these methods can be extremely expensive because they can take a long time to converge.
There are other iterative schemes which have magnitudes of order of speed up, e.g. Steepest Descent, Conjugate Gradient (CG), Preconditioned CG, Multi-Grid schemes, etc.
However, there is an alternative method to solving simple elliptical PDEs whereby we use the Discrete Sine Transform .
These are my notes about the formulation, the use cases and also my imagined pseudo-code.
Let’s take a field, u u u
u = u ( x ⃗ ) x ⃗ ∈ Ω ⊂ R D \begin{aligned}
u =\boldsymbol{u}(\vec{\mathbf{x}}) && && \vec{\mathbf{x}}\in\Omega\sub\mathbb{R}^D
\end{aligned} u = u ( x ) x ∈ Ω ⊂ R D and derive a time-independent solution which takes the form of a generic Elliptical PDE.
∇ 2 u = f ( x ⃗ ) x ⃗ ∈ Ω b ( x ⃗ ) = b u x ⃗ ∈ ∂ Ω \begin{aligned}
\nabla^2u &=
\boldsymbol{f}(\vec{\mathbf{x}})
&& && \vec{\mathbf{x}}\in\Omega \\
\boldsymbol{b}(\vec{\mathbf{x}})
&= b_u
&& && \vec{\mathbf{x}}\in\partial\Omega \\
\end{aligned} ∇ 2 u b ( x ) = f ( x ) = b u x ∈ Ω x ∈ ∂ Ω Technically, we can decompose the field, u u u
{ u = u 1 + u 2 x ⃗ ∈ Ω u = u 1 + u 2 x ⃗ ∈ ∂ Ω \begin{aligned}
\begin{cases}
u = u_1 + u_2 && && \vec{\mathbf{x}}\in\Omega \\
u = u_1 + u_2 && && \vec{\mathbf{x}}\in\partial\Omega
\end{cases}
\end{aligned} { u = u 1 + u 2 u = u 1 + u 2 x ∈ Ω x ∈ ∂ Ω We will propose one such decomposition here. We observe that all known boundary conditions can appear on the right hand side (RHS) of the PDE (5)
{ u = u I + u B x ⃗ ∈ Ω u = u I + u B x ⃗ ∈ ∂ Ω \begin{aligned}
\begin{cases}
u = u_\mathrm{I} + u_\mathrm{B} && && \vec{\mathbf{x}}\in\Omega \\
u = u_\mathrm{I} + u_\mathrm{B} && && \vec{\mathbf{x}}\in\partial\Omega
\end{cases}
\end{aligned} { u = u I + u B u = u I + u B x ∈ Ω x ∈ ∂ Ω We can say that the interior PDE, u I u_\mathrm{I} u I u B u_\mathrm{B} u B u B u_\mathrm{B} u B
∇ 2 u B = f B x ⃗ ∈ Ω u B = b u x ⃗ ∈ ∂ Ω \begin{aligned}
\nabla^2u_\mathrm{B} &= \boldsymbol{f}_\mathrm{B}
&& && \vec{\mathbf{x}}\in\Omega \\
u_\mathrm{B} &= b_u
&& && \vec{\mathbf{x}}\in\partial\Omega \\
\end{aligned} ∇ 2 u B u B = f B = b u x ∈ Ω x ∈ ∂ Ω where the values along the boundary ∂ Ω \partial\Omega ∂ Ω (5)
∇ 2 u I = − ∇ 2 u B + f x ⃗ ∈ Ω u I = 0 x ⃗ ∈ ∂ Ω \begin{aligned}
\nabla^2u_\mathrm{I} &= -\nabla^2 u_\mathrm{B} + \boldsymbol{f}
&& && \vec{\mathbf{x}}\in\Omega \\
u_\mathrm{I} &= 0
&& && \vec{\mathbf{x}}\in\partial\Omega \\
\end{aligned} ∇ 2 u I u I = − ∇ 2 u B + f = 0 x ∈ Ω x ∈ ∂ Ω where we see that we’re left with Dirichlet boundary conditions (BC) for the interior PDE. The extra term ∇ 2 u B \nabla^2 u_\mathrm{B} ∇ 2 u B
{ u = u I x ⃗ ∈ Ω u = u B x ⃗ ∈ ∂ Ω \begin{aligned}
\begin{cases}
u = u_\mathrm{I} && && \vec{\mathbf{x}}\in\Omega \\
u = u_\mathrm{B} && && \vec{\mathbf{x}}\in\partial\Omega
\end{cases}
\end{aligned} { u = u I u = u B x ∈ Ω x ∈ ∂ Ω whereby we simply need to solve the interior PDE given by equation (8) (5)
Now we have an elliptical PDE which is zero on all of the boundaries, i.e. Dirichlet boundary conditions.
This means we can use ultra-fast methods to solve this like the Discrete Sine Transform (DST).
This particular solver does work on linear Elliptical PDEs with Dirichlet boundary conditions.
The Laplacian operator in DST space is:
∇ ^ H 2 = 2 [ cos ( 2 π n M − 1 − 1 ) Δ x − 2 + cos ( 2 π m N − 1 − 1 ) Δ y − 2 ] \hat{\nabla}_H^2 =
2\left[
\cos\left(\frac{2\pi n}{M-1}-1\right)\Delta x^{-2} +
\cos\left(\frac{2\pi m}{N-1}-1\right)\Delta y^{-2}
\right] ∇ ^ H 2 = 2 [ cos ( M − 1 2 πn − 1 ) Δ x − 2 + cos ( N − 1 2 πm − 1 ) Δ y − 2 ] There are two DST types that we can use. The first one is the DST-I transform given by this equation
DST-I [ x ] k = ∑ n = 0 N − 1 x n sin [ π N − 1 ( n + 1 ) ( k + 1 ) ] k = 0 , … , N − 1 \begin{aligned}
\text{DST-I}[x]_k &= \sum_{n=0}^{N-1} x_n \sin\left[ \frac{\pi}{N-1}(n+1)(k+1)\right] && && k=0,\ldots, N-1
\end{aligned} DST-I [ x ] k = n = 0 ∑ N − 1 x n sin [ N − 1 π ( n + 1 ) ( k + 1 ) ] k = 0 , … , N − 1 This transformation is useful when we have an unknown, u u u r r r
DST-I [ x ] k = ∑ n = 0 N − 1 x n sin [ π N ( n + 1 2 ) ( k + 1 ) ] k = 0 , … , N − 1 \begin{aligned}
\text{DST-I}[x]_k &=
\sum_{n=0}^{N-1} x_n
\sin\left[
\frac{\pi}{N}(n+\frac{1}{2})(k+1)\right] && && k=0,\ldots, N-1
\end{aligned} DST-I [ x ] k = n = 0 ∑ N − 1 x n sin [ N π ( n + 2 1 ) ( k + 1 ) ] k = 0 , … , N − 1 This transformation is useful with u , r u,r u , r
So both the interior points and the boundary values are both satisfied.
Now, let’s say that the interior points, x ⃗ \vec{\mathbf{x}} x
{ u 1 = u x ⃗ ∈ Ω u 1 = 0 x ⃗ ∈ ∂ Ω { u 2 = 0 x ⃗ ∈ Ω u 2 = u x ⃗ ∈ ∂ Ω \begin{aligned}
\begin{cases}
u_1 = u && && \vec{\mathbf{x}}\in\Omega \\
u_1 = 0 && && \vec{\mathbf{x}}\in\partial\Omega
\end{cases} \\
\begin{cases}
u_2 = 0 && && \vec{\mathbf{x}}\in\Omega \\
u_2 = u && && \vec{\mathbf{x}}\in\partial\Omega
\end{cases}
\end{aligned} { u 1 = u u 1 = 0 x ∈ Ω x ∈ ∂ Ω { u 2 = 0 u 2 = u x ∈ Ω x ∈ ∂ Ω To rename them, let’s call PDE 1, u 1 u_1 u 1 u I u_I u I u 2 u_2 u 2 u B u_B u B
So going back to the original PDE (5)
Use Cases ¶ There are some PDEs where we define a PDE in terms of a derived variable however, we have access to some measurements of the original variable. For example, let’s say that we observe some variable, ξ, that has a relationship with a PDE field, u u u
N [ u ; θ ] ( x ⃗ ) = ξ \begin{aligned}
\mathcal{N}[u;\theta](\vec{\mathbf{x}}) &= \xi \\
\end{aligned} N [ u ; θ ] ( x ) = ξ ( α ∇ 2 − β ) u = ξ \begin{aligned}
(\alpha\nabla^2 - \beta)u &= \xi
\end{aligned} ( α ∇ 2 − β ) u = ξ where θ = { α , β } \theta=\{\alpha,\beta\} θ = { α , β }
Example: Quasi-Geostrophic Equations
We have many approximate models for sea surface height (SSH), η. One such model is called the Quasi-Geostrophic (QG) Equations.
In this model, we can directly relate SSH to the streamfunction, ψ, and subsequently the potential vorticity, q q q
u = u ( x ⃗ ) ψ = ψ ( x ⃗ ) q = q ( x ⃗ ) \begin{aligned}
u=\boldsymbol{u}(\vec{\mathbf{x}}) && &&
\psi=\boldsymbol{\psi}(\vec{\mathbf{x}})&& &&
q=\boldsymbol{q}(\vec{\mathbf{x}})
\end{aligned} u = u ( x ) ψ = ψ ( x ) q = q ( x ) where the input vector x ⃗ = [ x , y ] ⊤ \vec{\mathbf{x}}=[x,y]^\top x = [ x , y ] ⊤ x ⃗ ∈ Ω ⊂ R 2 \vec{\mathbf{x}}\in\Omega\sub\mathbb{R}^2 x ∈ Ω ⊂ R 2
q = ∇ 2 ψ − b ψ ψ = a η \begin{aligned}
q &= \nabla^2\psi - b\psi \\
\psi &= a\eta
\end{aligned} q ψ = ∇ 2 ψ − b ψ = a η where a = g f 0 , b = f 0 2 c 1 2 a=\frac{g}{f_0},b=\frac{f_0^2}{c_1^2} a = f 0 g , b = c 1 2 f 0 2
∂ q ∂ t + det J ( ψ , q ) = 0 x ⃗ ∈ Ω \begin{aligned}
\frac{\partial q}{\partial t} +
\det\boldsymbol{J}(\psi,q) &= 0 && &&
\vec{\mathbf{x}}\in\Omega \\
\end{aligned} ∂ t ∂ q + det J ( ψ , q ) = 0 x ∈ Ω To simplify the PDE, we can write everything in terms of SSH, η.
So substituting η into each of the derived variables ψ , q \psi,q ψ , q
q = a ∇ 2 η − a b η ψ = a η \begin{aligned}
q &= a\nabla^2\eta - ab\eta \\
\psi &= a\eta
\end{aligned} q ψ = a ∇ 2 η − ab η = a η Now, we can substitute these values in terms of SSH into the QG equation.
∂ t ( ∇ 2 − b ) η + a det J ( η , ∇ 2 η ) = 0 x ⃗ ∈ Ω B C [ η ] ( x ⃗ ) = b ( x ⃗ ) x ⃗ ∈ ∂ Ω \begin{aligned}
\partial_t \left(\nabla^2 - b \right)\eta &+
a\det J(\eta, \nabla^2\eta) = 0 && && \vec{\mathbf{x}}\in\Omega\\
\mathcal{BC}[\eta](\vec{\mathbf{x}}) &=
\boldsymbol{b}(\vec{\mathbf{x}}) && &&
\vec{\mathbf{x}} \in \partial\Omega
\end{aligned} ∂ t ( ∇ 2 − b ) η BC [ η ] ( x ) + a det J ( η , ∇ 2 η ) = 0 = b ( x ) x ∈ Ω x ∈ ∂ Ω However, the QG equations describe a PDE which time steps through a derived quantity of SSH, i.e. the Helmholtz operator, ( α ∇ 2 − β 2 ) (\alpha\nabla^2 - \beta^2) ( α ∇ 2 − β 2 )
Imagine we have to time step with η n \eta^n η n
∂ t ( ∇ 2 − b ) η = − a J ( η , ∇ 2 η ) \partial_t \left(\nabla^2 - b\right)\eta = -a\boldsymbol{J}(\eta,\nabla^2\eta) ∂ t ( ∇ 2 − b ) η = − a J ( η , ∇ 2 η ) but of course, now we have to solve the inverse Helmholtz operator to isolate η.
Let L H = ( ∇ 2 − β ) \mathbf{L}_\mathrm{H}=(\nabla^2 -\beta) L H = ( ∇ 2 − β ) F = − a det J ( η , ∇ 2 η ) \mathbf{F}=-a\det\boldsymbol{J}(\eta,\nabla^2\eta) F = − a det J ( η , ∇ 2 η )
L H η = F η = L H − 1 F \begin{aligned}
\mathbf{L}_\mathrm{H}\boldsymbol{\eta} &= \mathbf{F} \\
\boldsymbol{\eta} &= \mathbf{L}_\mathrm{H}^{-1}\mathbf{F}
\end{aligned} L H η η = F = L H − 1 F How we actually solve this equation is where we can use the generic algorithm that we outlined below.
Now we can step forward in time with our generic time stepper, g \boldsymbol{g} g t + 1 t+1 t + 1
η t + 1 = η t + g ( Δ t , η t , L H − 1 F ) \eta^{t+1} = \eta^t + \boldsymbol{g}(\Delta t, \eta^t, \mathbf{L}_\mathrm{H}^{-1}\mathbf{F}) η t + 1 = η t + g ( Δ t , η t , L H − 1 F ) Practical Algorithm ¶ Define field of interest
u = u ( x ) x = [ x , y ] ⊤ x ∈ Ω \begin{aligned}
u &= u(\mathbf{x}) && &&
\mathbf{x}= [x, y]^{\top}
&&\mathbf{x}\in\Omega
\end{aligned} u = u ( x ) x = [ x , y ] ⊤ x ∈ Ω domain: Domain = ...
u: Field = init_field(domain)
f: Field = init_source(domain)Solve PDE and BCs on the Full Field
N [ u ] ( x ) = f x ∈ Ω B C [ u ] ( x ) = b x ∈ ∂ Ω \begin{aligned}
\mathcal{N}[u](\mathbf{x}) &= f && && \mathbf{x}\in\Omega \\
\mathcal{BC}[u](\mathbf{x}) &= b && && \mathbf{x}\in\partial\Omega
\end{aligned} N [ u ] ( x ) BC [ u ] ( x ) = f = b x ∈ Ω x ∈ ∂ Ω u = PDESolver(f, rhs_fn, bcs_fn)Parse the field into a zero field but maintaining the boundaries
u B = 0 x ∈ Ω u B = b x ∈ ∂ Ω \begin{aligned}
u_\mathrm{B} &= 0 && && \mathbf{x}\in\Omega \\
u_\mathrm{B} &= b && && \mathbf{x}\in\partial\Omega
\end{aligned} u B u B = 0 = b x ∈ Ω x ∈ ∂ Ω u_b = zero_interior(u)
u_b = bcs_fn(u)Apply the same PDE and Boundary Conditions on the empty Field
N [ u B ] ( x ) = f x ∈ Ω B C [ u B ] ( x ) = b x ∈ ∂ Ω \begin{aligned}
\mathcal{N}[u_\mathrm{B}](\mathbf{x}) &= f && && \mathbf{x}\in\Omega \\
\mathcal{BC}[u_\mathrm{B}](\mathbf{x}) &= b && && \mathbf{x}\in\partial\Omega
\end{aligned} N [ u B ] ( x ) BC [ u B ] ( x ) = f = b x ∈ Ω x ∈ ∂ Ω f_b = PDESolver(u_b, rhs_fn, bcs_fn)Remove the Boundary Effects on the Solution
u I = u − u B x ∈ Ω \begin{aligned}
u_\mathrm{I} &= u - u_\mathrm{B} && && \mathbf{x}\in\Omega
\end{aligned} u I = u − u B x ∈ Ω f_I = f[1:-1,1:-1] - f_b[1:-1,1:-1]Do Inversion for Elliptical PDE
L u I = F u I = L − 1 F \begin{aligned}
\mathbf{Lu}_\mathrm{I} &= \mathbf{F} \\
\mathbf{u}_\mathrm{I} &= \mathbf{L}^{-1}\mathbf{F}
\end{aligned} Lu I u I = F = L − 1 F matvec_Ax: Callable = ...
u_I = LinearSolve(matvec_Ax, b=f_I)Add the boundaries to the field
u = u I x ∈ Ω u = u B x ∈ ∂ Ω \begin{aligned}
u &= u_\mathrm{I} && && \mathbf{x}\in\Omega \\
u &= u_\mathrm{B} && && \mathbf{x}\in\partial\Omega
\end{aligned} u u = u I = u B x ∈ Ω x ∈ ∂ Ω u = zeros_like(u)
u[1:-1,1:-1] = u_I
u[bc_indices] = u_BPseudoCde ¶ High Level ¶ # get derived variable
q, u_bc = load_state(...)
# define callable function
matvec_Ax: Callable = ...
# solve for variable
u: Field = LinearSolve(q, matvec_Ax)
# add boundary conditions
u: Field = boundary_conditions(u, u_bv)
# calculate RHS
rhs: Field = RHS(u, q)
# step once
u_new: Field = TimeStepper(dt, u, rhs)