add plot_analytical to plot recipes

Author: JoshuaLampert

docs/src/basic_example.md

@@ -127,15 +127,19 @@ nothing # hide
 
 ![shoaling solution](shoaling_solution.png)
 
-By default, this will plot the bathymetry, but not the initial (analytical) solution.
+By default, this will plot the bathymetry, but neither the initial nor the analytical solution.
 
-You can adjust this by passing the boolean values `plot_bathymetry` (if `true`, always plot bathymetry in the first subplot) and `plot_initial`. Note that `plot_initial = true` will evaluate and plot the initial condition function at the same time `t` as the numerical solution being displayed (the final time by default). This means if your initial condition function represents an analytical solution, setting `plot_initial = true` will plot the analytical solution at that specific time for comparison.
+You can adjust this by passing the boolean values `plot_bathymetry` (if `true`, always plot bathymetry in the first subplot),
+`plot_initial`, and `plot_analytical`. Note that `plot_analytical = true` will evaluate and plot the initial condition function
+at the same time `t` as the numerical solution being displayed (the final time by default). This means if your initial condition
+function represents an analytical solution, setting `plot_analytical = true` will plot the analytical solution at that specific
+time for comparison. On the other hand, `plot_initial = true` will always plot the initial condition at time `t = first(tspan)`.
 
 Plotting an animation over time can, e.g., be done by the following command, which uses `step` to plot the solution at a specific time step. Here `conversion = waterheight_total` makes it so that we only look at the total water height ``\eta`` and not also the velocity ``v``. More on tutorials for plotting can be found in the chapter [Plotting Simulation Results](@ref plotting).
 
 ```@example overview
 anim = @animate for step in 1:length(sol.u)
-    plot(semi => sol, plot_initial = true, conversion = waterheight_total, step = step, xlims = (-50, 20), ylims = (-0.8, 0.1),
+    plot(semi => sol, plot_analytical = true, conversion = waterheight_total, step = step, xlims = (-50, 20), ylims = (-0.8, 0.1),
          plot_title = @sprintf "BBM-BBM equations at t = %.2f" sol.t[step])
 end
 gif(anim, "shoaling_solution.gif", fps = 25)

docs/src/plotting.md

@@ -47,7 +47,7 @@ nothing # hide
 
 ## Variable Conversion and Visualization Options
 
-The plotting system supports different variable conversions and visualization options. By default, the physical variables (returned by [`prim2phys`](@ref)) are plotted. For hyperbolic approximations like [`HyperbolicSerreGreenNaghdiEquations1D`](@ref), this means auxiliary variables are not plotted by default - only the physical variables from the limit system.
+The plotting system supports different variable conversions and visualization options. By default, the physical variables (returned by [`prim2phys`](@ref)) are plotted. For hyperbolic approximations like [`HyperbolicSerreGreenNaghdiEquations1D`](@ref), this means auxiliary variables are not plotted by default – only the physical variables from the limit system.
 
 You can plot conservative variables, specific physical quantities, and control what additional information is displayed:
 
@@ -59,8 +59,8 @@ using Plots
 
 t = 13.37 # plot solution at (roughly) t = 13.37s
 step_idx = argmin(abs.(saveat .- t)) # get the closest point to 13.37
-p1 = plot(semi => sol, conversion = prim2prim, plot_bathymetry = false,
-          suptitle = "Primitive Variables", step = step_idx)
+p1 = plot(semi => sol, conversion = prim2prim, plot_bathymetry = false, plot_analytical = true,
+          suptitle = "Primitive Variables with Analytical Solution", step = step_idx)
 p2 = plot(semi => sol, conversion = prim2cons, plot_bathymetry = false,
           suptitle = "Conservative Variables", step = step_idx)
 p3 = plot(semi => sol, conversion = waterheight_total, plot_bathymetry = true,
@@ -75,8 +75,9 @@ nothing # hide
 
 ![variable conversions](variable_conversions.png)
 
-Note that the argument `plot_initial = true` plots the initial condition evaluated at the selected time step, which means that the analytical solution is plotted if the initial condition
-function describes an exact solution that varies with time.
+Note that the argument `plot_analytical = true` plots the initial condition evaluated at the selected time step, which means
+that the analytical solution is plotted if the initial condition function describes an exact solution that varies with time.
+Similarly, `plot_initial = true` always plots the initial condition at time `t = first(tspan)`.
 
 ## Time Series Analysis at Spatial Points
 

src/visualization.jl

@@ -1,6 +1,7 @@
 struct PlotData{Conversion}
     semisol::Pair{<:Semidiscretization, <:ODESolution}
     plot_initial::Bool
+    plot_analytical::Bool
     plot_bathymetry::Bool
     conversion::Conversion
     step::Integer
@@ -13,7 +14,7 @@ struct PlotDataOverTime{RealT, Conversion}
 end
 
 @recipe function f(plotdata::PlotData)
-    @unpack semisol, plot_initial, plot_bathymetry, conversion, step = plotdata
+    @unpack semisol, plot_initial, plot_analytical, plot_bathymetry, conversion, step = plotdata
     semi, sol = semisol
     equations = semi.equations
     names = varnames(conversion, equations)
@@ -47,7 +48,7 @@ end
         names[i] in ("D", "b") && continue
 
         subplot += 1
-        if plot_initial == true
+        if plot_analytical == true
             q_exact = compute_coefficients(initial_condition, t, semi)
             data_exact = zeros(nvars, nnodes(semi))
             for j in eachnode(semi)
@@ -57,11 +58,27 @@ end
             @series begin
                 subplot --> subplot
                 linestyle --> :solid
-                label --> "initial $(names[i])"
+                label --> "analytical $(names[i])"
                 grid(semi), data_exact[i, :]
             end
         end
 
+        if plot_initial == true
+            t0 = sol.t[1]
+            q_exact = compute_coefficients(initial_condition, t0, semi)
+            data_initial = zeros(nvars, nnodes(semi))
+            for j in eachnode(semi)
+                data_initial[:, j] .= conversion(get_node_vars(q_exact, equations, j),
+                                                 equations)
+            end
+            @series begin
+                subplot --> subplot
+                linestyle --> :solid
+                label --> "initial $(names[i])"
+                grid(semi), data_initial[i, :]
+            end
+        end
+
         @series begin
             subplot --> subplot
             label --> names[i]
@@ -139,13 +156,15 @@ end
 end
 
 @recipe function f(semisol::Pair{<:Semidiscretization, <:ODESolution}; plot_initial = false,
-                   plot_bathymetry = true, conversion = prim2phys, step = -1)
-    PlotData(semisol, plot_initial, plot_bathymetry, conversion, step)
+                   plot_analytical = false, plot_bathymetry = true, conversion = prim2phys,
+                   step = -1)
+    PlotData(semisol, plot_initial, plot_analytical, plot_bathymetry, conversion, step)
 end
 
 @recipe function f(semi::Semidiscretization, sol::ODESolution; plot_initial = false,
-                   plot_bathymetry = true, conversion = prim2phys, step = -1)
-    PlotData(semi => sol, plot_initial, plot_bathymetry, conversion, step)
+                   plot_analytical = false, plot_bathymetry = true, conversion = prim2phys,
+                   step = -1)
+    PlotData(semi => sol, plot_initial, plot_analytical, plot_bathymetry, conversion, step)
 end
 
 @recipe function f(semisol::Pair{<:Semidiscretization, <:ODESolution}, x_value;

test/test_visualization.jl

@@ -6,6 +6,7 @@
                   extra_analysis_integrals = (waterheight_total, custom_integral))
     @test_nowarn plot(semi => sol)
     @test_nowarn plot!(semi => sol, plot_initial = true)
+    @test_nowarn plot!(semi => sol, plot_analytical = true)
     @test_nowarn plot(semi, sol, conversion = prim2cons, plot_bathymetry = false)
     @test_nowarn plot(semi => sol, 0.0)
     @test_nowarn plot(semi, sol, 0.0, conversion = prim2cons)