make gifs for each element in the cvg series

plot/Manifest.toml

 # 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"]

plot/gif.jl

 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

plot/heateq_gaussian.toml → plot/heateq_sine_1.toml

 [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

plot/heateq_sine_2.toml

+[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

plot/heateq_sine_3.toml

+[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

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