Mercurial > repos > public > sbplib_julia

comparison examples/wave_eq.jl @ 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 wave_eq.jl@7e9ebd572deb
children
comparison
equal deleted inserted replaced
875:067a322e4f73 876:4f3924293894
1 using Sbplib.Grids, Sbplib.SbpOperators, Sbplib.LazyTensors, Sbplib.RegionIndices
2 using OrdinaryDiffEq, Plots, Printf, Base.Threads
3
4 function apply_tm!(f,u,tm,ind)
5 for I in ind
6 @inbounds f[I] = (tm*u)[I]
7 end
8 end
9
10 function apply_tm_all_regions!(f,u,tm,rinds)
11 apply_tm!(f,u,tm,rinds[1])
12 apply_tm!(f,u,tm,rinds[2])
13 apply_tm!(f,u,tm,rinds[3])
14 end
15
16 region_indices(L, N, ::Lower) = map(x->Index{Lower}(x),1:closure_size(L))
17 region_indices(L, N, ::Interior) = map(x->Index{Interior}(x),closure_size(L)+1:N-closure_size(L))
18 region_indices(L, N, ::Upper) = map(x->Index{Upper}(x),N-closure_size(L)+1:N)
19
20 function get_region_indices(L,N)
21 ind_lower = region_indices(L, N, Lower())
22 ind_interior = region_indices(L, N, Interior())
23 ind_upper = region_indices(L, N, Upper())
24 return (ind_lower, ind_interior, ind_upper)
25 end
26
27 function wave_eq_sim(alg,T,CFL)
28 # Domain
29 N = 101
30 g = EquidistantGrid(N,0.,1.)
31 dx = min(spacing(g)...)
32
33 # Spatial discretization
34 Δ = Laplace(g,sbp_operators_path()*"standard_diagonal.toml"; order=4)
35 (id_l, id_r) = boundary_identifiers(g)
36 SAT_l = boundary_condition(Δ,id_l,NeumannBC())
37 SAT_r = boundary_condition(Δ,id_r,NeumannBC())
38 tm = (Δ + SAT_l + SAT_r)
39
40 # RHS function
41 rinds = get_region_indices(Δ,N)
42 function f(du,u,p,t)
43 du[1:N] .= u[N+1:end]
44 apply_tm_all_regions!(view(du,N+1:2*N), view(u,1:N), tm, rinds)
45 end
46 # Initial condition
47 sigma = 0.1
48 ic_u = x->1/(sigma*sqrt(2*pi))*exp(-1/2*((x-0.5)^2/sigma^2))
49 ic_u_t = x->0
50 w0 = [evalOn(g,ic_u);
51 evalOn(g,ic_u_t)]
52 # Setup ODE and solve
53 tspan = (0.,T)
54 prob = ODEProblem(f,w0,tspan)
55 sol = solve(prob, alg, dt=CFL*dx, saveat=0.05)
56
57 # Plotting
58 x = [x[1] for x in points(g)]
59 anim = @animate for i ∈ eachindex(sol.t)
60 u_i = sol.u[i]
61 plot(x, u_i[1:N], ylims = (0,4), lw=3,ls=:dash,label="",title=@sprintf("u at t = %.3f", sol.t[i]))
62 end
63 gif(anim, "wave.gif", fps = 15)
64 end
65
66 wave_eq_sim(CarpenterKennedy2N54(),1.,0.25)
67