Mercurial > repos > public > sbplib_julia

changeset 876:4f3924293894 laplace_benchmarks

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add examples and benchmarks folders
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 27 Jan 2022 11:00:31 +0100
parents 067a322e4f73
children 0a856fb96db4
files benchmarks/laplace_benchmark.jl benchmarks/laplace_benchmark.m examples/wave_eq.jl laplace_benchmark.jl laplace_benchmark.m wave.gif wave_eq.jl
diffstat 7 files changed, 123 insertions(+), 154 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/benchmarks/laplace_benchmark.jl	Thu Jan 27 11:00:31 2022 +0100
@@ -0,0 +1,42 @@
+using Sbplib.Grids, Sbplib.SbpOperators, Sbplib.LazyTensors, Sbplib.RegionIndices
+using Profile, BenchmarkTools
+
+function apply_laplace!(f, u, L, inds)
+    for I in inds
+        f[I] = (L*u)[I]
+    end
+end
+
+apply_laplace!(f, u, L) = apply_laplace!(f, u, L, eachindex(u))
+
+region_indices(L, N, ::Lower) = map(x->Index{Lower}(x),1:closure_size(L))
+region_indices(L, N, ::Interior) = map(x->Index{Interior}(x),closure_size(L)+1:N-closure_size(L))
+region_indices(L, N, ::Upper) = map(x->Index{Upper}(x),N-closure_size(L)+1:N)
+
+function get_region_indices(L,N)
+    ind_lower = region_indices(L, N, Lower())
+    ind_interior = region_indices(L, N, Interior())
+    ind_upper = region_indices(L, N, Upper())
+    return (ind_lower, ind_interior, ind_upper)
+end
+
+function apply_laplace_regions!(f, u, L, region_inds)
+    apply_laplace!(f, u, L, region_inds[1])
+    apply_laplace!(f, u, L, region_inds[2])
+    apply_laplace!(f, u, L, region_inds[3])
+end
+
+# Domain
+N = 4001
+g = EquidistantGrid(N,0.,1.)
+
+# Operators
+L = Laplace(g,sbp_operators_path()*"standard_diagonal.toml"; order=4)
+u = evalOn(g,x->x^2)
+f = similar(u)
+
+apply_laplace!(f,u,L) #ensure compilation
+@btime apply_laplace!(f,u,L)
+rinds = get_region_indices(L,N)
+apply_laplace_regions!(f,u,L,rinds) #ensure compilation
+@btime apply_laplace_regions!(f,u,L,rinds)
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/benchmarks/laplace_benchmark.m	Thu Jan 27 11:00:31 2022 +0100
@@ -0,0 +1,14 @@
+m = 4001;
+ops = sbp.D2Standard(m,{0,1},4);
+D2 = ops.D2;
+u = linspace(0,1,m)';
+f = zeros(size(u));
+
+nsample = 10000;
+ts = zeros(nsample,1);
+
+for i = 1:nsample
+    tic; f = D2*u; t = toc;
+    ts(i) = t;
+end
+min(ts)
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/wave_eq.jl	Thu Jan 27 11:00:31 2022 +0100
@@ -0,0 +1,67 @@
+using Sbplib.Grids, Sbplib.SbpOperators, Sbplib.LazyTensors, Sbplib.RegionIndices
+using OrdinaryDiffEq, Plots, Printf, Base.Threads
+
+function apply_tm!(f,u,tm,ind)
+    for I in ind
+        @inbounds f[I] = (tm*u)[I]
+    end
+end
+
+function apply_tm_all_regions!(f,u,tm,rinds)
+    apply_tm!(f,u,tm,rinds[1])
+    apply_tm!(f,u,tm,rinds[2])
+    apply_tm!(f,u,tm,rinds[3])
+end
+
+region_indices(L, N, ::Lower) = map(x->Index{Lower}(x),1:closure_size(L))
+region_indices(L, N, ::Interior) = map(x->Index{Interior}(x),closure_size(L)+1:N-closure_size(L))
+region_indices(L, N, ::Upper) = map(x->Index{Upper}(x),N-closure_size(L)+1:N)
+
+function get_region_indices(L,N)
+    ind_lower = region_indices(L, N, Lower())
+    ind_interior = region_indices(L, N, Interior())
+    ind_upper = region_indices(L, N, Upper())
+    return (ind_lower, ind_interior, ind_upper)
+end
+
+function wave_eq_sim(alg,T,CFL)
+    # Domain
+    N = 101
+    g = EquidistantGrid(N,0.,1.)
+    dx = min(spacing(g)...)
+
+    # Spatial discretization
+    Δ = Laplace(g,sbp_operators_path()*"standard_diagonal.toml"; order=4)
+    (id_l, id_r) = boundary_identifiers(g)
+    SAT_l = boundary_condition(Δ,id_l,NeumannBC())
+    SAT_r = boundary_condition(Δ,id_r,NeumannBC())
+    tm = (Δ + SAT_l + SAT_r)
+
+    # RHS function
+    rinds = get_region_indices(Δ,N)
+    function f(du,u,p,t)
+        du[1:N] .= u[N+1:end]
+        apply_tm_all_regions!(view(du,N+1:2*N), view(u,1:N), tm, rinds)
+    end
+    # Initial condition
+    sigma = 0.1
+    ic_u = x->1/(sigma*sqrt(2*pi))*exp(-1/2*((x-0.5)^2/sigma^2))
+    ic_u_t = x->0
+    w0 = [evalOn(g,ic_u);
+          evalOn(g,ic_u_t)]
+    # Setup ODE and solve
+    tspan = (0.,T)
+    prob = ODEProblem(f,w0,tspan)
+    sol = solve(prob, alg, dt=CFL*dx, saveat=0.05)
+
+    # Plotting
+    x = [x[1] for x in points(g)]
+    anim = @animate for i ∈ eachindex(sol.t)
+        u_i = sol.u[i]
+        plot(x, u_i[1:N], ylims = (0,4), lw=3,ls=:dash,label="",title=@sprintf("u at t = %.3f", sol.t[i]))
+    end
+    gif(anim, "wave.gif", fps = 15)
+end
+
+wave_eq_sim(CarpenterKennedy2N54(),1.,0.25)
+
--- a/laplace_benchmark.jl	Thu Jan 27 10:55:08 2022 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,42 +0,0 @@
-using Sbplib.Grids, Sbplib.SbpOperators, Sbplib.LazyTensors, Sbplib.RegionIndices
-using Profile, BenchmarkTools
-
-function apply_laplace!(f, u, L, inds)
-    for I in inds
-        f[I] = (L*u)[I]
-    end
-end
-
-apply_laplace!(f, u, L) = apply_laplace!(f, u, L, eachindex(u))
-
-region_indices(L, N, ::Lower) = map(x->Index{Lower}(x),1:closure_size(L))
-region_indices(L, N, ::Interior) = map(x->Index{Interior}(x),closure_size(L)+1:N-closure_size(L))
-region_indices(L, N, ::Upper) = map(x->Index{Upper}(x),N-closure_size(L)+1:N)
-
-function get_region_indices(L,N)
-    ind_lower = region_indices(L, N, Lower())
-    ind_interior = region_indices(L, N, Interior())
-    ind_upper = region_indices(L, N, Upper())
-    return (ind_lower, ind_interior, ind_upper)
-end
-
-function apply_laplace_regions!(f, u, L, region_inds)
-    apply_laplace!(f, u, L, region_inds[1])
-    apply_laplace!(f, u, L, region_inds[2])
-    apply_laplace!(f, u, L, region_inds[3])
-end
-
-# Domain
-N = 4001
-g = EquidistantGrid(N,0.,1.)
-
-# Operators
-L = Laplace(g,sbp_operators_path()*"standard_diagonal.toml"; order=4)
-u = evalOn(g,x->x^2)
-f = similar(u)
-
-apply_laplace!(f,u,L) #ensure compilation
-@btime apply_laplace!(f,u,L)
-rinds = get_region_indices(L,N)
-apply_laplace_regions!(f,u,L,rinds) #ensure compilation
-@btime apply_laplace_regions!(f,u,L,rinds)
--- a/laplace_benchmark.m	Thu Jan 27 10:55:08 2022 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,14 +0,0 @@
-m = 4001;
-ops = sbp.D2Standard(m,{0,1},4);
-D2 = ops.D2;
-u = linspace(0,1,m)';
-f = zeros(size(u));
-
-nsample = 10000;
-ts = zeros(nsample,1);
-
-for i = 1:nsample
-    tic; f = D2*u; t = toc;
-    ts(i) = t;
-end
-min(ts)
\ No newline at end of file
Binary file wave.gif has changed
--- a/wave_eq.jl	Thu Jan 27 10:55:08 2022 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,98 +0,0 @@
-using Sbplib.Grids, Sbplib.SbpOperators, Sbplib.LazyTensors, Sbplib.RegionIndices
-using OrdinaryDiffEq, Plots, Printf, Base.Threads
-
-function apply_tm!(f,u,tm,ind)
-    for I in ind
-        @inbounds f[I] = (tm*u)[I]
-    end
-end
-
-function apply_tm_all_regions!(f,u,tm,rinds)
-    apply_tm!(f,u,tm,rinds[1])
-    apply_tm!(f,u,tm,rinds[2])
-    apply_tm!(f,u,tm,rinds[3])
-end
-
-region_indices(L, N, ::Lower) = map(x->Index{Lower}(x),1:closure_size(L))
-region_indices(L, N, ::Interior) = map(x->Index{Interior}(x),closure_size(L)+1:N-closure_size(L))
-region_indices(L, N, ::Upper) = map(x->Index{Upper}(x),N-closure_size(L)+1:N)
-
-function get_region_indices(L,N)
-    ind_lower = region_indices(L, N, Lower())
-    ind_interior = region_indices(L, N, Interior())
-    ind_upper = region_indices(L, N, Upper())
-    return (ind_lower, ind_interior, ind_upper)
-end
-
-function wave_eq_sim(alg,T,CFL)
-    # Domain
-    N = 101
-    g = EquidistantGrid(N,0.,1.)
-    dx = min(spacing(g)...)
-
-    # Spatial discretization
-    Δ = Laplace(g,sbp_operators_path()*"standard_diagonal.toml"; order=4)
-    (id_l, id_r) = boundary_identifiers(g)
-    SAT_l = boundary_condition(Δ,id_l,NeumannBC())
-    SAT_r = boundary_condition(Δ,id_r,NeumannBC())
-    tm = (Δ + SAT_l + SAT_r)
-
-    # RHS function
-    rinds = get_region_indices(Δ,N)
-    function f(du,u,p,t)
-        du[1:N] .= u[N+1:end]
-        apply_tm_all_regions!(view(du,N+1:2*N), view(u,1:N), tm, rinds)
-    end
-    # Initial condition
-    sigma = 0.1
-    ic_u = x->1/(sigma*sqrt(2*pi))*exp(-1/2*((x-0.5)^2/sigma^2))
-    ic_u_t = x->0
-    w0 = [evalOn(g,ic_u);
-          evalOn(g,ic_u_t)]
-    # Setup ODE and solve
-    tspan = (0.,T)
-    prob = ODEProblem(f,w0,tspan)
-    sol = solve(prob, alg, dt=CFL*dx, saveat=0.05)
-
-    # Plotting
-    x = [x[1] for x in points(g)]
-    anim = @animate for i ∈ eachindex(sol.t)
-        u_i = sol.u[i]
-        plot(x, u_i[1:N], ylims = (0,4), lw=3,ls=:dash,label="",title=@sprintf("u at t = %.3f", sol.t[i]))
-    end
-    gif(anim, "wave.gif", fps = 15)
-end
-
-wave_eq_sim(CarpenterKennedy2N54(),1.,0.25)
-
-# function boundary_condition(L,id; )
-#     neumann_closure
-#     neumann_penalty
-# end
-#
-# function neumann_bc(L,id)
-#     e = boundary_restriction.(L,ids)
-#     d = normal_derivative(L,ids)
-#     return (e'∘d,)
-# end
-#
-# function (closure, penalty) neumann_bc(L,id,g::Function)
-#     e = boundary_restriction.(L,ids)
-#     d = normal_derivative.(L,ids)
-#     return e'∘d
-# end
-#
-# function dirichlet_closure(L,id)
-#     e = boundary_restriction.(L,ids)
-#     d = normal_derivative.(L,ids)
-#     return ...
-# end
-#
-# function SAT(L,ids,BCs)
-#     BC = BCs(L,ids[1])
-#     for i = 2:length(ids)
-#         BC = BC+ BCs(L,ids[1])
-#     end
-#     H_inv = inverse_inner_product(L)
-#     return H_inv∘BC
-# end