Mercurial > repos > public > sbplib_julia

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/laplace_benchmark.jl	Thu Jan 27 10:55:08 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/laplace_benchmark.m	Thu Jan 27 10:55:08 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
--- a/src/SbpOperators/volumeops/laplace/laplace.jl	Tue Jan 25 10:36:13 2022 +0100
+++ b/src/SbpOperators/volumeops/laplace/laplace.jl	Thu Jan 27 10:55:08 2022 +0100
@@ -70,7 +70,7 @@
     d_pairs  = map(id -> Pair(id, normal_derivative(grid, d_closure_stencil, id)),                               ids)
     Hᵧ_pairs = map(id -> Pair(id, inner_product(boundary_grid(grid, id), H_inner_stencils, H_closure_stencils)), ids)

-    return Laplace(Δ, H, H⁻¹, StaticDict(e_pairs), StaticDict(d_pairs), StaticDict(Hᵧ_pairs))
+    return Laplace(Δ, H, H⁻¹, StaticDict(e_pairs), StaticDict(d_pairs), StaticDict(Hᵧ_pairs), order)
 end

 # TODO: Consider pretty printing of the following form
@@ -81,6 +81,14 @@
 LazyTensors.domain_size(L::Laplace) = LazyTensors.domain_size(L.D)
 LazyTensors.apply(L::Laplace, v::AbstractArray, I...) = LazyTensors.apply(L.D,v,I...)

+"""
+    closure_size(L::Lapalace)
+
+Returns the inner product operator associated with `L`
+
+"""
+closure_size(L::Laplace) = closure_size(L.D)
+export closure_size

 """
     inner_product(L::Laplace)
@@ -143,6 +151,25 @@
 boundary_quadrature(L::Laplace, ids::NTuple{N,BoundaryIdentifier}) where N = ntuple(i->L.H_boundary[ids[i]],N)
 boundary_quadrature(L::Laplace, ids::Vararg{BoundaryIdentifier,N}) where N = ntuple(i->L.H_boundary[ids[i]],N)

+abstract type BoundaryConditionType end
+struct NeumannBC <: BoundaryConditionType end
+struct DirichletBC <: BoundaryConditionType end
+export NeumannBC
+
+boundary_condition(L, id, ::NeumannBC) = neumann_bc(L, id)
+boundary_condition(L, id, ::DirichletBC) = dirichlet_bc(L, id)
+export boundary_condition
+
+function neumann_bc(L::Laplace, id::BoundaryIdentifier)
+    H_inv = inverse_inner_product(L)
+    e = boundary_restriction(L,id)
+    d = normal_derivative(L,id)
+    H_b = boundary_quadrature(L,id)
+    tau = H_inv∘e'∘H_b
+    closure = tau∘d
+    # TODO: Return penalty once we have implemented scalar scaling of the operators.
+    return closure
+end

 """
     laplace(grid::EquidistantGrid{Dim}, inner_stencil, closure_stencils)
Binary file wave.gif has changed
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/wave_eq.jl	Thu Jan 27 10:55:08 2022 +0100
@@ -0,0 +1,98 @@
+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