Mercurial > repos > public > sbplib_julia
comparison wave_eq.jl @ 874:7e9ebd572deb laplace_benchmarks
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 |
comparison
equal
deleted
inserted
replaced
| 873:9929c99754fb | 874:7e9ebd572deb |
|---|---|
| 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 | |
| 68 # function boundary_condition(L,id; ) | |
| 69 # neumann_closure | |
| 70 # neumann_penalty | |
| 71 # end | |
| 72 # | |
| 73 # function neumann_bc(L,id) | |
| 74 # e = boundary_restriction.(L,ids) | |
| 75 # d = normal_derivative(L,ids) | |
| 76 # return (e'∘d,) | |
| 77 # end | |
| 78 # | |
| 79 # function (closure, penalty) neumann_bc(L,id,g::Function) | |
| 80 # e = boundary_restriction.(L,ids) | |
| 81 # d = normal_derivative.(L,ids) | |
| 82 # return e'∘d | |
| 83 # end | |
| 84 # | |
| 85 # function dirichlet_closure(L,id) | |
| 86 # e = boundary_restriction.(L,ids) | |
| 87 # d = normal_derivative.(L,ids) | |
| 88 # return ... | |
| 89 # end | |
| 90 # | |
| 91 # function SAT(L,ids,BCs) | |
| 92 # BC = BCs(L,ids[1]) | |
| 93 # for i = 2:length(ids) | |
| 94 # BC = BC+ BCs(L,ids[1]) | |
| 95 # end | |
| 96 # H_inv = inverse_inner_product(L) | |
| 97 # return H_inv∘BC | |
| 98 # end |