Skip to content

Commit 90bfa8e

Browse files
GodotMisogiGodotMisogi
GodotMisogi
authored and
GodotMisogi
committed
Fixed typo + tested KenCarp and ARKODE solvers in KS pseudo-spectral
1 parent 8e03008 commit 90bfa8e

File tree

1 file changed

+85
-110
lines changed

1 file changed

+85
-110
lines changed

‎benchmarks/SimpleHandwrittenPDE/ks_spectral_wpd.jmd‎

Lines changed: 85 additions & 110 deletions
Original file line numberDiff line numberDiff line change
@@ -5,18 +5,18 @@ author: Arjit Seth
55

66
## Problem Description
77

8-
The Kuramoto-Sivashinsky partial differential equation is solved on the domain $[-L, L] \times [0, T] \in \mathbb R \times R,~L = 16,~T = 30$, with the following initial and periodic boundary conditions:
8+
The Kuramoto-Sivashinsky partial differential equation is solved on the domain $[-L, L] \times [0, T] \in \mathbb R \times R,~L = 16,~T = 10$, with the following initial and periodic boundary conditions:
99
```math
1010
\begin{align}
11-
\partial_t u & = -\kappa \partial_x \left[u(t,x) \partial_x u(t,x)\right] - \beta \partial_x^2 u(t,x) - \gamma \partial_x^4 u(t,x) \\
11+
\partial_t u & = -\kappa u(t,x) \partial_x u(t,x) - \beta \partial_x^2 u(t,x) - \gamma \partial_x^4 u(t,x) \\
1212
u(0,x) & = \cos\left(\frac{2\pi}{L} x\right), \\
1313
u(t,-L) & = u(t,L) = 1.
1414
\end{align}
1515
```
1616
The spatial derivative operators are represented via Fourier pseudospectral approximations. Here, the domain is discretized by projecting on an equispaced grid of points $x_s \in [-L, L]$; the solution is approximated on this grid via linear combinations of sinusoidal functions in space. The coefficients $\kappa = 1, \beta = 1/2,~\gamma = 1/8$ are chosen to produce `interesting' behavior as seen in the reference solution below.
1717
```math
1818
\begin{align}
19-
\frac{du}{dt} & = -\kappa D_x (u D_x u) - \beta D_x^2 u - \gamma D_x^4 u \\
19+
\frac{du}{dt} & = -\kappa u D_x u - \beta D_x^2 u - \gamma D_x^4 u \\
2020
u(0,x_s) & = \cos\left(\frac{2\pi}{L} x_s\right), \\
2121
u(t,-L) & = u(t,L) = 1.
2222
\end{align}
@@ -65,7 +65,7 @@ Using an adaptive timestepping method to solve the system of ordinary differenti
6565
```julia
6666
N = 128 # Number of grid points
6767
L = 16.0 # Domain length
68-
alpha = 30.0 # Time scaling factor
68+
alpha = 10.0 # Time scaling factor
6969
xs, prob = kuramoto_sivashinsky(N, L, alpha)
7070

7171
sol = solve(prob, Rodas5(autodiff=false); dt = 1e-4, reltol = 1e-12, abstol = 1e-12)
@@ -79,119 +79,13 @@ heatmap(xs, tslices, ys', xlabel = "x", ylabel = "t")
7979

8080
## Work-Precision Diagrams
8181

82-
### Low Tolerances
83-
84-
#### Implicit-Explicit Methods
85-
86-
Dense/banded linear solvers.
87-
88-
```julia
89-
abstols = 0.1 .^ (7:13)
90-
reltols = 0.1 .^ (4:10)
91-
setups = [
92-
# KenCarp methods take forever with adaptive timestepping for some reason
93-
# Dict(:alg => KenCarp3()),
94-
# Dict(:alg => KenCarp4()),
95-
# Dict(:alg => KenCarp5()),
96-
# ARKODE with banded linear solves fail with adaptive timestepping for some reason
97-
# Dict(:alg => ARKODE(Sundials.Implicit(), order=3, linear_solver=:Band, jac_upper=1, jac_lower=1)),
98-
# Dict(:alg => ARKODE(Sundials.Implicit(), order=4, linear_solver=:Band, jac_upper=1, jac_lower=1)),
99-
# Dict(:alg => ARKODE(Sundials.Implicit(), order=5, linear_solver=:Band, jac_upper=1, jac_lower=1))
100-
]
101-
labels = hcat(
102-
# "KenCarp3",
103-
# "KenCarp4",
104-
# "KenCarp5",
105-
# "ARKODE3",
106-
# "ARKODE4",
107-
# "ARKODE5",
108-
)
109-
@time wp = WorkPrecisionSet(prob, abstols, reltols, setups;
110-
print_names=true, names=labels, numruns=5, error_estimate=:l2,
111-
save_everystep=false, appxsol=test_sol, maxiters=Int(1e5));
112-
113-
plot(wp, label=labels, markershape=:auto, title="IMEX Methods, Band Linsolve, Low Tolerances")
114-
```
115-
116-
Krylov linear solvers.
117-
```julia
118-
abstols = 0.1 .^ (7:13)
119-
reltols = 0.1 .^ (4:10)
120-
setups = [
121-
# KenCarp methods take forever with adaptive timestepping for some reason
122-
# Dict(:alg => KenCarp3(linsolve=KrylovJL_GMRES())),
123-
# Dict(:alg => KenCarp4(linsolve=KrylovJL_GMRES())),
124-
# Dict(:alg => KenCarp5(linsolve=KrylovJL_GMRES())),
125-
Dict(:alg => ARKODE(Sundials.Implicit(), order=3, linear_solver=:GMRES)),
126-
Dict(:alg => ARKODE(Sundials.Implicit(), order=4, linear_solver=:GMRES)),
127-
Dict(:alg => ARKODE(Sundials.Implicit(), order=5, linear_solver=:GMRES)),
128-
]
129-
labels = hcat(
130-
# "KenCarp3",
131-
# "KenCarp4",
132-
# "KenCarp5",
133-
"ARKODE3",
134-
"ARKODE4",
135-
"ARKODE5",
136-
)
137-
@time wp = WorkPrecisionSet(prob, abstols, reltols, setups;
138-
print_names=true, names=labels, numruns=5, error_estimate=:l2,
139-
save_everystep=false, appxsol=test_sol, maxiters=Int(1e5));
140-
141-
plot(wp, label=labels, markershape=:auto, title="IMEX Methods, Krylov Linsolve, Low Tolerances")
142-
```
143-
144-
#### Exponential Integrators
145-
146-
```julia
147-
abstols = 0.1 .^ (7:11) # all fixed dt methods so these don't matter much
148-
reltols = 0.1 .^ (4:8)
149-
multipliers = 0.5 .^ (0:4)
150-
setups = [
151-
Dict(:alg => ETDRK3(), :dts => 1e-2 * multipliers),
152-
Dict(:alg => ETDRK4(), :dts => 1e-2 * multipliers),
153-
Dict(:alg => HochOst4(), :dts => 1e-2 * multipliers),
154-
]
155-
labels = hcat(
156-
"ETDRK3 (caching)",
157-
"ETDRK4 (caching)",
158-
"HochOst4 (caching)",
159-
)
160-
@time wp = WorkPrecisionSet(prob, abstols, reltols, setups;
161-
print_names=true, names=labels, numruns=5, error_estimate=:l2,
162-
save_everystep=false, appxsol=test_sol, maxiters=Int(1e5));
163-
164-
plot(wp, label=labels, markershape=:auto, title="ExpRK Methods, Low Tolerances")
165-
```
166-
167-
#### Comparisons Between Families
168-
169-
```julia
170-
abstols = 0.1 .^ (7:11)
171-
reltols = 0.1 .^ (4:8)
172-
multipliers = 0.5 .^ (0:4)
173-
setups = [
174-
Dict(:alg => ARKODE(Sundials.Implicit(), order=5, linear_solver=:GMRES)),
175-
Dict(:alg => ETDRK3(), :dts => 1e-2 * multipliers),
176-
Dict(:alg => ETDRK4(), :dts => 1e-2 * multipliers),
177-
]
178-
labels = hcat("ARKODE (Krylov linsolve)", "ETDRK3 ()", "ETDRK4 ()")
179-
@time wp = WorkPrecisionSet(prob, abstols, reltols, setups;
180-
print_names=true, names=labels, numruns=5, error_estimate=:l2,
181-
save_everystep=false, appxsol=test_sol, maxiters=Int(1e5));
182-
183-
plot(wp, label=labels, markershape=:auto, title="Between Families, Low Tolerances")
184-
```
185-
186-
18782
### High Tolerances
18883

18984
#### Implicit-Explicit Methods
19085

19186
```julia
19287
abstols = 0.1 .^ (5:8) # all fixed dt methods so these don't matter much
19388
reltols = 0.1 .^ (1:4)
194-
multipliers = 0.5 .^ (0:3)
19589
setups = [
19690
Dict(:alg => IMEXEuler(), :dts => 1e-4 * multipliers),
19791
Dict(:alg => CNAB2(), :dts => 1e-4 * multipliers),
@@ -264,6 +158,87 @@ labels = hcat(
264158
plot(wp, label=labels, markershape=:auto, title="Between Families, High Tolerances")
265159
```
266160

161+
162+
### Low Tolerances
163+
164+
#### Implicit-Explicit Methods
165+
166+
```julia
167+
abstols = 0.1 .^ (7:10)
168+
reltols = 0.1 .^ (4:7)
169+
setups = [
170+
Dict(:alg => KenCarp3()),
171+
Dict(:alg => KenCarp4()),
172+
Dict(:alg => KenCarp5()),
173+
Dict(:alg => KenCarp3(linsolve=KrylovJL_GMRES())),
174+
Dict(:alg => KenCarp4(linsolve=KrylovJL_GMRES())),
175+
Dict(:alg => KenCarp5(linsolve=KrylovJL_GMRES())),
176+
# ARKODE with banded linear solvers fail.
177+
Dict(:alg => ARKODE(Sundials.Implicit(), order=3, linear_solver=:GMRES)),
178+
Dict(:alg => ARKODE(Sundials.Implicit(), order=4, linear_solver=:GMRES)),
179+
Dict(:alg => ARKODE(Sundials.Implicit(), order=5, linear_solver=:GMRES)),
180+
]
181+
labels = hcat(
182+
"KenCarp3",
183+
"KenCarp4",
184+
"KenCarp5",
185+
"KenCarp3 (Krylov)",
186+
"KenCarp4 (Krylov)",
187+
"KenCarp5 (Krylov)",
188+
"ARKODE3 (Krylov)",
189+
"ARKODE4 (Krylov)",
190+
"ARKODE5 (Krylov)",
191+
)
192+
@time wp = WorkPrecisionSet(prob, abstols, reltols, setups;
193+
print_names=true, names=labels, numruns=5, error_estimate=:l2,
194+
save_everystep=false, appxsol=test_sol, maxiters=Int(1e5));
195+
196+
plot(wp, label=labels, markershape=:auto, title="IMEX Methods, Krylov Linsolve, Low Tolerances")
197+
```
198+
199+
#### Exponential Integrators
200+
201+
```julia
202+
abstols = 0.1 .^ (7:11) # all fixed dt methods so these don't matter much
203+
reltols = 0.1 .^ (4:8)
204+
multipliers = 0.5 .^ (0:4)
205+
setups = [
206+
Dict(:alg => ETDRK3(), :dts => 1e-2 * multipliers),
207+
Dict(:alg => ETDRK4(), :dts => 1e-2 * multipliers),
208+
Dict(:alg => HochOst4(), :dts => 1e-2 * multipliers),
209+
]
210+
labels = hcat(
211+
"ETDRK3 (caching)",
212+
"ETDRK4 (caching)",
213+
"HochOst4 (caching)",
214+
)
215+
@time wp = WorkPrecisionSet(prob, abstols, reltols, setups;
216+
print_names=true, names=labels, numruns=5, error_estimate=:l2,
217+
save_everystep=false, appxsol=test_sol, maxiters=Int(1e5));
218+
219+
plot(wp, label=labels, markershape=:auto, title="ExpRK Methods, Low Tolerances")
220+
```
221+
222+
#### Comparisons Between Families
223+
224+
```julia
225+
abstols = 0.1 .^ (7:11)
226+
reltols = 0.1 .^ (4:8)
227+
multipliers = 0.5 .^ (0:4)
228+
setups = [
229+
Dict(:alg => ARKODE(Sundials.Implicit(), order=5, linear_solver=:GMRES)),
230+
Dict(:alg => ETDRK3(), :dts => 1e-2 * multipliers),
231+
Dict(:alg => ETDRK4(), :dts => 1e-2 * multipliers),
232+
]
233+
labels = hcat("ARKODE (Krylov linsolve)", "ETDRK3 ()", "ETDRK4 ()")
234+
@time wp = WorkPrecisionSet(prob, abstols, reltols, setups;
235+
print_names=true, names=labels, numruns=5, error_estimate=:l2,
236+
save_everystep=false, appxsol=test_sol, maxiters=Int(1e5));
237+
238+
plot(wp, label=labels, markershape=:auto, title="Between Families, Low Tolerances")
239+
```
240+
241+
267242
```julia, echo = false
268243
using SciMLBenchmarks
269244
SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])

0 commit comments

Comments
 (0)