Skip to content
Draft
Show file tree
Hide file tree
Changes from 9 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
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
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
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
4 changes: 2 additions & 2 deletions src/DispersiveShallowWater.jl
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,7 @@ using RecursiveArrayTools: ArrayPartition
using Reexport: @reexport
using Roots: AlefeldPotraShi, find_zero

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

@reexport using StaticArrays: SVector
Expand Down Expand Up @@ -74,7 +74,7 @@ export LinearDispersionRelation, EulerEquations1D, wave_speed
export prim2prim, prim2cons, cons2prim, prim2phys,
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