Opportunities for improving efficiency

[ranocha]

This may not be crucial right now, but we should think about the following options:

In code like

Ariadne.jl/libs/Theseus/src/imex/imex.jl

Lines 30 to 39 in 8f0fff3

	function (::RKIMEX{N})(res, uₙ, Δt, f1!, du, du_tmp, u, p, t, stages_ex, stages_im, stage, RK) where {N}
	@. res = u
	for j in 1:(stage - 1)
	@. res = res - RK.a_ex[stage, j] * stages_ex[j] - RK.a_im[stage, j] * stages_im[j]
	end
	@. du = u * Δt + uₙ
	f1!(du_tmp, du, p, t + RK.c_im[stage] * Δt)
	@. res = res - RK.a_im[stage, stage] * du_tmp
	return res
	end

, adding the previous equations could be done once before calling the nonlinear solver.

We could use muladd for many of the updates, e.g.,

Ariadne.jl/libs/Theseus/src/imex/imex.jl

Lines 30 to 39 in 8f0fff3

	function (::RKIMEX{N})(res, uₙ, Δt, f1!, du, du_tmp, u, p, t, stages_ex, stages_im, stage, RK) where {N}
	@. res = u
	for j in 1:(stage - 1)
	@. res = res - RK.a_ex[stage, j] * stages_ex[j] - RK.a_im[stage, j] * stages_im[j]
	end
	@. du = u * Δt + uₙ
	f1!(du_tmp, du, p, t + RK.c_im[stage] * Δt)
	@. res = res - RK.a_im[stage, stage] * du_tmp
	return res
	end

.

Does something like

Ariadne.jl/libs/Theseus/src/imex/imex.jl

Line 53 in 8f0fff3

krylov_tol_abs::Float64

impact the use of other floating point types?

Is it better to use

Ariadne.jl/libs/Theseus/src/imex/imex.jl

Line 300 in 8f0fff3

fill!(integrator.u_tmp, zero(eltype(integrator.u_tmp)))

or

Ariadne.jl/libs/Theseus/src/imex/imex.jl

Line 218 in 8f0fff3

@. integrator.u_tmp = 0

, @vchuravy?

[MarcoArtiano]

Regarding the reset function, I was testing the following:

julia> b1 = @benchmark fill!($integrator.u_tmp, zero(eltype($integrator.u_tmp)))

BenchmarkTools.Trial: 10000 samples with 233 evaluations per sample.
 Range (min … max):  319.661 ns … 701.545 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     342.195 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   347.646 ns ±  19.659 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

        ▁▆▇▃▁  ▇█▆▃  ▆█▆▃▁ ▄▆▄▂   ▂       ▁▁    ▂▂▂   ▁▂▂▁      ▃
  ▇█▆▁▁▁█████▅████████████▇█████▆███▇▅▄▅▃▇███▇▅█████▇██████▆▇██ █
  320 ns        Histogram: log(frequency) by time        404 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> b2 = @benchmark @. $integrator.u_tmp = 0
BenchmarkTools.Trial: 10000 samples with 234 evaluations per sample.
 Range (min … max):  318.380 ns … 767.269 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     346.209 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   346.453 ns ±  19.462 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

         ▃      ▇▁    █▃    ▁▂                                   
  ▅▃▄▂▂▁▂██▆▃▂▂███▄▂▂▄██▄▂▂▂██▄▂▂▂▃▅▃▂▂▂▂▃▃▃▃▂▂▂▂▃▃▂▂▂▂▂▂▂▂▂▂▂▂ ▃
  318 ns           Histogram: frequency by time          403 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> b3 = @benchmark $integrator.u_tmp .= 0
BenchmarkTools.Trial: 10000 samples with 233 evaluations per sample.
 Range (min … max):  320.004 ns … 556.077 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     342.322 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   347.321 ns ±  17.328 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▁▁    ▄▇▆▄▂ ▃██▆▃ ▂█▇▆▃  ▅▆▄▁   ▂▁            ▂▂▂▁  ▁▂▂▁   ▁▁ ▃
  ███▄▃▁█████▆█████▇██████▇█████▆███▅▅▅▃▆▇███▇▇▇███████████████ █
  320 ns        Histogram: log(frequency) by time        404 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> b4 = @benchmark @. $integrator.u_tmp = zero(eltype($integrator.u_tmp))
BenchmarkTools.Trial: 10000 samples with 233 evaluations per sample.
 Range (min … max):  319.103 ns … 788.451 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     346.751 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   349.751 ns ±  24.011 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▃ ▄  ▅▆▇▅ ▆██▅ ▆█▇▃▁▅▇▅▂ ▃▄▂▁  ▁▁▂▁ ▁▂▂▂▁▁▃▃▁▁▁▁▁             ▃
  ████▁████████████████████████████████████████████████▆▇▅▅▅▆▅▅ █
  319 ns        Histogram: log(frequency) by time        427 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> b5 = @benchmark $integrator.u_tmp .= zero(eltype($integrator.u_tmp))

       # Visualizza i risultati
BenchmarkTools.Trial: 10000 samples with 234 evaluations per sample.
 Range (min … max):  318.338 ns … 656.547 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     342.594 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   346.715 ns ±  20.829 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▄▃▂   ▅▇▆▃  ▆█▇▄  ▇█▆▃  ▆▆▃▁  ▄▃     ▁▁▁   ▂▂▁  ▁▂▂▁   ▁▁     ▃
  ███▇▅▁█████▅█████▇█████▇████▇▇███▇▄▅▇███▇▆█████▇█████▇▇███▆▅▇ █
  318 ns        Histogram: log(frequency) by time        410 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

but when I run the code using TrixiBase timer, then I get that the last two seem to be the best choice:

filling zeros 3                 3.20k   2.03ms    0.0%   635ns     0.00B    0.0%    0.00B
filling zeros 1                 3.20k   2.00ms    0.0%   624ns     0.00B    0.0%    0.00B
filling zeros 2                 3.20k   1.52ms    0.0%   476ns     0.00B    0.0%    0.00B
filling zeros 5                 3.20k   1.43ms    0.0%   448ns     0.00B    0.0%    0.00B
filling zeros 4                 3.20k   1.41ms    0.0%   439ns     0.00B    0.0%    0.00B

where

   @trixi_timeit timer() "filling zeros 1"  fill!(integrator.u_tmp, zero(eltype(integrator.u_tmp)))
   @trixi_timeit timer() "filling zeros 2"  @. integrator.u_tmp = 0 
   @trixi_timeit timer() "filling zeros 3"   integrator.u_tmp .= 0 
   @trixi_timeit timer() "filling zeros 4"  @. integrator.u_tmp = zero(eltype(integrator.u_tmp)) 
   @trixi_timeit timer() "filling zeros 5"   integrator.u_tmp .= zero(eltype(integrator.u_tmp)) 

I'm not sure whether this is meaningful or not.

[vchuravy]

Eh, you might get order effects here. Try reordering your statements.
Personally I prefere fill! since it requires less of the compiler (broadcast is more code.)