adds plots

src/InitialData.jl

@@ -138,9 +138,9 @@ end
                                           36 * M^2 / r^4 * (mul(l) + 2 * M / r))
 end
 
-@inline function poschl_teller(l::Real, r::Real)
-    # return -(l * (l + 1) / 2) * sech(r)^2
-    return -sech(r)^2
+@inline function poschl_teller(l, r)
+    return -(l * (l + 1) / 2) * sech(r)^2
+    # return -sech(r)^2
 end
 
 end #end of module

src/ODE.jl

@@ -7,6 +7,7 @@ using ..InitialData
     N = params.N
     l = params.l
     M = params.M
+    L = params.L
     rstar = params.rcoord
     boundary_type = params.bc
     @fastmath @inbounds begin
@@ -29,8 +30,8 @@ using ..InitialData
                 drΨ = (Ψ[end] - Ψ[N1]) / h
                 drΠ = (Π[end] - Π[N1]) / h
             else
-                drΠ = (Π[i+1] - Π[i-1]) / h2
-                drΨ = (Ψ[i+1] - Ψ[i-1]) / h2
+                drΠ = (Π[i + 1] - Π[i - 1]) / h2
+                drΨ = (Ψ[i + 1] - Ψ[i - 1]) / h2
             end
             dtΦ[i] = Π[i]
             # Apply boundary conditions
@@ -77,10 +78,10 @@ using ..InitialData
                            $boundary_type")
                 end
             else
-                # ri = InitialData.r(rstar[i], M) # transform rstar to r for potential
-                # @assert ri > 2M
-                # dtΠ[i] = drΨ + InitialData.Vpot(l, rstar[i], M) * Φ[i]
-                dtΠ[i] = drΨ + InitialData.poschl_teller(l, rstar[i]) * Φ[i]
+                ri = InitialData.r(rstar[i], M) # transform rstar to r for potential
+                # dtΠ[i] = drΨ + InitialData.poschl_teller(l, rstar[i]) * Φ[i]
+                dtΠ[i] = drΨ - cos((π / 2L) * rstar[i]) * Φ[i]
+                # dtΠ[i] = drΨ + InitialData.poschl_teller(l, ri) * Φ[i]
                 dtΨ[i] = drΠ
             end
         end

src/Plot.jl

@@ -2,41 +2,78 @@ module Plot
 
 using Plots
 using GridFunctions
+using ..InitialData
 
-function plot_ti(sim_folder, ti)
-    rs = GridFunctions.coords(get_grid_from_yaml(sim_folder))
+function plot_ti(u, ti)
+    rs = GridFunctions.coords(u[2])
 
-    Φ = h5read(dataset, "Phi")
-    Π = h5read(dataset, "Pi")
-    Ψ = h5read(dataset, "Psi")
+    Φ = u[3][ti][:, 1]
+    Π = u[3][ti][:, 2]
+    Ψ = u[3][ti][:, 3]
 
-    p = plot(rs, Φ)
+    p = plot(rs, Φ; label="Φ")
+    plot!(rs, Π; label="Π")
+    plot!(rs, Ψ; label="Ψ")
     return p
 end
 
 function trapz(f, h)
     return (h / 2) * (f[1] + f[end] + 2sum(f[2:(end - 1)]))
 end
 
-function field_energy()
-    energy = @. 0.5 * (ppi^2 + psi^2 + InitialData.poschl_teller(l, rs) * phi)
+function field_energy(u, ti)
+    h = u[4].h
+    l = u[4].l
+    L = u[4].L
+    rs = GridFunctions.coords(u[2])
+    Φ = u[3][ti][:, 1]
+    Π = u[3][ti][:, 2]
+    Ψ = u[3][ti][:, 3]
+    # energy = @. 0.5 * (Π * Π + Ψ * Ψ + InitialData.poschl_teller((l), rs) * Φ)
+    energy = @. 0.5 * (Π * Π + Ψ * Ψ - cos((π / 2L) * rs) * Φ)
+
     return trapz(energy, h)
 end
 
-function plot_energy(sim_folder)
-    sims = readdir(string(sim_folder, "/sims"); join=true)
-    r = get_grid_from_yaml(sim_folder)
-    h = spacing(r)
-    n = length(sims)
+function plot_energy(u; show=true)
+    t = u[3].t
+    n = length(t)
     e = Vector{Float64}(undef, n)
-    t = GridFunctions.coords(get_time_grid_from_yaml(sim_folder))
-    for (i, dataset) in enumerate(sims)
-        Π = h5read(dataset, "Pi")
-        Ψ = h5read(dataset, "Psi")
-        e[i] = field_energy(h, Π, Ψ)
+
+    for i in 1:n
+        e[i] = field_energy(u, i)
     end
-    p = plot(t, e)
-    return p
+
+    show == true && (p = plot(t, e); display(p))
+    return t, e
+end
+
+function plot_dtE(u)
+    t, e = plot_energy(u; show=false)
+    dt = t[2] - t[1]
+    h = u[4].h
+    l = u[4].l
+    L = u[4].L
+    n = length(t)
+    r = GridFunctions.coords(u[2])
+    dedt = similar(e)
+    pt = similar(e)
+
+    dedt[1] = (e[2] - e[1]) / dt
+    dedt[end] = (e[end] - e[end - 1]) / dt
+    for i in 1:n
+        i > 1 && i < n && (dedt[i] = (e[i + 1] - e[i - 1]) / (2dt))
+        Π = u[3][i][:, 2]
+        Ψ = u[3][i][:, 3]
+        # term = @. 0.5 * Π * InitialData.poschl_teller.(l, r)
+        term = @. -0.5 * Π * cos((π / 2L) * r)
+        pt[i] = trapz(term, h) + Π[end] * Ψ[end] - Π[1] * Ψ[1]
+    end
+
+    p = plot(t, dedt; label="numerical")
+    plot!(p, t, pt; ls=:dash, label="semi-analytic")
+    display(p)
+    return nothing
 end
 
 end

src/Run.jl

@@ -8,10 +8,9 @@ using OrdinaryDiffEqTsit5
 
 export zerilli
 
-function zerilli(domain::Vector, ncells::Integer,
-                 tf::Real, cfl::Real, M::Real=1.0, l::Real=0;
-                 boundary_type::Symbol=:radiative,
-                 folder="", save_every=1)
+function zerilli(; domain::Vector, ncells::Integer,
+                 tf::Real, cfl::Real, M::Real=1.0, l::Real=0,
+                 boundary_type::Symbol=:radiative)
     @assert boundary_type === :radiative || boundary_type === :reflective
 
     grid = UniformGrid(domain, ncells)
@@ -24,19 +23,19 @@ function zerilli(domain::Vector, ncells::Integer,
 
     # Initial Data
 
-    # σ = 6
-    # r0 = 20# 20M / 2 + horizon
-    # Φ = InitialData.Gaussian1D.(0, coords(grid), σ, r0)
-    # Π = InitialData.dtGaussian1D.(0, coords(grid), σ, r0)
-    # Ψ = InitialData.drGaussian1D.(0, coords(grid), σ, r0)
+    σ = 6
+    r0 = 0 # 20# 20M / 2 + horizon
+    Φ = InitialData.Gaussian1D.(0.0, coords(grid), σ, r0)
+    Π = InitialData.dtGaussian1D.(0.0, coords(grid), σ, r0)
+    Ψ = InitialData.drGaussian1D.(0.0, coords(grid), σ, r0)
 
-    Φ = rand(N)
-    Π = rand(N)
-    Ψ = rand(N)
+    # Φ = rand(N)
+    # Π = rand(N)
+    # Ψ = rand(N)
 
     params = (h=spacing(grid), N=N, bc=boundary_type, t=t, ti=coords(t),
-              dt=spacing(t), save_every=save_every, M=M, grid=grid,
-              folder=folder, cfl=cfl, rcoord=coords(grid), l=l)
+              dt=spacing(t), M=M, grid=grid, L=domain[2],
+              cfl=cfl, rcoord=coords(grid), l=l)
     statevector = hcat(Φ, Π, Ψ)
 
     ode = ODEProblem{true}(ODE.rhs!, statevector, tspan, params)