Mercurial > repos > public > sbplib_julia

annotate wave_eq.jl @ 874:7e9ebd572deb laplace_benchmarks

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add file for wave equation
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 20 Jan 2022 21:51:53 +0100
parents
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
874
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
1 using Sbplib.Grids, Sbplib.SbpOperators, Sbplib.LazyTensors, Sbplib.RegionIndices
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
2 using OrdinaryDiffEq, Plots, Printf, Base.Threads
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
3
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
4 function apply_tm!(f,u,tm,ind)
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
5 for I in ind
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
6 @inbounds f[I] = (tm*u)[I]
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
7 end
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
8 end
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
9
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
10 function apply_tm_all_regions!(f,u,tm,rinds)
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
11 apply_tm!(f,u,tm,rinds[1])
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
12 apply_tm!(f,u,tm,rinds[2])
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
13 apply_tm!(f,u,tm,rinds[3])
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
14 end
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
15
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
16 region_indices(L, N, ::Lower) = map(x->Index{Lower}(x),1:closure_size(L))
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
17 region_indices(L, N, ::Interior) = map(x->Index{Interior}(x),closure_size(L)+1:N-closure_size(L))
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
18 region_indices(L, N, ::Upper) = map(x->Index{Upper}(x),N-closure_size(L)+1:N)
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
19
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
20 function get_region_indices(L,N)
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
21 ind_lower = region_indices(L, N, Lower())
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
22 ind_interior = region_indices(L, N, Interior())
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
23 ind_upper = region_indices(L, N, Upper())
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
24 return (ind_lower, ind_interior, ind_upper)
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
25 end
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
26
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
27 function wave_eq_sim(alg,T,CFL)
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
28 # Domain
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
29 N = 101
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
30 g = EquidistantGrid(N,0.,1.)
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
31 dx = min(spacing(g)...)
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
32
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
33 # Spatial discretization
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
34 Δ = Laplace(g,sbp_operators_path()*"standard_diagonal.toml"; order=4)
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
35 (id_l, id_r) = boundary_identifiers(g)
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
36 SAT_l = boundary_condition(Δ,id_l,NeumannBC())
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
37 SAT_r = boundary_condition(Δ,id_r,NeumannBC())
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
38 tm = (Δ + SAT_l + SAT_r)
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
39
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
40 # RHS function
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
41 rinds = get_region_indices(Δ,N)
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
42 function f(du,u,p,t)
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
43 du[1:N] .= u[N+1:end]
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
44 apply_tm_all_regions!(view(du,N+1:2*N), view(u,1:N), tm, rinds)
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
45 end
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
46 # Initial condition
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
47 sigma = 0.1
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
48 ic_u = x->1/(sigma*sqrt(2*pi))*exp(-1/2*((x-0.5)^2/sigma^2))
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
49 ic_u_t = x->0
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
50 w0 = [evalOn(g,ic_u);
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
51 evalOn(g,ic_u_t)]
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
52 # Setup ODE and solve
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
53 tspan = (0.,T)
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
54 prob = ODEProblem(f,w0,tspan)
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
55 sol = solve(prob, alg, dt=CFL*dx, saveat=0.05)
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
56
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
57 # Plotting
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
58 x = [x[1] for x in points(g)]
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
59 anim = @animate for i ∈ eachindex(sol.t)
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
60 u_i = sol.u[i]
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
61 plot(x, u_i[1:N], ylims = (0,4), lw=3,ls=:dash,label="",title=@sprintf("u at t = %.3f", sol.t[i]))
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
62 end
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
63 gif(anim, "wave.gif", fps = 15)
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
64 end
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
65
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
66 wave_eq_sim(CarpenterKennedy2N54(),1.,0.25)
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
67
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
68 # function boundary_condition(L,id; )
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
69 # neumann_closure
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
70 # neumann_penalty
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
71 # end
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
72 #
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
73 # function neumann_bc(L,id)
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
74 # e = boundary_restriction.(L,ids)
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
75 # d = normal_derivative(L,ids)
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
76 # return (e'∘d,)
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
77 # end
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
78 #
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
79 # function (closure, penalty) neumann_bc(L,id,g::Function)
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
80 # e = boundary_restriction.(L,ids)
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
81 # d = normal_derivative.(L,ids)
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
82 # return e'∘d
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
83 # end
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
84 #
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
85 # function dirichlet_closure(L,id)
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
86 # e = boundary_restriction.(L,ids)
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
87 # d = normal_derivative.(L,ids)
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
88 # return ...
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
89 # end
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
90 #
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
91 # function SAT(L,ids,BCs)
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
92 # BC = BCs(L,ids[1])
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
93 # for i = 2:length(ids)
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
94 # BC = BC+ BCs(L,ids[1])
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
95 # end
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
96 # H_inv = inverse_inner_product(L)
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
97 # return H_inv∘BC
7e9ebd572deb Add file for wave equation
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
98 # end