improved convergencehistory! plotter, export of some low level stuff, some other small improvements

Author: chmerdon

Project.toml

@@ -27,8 +27,8 @@ DiffResults = "1"
 CommonSolve = "0.2"
 DocStringExtensions = "0.8,0.9"
 ExtendableGrids = "0.9.16"
-ExtendableSparse = "1.0.1"
-ForwardDiff = "^0.10"
+ExtendableSparse = "1.2"
+ForwardDiff = "^0.10.35"
 GridVisualize = "0.5.2,0.6.1,1"
 LinearSolve = "1.37"
 SparseDiffTools = "^1.19,2"

examples/Example203_ReactionConvectionDiffusion2D.jl

@@ -127,9 +127,9 @@ function main(; verbosity = 0, Plotter = nothing, ν = 1e-5, τ = 1e-2, nlevels
     ## plot
     p = GridVisualizer(; Plotter = Plotter, layout = (1,3), clear = true, resolution = (1500,500))
     scalarplot!(p[1,1], xgrid, view(nodevalues(Solution[1]),1,:), levels = 7, title = "u_h")
-    scalarplot!(p[1,2], xgrid, view(nodevalues(Solution[1], Gradient; abs = true),1,:), levels = 0, colorbarticks = 9, title = "∇u_h (abs + quiver)")
+    scalarplot!(p[1,2], xgrid, view(nodevalues(Solution[1], Gradient; abs = true),1,:), levels = 7, colorbarticks = 9, title = "∇u_h (abs + quiver)")
     vectorplot!(p[1,2], xgrid, evaluate(PointEvaluator(Solution[1], Gradient)), vscale = 0.8, clear = false)
-    convergencehistory!(p[1,3], NDofs, Results; add_h_powers = [order,order+1], X_to_h = X -> X.^(-1/2), legend = :lb, fontsize = 20, ylabels = ["|| u - u_h ||", "|| u - Iu ||", "|| ∇(u - u_h) ||", "|| ∇(u - Iu) ||"])
+    convergencehistory!(p[1,3], NDofs, Results; add_h_powers = [order,order+1], X_to_h = X -> X.^(-1/2), legend = :lb, fontsize = 20, ylabels = ["|| u - u_h ||", "|| u - Iu ||", "|| ∇(u - u_h) ||", "|| ∇(u - Iu) ||"], limits = (1e-8,1e-1))
 
     ## print convergence history
     print_convergencehistory(NDofs, Results; X_to_h = X -> X.^(-1/2), ylabels = ["|| u - u_h ||", "|| u - Iu ||", "|| ∇(u - u_h) ||", "|| ∇(u - Iu) ||"])

examples/Example219_AxisymmetricStokes25D.jl

@@ -41,7 +41,8 @@ a(u,v) := \int_{\Omega} \left( \nabla u : \nabla v + r^{-2} u_r v_r \right) r dr
 b(q,v) := \int_{\Omega} q \left( \mathrm{div}(v) + r^{-1} u_r \right) r dr dz
 \end{aligned}
 ```
-where the usual Cartesian differential operators can be used.
+where the usual Cartesian differential operators can be used. The factor ``2\pi`` from the
+integral over the rotation angle drops out on both sides.
 
 
 =#
@@ -52,16 +53,6 @@ using GradientRobustMultiPhysics
 using ExtendableGrids
 using GridVisualize
 
-## data for axisymmetric Hagen-Poiseuille flow
-function exact_pressure!(μ)
-    function closure(result,x)
-        result[1] = 4*μ*(1-x[1])
-    end
-end
-function exact_velocity!(result,x)
-    result[1] = (1.0-x[2]^2);
-    result[2] = 0.0;
-end
 
 ## custom bilinearform a for the axisymmetric Stokes problem
 function ASLaplaceOperator(μ)
@@ -108,14 +99,36 @@ function ASPressureOperatorTransposed()
 end
 
 ## everything is wrapped in a main function
-function main(; verbosity = 0, Plotter = nothing, μ = 1)
+## problem = 1 is Hagen-Poiseuille
+## problem = 2 is stagnation point flow without solid surface
+## μ = viscosity
+## k = parameter for problem = 2
+function main(; problem = 1, k = 1, verbosity = 0, Plotter = nothing, μ = 1)
 
     ## set log level
     set_verbosity(verbosity)
 
+    ## boundary/comparison data for axisymmetric Hagen-Poiseuille flow
+    function exact_pressure!(result, x)
+        if problem == 1 # Hagen-Poiseuille
+            result[1] = 4*μ*(1-x[1])
+        elseif problem == 2 # stagnation point flow without solid surface
+            result[1] = k
+        end
+    end
+    function exact_velocity!(result, x)
+        if problem == 1 # Hagen-Poiseuille
+            result[1] = (1.0-x[2]^2);
+            result[2] = 0.0;
+        elseif problem == 2 # stagnation point flow without solid surface
+            result[1] = -2*k*x[1]
+            result[2] = k*x[2]
+        end
+    end
+
     ## negotiate data functions to the package
     u = DataFunction(exact_velocity!, [2,2]; name = "u", dependencies = "X", bonus_quadorder = 2)
-    p = DataFunction(exact_pressure!(μ), [1,2]; name = "p", dependencies = "X", bonus_quadorder = 1)
+    p = DataFunction(exact_pressure!, [1,2]; name = "p", dependencies = "X", bonus_quadorder = 1)
 
     ## grid
     xgrid = uniform_refine(grid_unitsquare(Triangle2D), 5);
@@ -132,9 +145,14 @@ function main(; verbosity = 0, Plotter = nothing, μ = 1)
     add_operator!(Problem, [2,1], ASPressureOperatorTransposed())
 
     ## boundary conditions
-    add_boundarydata!(Problem, 1, [4], BestapproxDirichletBoundary; data = u)
-    add_boundarydata!(Problem, 1, [3], HomogeneousDirichletBoundary)
-    add_boundarydata!(Problem, 1, [1], HomogeneousDirichletBoundary; mask = (0,1)) 
+    if problem == 1
+        add_boundarydata!(Problem, 1, [3,4], InterpolateDirichletBoundary; data = u)
+        add_boundarydata!(Problem, 1, [1], HomogeneousDirichletBoundary; mask = (0,1)) 
+    elseif problem == 2
+        add_boundarydata!(Problem, 1, [2], InterpolateDirichletBoundary; data = u)
+        add_boundarydata!(Problem, 1, [4], HomogeneousDirichletBoundary; mask = (1,0))
+        add_boundarydata!(Problem, 1, [1], HomogeneousDirichletBoundary; mask = (0,1))
+    end
     @show Problem
 
     ## generate FESpaces
@@ -152,8 +170,9 @@ function main(; verbosity = 0, Plotter = nothing, μ = 1)
 
     ## plot
     p = GridVisualizer(; Plotter = Plotter, layout = (1,2), clear = true, resolution = (1000,500))
-    scalarplot!(p[1,1],xgrid,view(nodevalues(Solution[1]; abs = true),1,:), levels = 9)
-    scalarplot!(p[1,2],xgrid,view(nodevalues(Solution[2]),1,:), levels = 11, title = "p_h")
+    scalarplot!(p[1,1],xgrid,view(nodevalues(Solution[1]; abs = true),1,:), levels = 6, xlabel = "z", ylabel = "r")
+    vectorplot!(p[1,1],xgrid,evaluate(PointEvaluator(Solution[1], Identity)), spacing = 0.1, clear = false)
+    scalarplot!(p[1,2],xgrid,view(nodevalues(Solution[2]),1,:), levels = 5, xlabel = "z", ylabel = "r", title = "p_h")
 end
 
 end

src/GradientRobustMultiPhysics.jl

@@ -92,7 +92,7 @@ export nodevalues_view
 export FEVectorBlock, FEVector
 export dot, norm, norms
 export FEMatrixBlock, FEMatrix, _addnz
-export fill!, addblock!, addblock_matmul!, lrmatmul, add!, apply_penalties!
+export fill!, addblock!, addblock_matmul!, lrmatmul, mul!, add!, apply_penalties!
 
 export get_reconstruction_matrix
 
@@ -120,9 +120,10 @@ include("assemblypatterns.jl")
 export AssemblyPatternType, AssemblyPreparations
 export ItemIntegrator, L2ErrorIntegrator, L2NormIntegrator, L2DifferenceIntegrator
 export DiscreteLinearForm
-export APT_BilinearForm, APT_SymmetricBilinearForm, APT_LumpedBilinearForm
+export AssemblyPattern, APT_BilinearForm, APT_SymmetricBilinearForm, APT_LumpedBilinearForm
 export DiscreteBilinearForm, DiscreteSymmetricBilinearForm, DiscreteLumpedBilinearForm
 export DiscreteNonlinearForm
+export assemble!, full_assemble!
 export prepare_assembly!
 export assemble!, evaluate!, evaluate
 export AssemblyManager, update_assembly!

src/plots.jl

@@ -1,4 +1,4 @@
-function convergencehistory!(target, X, Y; add_h_powers = [], X_to_h = X -> X, colors = [:blue, :green, :red, :magenta, :lightblue], title = "convergence history", legend = :best, ylabel = "", ylabels = [], xlabel = "ndofs", clear = true, fontsize = 14)
+function convergencehistory!(target, X, Y; add_h_powers = [], X_to_h = X -> X, colors = [:blue, :green, :red, :magenta, :lightblue], title = "convergence history", legend = :best, ylabel = "", ylabels = [], xlabel = "ndofs", markershape = :circle, markevery = 1, clear = true, args...)
     for j = 1 : size(Y,2)
         Xk = []
         Yk = []
@@ -13,15 +13,15 @@ function convergencehistory!(target, X, Y; add_h_powers = [], X_to_h = X -> X, c
         else
             label = "Data $j"
         end
-        scalarplot!(target, simplexgrid(Xk), Yk; fontsize = fontsize, xlabel = xlabel, ylabel = ylabel, color = length(colors) >= j ? colors[j] : :black, clear = j == 1 ? clear : false, xscale = :log, yscale = :log, markershape = :circle, markevery = 1, label = label, legend = legend, title = title)
+        scalarplot!(target, simplexgrid(Xk), Yk; xlabel = xlabel, ylabel = ylabel, color = length(colors) >= j ? colors[j] : :black, clear = j == 1 ? clear : false, markershape = markershape, markevery = markevery, xscale = :log, yscale = :log, label = label, legend = legend, title = title, args...)
     end
     for p in add_h_powers
         label = "h^$p"
-        scalarplot!(target, simplexgrid(X), X_to_h(X).^p; fontsize = fontsize, linestyle = :dot, xlabel = xlabel, ylabel = ylabel, color = :gray, clear = false, markershape = :none, xscale = :log, yscale = :log, label = label, legend = legend, title = title)
+        scalarplot!(target, simplexgrid(X), X_to_h(X).^p; linestyle = :dot, xlabel = xlabel, ylabel = ylabel, color = :gray, clear = false, markershape = :none, xscale = :log, yscale = :log, label = label, legend = legend, title = title, args...)
     end
 end
 
-function plot_convergencehistory(X, Y; Plotter = nothing, resolution = (800,600), add_h_powers = [], X_to_h = X -> X, colors = [:blue, :green, :red, :magenta, :lightblue], title = title, legend = :best, ylabel = "", ylabels = [], xlabel = "ndofs", clear = true, fontsize = 14)
+function plot_convergencehistory(X, Y; Plotter = nothing, resolution = (800,600), add_h_powers = [], X_to_h = X -> X, colors = [:blue, :green, :red, :magenta, :lightblue], legend = :best, ylabel = "", ylabels = [], xlabel = "ndofs", clear = true, args...)
     p=GridVisualizer(; Plotter = Plotter, layout = (1,1), clear = true, resolution = resolution)
-    plot_convergencehistory!(p[1,1], X, Y; add_h_powers = add_h_powers, X_to_h = X_to_h, colors = colors, title = title, legend = legend, ylabel = ylabel, ylabels = ylabels, xlabel = xlabel, clear = clear, fontsize = fontsize)
+    plot_convergencehistory!(p[1,1], X, Y; add_h_powers = add_h_powers, X_to_h = X_to_h, colors = colors, legend = legend, ylabel = ylabel, ylabels = ylabels, xlabel = xlabel, clear = clear, args...)
 end