Skip to content
44 changes: 44 additions & 0 deletions examples/kdv_1d/kdv_1d_IMEX.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,44 @@
using OrdinaryDiffEqSDIRK
using DispersiveShallowWater
using SummationByPartsOperators: upwind_operators, periodic_derivative_operator

###############################################################################
# Semidiscretization of the KdV equation

equations = KdVEquation1D(gravity = 9.81, D = 1.0)
initial_condition = initial_condition_convergence_test
boundary_conditions = boundary_condition_periodic

# create homogeneous mesh
coordinates_min = -50.0
coordinates_max = 50.0
N = 512
mesh = Mesh1D(coordinates_min, coordinates_max, N)

# Create solver with periodic SBP operators of accuracy order 3,
# which results in a 4th order accurate semi discretizations.
# We can set the accuracy order of the upwind operators to 3 since
# we only use central versions/combinations of the upwind operators.
D1_upwind = upwind_operators(periodic_derivative_operator;
derivative_order = 1, accuracy_order = 3,
xmin = xmin(mesh), xmax = xmax(mesh),
N = nnodes(mesh))
solver = Solver(D1_upwind)

semi = Semidiscretization(mesh, equations, initial_condition, solver,
boundary_conditions = boundary_conditions)

tspan = (0.0, 5.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()
analysis_callback = AnalysisCallback(semi; interval = 100,
extra_analysis_errors = (:conservation_error,),
extra_analysis_integrals = (waterheight_total,
waterheight))
callbacks = CallbackSet(analysis_callback, summary_callback)
saveat = range(tspan..., length = 100)

alg = KenCarp4() # use an IMEX method
sol = solve(ode, alg, abstol = 1e-7, reltol = 1e-7,
save_everystep = false, callback = callbacks, saveat = saveat)
2 changes: 1 addition & 1 deletion examples/kdv_1d/kdv_1d_basic.jl
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@ semi = Semidiscretization(mesh, equations, initial_condition, solver,
boundary_conditions = boundary_conditions)

tspan = (0.0, 5.0)
ode = semidiscretize(semi, tspan)
ode = semidiscretize(semi, tspan, split_ode = Val{false}()) # no IMEX for now

summary_callback = SummaryCallback()
analysis_callback = AnalysisCallback(semi; interval = 100,
Expand Down
2 changes: 1 addition & 1 deletion examples/kdv_1d/kdv_1d_fourier.jl
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,7 @@ semi = Semidiscretization(mesh, equations, initial_condition, solver,
boundary_conditions = boundary_conditions)

tspan = (0.0, 5.0)
ode = semidiscretize(semi, tspan)
ode = semidiscretize(semi, tspan, split_ode = Val{false}()) # no IMEX for now

summary_callback = SummaryCallback()
analysis_callback = AnalysisCallback(semi; interval = 100,
Expand Down
2 changes: 1 addition & 1 deletion examples/kdv_1d/kdv_1d_implicit.jl
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@ semi = Semidiscretization(mesh, equations, initial_condition, solver,
boundary_conditions = boundary_conditions)

tspan = (0.0, 5.0)
ode = semidiscretize(semi, tspan)
ode = semidiscretize(semi, tspan, split_ode = Val{false}())

summary_callback = SummaryCallback()
analysis_callback = AnalysisCallback(semi; interval = 100,
Expand Down
123 changes: 122 additions & 1 deletion examples/kdv_1d/kdv_1d_manufactured.jl
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@ semi = Semidiscretization(mesh, equations, initial_condition, solver,
source_terms = source_terms)

tspan = (0.0, 1.0)
ode = semidiscretize(semi, tspan)
ode = semidiscretize(semi, tspan, split_ode = Val{false}())

summary_callback = SummaryCallback()
analysis_callback = AnalysisCallback(semi; interval = 100,
Expand All @@ -46,3 +46,124 @@ saveat = range(tspan..., length = 100)
alg = Rodas5()
sol = solve(ode, alg, abstol = 1e-12, reltol = 1e-12,
save_everystep = false, callback = callbacks, saveat = saveat)

"""
using alg = KenCarp4() I get the following Benchmarks:
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could you please also check the performance for a physically relevant setup, i.e., not a manufactured solution? Manufactured setups can be unusual and behave differently from what we observe in physically relevant setups.

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here are some benchmarks using initial_condition_convergence_test and alg = KenCarp4() and alg = Tsit5().

alg = KenCarp4()
IMEX:
# @benchmark, with callbacks = nothing: 
BenchmarkTools.Trial: 23 samples with 1 evaluation per sample.
 Range (min  max):  211.648 ms  225.053 ms  ┊ GC (min  max): 0.00%  3.14%
 Time  (median):     216.464 ms               ┊ GC (median):    1.19%
 Time  (mean ± σ):   217.090 ms ±   3.438 ms  ┊ GC (mean ± σ):  1.16% ± 1.06%

       ▃        ▃  ▃           █      ▃
  ▇▁▁▁▁█▁▇▁▁▁▇▁▁█▁▁█▁▇▇▇▁▁▁▁▁▇▇█▇▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▇▁▁▁▁▁▁▇▁▁▁▁▇ ▁
  212 ms           Histogram: frequency by time          225 ms <

 Memory estimate: 27.45 MiB, allocs estimate: 49627.
#  Callback (for ncalls etc:)
──────────────────────────────────────────────────────────────────────────────────────
           DispersiveSWE                     Time                    Allocations
                                    ───────────────────────   ────────────────────────
         Tot / % measured:               225ms /  12.6%           23.3MiB /   0.2%

Section                     ncalls     time    %tot     avg     alloc    %tot      avg
──────────────────────────────────────────────────────────────────────────────────────
rhs_split_stiff!             3.51k   22.1ms   77.6%  6.29μs      672B    1.4%    0.19B
  third-order derivatives    3.51k   21.5ms   75.7%  6.14μs     0.00B    0.0%    0.00B
  ~rhs_split_stiff!~         3.51k    528μs    1.9%   150ns      672B    1.4%    0.19B
analyze solution                 2   5.75ms   20.2%  2.87ms   43.7KiB   96.4%  21.9KiB
rhs_split_nonstiff!            771    619μs    2.2%   803ns      976B    2.1%    1.27B
  hyperbolic                   771    364μs    1.3%   473ns     0.00B    0.0%    0.00B
  ~rhs_split_nonstiff!~        771    238μs    0.8%   308ns      976B    2.1%    1.27B
  source terms                 771   16.9μs    0.1%  21.9ns     0.00B    0.0%    0.00B
──────────────────────────────────────────────────────────────────────────────────────

No IMEX:
# @benchmark, with callbacks = nothing: 
BenchmarkTools.Trial: 19 samples with 1 evaluation per sample.
 Range (min  max):  250.372 ms  287.337 ms  ┊ GC (min  max): 0.00%  2.31%
 Time  (median):     266.741 ms               ┊ GC (median):    1.01%
 Time  (mean ± σ):   265.551 ms ±   9.897 ms  ┊ GC (mean ± σ):  1.04% ± 0.96%

  █  ▁ ▁    ▁     ▁   ▁▁  ▁  ▁ █   ▁ ▁▁▁▁      ▁              ▁
  █▁▁█▁█▁▁▁▁█▁▁▁▁▁█▁▁▁██▁▁█▁▁█▁█▁▁▁█▁████▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
  250 ms           Histogram: frequency by time          287 ms <

 Memory estimate: 29.41 MiB, allocs estimate: 48408.
# Callback (for ncalls etc:)
──────────────────────────────────────────────────────────────────────────────────────
           DispersiveSWE                     Time                    Allocations
                                    ───────────────────────   ────────────────────────
         Tot / % measured:               263ms /  24.3%           25.2MiB /   0.1%

Section                     ncalls     time    %tot     avg     alloc    %tot      avg
──────────────────────────────────────────────────────────────────────────────────────
rhs!                         3.33k   59.4ms   92.7%  17.8μs   1.25KiB    5.4%    0.38B
  third-order derivatives    3.33k   29.9ms   46.6%  8.96μs     0.00B    0.0%    0.00B
  hyperbolic                 3.33k   28.1ms   43.8%  8.42μs     0.00B    0.0%    0.00B
  ~rhs!~                     3.33k   1.39ms    2.2%   417ns   1.25KiB    5.4%    0.38B
  source terms               3.33k   83.6μs    0.1%  25.1ns     0.00B    0.0%    0.00B
analyze solution                 1   4.69ms    7.3%  4.69ms   21.8KiB   94.6%  21.8KiB
──────────────────────────────────────────────────────────────────────────────────────


alg = Tsit5() # (IMEX or No IMEX should make almost no difference)
IMEX:
# @benchmark, with callbacks = nothing: 
BenchmarkTools.Trial: 608 samples with 1 evaluation per sample.
 Range (min  max):  7.119 ms  25.909 ms  ┊ GC (min  max): 0.00%  68.83%
 Time  (median):     7.961 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   8.203 ms ±  1.102 ms  ┊ GC (mean ± σ):  0.68% ±  3.79%

       ▄▅█▅▂▂▁
  ▃▃▅▇▇███████▇▇▇▆▄▃▄▃▄▃▃▃▃▃▂▃▁▃▂▁▂▂▂▂▃▂▁▁▁▂▁▁▁▁▁▂▁▁▁▁▁▁▂▁▁▂ ▃
  7.12 ms        Histogram: frequency by time        12.5 ms <

 Memory estimate: 559.71 KiB, allocs estimate: 4637.

# Callback (for ncalls etc:)
──────────────────────────────────────────────────────────────────────────────────────
           DispersiveSWE                     Time                    Allocations
                                    ───────────────────────   ────────────────────────
         Tot / % measured:              22.9ms /  78.0%            624KiB /  21.3%

Section                     ncalls     time    %tot     avg     alloc    %tot      avg
──────────────────────────────────────────────────────────────────────────────────────
analyze solution                 6   13.4ms   74.8%  2.23ms    131KiB   98.8%  21.9KiB
rhs_split_nonstiff!          3.56k   2.54ms   14.2%   714ns      976B    0.7%    0.27B
  hyperbolic                 3.56k   1.51ms    8.5%   424ns     0.00B    0.0%    0.00B
  ~rhs_split_nonstiff!~      3.56k    944μs    5.3%   266ns      976B    0.7%    0.27B
  source terms               3.56k   87.2μs    0.5%  24.5ns     0.00B    0.0%    0.00B
rhs_split_stiff!             3.56k   1.95ms   10.9%   548ns      672B    0.5%    0.19B
  third-order derivatives    3.56k   1.47ms    8.2%   412ns     0.00B    0.0%    0.00B
  ~rhs_split_stiff!~         3.56k    482μs    2.7%   136ns      672B    0.5%    0.19B
──────────────────────────────────────────────────────────────────────────────────────

No IMEX:
# @benchmark, with callbacks = nothing: 
BenchmarkTools.Trial: 641 samples with 1 evaluation per sample.
 Range (min  max):  6.771 ms  41.173 ms  ┊ GC (min  max): 0.00%  0.00%
 Time  (median):     7.613 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   7.806 ms ±  1.679 ms  ┊ GC (mean ± σ):  0.69% ± 3.97%

        ▃▃▅█▅▆█▂
  ▃▃▅█▇█████████▇▇▇▅▆▅▅▅▃▄▂▃▂▃▂▁▁▂▂▂▂▁▁▂▂▁▁▂▂▂▂▂▁▂▂▁▁▁▁▁▁▁▁▂ ▃
  6.77 ms        Histogram: frequency by time        11.1 ms <

 Memory estimate: 554.96 KiB, allocs estimate: 4632.
# Callback (for ncalls etc:)
──────────────────────────────────────────────────────────────────────────────────────
           DispersiveSWE                     Time                    Allocations
                                    ───────────────────────   ────────────────────────
         Tot / % measured:              21.4ms /  79.6%            622KiB /  21.2%

Section                     ncalls     time    %tot     avg     alloc    %tot      avg
──────────────────────────────────────────────────────────────────────────────────────
analyze solution                 6   12.5ms   73.3%  2.08ms    130KiB   99.1%  21.7KiB
rhs!                         3.56k   4.55ms   26.7%  1.28μs   1.25KiB    0.9%    0.36B
  hyperbolic                 3.56k   1.55ms    9.1%   437ns     0.00B    0.0%    0.00B
  third-order derivatives    3.56k   1.48ms    8.7%   417ns     0.00B    0.0%    0.00B
  ~rhs!~                     3.56k   1.42ms    8.3%   399ns   1.25KiB    0.9%    0.36B
  source terms               3.56k   91.7μs    0.5%  25.8ns     0.00B    0.0%    0.00B
──────────────────────────────────────────────────────────────────────────────────────

For some reason IMEX seems to perform worse - doing more steps
no IMEX: (split_ode = Val{false}())
──────────────────────────────────────────────────────────────────────────────────────
DispersiveSWE Time Allocations
─────────────────────── ────────────────────────
Tot / % measured: 611ms / 34.2% 40.0MiB / 10.0%

Section ncalls time %tot avg alloc %tot avg
──────────────────────────────────────────────────────────────────────────────────────
rhs! 6.62k 196ms 93.6% 29.6μs 3.94MiB 98.4% 624B
source terms 6.62k 108ms 51.5% 16.3μs 3.94MiB 98.4% 624B
third-order derivatives 6.62k 44.1ms 21.1% 6.66μs 0.00B 0.0% 0.00B
hyperbolic 6.62k 41.1ms 19.7% 6.22μs 0.00B 0.0% 0.00B
~rhs!~ 6.62k 2.97ms 1.4% 450ns 1.25KiB 0.0% 0.19B
analyze solution 3 13.3ms 6.4% 4.45ms 64.8KiB 1.6% 21.6KiB
──────────────────────────────────────────────────────────────────────────────────────

IMEX with sources in nonstiff part:
──────────────────────────────────────────────────────────────────────────────────────
DispersiveSWE Time Allocations
─────────────────────── ────────────────────────
Tot / % measured: 2.85s / 12.3% 98.6MiB / 5.7%

Section ncalls time %tot avg alloc %tot avg
──────────────────────────────────────────────────────────────────────────────────────
rhs_split_nonstiff! 8.93k 169ms 48.3% 18.9μs 5.31MiB 93.8% 624B
source terms 8.93k 160ms 45.6% 17.9μs 5.31MiB 93.8% 624B
hyperbolic 8.93k 5.67ms 1.6% 635ns 0.00B 0.0% 0.00B
~rhs_split_nonstiff!~ 8.93k 3.44ms 1.0% 386ns 976B 0.0% 0.11B
rhs_split_stiff! 45.4k 103ms 29.4% 2.27μs 672B 0.0% 0.01B
third-order derivatives 45.4k 95.9ms 27.4% 2.11μs 0.00B 0.0% 0.00B
~rhs_split_stiff!~ 45.4k 7.07ms 2.0% 156ns 672B 0.0% 0.01B
analyze solution 15 78.2ms 22.3% 5.21ms 360KiB 6.2% 24.0KiB
──────────────────────────────────────────────────────────────────────────────────────

IMEX with sources in stiff part:
──────────────────────────────────────────────────────────────────────────────────────
DispersiveSWE Time Allocations
─────────────────────── ────────────────────────
Tot / % measured: 4.89s / 28.2% 88.7MiB / 49.0%

Section ncalls time %tot avg alloc %tot avg
──────────────────────────────────────────────────────────────────────────────────────
rhs_split_stiff! 72.2k 1.27s 92.0% 17.6μs 43.0MiB 98.9% 624B
source terms 72.2k 1.17s 85.2% 16.3μs 43.0MiB 98.9% 624B
third-order derivatives 72.2k 70.8ms 5.1% 981ns 0.00B 0.0% 0.00B
~rhs_split_stiff!~ 72.2k 22.1ms 1.6% 306ns 976B 0.0% 0.01B
analyze solution 23 102ms 7.4% 4.42ms 501KiB 1.1% 21.8KiB
rhs_split_nonstiff! 13.3k 8.86ms 0.6% 665ns 672B 0.0% 0.05B
hyperbolic 13.3k 6.82ms 0.5% 512ns 0.00B 0.0% 0.00B
~rhs_split_nonstiff!~ 13.3k 2.03ms 0.1% 153ns 672B 0.0% 0.05B
──────────────────────────────────────────────────────────────────────────────────────






No with fixed time step dt = 1e-2:



no IMEX
──────────────────────────────────────────────────────────────────────────────────────
DispersiveSWE Time Allocations
─────────────────────── ────────────────────────
Tot / % measured: 6.33s / 57.1% 0.99GiB / 1.3%

Section ncalls time %tot avg alloc %tot avg
──────────────────────────────────────────────────────────────────────────────────────
rhs! 22.0k 3.46s 96.0% 157μs 13.1MiB 99.8% 624B
third-order derivatives 22.0k 1.61s 44.6% 73.1μs 0.00B 0.0% 0.00B
hyperbolic 22.0k 1.49s 41.2% 67.5μs 0.00B 0.0% 0.00B
source terms 22.0k 355ms 9.8% 16.1μs 13.1MiB 99.8% 624B
~rhs!~ 22.0k 13.5ms 0.4% 615ns 1.25KiB 0.0% 0.06B
analyze solution 1 146ms 4.0% 146ms 21.9KiB 0.2% 21.9KiB
──────────────────────────────────────────────────────────────────────────────────────

IMEX with sources in nonstiff part:
──────────────────────────────────────────────────────────────────────────────────────
DispersiveSWE Time Allocations
─────────────────────── ────────────────────────
Tot / % measured: 4.34s / 39.4% 0.98GiB / 0.0%

Section ncalls time %tot avg alloc %tot avg
──────────────────────────────────────────────────────────────────────────────────────
rhs_split_stiff! 22.2k 1.70s 99.0% 76.3μs 672B 0.2% 0.03B
third-order derivatives 22.2k 1.69s 98.7% 76.1μs 0.00B 0.0% 0.00B
~rhs_split_stiff!~ 22.2k 5.55ms 0.3% 250ns 672B 0.2% 0.03B
rhs_split_nonstiff! 601 15.1ms 0.9% 25.0μs 367KiB 89.4% 626B
source terms 601 12.9ms 0.8% 21.5μs 366KiB 89.2% 624B
hyperbolic 601 1.41ms 0.1% 2.35μs 0.00B 0.0% 0.00B
~rhs_split_nonstiff!~ 601 713μs 0.0% 1.19μs 976B 0.2% 1.62B
analyze solution 1 1.86ms 0.1% 1.86ms 42.7KiB 10.4% 42.7KiB
──────────────────────────────────────────────────────────────────────────────────────

IMEX with sources in stiff part:
──────────────────────────────────────────────────────────────────────────────────────
DispersiveSWE Time Allocations
─────────────────────── ────────────────────────
Tot / % measured: 4.29s / 48.3% 0.99GiB / 1.3%

Section ncalls time %tot avg alloc %tot avg
──────────────────────────────────────────────────────────────────────────────────────
rhs_split_stiff! 22.2k 2.06s 99.7% 92.9μs 13.2MiB 99.8% 624B
third-order derivatives 22.2k 1.70s 82.1% 76.4μs 0.00B 0.0% 0.00B
source terms 22.2k 356ms 17.2% 16.0μs 13.2MiB 99.8% 624B
~rhs_split_stiff!~ 22.2k 8.90ms 0.4% 401ns 976B 0.0% 0.04B
analyze solution 1 3.62ms 0.2% 3.62ms 21.9KiB 0.2% 21.9KiB
rhs_split_nonstiff! 601 1.61ms 0.1% 2.68μs 672B 0.0% 1.12B
hyperbolic 601 1.16ms 0.1% 1.93μs 0.00B 0.0% 0.00B
~rhs_split_nonstiff!~ 601 448μs 0.0% 746ns 672B 0.0% 1.12B
──────────────────────────────────────────────────────────────────────────────────────


Also for dt = 1e-1 the solution for IMEX already looks bad(ish),
will for split_ode = Val{false}() it looks still good.
"""
2 changes: 1 addition & 1 deletion examples/kdv_1d/kdv_1d_narrow_stencil.jl
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ semi = Semidiscretization(mesh, equations, initial_condition, solver,
boundary_conditions = boundary_conditions)

tspan = (0.0, 5.0)
ode = semidiscretize(semi, tspan)
ode = semidiscretize(semi, tspan, split_ode = Val{false}()) # no IMEX for now

summary_callback = SummaryCallback()
analysis_callback = AnalysisCallback(semi; interval = 100,
Expand Down
5 changes: 3 additions & 2 deletions src/DispersiveShallowWater.jl
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,8 @@ using RecursiveArrayTools: ArrayPartition
using Reexport: @reexport
using Roots: AlefeldPotraShi, find_zero

using SciMLBase: SciMLBase, DiscreteCallback, ODEProblem, ODESolution
using SciMLBase: SciMLBase, DiscreteCallback, ODEProblem, ODESolution, SplitFunction,
SplitODEProblem
import SciMLBase: u_modified!

@reexport using StaticArrays: SVector
Expand Down Expand Up @@ -84,7 +85,7 @@ export prim2prim, prim2cons, cons2prim, prim2phys,
prim2nondim, nondim2prim,
waterheight_total, waterheight,
velocity, momentum, discharge,
gravity,
gravity, have_stiff_terms,
bathymetry, still_water_surface,
energy_total, entropy, lake_at_rest_error,
energy_total_modified, entropy_modified,
Expand Down
15 changes: 15 additions & 0 deletions src/equations/equations.jl
Original file line number Diff line number Diff line change
Expand Up @@ -265,6 +265,21 @@ Return the gravitational acceleration ``g`` for a given set of `equations`.
return equations.gravity
end

"""
DispersiveShallowWater.have_stiff_terms(equations)

Returns `Val{true}()` if the equations have stiff terms that benefit from
implicit time integration methods and `Val{false}()` otherwise (default).

This trait is used to determine whether to create a `SplitFunction` in
[`semidiscretize`](@ref) for IMEX time integration methods.

!!! note "Implementation details"
This function is used for internal dispatch to determine the appropriate
ODE problem formulation.
"""
have_stiff_terms(::AbstractEquations) = Val{false}()

"""
energy_total(q, equations)

Expand Down
Loading
Loading