add blastwave1 parfiles and plot script

srhd_blastwave_1/avmda1_approx_covariant/avmda1_approx_covariant.toml

+[EquationOfState]
+eos = "idealgas"
+idealgas_gamma = 1.6666666666666666666667
+
+[SRHD]
+id = "blastwave1"
+id_smooth = false
+bc = "from_id"
+av_regularization = "mono_approx_covariant"
+
+[Mesh]
+range = [ -0.5, 0.5 ]
+n = 5
+k = 54
+basis = "lgl"
+periodic = false
+
+[Output]
+variables1d     = [ "D", "S", "tau", "p", "rho", "eps", "v", "smoothed_mu" ]
+aligned_ts      = "$(collect(range(0.01,0.4,step=0.01)))"
+
+[Evolution]
+cfl = 0.1
+tend = 0.4
+
+[HRSC]
+method = "av"
+av_method = "mda"

srhd_blastwave_1/avmda1_approx_covariant/dirty

+git commit hash: 4ad53e2aa88b3259d77af94910bdad4c76d112ea
+====================================================================================================
+

srhd_blastwave_1/avmda2_approx_covariant/avmda2_approx_covariant.toml

+[EquationOfState]
+eos = "idealgas"
+idealgas_gamma = 1.6666666666666666666667
+
+[SRHD]
+id = "blastwave1"
+id_smooth = false
+bc = "from_id"
+av_regularization = "mono_approx_covariant"
+
+[Mesh]
+range = [ -0.5, 0.5 ]
+n = 5
+k = 104
+basis = "lgl"
+periodic = false
+
+[Output]
+variables1d     = [ "D", "S", "tau", "p", "rho", "eps", "v", "smoothed_mu" ]
+aligned_ts      = "$(collect(range(0.01,0.4,step=0.01)))"
+
+[Evolution]
+cfl = 0.1
+tend = 0.4
+
+[HRSC]
+method = "av"
+av_method = "mda"

srhd_blastwave_1/avmda2_approx_covariant/dirty

+git commit hash: 4ad53e2aa88b3259d77af94910bdad4c76d112ea
+====================================================================================================
+

srhd_blastwave_1/avmda3_approx_covariant/avmda3_approx_covariant.toml

+[EquationOfState]
+eos = "idealgas"
+idealgas_gamma = 1.6666666666666666666667
+
+[SRHD]
+id = "blastwave1"
+id_smooth = false
+bc = "from_id"
+av_regularization = "mono_approx_covariant"
+
+[Mesh]
+range = [ -0.5, 0.5 ]
+n = 5
+k = 204
+basis = "lgl"
+periodic = false
+
+[Output]
+variables1d     = [ "D", "S", "tau", "p", "rho", "eps", "v", "smoothed_mu" ]
+aligned_ts      = "$(collect(range(0.01,0.4,step=0.01)))"
+
+[Evolution]
+cfl = 0.1
+tend = 0.4
+
+[HRSC]
+method = "av"
+av_method = "mda"

srhd_blastwave_1/avmda3_approx_covariant/dirty

+git commit hash: 4ad53e2aa88b3259d77af94910bdad4c76d112ea
+====================================================================================================
+

srhd_blastwave_1/plot_profiles.jl

+using CairoMakie
+using HDF5
+using FFTW
+using TOML
+
+
+include(joinpath(@__DIR__,"..","theme.jl"))
+include(joinpath(@__DIR__,"..","utils.jl"))
+
+
+function make_plot(rootdir, dirs)
+
+    for dir in dirs
+        @assert isdir(joinpath(rootdir,dir)) "missing output directory $(joinpath(rootdir,dir))"
+    end
+
+    ### load
+    data = [ readdata(joinpath(rootdir, dir), ["rho","p","v","smoothed_mu"],
+                      output="output1d.h5") for dir in dirs ]
+    pars = [ TOML.parsefile(joinpath(rootdir,dir,"$(replace(dir,".previous"=>"")).toml")) for dir in dirs ]
+    ns = [ par["Mesh"]["n"] for par in pars ]
+    ks = [ par["Mesh"]["k"] for par in pars ]
+
+    ### sort wrt k
+    P          = sortperm(ks)
+    data       = data[P]
+    pars       = pars[P]
+    ns         = ns[P]
+    ks         = ks[P]
+    dirs       = dirs[P]
+
+    ts = [ d.t for d in data ]
+    xs = [ d.x for d in data ]
+    ρs = [ d.rho for d in data ]
+    ps = [ d.p for d in data ]
+    vs = [ d.v for d in data ]
+    μs = [ d.smoothed_mu for d in data ]
+
+    stack_μs   = stack.(μs)
+    min_μ = minimum(minimum, stack_μs)
+    max_μ = maximum(maximum, stack_μs)
+
+    #######################################################################
+    #                              plotting                               #
+    #######################################################################
+
+    fig = Figure(size=(1200,800))
+    gl_ρ = fig[1,1:2] = GridLayout()
+    gl_μs = fig[2,1] = GridLayout()
+    gl_cb = fig[2,2] = GridLayout(; alignmode=Outside())
+
+    ax_ρ = Axis(gl_ρ[1,1], xlabel=L"x")
+    ax_μs = [ Axis(gl_μs[1,i], xlabel=L"x") for i in 1:length(dirs) ]
+    inset_ax_ρ = place_inset_axis!(gl_ρ[1,1],
+                                   inset_halign=0.9, inset_valign=0.9,
+                                   inset_width=Relative(0.5), inset_height=Relative(0.6))
+
+    # plot density and inset
+    # ylabel_ρ=L"\rho(x), p(x), v(x)"
+    # Label(fig[1,3], ylabel_ρ,
+    #       tellwidth=true, tellheight=false, rotation=pi/2)
+    for (i,dir) in enumerate(dirs)
+        lbl = "K = $(ks[i])"
+        for ax in (ax_ρ, inset_ax_ρ)
+            lines!(ax, xs[i], ρs[i][end #= final time =#], label=lbl,
+                   color=colors[periodic_index(i,length(colors))], linestyle=:solid)
+            lines!(ax, xs[i], ps[i][end #= final time =#],
+                   color=colors[periodic_index(i,length(colors))], linestyle=:dash)
+            lines!(ax, xs[i], vs[i][end #= final time =#],
+                   color=colors[periodic_index(i,length(colors))], linestyle=:dashdot)
+        end
+    end
+
+    # place a label to avoid on top of data inside plot
+    x0_text = first(xs)[1]
+    ρ0_text = first(ρs)[end][1]
+    p0_text = first(ps)[end][1]
+    v0_text = first(vs)[end][1]
+    @show x0_text, ρ0_text, v0_text
+    text!(ax_ρ, x0_text, 1.02*ρ0_text, text=L"ρ(x)")
+    text!(ax_ρ, x0_text, 1.02*p0_text, text=L"p(x)")
+    text!(ax_ρ, x0_text, 1.02*v0_text, text=L"v(x)")
+
+    xlims!(ax_ρ, (-0.52,0.52))
+    # ylims!(ax_ρ, (0.6,13.6))
+    ylims!(ax_ρ, (-0.8,16.6))
+    ax_ρ.xticks[] = range(-0.5,0.5,step=0.1)
+
+    # draw the frame for inset
+    xlims!(inset_ax_ρ, (0.06,0.38))
+    ylims!(inset_ax_ρ, (0.2,5.4))
+    inset_ax_ρ.yticks[] = WilkinsonTicks(3)
+    reset_limits!(inset_ax_ρ, xauto=false, yauto=false)
+    draw_inset_frame!(ax_ρ, inset_ax_ρ)
+
+
+    # plot the heatmaps for the viscosity
+    ylabel_μ=L"\mu(x)"
+    Label(fig[2,3], ylabel_μ,
+          tellwidth=true, tellheight=false, rotation=pi/2)
+    Label(fig[2,0], L"t",
+          tellwidth=true, tellheight=false, rotation=pi/2)
+    heatmaps = []
+    for (i,dir) in enumerate(dirs)
+        h = heatmap!(ax_μs[i], xs[i], ts[i], stack_μs[i], levels=30, colormap=:jet1,
+                     colorrange=(min_μ,max_μ))
+        push!(heatmaps, h)
+    end
+    for ax in ax_μs
+        xlims!(ax, (-0.48,0.48))
+        ax.xticks[] = [-0.4,-0.2,0.0,0.2,0.4]
+    end
+    for ax in ax_μs[2:end]
+        hideydecorations!(ax)
+        linkaxes!(ax_μs[1], ax)
+    end
+
+    Colorbar(gl_cb[1,1], heatmaps[end])
+    Legend(gl_ρ[0,1], ax_ρ, colgap=40, orientation=:horizontal, width=nothing, tellwidth=false)
+    rowgap!(gl_ρ, 1, 0)
+
+    return fig
+
+end
+
+#######################################################################
+#                             parameters                              #
+#######################################################################
+
+let
+
+rootdir = @__DIR__
+
+dirs = [
+         "avmda1_approx_covariant",
+         "avmda2_approx_covariant",
+         "avmda3_approx_covariant",
+         ]
+fig = with_theme(theme; fontsize=26) do
+    make_plot(rootdir, dirs)
+end
+display(fig)
+save(plotname(@__FILE__, "rho_mu.pdf"), fig)
+
+end