Skip to content

Commit df51ae4

Browse files
AdelekeBankolejedbrownjeremylt
authored
SUPG Advection (#56)
* add supg advection example * add supg advection jupyter notebook example * append supg advection convenience operator * add supg advection example to contents * add supgadvection operator to docs * add supg advection to GalleryOperator * update references in docs * supg advection example * supg advection jupyter demo * work on size mismatch * convenience constructor for supgadvection * add eigenvalue tests for supgadvection example * add supg advection formulation in docs * SUPG demo and fixes for element matrix with multiple input/output fields * SUPG tweaks * update references * update problem formulation * add supgmass weak form * shape mismatch issues * remove some @show * fix style * update Plots package * add supg mass operator * clean up example * clean up * add convenience constructor for supgmass * update min and max eigenvalues in example * define wind in convenience constructors * wip: supg advection demo * test passing issues with constructors * Update src/Operator/Constructors.jl Co-authored-by: Jed Brown <[email protected]> * Update src/Operator/Constructors.jl Co-authored-by: Jed Brown <[email protected]> * remove plots in Project.toml * Update docs/src/examples/advection_supg.md Co-authored-by: Jed Brown <[email protected]> * Update docs/src/examples/advection_supg.md Co-authored-by: Jed Brown <[email protected]> * Update docs/src/public/operatorgallery.md Co-authored-by: Jed Brown <[email protected]> * Update src/Operator/Base.jl Co-authored-by: Jed Brown <[email protected]> * Update docs/src/examples/advection_supg.md Co-authored-by: Jed Brown <[email protected]> * supress mapping * supg advection demo * remove wind as keyword argument * Fix style * remove wind as keyword argument * update on jupyter demo figure * fix shape issue in Constructors.jl * update GalleryOperator * supg advection demo notebook * Update test/runtests.jl Co-authored-by: Jeremy L Thompson <[email protected]> * Update src/Operator/Constructors.jl Co-authored-by: Jeremy L Thompson <[email protected]> * Update test/runtests.jl Co-authored-by: Jeremy L Thompson <[email protected]> * replace wind as argument * GalleryOperator interface * parameter as kwargs * fix style * Update src/Operator/Constructors.jl Co-authored-by: Jeremy L Thompson <[email protected]> * Update src/Operator/Constructors.jl Co-authored-by: Jeremy L Thompson <[email protected]> * Update src/Operator/Constructors.jl Co-authored-by: Jeremy L Thompson <[email protected]> * Update src/Operator/Constructors.jl Co-authored-by: Jeremy L Thompson <[email protected]> * fix GalleryOperator interface * documentation update * Update src/Operator/Constructors.jl Co-authored-by: Jeremy L Thompson <[email protected]> * Update src/Operator/Constructors.jl Co-authored-by: Jeremy L Thompson <[email protected]> * Update src/Operator/Constructors.jl Co-authored-by: Jeremy L Thompson <[email protected]> * Update src/Operator/Constructors.jl Co-authored-by: Jeremy L Thompson <[email protected]> * Update src/Operator/Constructors.jl Co-authored-by: Jeremy L Thompson <[email protected]> * Update src/Operator/Constructors.jl Co-authored-by: Jeremy L Thompson <[email protected]> * Update src/Operator/Constructors.jl Co-authored-by: Jeremy L Thompson <[email protected]> * Update src/Operator/Constructors.jl Co-authored-by: Jeremy L Thompson <[email protected]> * Update src/Operator/Constructors.jl Co-authored-by: Jeremy L Thompson <[email protected]> * include parameters in advection constructor * Update src/Operator/Constructors.jl Co-authored-by: Jeremy L Thompson <[email protected]> * parameter fix * warning for numbercomponents * revert to previous state * revert to previous state * Update src/Operator/Constructors.jl Co-authored-by: Jeremy L Thompson <[email protected]> * update basis components and check tests * basis components and test issues * rm advection example from GalleryVectorOperator Co-authored-by: Jed Brown <[email protected]> Co-authored-by: Jeremy L Thompson <[email protected]>
1 parent 7dd4629 commit df51ae4

File tree

12 files changed

+749
-94
lines changed

12 files changed

+749
-94
lines changed

‎Project.toml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -4,8 +4,8 @@ authors = ["Jeremy L Thompson <[email protected]>"]
44
version = "0.5.0"
55

66
[deps]
7-
Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
87
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
8+
Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
99
Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
1010
SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
1111

‎docs/src/examples.md‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,7 @@ This section documents the LFAToolkit examples.
77
```@contents
88
Pages = [
99
"examples/advection.md",
10+
"examples/advection_supg.md",
1011
"examples/diffusion.md",
1112
"examples/linear_elasticity.md",
1213
"examples/hyperelasticity.md",

‎docs/src/examples/advection_supg.md‎

Lines changed: 43 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,43 @@
1+
## Streamline-Upwind/Petrov-Galerkin (SUPG) advection operator
2+
3+
This is an example of the 1D scalar SUPG advection problem.
4+
5+
### Problem formulation
6+
7+
The advection problem is given by
8+
9+
```math
10+
\dfrac{\partial u}{\partial t} + c \dfrac{\partial u}{\partial x} = 0
11+
```
12+
where ``u = u(x,t)`` and ``c`` is its associated phase speed.
13+
14+
The SUPG advection weak form is given by
15+
16+
```math
17+
\int_\Omega v \left(\dfrac{\partial u}{\partial t} \right) dV - \int_\Omega \dfrac{\partial v}{\partial x} \left(c u - c τ \left( \dfrac{\partial u}{\partial t} + c \dfrac{\partial u}{\partial x} \right) \right) dV = 0, \forall v \in V
18+
```
19+
for an appropriate test space ``V \subseteq H^1 \left( \Omega \right)`` on the domain.
20+
In this weak formulation, boundary terms have been omitted, as they are not present on the infinite grid for Local Fourier Analysis.
21+
The SUPG stabilization is controlled by the parameter ``τ``, where ``τ = 0`` gives the classical Galerkin formulation and ``τ = \dfrac{h}{2}`` gives a nodally exact solution to the steady advection equation with source when using linear elements (this can be extended to advection-diffusion with a further scaling that depends on the cell Péclet number).
22+
23+
### LFAToolkit code
24+
25+
The symbol of the continuous advection operator ``u_t + c u_x = 0`` applied to the Fourier mode ``e^{i\theta x}`` is ``i\theta``.
26+
The finite element discretization yields ``M u_t + A u = 0``, and thus we are interested in the symbol of ``-M^{-1} A``.
27+
One may compare the continuous spectrum with the discrete symbol, which is necessarily periodic on ``[-\pi, \pi)``, to understand the behavior for high wave numbers, including the high frequencies that will limit stable time steps.
28+
To understand dispersion within the resolved frequencies, we instead plot the phase speed ``\lambda/\theta``, which should be very close to ``c`` through the resolved frequencies.
29+
Here we show the SUPG advection operator on ``H^1`` Lagrange basis.
30+
31+
For understanding about SUPG in this work, see papers by Hughes TJR, Brooks AN (1979, 1982) and C.H. Whiting
32+
A multi-dimensional upwind scheme with no crosswind diffusion. In: Hughes TJR, editor. Finite element methods for convection dominated flows, AMD-vol. 34. New York: ASME, (1979), pp. 19-35.
33+
Streamline upwind/Petrov–Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier–Stokes equations. Comput Meth Appl Mech Eng, 32 (1982), pp. 199-259.
34+
Hierarchical basis for stabilized finite element methods for compressible flows. Comput. Methods Appl. Mech. Engrg. 192, (2003), pp. 5167-5185.
35+
36+
````@eval
37+
using Markdown
38+
Markdown.parse("""
39+
```julia
40+
$(read("../../../examples/ex011_advection_supg.jl", String))
41+
```
42+
""")
43+
````

‎docs/src/private/operator.md‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -12,4 +12,6 @@ LFAToolkit.getnodecoordinatedifferences(::Operator)
1212
LFAToolkit.massoperator
1313
LFAToolkit.diffusionoperator
1414
LFAToolkit.advectionoperator
15+
LFAToolkit.supgadvectionoperator
16+
LFAToolkit.supgmassoperator
1517
```

‎docs/src/public/operatorgallery.md‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,7 @@
22

33
The following PDEs are available in the operator gallery.
44

5-
### Scalar Mass, Diffusion, and Advection Operators
5+
### Scalar Mass, Diffusion, Advection, and SUPG Advection Operators
66

77
```@docs
88
GalleryOperator

‎docs/src/references.md‎

Lines changed: 7 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -2,10 +2,16 @@
22

33
A. Brandt, *Multi-level adaptive solutions to boundary-value problems*, Math. Comp., 31(138) (1977), pp. 33-390.
44

5+
A. N. Brooks and T. J. R. Hughes *Streamline upwind/Petrov–Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier–Stokes equations*, Comput Meth Appl Mech Eng, 32 (1982), pp. 199-259.
6+
57
J. Brown, *Efficient nonlinear solvers for nodal high-order finite elements in 3D*, Journal of Scientific Computing, 45 (2010), pp. 48-63.
68

79
M. Gutknecht and S. Röllin, *The Chebyshev iteration revisited*, Parallel Computing, 28 (2002), pp. 263-283.
810

911
N. Hale and L. N. Trefethen, [New quadrature formulas from conformal maps](https://doi.org/10.1137/07068607X), SIAM Journal on Numerical Analysis. Vol. 46, No. 2, (2008), pp. 930-948.
1012

11-
T. Melvin, A. Staniforth and J. Thuburn, *Dispersion analysis of the spectral element method*, Q.J.R. Meteorol. Soc. 138, (2012), pp. 1934-1947.
13+
T. J. R. Hughes and A. N. Brooks *A multi-dimensional upwind scheme with no crosswind diffusion*, In: Hughes TJR, editor. Finite element methods for convection dominated flows, AMD-vol. 34. New York: ASME, (1979), pp. 19-35.
14+
15+
T. Melvin, A. Staniforth and J. Thuburn, *Dispersion analysis of the spectral element method*, Q.J.R. Meteorol. Soc. 138, (2012), pp. 1934-1947.
16+
17+
C. H. Whiting, K. E. Jansen and S. Dey, *Hierarchical basis for stabilized finite element methods for compressible flows*, Comput. Methods Appl. Mech. Engrg. 192, (2003), pp. 5167-5185.

‎examples/ex011_advection_supg.jl‎

Lines changed: 81 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,81 @@
1+
# ---------------------------------------------------------------
2+
# SUPG advection example following Melvin, Staniforth and Thuburn
3+
# Q.J:R. Meteorol. Soc. 138: 1934-1947, Oct (2012)
4+
# non polynomial advection example using transplanted (transformed)
5+
# quadrature formulas from conformal maps following Hale and
6+
# Trefethen (2008) SIAM J. NUMER. ANAL. Vol. 46, No. 2
7+
# In this example, we can call other mapping options available
8+
# i.e, sausage_transformation(9), kosloff_talezer_transformation(0.98)
9+
# otherwise, run without mapping
10+
# ---------------------------------------------------------------
11+
using LFAToolkit
12+
using LinearAlgebra
13+
14+
# setup
15+
mesh = Mesh1D(1.0)
16+
P = 2;
17+
Q = P;
18+
collocate = false
19+
mapping = nothing
20+
basis =
21+
TensorH1LagrangeBasis(P, Q, 1, 1, collocatedquadrature = collocate, mapping = mapping)
22+
23+
# frequency set up
24+
numbersteps = 100
25+
θ_min = 0
26+
θ_max = (P - 1) * π
27+
θ = LinRange(θ_min, θ_max, numbersteps)
28+
29+
# associated phase speed
30+
c = 1.0
31+
32+
# Tau scaling for SUPG
33+
τ = 0.5 / (P - 1) # 0 returns Galerkin method; empirically, 0.25 is too big for P=3 (quadratic)
34+
35+
# weak form for SUPG advection
36+
function supgadvectionweakform(U::Matrix{Float64}, w::Array{Float64})
37+
u = U[1, :]
38+
du = U[2, :]
39+
dv = (c * u - c * τ * (c * du)) * w[1]
40+
return [dv]
41+
end
42+
43+
# supg advection operator
44+
inputs = [
45+
OperatorField(
46+
basis,
47+
[EvaluationMode.interpolation, EvaluationMode.gradient],
48+
"advected field",
49+
),
50+
OperatorField(basis, [EvaluationMode.quadratureweights], "quadrature weights"),
51+
]
52+
outputs = [OperatorField(basis, [EvaluationMode.gradient])]
53+
54+
supgadvection = Operator(supgadvectionweakform, mesh, inputs, outputs)
55+
56+
# supg mass operator
57+
function supgmassweakform(udot::Array{Float64}, w::Array{Float64})
58+
v = udot * w[1]
59+
dv = c * τ * udot * w[1]
60+
return ([v; dv],)
61+
end
62+
supgmass = Operator(
63+
supgmassweakform,
64+
mesh,
65+
[
66+
OperatorField(basis, [EvaluationMode.interpolation], "u_t"),
67+
OperatorField(basis, [EvaluationMode.quadratureweights], "quadrature weights"),
68+
],
69+
[OperatorField(basis, [EvaluationMode.interpolation, EvaluationMode.gradient])],
70+
)
71+
72+
# compute operator symbols
73+
function advection_supg_symbol(θ)
74+
A = computesymbols(supgadvection, [θ]) * 2 # transform from reference to physical on dx=1 grid
75+
M = computesymbols(supgmass, [θ]) # mass matrix
76+
return sort(imag.(eigvals(-M \ A)))
77+
end
78+
79+
eigenvalues = hcat(advection_supg_symbol.(θ)...)'
80+
81+
# ---------------------------------------------------------------

‎examples/jupyter/demo011_advection_supg.ipynb‎

Lines changed: 158 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,158 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "markdown",
5+
"metadata": {},
6+
"source": [
7+
"# SUPG Advection Example: Dispersion and Phase Speed Diagrams"
8+
]
9+
},
10+
{
11+
"cell_type": "code",
12+
"execution_count": null,
13+
"metadata": {},
14+
"outputs": [],
15+
"source": [
16+
"# dependencies\n",
17+
"using Pkg\n",
18+
"Pkg.activate(\".\")\n",
19+
"Pkg.instantiate()\n",
20+
"Pkg.develop(path=\"../..\")\n",
21+
"using Plots\n",
22+
"using LFAToolkit\n",
23+
"using LinearAlgebra "
24+
]
25+
},
26+
{
27+
"cell_type": "code",
28+
"execution_count": null,
29+
"metadata": {},
30+
"outputs": [],
31+
"source": [
32+
"# setup\n",
33+
"p = 2\n",
34+
"q = p\n",
35+
"dimension = 1\n",
36+
"mapping = nothing\n",
37+
"collocate = false\n",
38+
"basis = TensorH1LagrangeBasis(p, q, 1, 1, collocatedquadrature = collocate, mapping = mapping)\n",
39+
"\n",
40+
"mesh = []\n",
41+
"if dimension == 1\n",
42+
" mesh = Mesh1D(1.0)\n",
43+
"elseif dimension == 2\n",
44+
" mesh = Mesh2D(1.0, 1.0)\n",
45+
"end\n",
46+
"\n",
47+
"# advective speed\n",
48+
"c = 1.0\n",
49+
"# Tau scaling for SUPG, 0 returns Galerkin method\n",
50+
"τ = 0.5 / (p - 1)"
51+
]
52+
},
53+
{
54+
"cell_type": "code",
55+
"execution_count": null,
56+
"metadata": {},
57+
"outputs": [],
58+
"source": [
59+
"# SUPG advection weakform \n",
60+
"function supgadvectionweakform(U::Matrix{Float64}, w::Array{Float64})\n",
61+
" u = U[1, :]\n",
62+
" du = U[2, :]\n",
63+
" dv = (c * u - c * τ * (c * du)) * w[1]\n",
64+
" return [dv]\n",
65+
"end\n",
66+
"\n",
67+
"# SUPG advection operator\n",
68+
"inputs = [\n",
69+
" OperatorField(\n",
70+
" basis,\n",
71+
" [EvaluationMode.interpolation, EvaluationMode.gradient],\n",
72+
" \"advected field\",\n",
73+
" ),\n",
74+
" OperatorField(basis, [EvaluationMode.quadratureweights], \"quadrature weights\"),\n",
75+
"]\n",
76+
"outputs = [OperatorField(basis, [EvaluationMode.gradient])]\n",
77+
"supgadvection = Operator(supgadvectionweakform, mesh, inputs, outputs)\n",
78+
"\n",
79+
"# SUPG mass weakform and mass operator\n",
80+
"function supgmassweakform(udot::Array{Float64}, w::Array{Float64})\n",
81+
" v = udot * w[1]\n",
82+
" dv = c * τ * udot * w[1]\n",
83+
" return ([v; dv],)\n",
84+
"end\n",
85+
"supgmass = Operator(\n",
86+
" supgmassweakform,\n",
87+
" mesh,\n",
88+
" [\n",
89+
" OperatorField(basis, [EvaluationMode.interpolation], \"u_t\"),\n",
90+
" OperatorField(basis, [EvaluationMode.quadratureweights], \"quadrature weights\"),\n",
91+
" ],\n",
92+
" [OperatorField(basis, [EvaluationMode.interpolation, EvaluationMode.gradient])],\n",
93+
")"
94+
]
95+
},
96+
{
97+
"cell_type": "code",
98+
"execution_count": null,
99+
"metadata": {},
100+
"outputs": [],
101+
"source": [
102+
"# compute full operator symbols\n",
103+
"numbersteps = 100\n",
104+
"maxeigenvalue = 0\n",
105+
"θ_min = 0\n",
106+
"θ_max = (p-1) * π\n",
107+
"θ_range = LinRange(θ_min, θ_max, numbersteps)\n",
108+
"\n",
109+
"# compute and plot dispersion and phase speed\n",
110+
"# -- 1D --\n",
111+
"function advection_supg_symbol(θ_range)\n",
112+
" A = computesymbols(supgadvection, [θ_range]) * 2 # transform from reference to physical on dx=1 grid\n",
113+
" M = computesymbols(supgmass, [θ_range]) # mass matrix\n",
114+
" return sort(imag.(eigvals(-M \\ A)))\n",
115+
"end\n",
116+
"\n",
117+
"eigenvalues = hcat(advection_supg_symbol.(θ_range)...)'\n",
118+
"\n",
119+
"plot1 = plot(θ_range / π, eigenvalues ./ π, linewidth = 3) # Dispersion diagram\n",
120+
" plot!(\n",
121+
" identity,\n",
122+
" xlabel = \"θ/π\",\n",
123+
" ylabel = \"Eigenvalues\",\n",
124+
" label = \"exact\",\n",
125+
" legend = :none,\n",
126+
" color = :black,\n",
127+
" title = \"Dispersion relation, p=$p, collocate=$collocate, τ=$τ\",\n",
128+
")\n",
129+
"plot2 = plot(θ_range / π, eigenvalues ./ θ_range, linewidth = 3) # Phase speed diagram\n",
130+
" plot!(\n",
131+
" one,\n",
132+
" xlabel = \"θ/π\",\n",
133+
" ylabel = \"Eigvalues/θ\",\n",
134+
" legend = :none,\n",
135+
" color = :black,\n",
136+
" title = \"Phase speed, p=$p, collocate=$collocate, τ=$τ\",\n",
137+
")\n",
138+
"plot!(plot1, plot2, layout = (1, 2), size = (1050, 450))"
139+
]
140+
}
141+
],
142+
"metadata": {
143+
"kernelspec": {
144+
"display_name": "Julia 1.7.0",
145+
"language": "julia",
146+
"name": "julia-1.7"
147+
},
148+
"language_info": {
149+
"file_extension": ".jl",
150+
"mimetype": "application/julia",
151+
"name": "julia",
152+
"version": "1.7.0"
153+
},
154+
"orig_nbformat": 4
155+
},
156+
"nbformat": 4,
157+
"nbformat_minor": 2
158+
}

0 commit comments

Comments
 (0)