diff --git a/plot/Manifest.toml b/plot/Manifest.toml
index d5783786d5daa0886a51413e8688a9f8c7550036..e247526443beb434175359a03336928816588550 100644
--- a/plot/Manifest.toml
+++ b/plot/Manifest.toml
@@ -1,6 +1,6 @@
# This file is machine-generated - editing it directly is not advised
-julia_version = "1.10.2"
+julia_version = "1.10.5"
manifest_format = "2.0"
project_hash = "43a80ecfaaf92590af56c030ba6fac5bbf756dc8"
@@ -243,7 +243,7 @@ weakdeps = ["Dates", "LinearAlgebra"]
[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
-version = "1.1.0+0"
+version = "1.1.1+0"
[[deps.CompositionsBase]]
git-tree-sha1 = "802bb88cd69dfd1509f6670416bd4434015693ad"
@@ -1745,7 +1745,7 @@ version = "0.15.1+0"
[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
-version = "5.8.0+1"
+version = "5.11.0+0"
[[deps.libevent_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "OpenSSL_jll"]
diff --git a/plot/gif.jl b/plot/gif.jl
index 71454f2fb6b2ef5b742e51a34c48b1dba2352151..03c3e15325e0c1a82d669c900f7868b9eff33ec9 100644
--- a/plot/gif.jl
+++ b/plot/gif.jl
@@ -1,13 +1,13 @@
using GLMakie
+using Printf
GLMakie.activate!()
using dg1d
-let
MT = GLMakie.Makie.MathTeXEngine
mt_fonts_dir = joinpath(dirname(pathof(MT)), "..", "assets", "fonts", "NewComputerModern")
-theme = Theme(
+mytheme = Theme(
fonts = ( regular = joinpath(mt_fonts_dir, "NewCM10-Regular.otf"),
bold = joinpath(mt_fonts_dir, "NewCM10-Bold.otf"),
bold_italic = joinpath(mt_fonts_dir, "NewCM10-BoldItalic.otf"),
@@ -17,26 +17,87 @@ theme = Theme(
)
-dir = joinpath(@__DIR__, "heateq_gaussian")
+function read_data(dir::String)
+ output = Output(joinpath(dir, "output1d.h5"))
+ variables = output.variables
+ varnames = [ last(split(var, ".")) for var in variables ]
+ nvars = length(variables)
+ ts = output.times
+ Nts = length(ts)
+ u(t) = read_variable(output, "u", t)
+ return ts, output.x, u
+end
+
+function make_gif(filename, ts, x, u, N, K)
+ framerate=10
+ buf_x = zeros((N+2)*K)
+ buf_y = zeros((N+2)*K)
+ # separate cell data by NaN rows to make lines discontinuous
+ fill!(buf_x, NaN); fill!(buf_y, NaN)
+ # idxs = [ i for i in 1:length(buf_x) if i % (N+2) != 0 ]
+ idxs = [ i%(N+2)!=0 for i in 1:length(buf_x) ]
+ buf_x[idxs] .= x
+ buf_y[idxs] .= u(0.0)
+ obs_y = Observable(buf_y)
+ obs_t = Observable(ts[1])
+ lbl_t = lift(obs_t) do t
+ fmt_t = @sprintf "%.02f" t
+ L"t = %$fmt_t"
+ end
+ with_theme(mytheme) do
+ fig = Figure()
+ ax = Axis(fig[1,1], xlabel=L"x", ylabel="u")
+ lines!(ax, buf_x, obs_y, label=L"K=%$K")
+ text!(ax, -1.0, 0.5; text=lbl_t)
+ axislegend(ax)
+ record(fig, joinpath(@__DIR__,filename), ts[2:end]; framerate) do t
+ obs_y.val[idxs] .= u(t); notify(obs_y)
+ obs_t[] = t
+ end
+ end
+end
-# read data
-output = Output(joinpath(dir, "output1d.h5"))
-variables = output.variables
-varnames = [ last(split(var, ".")) for var in variables ]
-nvars = length(variables)
-ts = output.times
-Nts = length(ts)
-q(t) = read_variable(output, "q", t)
-obs_y = Observable(q(0.0))
+let
-with_theme(theme) do
+
+dir = joinpath(@__DIR__, "heateq_sine_1")
+ts, x, u = read_data(dir)
+N, K = 1, 10
+make_gif("heateq_sine_1.gif", ts, x, u, N, K)
+
+dir = joinpath(@__DIR__, "heateq_sine_2")
+ts, x, u = read_data(dir)
+N, K = 1, 40
+make_gif("heateq_sine_2.gif", ts, x, u, N, K)
+
+dir = joinpath(@__DIR__, "heateq_sine_3")
+ts, x, u = read_data(dir)
+N, K = 1, 80
+make_gif("heateq_sine_3.gif", ts, x, u, N, K)
+
+# analytical solution
+L = 4.0
+x = range(extrema(x)..., length=500)
+u(t, x) = exp(-4*pi^2/L^2 * t) * sin(2*pi*x/L)
+framerate=10
+obs_y = Observable(u.(0.0,x))
+obs_t = Observable(ts[1])
+lbl_t = lift(obs_t) do t
+ fmt_t = @sprintf "%.02f" t
+ L"t = %$fmt_t"
+end
+with_theme(mytheme) do
fig = Figure()
- ax = Axis(fig[1,1], xlabel=L"x", ylabel="q")
- lines!(ax, output.x, obs_y)
- record(fig, "heateq.gif", ts[2:end]) do t
- obs_y[] = q(t)
+ ax = Axis(fig[1,1], xlabel=L"x", ylabel="u")
+ lines!(ax, x, obs_y, label="exact")
+ text!(ax, -1.0, 0.5; text=lbl_t)
+ axislegend(ax)
+ record(fig, joinpath(@__DIR__,"heateq_sine_exact.gif"), ts[2:end]; framerate) do t
+ obs_y[] = u.(t,x)
+ obs_t[] = t
end
end
+
end
diff --git a/plot/heateq_gaussian.toml b/plot/heateq_sine_1.toml
similarity index 63%
rename from plot/heateq_gaussian.toml
rename to plot/heateq_sine_1.toml
index 5b3669906502731b49347286dae25460b8dfaf70..34fd4c6949de906ff7530961d381ae0d9a84ceb1 100644
--- a/plot/heateq_gaussian.toml
+++ b/plot/heateq_sine_1.toml
@@ -1,6 +1,6 @@
[Evolution]
-cfl = 0.8
-tend = 2.4
+cfl = 0.1
+tend = 1.2
[HeatEq]
id = "sine"
@@ -8,11 +8,11 @@ id_sine_width = 4.0
[Mesh]
range = [-2.0, 2.0]
-k = 20
+k = 10
basis = "lgl"
-n = 5
+n = 1
[Output]
variables1d = ["u", "q"]
-aligned_ts = "$(collect(range(0.02,10.4,step=0.02)))"
+aligned_ts = "$(collect(range(0.01,10.4,step=0.01)))"
every_iteration = 0
diff --git a/plot/heateq_sine_2.toml b/plot/heateq_sine_2.toml
new file mode 100644
index 0000000000000000000000000000000000000000..66e166233eda1c7693b36284fa03b7ec482368b0
--- /dev/null
+++ b/plot/heateq_sine_2.toml
@@ -0,0 +1,18 @@
+[Evolution]
+cfl = 0.1
+tend = 1.2
+
+[HeatEq]
+id = "sine"
+id_sine_width = 4.0
+
+[Mesh]
+range = [-2.0, 2.0]
+k = 40
+basis = "lgl"
+n = 1
+
+[Output]
+variables1d = ["u", "q"]
+aligned_ts = "$(collect(range(0.01,10.4,step=0.01)))"
+every_iteration = 0
diff --git a/plot/heateq_sine_3.toml b/plot/heateq_sine_3.toml
new file mode 100644
index 0000000000000000000000000000000000000000..2cef74cc4893cc3450d91a57b51548907813f93d
--- /dev/null
+++ b/plot/heateq_sine_3.toml
@@ -0,0 +1,18 @@
+[Evolution]
+cfl = 0.1
+tend = 1.2
+
+[HeatEq]
+id = "sine"
+id_sine_width = 4.0
+
+[Mesh]
+range = [-2.0, 2.0]
+k = 80
+basis = "lgl"
+n = 1
+
+[Output]
+variables1d = ["u", "q"]
+aligned_ts = "$(collect(range(0.01,10.4,step=0.01)))"
+every_iteration = 0
diff --git a/src/HeatEq/rhs.jl b/src/HeatEq/rhs.jl
index 878cec8fc9db7db41f56d4d9409be6e8b5da383a..e13ff7560d2efc62e0015e530e6e93ee556092cb 100644
--- a/src/HeatEq/rhs.jl
+++ b/src/HeatEq/rhs.jl
@@ -17,11 +17,9 @@ function timestep(env, P::Project, ::Mesh1d)
@unpack element = mesh
@unpack N, z = element
- dz = z[2]-z[1]
- dx, = minimum(dg1d.widths(mesh.boxes[1]))
- dl = dx * dz
+ dl = min_grid_spacing(mesh)
- dt = dl^2 / (N^2 * vmax) # just a guess to use 1/N^2 here since we have 2nd derivatives
+ dt = dl^2 / (4*vmax)
return dt
end