new Example295 on a sliding droplet + improvements in other examples

CHANGELOG.md

@@ -1,5 +1,12 @@
 # CHANGES
 
+## v1.8.0
+
+### Added
+  - new keyword streamlines that can be used in the plotfunction to trigger GridVisualize.streamplot, now used in Example252
+  - new Example295 on a sliding droplet
+
+
 ## v1.7.2
 
 ### Fixed

Project.toml

@@ -1,6 +1,6 @@
 name = "ExtendableFEM"
 uuid = "a722555e-65e0-4074-a036-ca7ce79a4aed"
-version = "1.7.2"
+version = "1.8.0"
 authors = ["Christian Merdon <[email protected]>", "Patrick Jaap <[email protected]>"]
 
 [deps]

docs/make.jl

@@ -48,6 +48,7 @@ function make_all(; with_examples::Bool = true, modules = :all, run_examples::Bo
                 "Example284_LevelSetMethod.jl",
                 "Example285_CahnHilliard.jl",
                 "Example290_PoroElasticity.jl",
+                "Example295_SlidingDroplet.jl",
                 "Example301_PoissonProblem.jl",
                 "Example310_DivFreeBasis.jl",
                 "Example312_PeriodicBoundary3D.jl",

examples/Example211_LshapeAdaptiveEQPoissonProblem.jl

@@ -118,6 +118,7 @@ function main(; maxdofs = 4000, μ = 1, order = 2, nlevels = 16, θ = 0.5, Plott
     ResultsH1 = zeros(Float64, 0)
     Resultsη = zeros(Float64, 0)
     sol = nothing
+    η4cell = nothing
     level = 0
     while (true)
         level += 1
@@ -132,6 +133,7 @@ function main(; maxdofs = 4000, μ = 1, order = 2, nlevels = 16, θ = 0.5, Plott
         ## evaluate eqilibration error estimator and append it to sol vector (for plotting etc.)
         local_equilibration_estimator!(sol, FETypeDual)
         η4cell = evaluate(EQIntegrator, sol)
+        @info η4cell
         push!(Resultsη, sqrt(sum(view(η4cell, 1, :))))
 
         ## calculate L2 error, H1 error, estimator, dual L2 error and write to results
@@ -168,7 +170,7 @@ function main(; maxdofs = 4000, μ = 1, order = 2, nlevels = 16, θ = 0.5, Plott
     ## print/plot convergence history
     print_convergencehistory(NDofs, [ResultsL2 ResultsH1 Resultsη]; X_to_h = X -> X .^ (-1 / 2), ylabels = ["|| u - u_h ||", "|| ∇(u - u_h) ||", "η"])
 
-    return sol, plt
+    return [sol, η4cell], plt
 end
 
 ## this function computes the local equilibrated fluxes
@@ -363,8 +365,9 @@ end
 
 generateplots = ExtendableFEM.default_generateplots(Example211_LshapeAdaptiveEQPoissonProblem, "example211.png") #hide
 function runtests() #hide
-    sol, plt = main(; maxdofs = 1000, order = 2) #hide
-    @test length(sol.entries) == 8641 #hide
+    results, plt = main(; maxdofs = 21, order = 2) #hide
+    @test length(results[1].entries) == 96 #hide
+    @test all(results[2] .≈ [0.0005948849237161765 0.02155094069716629 0.0012248396092485094 0.01840063545660059 0.02155094068517524 0.0005948849232996949])
     return nothing #hide
 end #hide
 end

examples/Example250_NSELidDrivenCavity.jl

@@ -115,9 +115,6 @@ function main(;
         FESpace{H1Pk{1, 2, order - 1}}(xgrid),
     ]
 
-    ## prepare plots
-    plt = GridVisualizer(; Plotter = Plotter, layout = (1, 2), clear = true, size = (1600, 800))
-
     ## solve by μ embedding
     step = 0
     sol = nothing
@@ -139,9 +136,6 @@ function main(;
         if step == 1
             initialize!(PE, sol)
         end
-        scalarplot!(plt[1, 1], xgrid, nodevalues(sol[1]; abs = true)[1, :]; title = "velocity (μ = $(extra_params[1]))", Plotter = Plotter)
-        vectorplot!(plt[1, 1], xgrid, eval_func_bary(PE), rasterpoints = 20, clear = false)
-        streamplot!(plt[1, 2], xgrid, eval_func_bary(PE), rasterpoints = 50, density = 2, title = "streamlines")
 
         if extra_params[1] <= μ_final
             break
@@ -150,9 +144,7 @@ function main(;
         end
     end
 
-    scalarplot!(plt[1, 1], xgrid, nodevalues(sol[1]; abs = true)[1, :]; title = "velocity (μ = $(extra_params[1]))", Plotter = Plotter)
-    vectorplot!(plt[1, 1], xgrid, eval_func_bary(PE), rasterpoints = 20, clear = false)
-    streamplot!(plt[1, 2], xgrid, eval_func_bary(PE), rasterpoints = 50, density = 2, title = "streamlines")
+    plt = plot([id(u), streamlines(u)], sol; Plotter = Plotter, rasterpoints = 30, width = 800, height = 400, title_add = " (μ = $(μ_final))")
 
     return sol, plt
 end

examples/Example284_LevelSetMethod.jl

@@ -102,7 +102,8 @@ function main(;
         t = 0
         for it in 1:Int(floor(T / τ))
             t += τ
-            ExtendableFEM.solve(PD, FES, SC; time = t)
+            @info "t = $t"
+            ExtendableFEM.solve(PD, FES, SC; time = t, timeroutputs = :hide)
             #scalarplot!(plt[1, 2], id(ϕ), sol; levels = [0.5], flimits = [-0.05, 1.05], colorbarticks = [0, 0.25, 0.5, 0.75, 1], title = "ϕ (t = $t)")
         end
     end

examples/Example290_PoroElasticity.jl

@@ -191,7 +191,7 @@ function main(; α = 0.93, E = 1.0e5, ν = 0.4, K = 1.0e-7, nrefs = 6, T = 0.5,
     for it in 1:Int(floor(T / τ))
         t += τ
         @info "t = $t"
-        ExtendableFEM.solve(PD, FES, SC; time = t)
+        ExtendableFEM.solve(PD, FES, SC; time = t, timeroutputs = :hide)
     end
 
     ## error calculation

examples/Example295_SlidingDroplet.jl

@@ -0,0 +1,235 @@
+#=
+
+# 295 : Sliding Droplet
+([source code](@__SOURCE_URL__))
+
+This problem describes the  motion of a liquid droplet sliding down a solid surface under gravity,
+while maintaining its shape due to  surface tension and experiencing  slip at the solid boundary.
+The problem is solved using an ALE (Arbitrary Lagrangian-Eulerian) approach where the mesh is moved
+according to the extension velocity ``\mathbf{w}``.
+
+It seeks a velocity ``\mathbf{u}``, a pressure ``p`` such that
+```math
+\begin{aligned}
+	- 2\mu \varepsilon(\mathbf{u}) + \nabla p & = \mathbf{g} \\
+	\mathrm{div} \mathbf{u} & = 0
+\end{aligned}
+```
+on a moving domain ``\Omega(t)`` (a sliding droplet) with boundary conditions
+```math
+\begin{aligned}
+	\mathbf{u} \cdot \mathbf{n} & = 0 && \quad \text{no penetration along liquid-solid interface} \\
+	\mathbf{u} \cdot \mathbf{t} & = \mu_{sl}^{-1} \sigma_{nt} && \quad \text{Navier slip condition along liquid-solid interface} \\
+	\sigma \mathbf{n} & = \gamma \kappa \mathbf{n} && \quad \text{surface tension at liquid-air interface} \\
+	\sigma \mathbf{n} & = \gamma_{sl} \kappa \mathbf{n} && \quad \text{surface tension at liquid-solid interface}
+\end{aligned}
+```
+where ``\sigma = -p I + 2\mu \varepsilon(\mathbf{u})`` is the stress tensor, ``\kappa`` is the curvature,
+``\mathbf{n}`` is the normal vector and ``\mathbf{t}`` is the tangent vector.
+
+Then, it computes the displacement ``\mathbf{w}`` from ``\mathbf{u}`` by Laplace smoothing to prevent mesh folding.
+
+The weak formulation chracterizes ``(\mathbf{u},p,\mathbf{w}) \in V \times Q \times W`` by
+```math
+\begin{aligned}
+	(2\mu \varepsilon(\mathbf{u}), \varepsilon(\mathbf{v})) - (\mathrm{div} \mathbf{v}, p)
+	+ \gamma \tau (\nabla_s \mathbf{u}, \nabla_s \mathbf{v})_{\Gamma_{air}} & \\
+	+ \gamma_{sl} \tau (\nabla_s \mathbf{u}, \nabla_s \mathbf{v})_{\Gamma_{solid}}
+	+ \mu_{sl} (\mathbf{u} \cdot \mathbf{t}, \mathbf{v} \cdot \mathbf{t})_{\Gamma_{solid}}
+	& = (\mathbf{g}, \mathbf{v}) - \gamma (\nabla_s \mathbf{x}, \nabla_s \mathbf{v})_{\Gamma_{air}}
+	- \gamma_{sl} (\nabla_s \mathbf{x}, \nabla_s \mathbf{v})_{\Gamma_{solid}}
+	&& \quad \text{for all } \mathbf{v} \in V,\\
+(\mathrm{div} \mathbf{u}, q) & = 0 && \quad \text{for all } q \in Q,\\
+(\nabla \mathbf{w}, \nabla \mathbf{z}) + \epsilon^{-1} (\mathbf{w} \cdot \mathbf{n}, \mathbf{z} \cdot \mathbf{n})_{\partial\Omega}
+	& = \epsilon^{-1} (\mathbf{u} \cdot \mathbf{n}, \mathbf{z} \cdot \mathbf{n})_{\partial\Omega}
+	&& \quad \text{for all } \mathbf{z} \in W.
+\end{aligned}
+```
+
+When the droplet speed becomes constant im time, a travelling solution is reached.
+For the default parameters the result looks like this:
+
+![](example295.png)
+
+
+!!! reference
+
+	"Resolving the microscopic hydrodynamics at the moving contact line",\
+	A. K. Giri, P. Malgaretti, D. Peschka, and M. Sega,\
+	Phys. Rev. Fluids 7 (2022),\
+	[>Journal-Link<](https://doi.org/10.1103/PhysRevFluids.7.L102001)
+=#
+
+module Example295_SlidingDroplet
+
+using ExtendableFEM
+using ExtendableGrids
+using ExtendableSparse
+using LinearAlgebra
+using Triangulate
+using SimplexGridFactory
+using GridVisualize
+using Test #hide
+
+## gavity
+function g!(result, qpinfo)
+    result[1] = 1.0
+    return result[2] = 0.0
+end
+
+function surface_tension!(result, input, qpinfo)
+    t1, t2 = qpinfo.normal[2], -qpinfo.normal[1] # tangent
+    p1 = (input[1] * t1 + input[2] * t2)
+    p2 = (input[3] * t1 + input[4] * t2)
+    result[1] = p1 * t1
+    result[2] = p1 * t2
+    result[3] = p2 * t1
+    return result[4] = p2 * t2
+end
+
+function initial_grid(nref; radius = 1)
+    builder = SimplexGridBuilder(Generator = Triangulate)
+    n = 2^(nref + 3)
+    maxvol = 2.0^(-nref - 3)
+    points = [point!(builder, radius * sin(t), radius * cos(t)) for t in range(-π / 2, π / 2, length = n)]
+
+    facetregion!(builder, 1)
+    for i in 1:(n - 1)
+        facet!(builder, points[i], points[i + 1])
+    end
+    facetregion!(builder, 2)
+    facet!(builder, points[end], points[1])
+
+    return simplexgrid(builder, maxvolume = maxvol), [1, length(points)]
+end
+
+function main(;
+        order = 2,
+        g = 2,             ## gravity factor (in right direction)
+        μ = 0.7,              ## bulk viscosity
+        μ_sl = 0.1,         ## coefficient for Navier slip condition
+        μ_dyn = 0,          ## dynamic contact angle
+        γ_la = 1,           ## surface tension coefficient at liquid <> air interface
+        γ_sl = 0,           ## surface tension coefficient at liquid <> solid interface
+        nsteps = 1000,        ## ALE steps
+        τ = 0.1,         ## ALE stepsize
+
+        nrefs = 4,
+        Plotter = nothing, kwargs...
+    )
+
+    ## grid
+    xgrid, triple_nodes = initial_grid(nrefs)
+
+    ## Stokes problem description
+    PD = ProblemDescription("Stokes problem")
+    u = Unknown("u"; name = "velocity")
+    p = Unknown("p"; name = "pressure")
+    x = Unknown("x"; name = "x")
+    w = Unknown("w"; name = "extension")
+    q = Unknown("q"; name = "Lagrange multiplier for normal flux")
+    v = Unknown("v"; name = "comoving velocity")
+
+    assign_unknown!(PD, u)
+    assign_unknown!(PD, p)
+
+    ## BLFs a(u,v), b(u,p)
+    assign_operator!(PD, BilinearOperator([εV(u, 1.0)]; factor = 2 * μ, kwargs...))
+    assign_operator!(PD, BilinearOperatorDG(surface_tension!, [grad(u)]; entities = ON_BFACES, factor = γ_la * τ, regions = [1], kwargs...))
+    assign_operator!(PD, BilinearOperatorDG(surface_tension!, [grad(u)]; entities = ON_BFACES, factor = γ_sl * τ, regions = [2], kwargs...))
+    assign_operator!(PD, BilinearOperator([id(p)], [div(u)]; transposed_copy = 1, factor = -1, kwargs...))
+    assign_operator!(PD, BilinearOperator([id(u)]; entities = ON_BFACES, regions = [2], factor = μ_sl, kwargs...))
+    if abs(μ_dyn) > 0
+        assign_operator!(PD, CallbackOperator(ctriple_junction_kernel!(triple_nodes, μ_dyn), [u]; kwargs...))
+    end
+
+    ## RHS
+    assign_operator!(PD, LinearOperator(g!, [id(u)]; factor = g, kwargs...))
+    assign_operator!(PD, LinearOperatorDG(surface_tension!, [grad(u)], [grad(x)]; entities = ON_BFACES, quadorder = 2, factor = -γ_la, regions = [1], kwargs...))
+    assign_operator!(PD, LinearOperatorDG(surface_tension!, [grad(u)], [grad(x)]; entities = ON_BFACES, quadorder = 2, factor = -γ_sl, regions = [2], kwargs...))
+
+    ## boundary conditions
+    assign_operator!(PD, HomogeneousBoundaryData(u; regions = [2], mask = [0, 1, 1]))
+
+    ## ALE problem description
+    PDALE = ProblemDescription("ALE problem")
+    assign_unknown!(PDALE, w)
+    assign_unknown!(PDALE, q)
+    assign_operator!(PDALE, BilinearOperator([grad(w)]; kwargs...))
+    assign_operator!(PDALE, BilinearOperator([normalflux(w)], [id(q)]; transposed_copy = 1, regions = [1, 2], entities = ON_BFACES))
+    assign_operator!(PDALE, LinearOperator([id(q)], [normalflux(u)]; regions = [1, 2], entities = ON_BFACES))
+    assign_operator!(PDALE, HomogeneousBoundaryData(w; regions = [2], mask = [0, 1, 1]))
+
+    ## prepare FESpace and solution vector
+    FES = [FESpace{H1Pk{2, 2, order}}(xgrid), FESpace{H1Pk{1, 2, order - 1}}(xgrid), FESpace{H1Pk{1, 1, order - 1}, ON_BFACES}(xgrid), FESpace{H1P1{2}}(xgrid)]
+    sol = FEVector([FES[1], FES[2]]; tags = [u, p])
+    append!(sol, FES[1]; tag = w)
+    append!(sol, FES[3]; tag = q)
+    append!(sol, FES[4]; tag = x)
+
+    SC = SolverConfiguration(PD, FES[[1, 2]]; init = sol, maxiterations = 1, verbosity = -1, timeroutputs = :none, kwargs...)
+    SCALE = SolverConfiguration(PDALE, FES[[1, 3]]; init = sol, maxiterations = 1, verbosity = -1, timeroutputs = :none, kwargs...)
+
+    ## prepare plot
+
+    ## time loop
+    time = 0.0
+    v0 = nothing
+    v0_old = 0
+    for step in 1:nsteps
+        @info "STEP $step, τ = $τ time = $(Float16(time))"
+
+        ## compute x for computation of the tangential identity grad(x)
+        interpolate!(sol[x], (result, qpinfo) -> (result .= qpinfo.x))
+
+        ## solve Stokes problem
+        solve(PD, FES[[1, 2]], SC; kwargs...)
+
+        ## solve ALE problem
+        solve(PDALE, FES[[1, 3]], SCALE; kwargs...)
+
+        ## displace mesh
+        time += τ
+        displace_mesh!(xgrid, sol[w]; magnify = τ)
+
+        ## calculate final droplet speed
+        v0_old = v0
+        v0 = sum(sol.entries[1:FES[1].coffset]) / FES[1].coffset
+
+        #if step > 1
+        #    @info 1e-10/abs(v0 - v0_old)
+        #    τ = min(10, max(1e-3, 1e-4/abs(v0 - v0_old)))
+        #end
+        @info "droplet_speed = $v0"
+    end
+
+    ## compute comoving velocity (= u - v0)
+    append!(sol, FES[1]; tag = v)
+    view(sol[v]) .= view(sol[u])
+    sol[v][1:FES[1].coffset] .-= v0
+
+    plt = plot([grid(u), id(u), id(p), streamlines(v)], sol; Plotter = Plotter, levels = 0, colorbarticks = 7, rasterpoints = 30)
+
+    return sol, plt
+end
+
+## slip condition for triplet junctions (air-surface-solid)
+function triple_junction_kernel!(triple_nodes, μ_dyn)
+    factors = [1, 1, 0, 0]
+    return function closure(A, b, args; assemble_matrix = true, kwargs...)
+        FES = args[1].FES
+        triple_dofs = copy(triple_nodes)
+        append!(triple_dofs, triple_nodes .+ FES.coffset)
+        if assemble_matrix
+            for j in 1:length(triple_dofs)
+                dof = triple_dofs[j]
+                A[dof, dof] += factors[j] * μ_dyn
+            end
+            flush!(A)
+        end
+        return nothing
+    end
+end
+
+generateplots = ExtendableFEM.default_generateplots(Example295_SlidingDroplet, "example295.png") #hide
+end # module

src/ExtendableFEM.jl

@@ -66,7 +66,7 @@ using ExtendableSparse: ExtendableSparse, ExtendableSparseMatrix, flush!,
     rawupdateindex!
 using ForwardDiff: ForwardDiff
 using GridVisualize: GridVisualize, GridVisualizer, gridplot!, reveal, save,
-    scalarplot!, vectorplot!
+    scalarplot!, vectorplot!, streamplot!
 using LinearAlgebra: LinearAlgebra, copyto!, isposdef, mul!, norm
 using LinearSolve: LinearSolve, LinearProblem, UMFPACKFactorization, deleteat!,
     init, solve
@@ -114,7 +114,7 @@ export print_table
 
 include("unknowns.jl")
 export Unknown
-export grid, dofgrid
+export grid, dofgrid, streamlines
 export id, grad, hessian, div, normalflux, tangentialflux, Δ, apply, curl1, curl2, curl3, laplace, tangentialgrad, symgrad_voigt, εV
 
 include("operators.jl")

src/plots.jl

@@ -62,6 +62,14 @@ function plot!(
             gridplot!(p[row, col], sol[op[1]].FES.xgrid; linewidth = linewidth, kwargs...)
         elseif op[2] == "dofgrid"
             gridplot!(p[row, col], sol[op[1]].FES.dofgrid; linewidth = linewidth, kwargs...)
+        elseif op[2] == "streamlines"
+            if typeof(op[1]) <: Unknown
+                title = String(op[1].identifier)
+            else
+                title = "$(sol[op[1]].name)"
+            end
+            PE = PointEvaluator([apply(op[1], Identity)], sol)
+            streamplot!(p[row, col], sol[op[1]].FES.dofgrid, eval_func_bary(PE); rasterpoints = rasterpoints, title = title * " (streamlines)" * title_add, kwargs...)
         else
             ncomponents = get_ncomponents(sol[op[1]])
             edim = size(sol[op[1]].FES.xgrid[Coordinates], 1)
@@ -138,7 +146,7 @@ function plot(ops, sol; add = 0, Plotter = nothing, ncols = min(2, length(ops) +
     for op in ops
         ncomponents = get_ncomponents(sol[op[1]])
         edim = size(sol[op[1]].FES.xgrid[Coordinates], 1)
-        if !(op[2] in ["grid", "dofgrid"])
+        if !(op[2] in ["grid", "dofgrid", "streamlines"])
             resultdim = Length4Operator(op[2], edim, ncomponents)
             if resultdim > 1 && do_abs == false
                 nplots += resultdim - 1