diff --git a/CHANGELOG.md b/CHANGELOG.md index 59c09649..229e365a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,14 @@ # CHANGES +## v.1.11.0 + +### Changed + - `UnicodePlots` is now the default Plotter in most 1d and 2d examples + - functions for plotting of FEVectorBlocks moved to ExtendableFEMBase#v1.6 + +### Added + - manual assembly! calls for LinearOperatorFromVector and LinearOperatorFromMatrix + ## v.1.10.5 diff --git a/Project.toml b/Project.toml index 530ad6b2..3deaeb01 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "ExtendableFEM" uuid = "a722555e-65e0-4074-a036-ca7ce79a4aed" -version = "1.10.5" +version = "1.11.0" authors = ["Christian Merdon ", "Patrick Jaap "] [deps] @@ -36,12 +36,13 @@ DifferentiationInterface = "0.7.4" DocStringExtensions = "0.8,0.9" ExampleJuggler = "2.2.1" ExplicitImports = "1" -ExtendableFEMBase = "1.3.0" -ExtendableGrids = "1.10.3" +ExtendableFEMBase = "1.6.0" +ExtendableGrids = "1.17.0" ExtendableSparse = "1.5.3, 2" ForwardDiff = "0.10.35,1" -GridVisualize = "1.8.1" +GridVisualize = "1.20" IncompleteLU = "0.2.1" +KrylovKit = "0.10.2" LinearAlgebra = "1.9" LinearSolve = "2, 3" Metis = "1.5.0" @@ -55,6 +56,7 @@ SparseConnectivityTracer = "1.0.0" SparseMatrixColorings = "0.4.21" StaticArrays = "1.9.13" Symbolics = "4.2,5,6, 7" +Term = "2" Test = "1" TetGen = "2.0.0" TimerOutputs = "0.5.29" @@ -67,15 +69,18 @@ Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" ExampleJuggler = "3bbe58f8-ed81-4c4e-a134-03e85fcf4a1a" ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7" IncompleteLU = "40713840-3770-5561-ab4c-a76e7d0d7895" +KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" Metis = "2679e427-3c69-5b7f-982b-ece356f1e94b" OrdinaryDiffEqRosenbrock = "43230ef6-c299-4910-a778-202eb28ce4ce" OrdinaryDiffEqSDIRK = "2d112036-d095-4a1e-ab9a-08536f3ecdbf" SimplexGridFactory = "57bfcd06-606e-45d6-baf4-4ba06da0efd5" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" +Term = "22787eb5-b846-44ae-b979-8e399b8463ab" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" TetGen = "c5d3f3f7-f850-59f6-8a2e-ffc6dc1317ea" Triangulate = "f7e6ffb2-c36d-4f8f-a77e-16e897189344" +UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228" [targets] -test = ["Aqua", "ExampleJuggler", "ExplicitImports", "IncompleteLU", "Metis", "OrdinaryDiffEqRosenbrock", "OrdinaryDiffEqSDIRK", "SimplexGridFactory", "StaticArrays", "Symbolics", "Test", "TetGen", "Triangulate"] +test = ["Aqua", "ExampleJuggler", "ExplicitImports", "IncompleteLU", "KrylovKit", "Metis", "OrdinaryDiffEqRosenbrock", "OrdinaryDiffEqSDIRK", "SimplexGridFactory", "StaticArrays", "Symbolics", "Test", "Term", "TetGen", "Triangulate", "UnicodePlots"] diff --git a/README.md b/README.md index 1416db31..f5adb0dd 100644 --- a/README.md +++ b/README.md @@ -54,8 +54,8 @@ FES = FESpace{FEType}(xgrid) sol = solve(PD, [FES]; plot = true, timeroutputs = :hide) # Plot the solution -import PythonPlot -plot([id(u), grid(u)], sol; Plotter = PythonPlot) +import UnicodePlots, Term +plot([id(u), grid(u)], sol; Plotter = UnicodePlots) ``` ## Running examples from documentation @@ -64,7 +64,7 @@ In the [documentation](https://wias-pdelib.github.io/ExtendableFEM.jl/stable/ind Each of these examples is implemented as a module that needs to be included first. Afterwards the main file of the module can be run. Usually the main function has a Plotter argument that can be used to toggle a plot of the solution with an -(already installed) backend of your choice (e.g. `Plotter=PythonPlot`, `GLMakie`, `Plots` or others supported by [GridVisualize.jl](https://github.com/WIAS-PDELib/GridVisualize.jl)). +(already installed) backend of your choice (e.g. `Plotter=UnicodePlots`, `PythonPlot`, `GLMakie`, `Plots` or others supported by [GridVisualize.jl](https://github.com/WIAS-PDELib/GridVisualize.jl)). Some examples need several further dependencies. To ensure an environment were every dependency is installed, one can use the test environment via the package [TestEnv](https://github.com/JuliaTesting/TestEnv.jl). The following script runs Example201: ```julia diff --git a/docs/Project.toml b/docs/Project.toml index 817f319f..556b795a 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -23,6 +23,7 @@ SimplexGridFactory = "57bfcd06-606e-45d6-baf4-4ba06da0efd5" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" +Term = "22787eb5-b846-44ae-b979-8e399b8463ab" TetGen = "c5d3f3f7-f850-59f6-8a2e-ffc6dc1317ea" Triangulate = "f7e6ffb2-c36d-4f8f-a77e-16e897189344" UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228" diff --git a/examples/Example103_BurgersEquation.jl b/examples/Example103_BurgersEquation.jl deleted file mode 100644 index f9e522fa..00000000 --- a/examples/Example103_BurgersEquation.jl +++ /dev/null @@ -1,118 +0,0 @@ -#= - -# 103 : Burger's Equation -([source code](@__SOURCE_URL__)) - -This example solves the Burger's equation -```math -\begin{aligned} -u_t - \mu \Delta u + \mathrm{div} f(u) & = 0 -\end{aligned} -``` -with periodic boundary conditions. -This script demonstrates how a time-dependent PDE can be solved -with DifferentialEquations or by a manual implicit Euler scheme. - -The initial condition and the final solution for the default parameters looks like this: - -![](example103.png) - -=# - -module Example103_BurgersEquation - -using ExtendableFEM -using ExtendableGrids -using OrdinaryDiffEqRosenbrock -using Test #hide - -## nonlinear kernel, i.e. f(u) -function f!(result, input, qpinfo) - result[1] = input[1]^2 / 2 - return nothing -end - -## initial condition -function uinit!(result, qpinfo) - result[1] = abs(qpinfo.x[1]) < 0.5 ? 1 : 0 - return nothing -end - -## everything is wrapped in a main function -function main(; - ν = 0.01, - h = 0.005, - T = 2, - order = 2, - τ = 0.001, - Plotter = nothing, - use_diffeq = true, - solver = Rosenbrock23(), - kwargs... - ) - - ## load mesh and exact solution - xgrid = simplexgrid(-2:h:2) - - ## generate empty PDEDescription for three unknowns (h, u) - PD = ProblemDescription("Burger's Equation") - u = Unknown("u"; name = "u") - assign_unknown!(PD, u) - assign_operator!(PD, NonlinearOperator(f!, [grad(u)], [id(u)]; bonus_quadorder = 2)) - assign_operator!(PD, BilinearOperator([grad(u)]; store = true, factor = ν)) - assign_operator!(PD, CombineDofs(u, u, [1], [num_nodes(xgrid)], [1.0]; kwargs...)) - - ## prepare solution vector and initial data - FES = FESpace{H1Pk{1, 1, order}}(xgrid) - sol = FEVector(FES; tags = PD.unknowns) - interpolate!(sol[u], uinit!) - - ## init plotter and plot u0 - plt = plot([id(u), id(u)], sol; Plotter = Plotter, title_add = " (t = 0)") - - ## generate mass matrix - M = FEMatrix(FES) - assemble!(M, BilinearOperator([id(1)]; lump = 2)) - - if (use_diffeq) - ## generate DifferentialEquations.ODEProblem - prob = ExtendableFEM.generate_ODEProblem(PD, FES, (0.0, T); init = sol, mass_matrix = M) - - ## solve ODE problem - de_sol = solve(prob, solver, abstol = 1.0e-6, reltol = 1.0e-3, dt = τ, dtmin = 1.0e-6, adaptive = true) - @info "#tsteps = $(length(de_sol.u))" - - ## extract final solution - sol.entries .= de_sol.u[end] - else - ## add backward Euler time derivative - assign_operator!(PD, BilinearOperator(M, [u]; factor = 1 / τ, kwargs...)) - assign_operator!(PD, LinearOperator(M, [u], [u]; factor = 1 / τ, kwargs...)) - - ## generate solver configuration - SC = SolverConfiguration(PD, FES; init = sol, maxiterations = 1, kwargs...) - - ## iterate tspan - t = 0 - for it in 1:Int(floor(T / τ)) - t += τ - ExtendableFEM.solve(PD, FES, SC; time = t, verbosity = -1, timeroutputs = :hide) - end - show(timeroutputs(SC)) - end - - ## plot final state - plot!(plt, [id(u)], sol; keep = 1, title_add = " (t = $T)") - - return sol, plt -end - -generateplots = ExtendableFEM.default_generateplots(Example103_BurgersEquation, "example103.png") #hide -function runtests() #hide - sol, plt = main(; h = 0.01, τ = 0.01, T = 1, use_diffeq = false) #hide - @test maximum(sol.entries) ≈ 0.991626266834898 #hide - sol, plt = main(; h = 0.01, τ = 0.01, T = 1, use_diffeq = true) #hide - @test maximum(sol.entries) ≈ 0.9974553634892214 #hide - return nothing #hide -end #hide -end diff --git a/examples/Example105_NonlinearPoissonEquation.jl b/examples/Example105_NonlinearPoissonEquation.jl deleted file mode 100644 index a7fef56d..00000000 --- a/examples/Example105_NonlinearPoissonEquation.jl +++ /dev/null @@ -1,82 +0,0 @@ -#= - -# 105 : Nonlinear Poisson Equation -([source code](@__SOURCE_URL__)) - -This examples solves the nonlinear Poisson problem -```math -\begin{aligned} -- \epsilon \partial^2 u / \partial x^2 + e^u - e^{-u} & = f && \text{in } \Omega -\end{aligned} -``` -where -```math -f(x) = \begin{cases} -1 & x \geq 0.5, --1 & x < 0.5. -\end{cases} -``` -on the domain ``\Omega := (0,1)`` with Dirichlet boundary conditions ``u(0) = 0`` and ``u(1) = 1``. - -The solution looks like this: - -![](example105.png) - -=# - -module Example105_NonlinearPoissonEquation - -using ExtendableFEM -using ExtendableGrids -using Test #hide - -## rigt-hand side data -function f!(result, qpinfo) - result[1] = qpinfo.x[1] < 0.5 ? -1 : 1 - return nothing -end -## boundary data -function boundary_data!(result, qpinfo) - result[1] = qpinfo.x[1] - return nothing -end - -## kernel for the (nonlinear) reaction-convection-diffusion operator -function nonlinear_kernel!(result, input, qpinfo) - u, ∇u, ϵ = input[1], input[2], qpinfo.params[1] - result[1] = exp(u) - exp(-u) - result[2] = ϵ * ∇u - return nothing -end - -## everything is wrapped in a main function -function main(; Plotter = nothing, h = 1.0e-2, ϵ = 1.0e-3, order = 2, kwargs...) - - ## problem description - PD = ProblemDescription("Nonlinear Poisson Equation") - u = Unknown("u"; name = "u") - assign_unknown!(PD, u) - assign_operator!(PD, NonlinearOperator(nonlinear_kernel!, [id(u), grad(u)]; params = [ϵ], kwargs...)) - assign_operator!(PD, LinearOperator(f!, [id(u)]; store = true, kwargs...)) - assign_operator!(PD, InterpolateBoundaryData(u, boundary_data!; kwargs...)) - - ## discretize: grid + FE space - xgrid = simplexgrid(0:h:1) - FES = FESpace{H1Pk{1, 1, order}}(xgrid) - - ## generate a solution vector and solve - sol = solve(PD, FES; kwargs...) - - ## plot discrete and exact solution (on finer grid) - plt = plot([id(u)], sol; Plotter = Plotter) - - return sol, plt -end - -generateplots = ExtendableFEM.default_generateplots(Example105_NonlinearPoissonEquation, "example105.png") #hide -function runtests() #hide - sol, plt = main(; h = 0.01, τ = 0.1, T = 1, use_diffeq = false) #hide - @test maximum(sol.entries) ≈ 0.4812118250102083 #hide - return nothing #hide -end #hide -end diff --git a/examples/Example106_NonlinearDiffusion.jl b/examples/Example106_NonlinearDiffusion.jl deleted file mode 100644 index 5250bd18..00000000 --- a/examples/Example106_NonlinearDiffusion.jl +++ /dev/null @@ -1,133 +0,0 @@ -#= - -# 106 : Nonlinear Diffusion -([source code](@__SOURCE_URL__)) - -This example solves the nonlinear diffusion equation -```math -\begin{aligned} -u_t - \Delta u^m & = 0 -\end{aligned} -``` -in ``\Omega := (-1,1)`` -with homogeneous Neumann boundary conditions. - -The solution looks like this: - -![](example106.png) - -=# - -module Example106_NonlinearDiffusion - -using ExtendableFEM -using ExtendableGrids -using OrdinaryDiffEqSDIRK -using GridVisualize -using Test #hide - -## Barenblatt solution -## (see Barenblatt, G. I. "On nonsteady motions of gas and fluid in porous medium." Appl. Math. and Mech.(PMM) 16.1 (1952): 67-78.) -function u_exact!(result, qpinfo) - t = qpinfo.time - x = qpinfo.x[1] - m = qpinfo.params[1] - tx = t^(-1.0 / (m + 1.0)) - xx = x * tx - xx = xx * xx - xx = 1 - xx * (m - 1) / (2.0 * m * (m + 1)) - if xx < 0.0 - xx = 0.0 - end - result[1] = tx * xx^(1.0 / (m - 1.0)) - return nothing -end - -function kernel_nonlinear!(result, input, qpinfo) - u, ∇u = input[1], input[2] - m = qpinfo.params[1] - result[1] = m * u^(m - 1) * ∇u - return nothing -end - -## everything is wrapped in a main function -function main(; - m = 2, - h = 0.05, - t0 = 0.001, - T = 0.01, - order = 1, - τ = 0.0001, - Plotter = nothing, - use_diffeq = true, - use_masslumping = true, - solver = ImplicitEuler(), - kwargs... - ) - - ## load mesh and exact solution - xgrid = simplexgrid(-1:h:1) - - ## set finite element types [surface height, velocity] - FEType = H1Pk{1, 1, order} - - ## generate empty PDEDescription for three unknowns (h, u) - PD = ProblemDescription("Nonlinear Diffusion Equation") - u = Unknown("u"; name = "u") - assign_unknown!(PD, u) - assign_operator!(PD, NonlinearOperator(kernel_nonlinear!, [grad(u)], [id(u), grad(u)]; params = [m], bonus_quadorder = 2)) - - ## prepare solution vector and initial data - FES = FESpace{FEType}(xgrid) - sol = FEVector(FES; tags = PD.unknowns) - interpolate!(sol[u], u_exact!; time = t0, params = [m]) - - ## init plotter and plot u0 - plt = GridVisualizer(; Plotter = Plotter, layout = (1, 2), size = (800, 400)) - scalarplot!(plt[1, 1], id(u), sol; label = "u_h", markershape = :circle, markevery = 1, title = "t = $t0") - - ## generate mass matrix (with mass lumping) - M = FEMatrix(FES) - assemble!(M, BilinearOperator([id(1)]; lump = 2 * use_masslumping)) - - if (use_diffeq) - ## generate ODE problem - prob = ExtendableFEM.generate_ODEProblem(PD, FES, (t0, T); init = sol, mass_matrix = M.entries.cscmatrix) - - ## solve ODE problem - de_sol = solve(prob, solver, abstol = 1.0e-6, reltol = 1.0e-3, dt = τ, dtmin = 1.0e-8, adaptive = true) - @info "#tsteps = $(length(de_sol.u))" - - ## get final solution - sol.entries .= de_sol.u[end] - else - ## add backward Euler time derivative - assign_operator!(PD, BilinearOperator(M, [u]; factor = 1 / τ, kwargs...)) - assign_operator!(PD, LinearOperator(M, [u], [u]; factor = 1 / τ, kwargs...)) - - ## generate solver configuration - SC = SolverConfiguration(PD, FES; init = sol, maxiterations = 1, kwargs...) - - ## iterate tspan - t = 0 - for it in 1:Int(floor((T - t0) / τ)) - t += τ - ExtendableFEM.solve(PD, FES, SC; time = t) - end - end - - ## plot final state and exact solution for comparison - scalarplot!(plt[1, 2], id(u), sol; label = "u_h", markershape = :circle, markevery = 1) - interpolate!(sol[1], u_exact!; time = T, params = [m]) - scalarplot!(plt[1, 2], id(u), sol; clear = false, color = :green, label = "u", title = "t = $T", legend = :best) - - return sol, plt -end - -generateplots = ExtendableFEM.default_generateplots(Example106_NonlinearDiffusion, "example106.png") #hide -function runtests(; T = 0.01, m = 2, kwargs...) #hide - sol, plt = main(; T = T, m = m, use_diffeq = false, kwargs...) #hide - @test maximum(sol.entries) ≈ 4.641588833612778 #hide - return nothing #hide -end #hide -end diff --git a/examples/Example108_RobinBoundaryCondition.jl b/examples/Example108_RobinBoundaryCondition.jl deleted file mode 100644 index 662d7a35..00000000 --- a/examples/Example108_RobinBoundaryCondition.jl +++ /dev/null @@ -1,96 +0,0 @@ -#= - -# 108 : Robin Boundary Condition -([source code](@__SOURCE_URL__)) - -This demonstrates the assignment of a mixed Robin boundary condition for a nonlinear 1D convection-diffusion-reaction PDE on the unit interval, i.e. -```math -\begin{aligned} --\partial^2 u / \partial x^2 + u \partial u / \partial x + u & = f && \text{in } \Omega\\ -u + \partial u / \partial_x & = g && \text{at } \Gamma_1 = \{ 0 \}\\ -u & = u_D && \text{at } \Gamma_2 = \{ 1 \} -\end{aligned} -``` -tested with data ``f(x) = e^{2x}``, ``g = 2`` and ``u_D = e`` such that ``u(x) = e^x`` is the exact solution. - -The solution looks like this: - -![](example108.png) -=# - -module Example108_RobinBoundaryCondition - -using ExtendableFEM -using ExtendableGrids -using GridVisualize -using Test #hide - -## data and exact solution -function f!(result, qpinfo) - result[1] = exp(2 * qpinfo.x[1]) - return nothing -end -function u!(result, qpinfo) - result[1] = exp(qpinfo.x[1]) - return nothing -end - -## kernel for the (nonlinear) reaction-convection-diffusion operator -function nonlinear_kernel!(result, input, qpinfo) - u, ∇u = input[1], input[2] - result[1] = u * ∇u + u # convection + reaction (will be multiplied with v) - result[2] = ∇u # diffusion (will be multiplied with ∇v) - return nothing -end - -## kernel for Robin boundary condition -function robin_kernel!(result, input, qpinfo) - result[1] = 2 - input[1] # = g - u (will be multiplied with v) - return nothing -end - -## everything is wrapped in a main function -function main(; Plotter = nothing, h = 1.0e-1, h_fine = 1.0e-3, order = 2, kwargs...) - - ## problem description - PD = ProblemDescription() - u = Unknown("u"; name = "u") - assign_unknown!(PD, u) - assign_operator!(PD, NonlinearOperator(nonlinear_kernel!, [id(u), grad(u)]; kwargs...)) - assign_operator!(PD, BilinearOperator(robin_kernel!, [id(u)]; entities = ON_BFACES, regions = [1], kwargs...)) - assign_operator!(PD, LinearOperator(f!, [id(u)]; kwargs...)) - assign_operator!(PD, InterpolateBoundaryData(u, u!; regions = [2], kwargs...)) - - ## generate coarse and fine mesh - xgrid = simplexgrid(0:h:1) - - ## choose finite element type and generate FESpace - FEType = H1Pk{1, 1, order} - FES = FESpace{FEType}(xgrid) - - ## generate a solution vector and solve - sol = solve(PD, FES; kwargs...) - - ## plot discrete and exact solution (on finer grid) - plt = GridVisualizer(Plotter = Plotter, layout = (1, 1)) - scalarplot!(plt[1, 1], id(u), sol; color = :black, label = "u_h", markershape = :circle, markersize = 10, markevery = 1) - xgrid_fine = simplexgrid(0:h_fine:1) - scalarplot!(plt[1, 1], xgrid_fine, view(nodevalues(xgrid_fine, u!), 1, :), clear = false, color = :red, label = "u", legend = :rb, markershape = :none) - - return sol, plt -end - -generateplots = ExtendableFEM.default_generateplots(Example108_RobinBoundaryCondition, "example108.png") #hide -function exact_error!(result, u, qpinfo) #hide - u!(result, qpinfo) #hide - result .= (result .- u) .^ 2 #hide - return nothing #hide -end #hide -function runtests(; kwargs...) #hide - sol, plt = main(; order = 2, kwargs...) #hide - L2error = ItemIntegrator(exact_error!, [id(1)]; quadorder = 4, kwargs...) #hide - error = sqrt(sum(evaluate(L2error, sol))) #hide - @test error ≈ 9.062544216508815e-6 #hide - return nothing #hide -end #hide -end diff --git a/examples/Example201_PoissonProblem.jl b/examples/Example201_PoissonProblem.jl deleted file mode 100644 index 1691cf70..00000000 --- a/examples/Example201_PoissonProblem.jl +++ /dev/null @@ -1,86 +0,0 @@ -#= - -# 201 : Poisson-Problem -([source code](@__SOURCE_URL__)) - -This example computes the solution ``u`` of the two-dimensional Poisson problem -```math -\begin{aligned} --\Delta u & = f \quad \text{in } \Omega -\end{aligned} -``` -with right-hand side ``f(x,y) \equiv xy`` and homogeneous Dirichlet boundary conditions -on the unit square domain ``\Omega`` on a given grid. - -The computed solution for the default parameters looks like this: - -![](example201.png) - -=# - -module Example201_PoissonProblem - -using ExtendableFEM -using ExtendableGrids -using GridVisualize -using Metis -using Test #hide - -## define variables -u = Unknown("u"; name = "potential") - -## default right-hand side f(x,y) = xy -function default_f!(fval, qpinfo) - fval[1] = qpinfo.x[1] * qpinfo.x[2] - return nothing -end - -function main(; - μ = 1.0, - nrefs = 4, - order = 2, - f! = default_f!, - use_restriction = true, - parallel = false, - npart = parallel ? 8 : 1, - Plotter = nothing, - kwargs... - ) - ## problem description - PD = ProblemDescription() - assign_unknown!(PD, u) - assign_operator!(PD, BilinearOperator([grad(u)]; parallel = parallel, factor = μ, kwargs...)) - assign_operator!(PD, LinearOperator(f!, [id(u)]; parallel = parallel, kwargs...)) - if use_restriction - assign_restriction!(PD, BoundaryDataRestriction(u; regions = 1:4, value = 0)) - else - assign_operator!(PD, HomogeneousBoundaryData(u; regions = 1:4)) - end - - ## discretize - xgrid = uniform_refine(grid_unitsquare(Triangle2D), nrefs) - if npart > 1 - xgrid = partition(xgrid, PlainMetisPartitioning(npart = npart)) - end - FES = FESpace{H1Pk{1, 2, order}}(xgrid) - - ## solve - sol = solve(PD, FES; kwargs...) - - ## plot - plt = plot([id(u), grad(u)], sol; Plotter = Plotter) - - return sol, plt -end - -generateplots = ExtendableFEM.default_generateplots(Example201_PoissonProblem, "example201.png") #hide -function runtests() #hide - sol, plt = main(; μ = 1.0, nrefs = 2, parallel = false, order = 2, npart = 2) #hide - @test sum(sol.entries) ≈ 1.1140313632246377 #hide - sol_parallel, plt = main(; μ = 1.0, nrefs = 2, order = 2, parallel = true, npart = 2) #hide - @assert sum((sol_parallel.entries .- sol.entries) .^ 2) ≈ 0.0 #hide - sol_restrict, plt = main(; μ = 1.0, nrefs = 2, use_restriction = false, parallel = false, order = 2, npart = 2) #hide - @assert sqrt(sum((sol_restrict.entries .- sol.entries) .^ 2)) < 1.0e-15 #hide - return nothing #hide -end #hide -end # module diff --git a/examples/Example202_MixedPoissonProblem.jl b/examples/Example202_MixedPoissonProblem.jl deleted file mode 100644 index 4948173e..00000000 --- a/examples/Example202_MixedPoissonProblem.jl +++ /dev/null @@ -1,86 +0,0 @@ -#= - -# 202 : Poisson-Problem (Mixed) -([source code](@__SOURCE_URL__)) - -This example computes the solution ``u`` and its stress ``\mathbf{\sigma} := - \mu \nabla u`` -of the two-dimensional Poisson problem in the mixed form -```math -\begin{aligned} -\mathbf{\sigma} + \mu \nabla u &= 0\\ -\mathrm{div} \mathbf{\sigma} & = f \quad \text{in } \Omega -\end{aligned} -``` -with right-hand side ``f(x,y) \equiv xy`` and homogeneous Dirichlet boundary conditions -on the unit square domain ``\Omega`` on a given grid. - -The computed solution looks like this: - -![](example202.png) -=# - -module Example202_MixedPoissonProblem - -using ExtendableFEM -using ExtendableGrids -using Test #hide - -## define unknowns -σ = Unknown("σ"; name = "pseudostress") -u = Unknown("u"; name = "potential") - -## bilinearform kernel for mixed Poisson problem -function blf!(result, u_ops, qpinfo) - σ, divσ, u = view(u_ops, 1:2), view(u_ops, 3), view(u_ops, 4) - μ = qpinfo.params[1] - result[1] = σ[1] / μ - result[2] = σ[2] / μ - result[3] = -u[1] - result[4] = divσ[1] - return nothing -end -## right-hand side data -function f!(fval, qpinfo) - fval[1] = qpinfo.x[1] * qpinfo.x[2] - return nothing -end -## boundary data -function boundarydata!(result, qpinfo) - result[1] = 0 - return nothing -end - -function main(; nrefs = 5, μ = 0.25, order = 0, Plotter = nothing, kwargs...) - - ## problem description - PD = ProblemDescription() - assign_unknown!(PD, u) - assign_unknown!(PD, σ) - assign_operator!(PD, BilinearOperator(blf!, [id(σ), div(σ), id(u)]; params = [μ], kwargs...)) - assign_operator!(PD, LinearOperator(boundarydata!, [normalflux(σ)]; entities = ON_BFACES, regions = 1:4, kwargs...)) - assign_operator!(PD, LinearOperator(f!, [id(u)]; kwargs...)) - assign_operator!(PD, FixDofs(u; dofs = [1], vals = [0])) - - ## discretize - xgrid = uniform_refine(grid_unitsquare(Triangle2D), nrefs) - FES = Dict( - u => FESpace{order == 0 ? L2P0{1} : H1Pk{1, 2, order}}(xgrid; broken = true), - σ => FESpace{HDIVRTk{2, order}}(xgrid) - ) - - ## solve - sol = ExtendableFEM.solve(PD, FES; kwargs...) - - ## plot - plt = plot([id(u), id(σ)], sol; Plotter = Plotter) - - return sol, plt -end - -generateplots = ExtendableFEM.default_generateplots(Example202_MixedPoissonProblem, "example202.png") #hide -function runtests() #hide - sol, plt = main(; μ = 0.25, order = 0, nrefs = 2) #hide - @test maximum(view(sol[1])) ≈ 0.08463539106946043 #hide - return nothing #hide -end #hide -end # module diff --git a/examples/Example203_PoissonProblemDG.jl b/examples/Example203_PoissonProblemDG.jl deleted file mode 100644 index 786fca48..00000000 --- a/examples/Example203_PoissonProblemDG.jl +++ /dev/null @@ -1,142 +0,0 @@ -#= - -# 203 : Poisson-Problem DG -([source code](@__SOURCE_URL__)) - -This example computes the solution ``u`` of the two-dimensional Poisson problem -```math -\begin{aligned} --\Delta u & = f \quad \text{in } \Omega -\end{aligned} -``` -with right-hand side ``f`` and inhomogeneous Dirichlet boundary conditions -chosen such that ``u(x,y) = x^3 - 3xy^2``. -This time the problem is solved on a given grid via the discontinuous Galerkin method. - -The computed solution looks like this: - -![](example203.png) -=# - -module Example203_PoissonProblemDG - -using ExtendableFEM -using ExtendableGrids -using LinearAlgebra -using Metis -using Symbolics -using Test #hide - -## exact data for problem by Symbolics -function prepare_data(; μ = 1) - - @variables x y - - ## exact solution - u = x^3 - 3 * x * y^2 - ∇u = Symbolics.gradient(u, [x, y]) - - ## right-hand side - Δu = Symbolics.gradient(∇u[1], [x]) + Symbolics.gradient(∇u[2], [y]) - f = -μ * Δu[1] - - ## build functions - u_eval = build_function(u, x, y, expression = Val{false}) - ∇u_eval = build_function(∇u, x, y, expression = Val{false}) - f_eval = build_function(f, x, y, expression = Val{false}) - return f_eval, u_eval, ∇u_eval[2] -end - -function main(; dg = true, μ = 1.0, τ = 10.0, nrefs = 4, order = 2, bonus_quadorder = 2, parallel = false, npart = parallel ? 8 : 1, Plotter = nothing, kwargs...) - - ## prepare problem data - f_eval, u_eval, ∇u_eval = prepare_data(; μ = μ) - rhs!(result, qpinfo) = (result[1] = f_eval(qpinfo.x[1], qpinfo.x[2])) - exact_u!(result, qpinfo) = (result[1] = u_eval(qpinfo.x[1], qpinfo.x[2])) - exact_∇u!(result, qpinfo) = (∇u_eval(result, qpinfo.x[1], qpinfo.x[2])) - - ## problem description - PD = ProblemDescription("Poisson problem") - u = Unknown("u"; name = "potential") - assign_unknown!(PD, u) - assign_operator!(PD, BilinearOperator([grad(u)]; factor = μ, parallel = parallel, kwargs...)) - assign_operator!(PD, LinearOperator(rhs!, [id(u)]; bonus_quadorder = bonus_quadorder, parallel = parallel, kwargs...)) - assign_operator!(PD, InterpolateBoundaryData(u, exact_u!; bonus_quadorder = bonus_quadorder, regions = 1:4)) - - ## discretize - xgrid = uniform_refine(grid_unitsquare(Triangle2D), nrefs) - if npart > 1 - xgrid = partition(xgrid, PlainMetisPartitioning(npart = npart); edges = true) - end - FES = FESpace{order == 0 ? L2P0{1} : H1Pk{1, 2, order}}(xgrid; broken = dg) - - ## add DG terms - assign_operator!(PD, BilinearOperatorDG(dg_kernel, [jump(id(u))], [average(grad(u))]; entities = ON_FACES, factor = -μ, transposed_copy = 1, parallel = parallel, kwargs...)) - assign_operator!(PD, LinearOperatorDG(dg_kernel_bnd(exact_u!), [average(grad(u))]; entities = ON_BFACES, factor = -μ, bonus_quadorder = bonus_quadorder, parallel = parallel, kwargs...)) - assign_operator!(PD, BilinearOperatorDG(dg_kernel2, [jump(id(u))]; entities = ON_FACES, factor = μ * τ, parallel = parallel, kwargs...)) - assign_operator!(PD, LinearOperatorDG(dg_kernel2_bnd(exact_u!), [id(u)]; entities = ON_BFACES, regions = 1:4, factor = μ * τ, bonus_quadorder = bonus_quadorder, parallel = parallel, kwargs...)) - - - ## solve - sol = solve(PD, FES; kwargs...) - - ## prepare error calculation - function exact_error!(result, u, qpinfo) - exact_u!(result, qpinfo) - exact_∇u!(view(result, 2:3), qpinfo) - result .-= u - result .= result .^ 2 - return nothing - end - function dgjumps!(result, u, qpinfo) - result .= u[1]^2 / qpinfo.volume - return nothing - end - ErrorIntegratorExact = ItemIntegrator(exact_error!, [id(u), grad(u)]; quadorder = 2 * (order + 1), params = [μ], kwargs...) - DGJumpsIntegrator = ItemIntegratorDG(dgjumps!, [jump(id(u))]; entities = ON_IFACES, kwargs...) - - ## calculate error - error = evaluate(ErrorIntegratorExact, sol) - dgjumps = sqrt(sum(evaluate(DGJumpsIntegrator, sol))) - L2error = sqrt(sum(view(error, 1, :))) - H1error = sqrt(sum(view(error, 2, :)) + sum(view(error, 3, :))) - @info "L2 error = $L2error" - @info "H1 error = $H1error" - @info "dgjumps = $dgjumps" - - ## plot - plt = plot([id(u), grad(u)], sol; Plotter = Plotter) - - return L2error, plt -end - -function dg_kernel(result, input, qpinfo) - result[1] = dot(input, qpinfo.normal) - return nothing -end -function dg_kernel_bnd(uDb! = nothing) - return function closure(result, qpinfo) - uDb!(result, qpinfo) - result[1:2] = result[1] .* qpinfo.normal - return nothing - end -end -function dg_kernel2(result, input, qpinfo) - result .= input / qpinfo.volume - return nothing -end -function dg_kernel2_bnd(uDb! = nothing) - return function closure(result, qpinfo) - uDb!(result, qpinfo) - result /= qpinfo.volume - return nothing - end -end - -generateplots = ExtendableFEM.default_generateplots(Example203_PoissonProblemDG, "example203.png") #hide -function runtests(; kwargs...) #hide - L2error, ~ = main(; μ = 0.25, nrefs = 2, order = 2, kwargs...) #hide - @test L2error ≈ 0.00020400470505497443 #hide - return nothing #hide -end #hide -end # module diff --git a/examples/Example204_LaplaceEVProblem.jl b/examples/Example204_LaplaceEVProblem.jl deleted file mode 100644 index 8a0b5b27..00000000 --- a/examples/Example204_LaplaceEVProblem.jl +++ /dev/null @@ -1,69 +0,0 @@ -#= - -# 204 : Eigenvalue problem for the Laplacian -([source code](@__SOURCE_URL__)) - -This example computes the pairs of eigenvalues and eigenvectors -``(\lambda,u) \in \mathbb{R} \times H^1_0(\Omega)`` -of the Laplacian, i.e, -```math -\begin{aligned} --\Delta u & = \lambda u \quad \text{in } \Omega -\end{aligned} -``` -on a two-dimensional L-shaped domain with homogeneous boundary conditions -with the help of an iterative solver from [KrylovKit.jl](https://github.com/Jutho/KrylovKit.jl). -The first twelve computed eigenvectors look like this: - -![](example204.png) -=# - -module Example204_LaplaceEVProblem - -using ExtendableFEM -using ExtendableGrids -using ExtendableSparse -using LinearAlgebra -using GridVisualize -using KrylovKit - -function main(; which = 1:12, ncols = 3, nrefs = 4, order = 1, Plotter = nothing, kwargs...) - - ## discretize - xgrid = uniform_refine(grid_lshape(Triangle2D), nrefs) - FES = FESpace{H1Pk{1, 2, order}}(xgrid) - - ## assemble operators - A = FEMatrix(FES) - B = FEMatrix(FES) - u = FEVector(FES; name = "u") - assemble!(A, BilinearOperator([grad(1)]; kwargs...)) - assemble!(A, BilinearOperator([id(1)]; entities = ON_BFACES, factor = 1.0e4, kwargs...)) - assemble!(B, BilinearOperator([id(1)]; kwargs...)) - - ## solver generalized eigenvalue problem iteratively with KrylovKit - λs, x, info = geneigsolve((A.entries, B.entries), maximum(which), :SR; maxiter = 4000, issymmetric = true, tol = 1.0e-8) - @show info - @assert info.converged >= maximum(which) - - ## plot requested eigenvalue pairs - nEVs = length(which) - nrows = Int(ceil(nEVs / ncols)) - plt = GridVisualizer(; Plotter = Plotter, layout = (nrows, ncols), clear = true, resolution = (900, 900 / ncols * nrows)) - col, row = 0, 1 - for j in which - col += 1 - if col == ncols + 1 - col, row = 1, row + 1 - end - λ = λs[j] - @info "λ[$j] = $λ, residual = $(sum(info.residual[j]))" - u.entries .= Real.(x[j]) - scalarplot!(plt[row, col], id(1), u; Plotter = Plotter, title = "λ[$j] = $(Float16(λ))") - end - - return u, plt -end - -generateplots = ExtendableFEM.default_generateplots(Example204_LaplaceEVProblem, "example204.png") #hide -end # module diff --git a/examples/Example205_HeatEquation.jl b/examples/Example205_HeatEquation.jl deleted file mode 100644 index 0ef9ad3d..00000000 --- a/examples/Example205_HeatEquation.jl +++ /dev/null @@ -1,99 +0,0 @@ -#= - -# 205 : Heat equation -([source code](@__SOURCE_URL__)) - -This example computes the solution ``u`` of the two-dimensional heat equation -```math -\begin{aligned} -u_t - \Delta u & = 0 \quad \text{in } \Omega -\end{aligned} -``` -for homogeneous Dirichlet boundary conditions and some given initial state -on the unit square domain ``\Omega`` on a given grid. - -The initial condition and the final solution for the default parameters looks like this: - -![](example205.png) - -=# - -module Example205_HeatEquation - -using ExtendableFEM -using ExtendableGrids -using OrdinaryDiffEqSDIRK -using Test #hide - -## initial state u at time t0 -function initial_data!(result, qpinfo) - x = qpinfo.x - result[1] = exp(-5 * x[1]^2 - 5 * x[2]^2) - return nothing -end - -function main(; - nrefs = 4, T = 2.0, τ = 1.0e-3, order = 2, use_diffeq = true, - solver = ImplicitEuler(), Plotter = nothing, kwargs... - ) - - ## problem description - PD = ProblemDescription("Heat Equation") - u = Unknown("u") - assign_unknown!(PD, u) - assign_operator!(PD, BilinearOperator([grad(u)]; store = true, kwargs...)) - assign_operator!(PD, HomogeneousBoundaryData(u; regions = 1:4)) - - ## grid - xgrid = uniform_refine(grid_unitsquare(Triangle2D; scale = [4, 4], shift = [-0.5, -0.5]), nrefs) - - ## prepare solution vector and initial data u0 - FES = FESpace{H1Pk{1, 2, order}}(xgrid) - sol = FEVector(FES; tags = PD.unknowns) - interpolate!(sol[u], initial_data!; bonus_quadorder = 5) - - ## init plotter and plot u0 - plt = plot([id(u)], sol; add = 1, Plotter = Plotter, title_add = " (t = 0)") - - if (use_diffeq) - ## generate DifferentialEquations.ODEProblem - prob = generate_ODEProblem(PD, FES, (0.0, T); init = sol, constant_matrix = true) - - ## solve ODE problem - de_sol = solve(prob, solver, abstol = 1.0e-6, reltol = 1.0e-3, dt = τ, dtmin = 1.0e-6, adaptive = true) - @info "#tsteps = $(length(de_sol.u))" - - ## get final solution - sol.entries .= de_sol.u[end] - else - ## add backward Euler time derivative - M = FEMatrix(FES) - assemble!(M, BilinearOperator([id(1)])) - assign_operator!(PD, BilinearOperator(M, [u]; factor = 1 / τ, kwargs...)) - assign_operator!(PD, LinearOperator(M, [u], [u]; factor = 1 / τ, kwargs...)) - - ## generate solver configuration - SC = SolverConfiguration(PD, FES; init = sol, maxiterations = 1, constant_matrix = true, kwargs...) - - ## iterate tspan - t = 0 - for it in 1:Int(floor(T / τ)) - t += τ - ExtendableFEM.solve(PD, FES, SC; time = t, verbosity = -1, timeroutputs = :hide) - end - @show timeroutputs(SC) - end - - ## plot final state - plot!(plt, [id(u)], sol; keep = 1, title_add = " (t = $T)") - - return sol, plt -end - -generateplots = ExtendableFEM.default_generateplots(Example205_HeatEquation, "example205.png") #hide -function runtests(; verbosity = -1, kwargs...) #hide - sol, plt = main(; nrefs = 2, T = 1, use_diffeq = false, kwargs...) #hide - @test maximum(sol.entries) ≈ 0.041490419236077006 #hide - return nothing #hide -end #hide -end # module diff --git a/examples/Example206_CoupledSubGridProblems.jl b/examples/Example206_CoupledSubGridProblems.jl deleted file mode 100644 index 0e2637eb..00000000 --- a/examples/Example206_CoupledSubGridProblems.jl +++ /dev/null @@ -1,129 +0,0 @@ -#= - -# 206 : CoupledSubGridProblems -([source code](@__SOURCE_URL__)) - -This example demonstrates how to solve a coupled problem where two -variables only live on a sub-domain and are coupled through an interface condition. -Consider the unit square domain cut in half through on of its diagonals. -On each subdomain a solutiong ``u_j`` of the two-dimensional Poisson problem -```math -\begin{aligned} --\Delta u & = f \quad \text{in } \Omega -\end{aligned} -``` -with inhomogeneous boundary conditions on the former boundaries of the full square is searched. -Along the common boundary between the two subdomains a new interface region is assigned (appended to BFaceNodes) -and an interface condition is assembled that couples the two solutions ``u_1`` and ``u_2`` -to each other. -In this toy example, this interface conditions penalizes the jump between the two solutions on each side -of the diagonal. Oberserve, that if the penalization factor ``\tau`` is large, the two solutions are -almost equal along the interface. - -The computed solution(s) looks like this: - -![](example206.png) - -Each column of the plot shows the solution, the subgrid it lives on. The last row shows the full grid. - -=# - -module Example206_CoupledSubGridProblems - -using ExtendableFEM -using ExtendableGrids -using Test # - - -function boundary_conditions!(result, qpinfo) - result[1] = 1 - qpinfo.x[1] - qpinfo.x[2] # used for both subsolutions - return nothing -end - -function interface_condition!(result, u, qpinfo) - result[1] = u[1] - u[2] - result[2] = -result[1] - return nothing -end - -function interface_condition_LM!(result, u, qpinfo) - result[1] = (u[1] - u[2]) - return nothing -end - - -function main(; μ = [1.0, 1.0], f = [10, -10], τ = 1, use_LM = true, nref = 4, order = 2, Plotter = nothing, kwargs...) - - ## Finite element type - FEType = H1Pk{1, 2, order} - FETypeLM = H1Pk{1, 1, order} - - ## generate mesh - xgrid = grid_unitsquare(Triangle2D) - - ## define regions - xgrid[CellRegions] = Int32[1, 2, 2, 1] - - ## add an interface between region 1 and 2 - ## (one can use the BFace storages for that) - xgrid[BFaceNodes] = Int32[xgrid[BFaceNodes] [2 5; 5 4]] - append!(xgrid[BFaceRegions], [5, 5]) - xgrid[FaceRegions][xgrid[BFaceFaces][(end - 1):end]] .= 5 - xgrid[BFaceGeometries] = VectorOfConstants{ElementGeometries, Int}(Edge1D, 6) - - ## refine - xgrid = uniform_refine(xgrid, nref) - - ## define an FESpace just on region 1 and one just on region 2 - FES1 = FESpace{FEType}(xgrid; regions = [1]) - FES2 = FESpace{FEType}(xgrid; regions = [2]) - if use_LM - FES3 = FESpace{FETypeLM, ON_FACES}(xgrid; regions = [5]) - @show FES3.xgrid FES3.dofgrid - end - - ## define variables - u1 = Unknown("u1"; name = "potential in region 1") - u2 = Unknown("u2"; name = "potential in region 2") - p = Unknown("p"; name = "LM for interface condition") - - ## problem description - PD = ProblemDescription() - assign_unknown!(PD, u1) - assign_unknown!(PD, u2) - assign_operator!(PD, BilinearOperator([grad(u1)]; regions = [1], factor = μ[1], kwargs...)) - assign_operator!(PD, BilinearOperator([grad(u2)]; regions = [2], factor = μ[2], kwargs...)) - assign_operator!(PD, LinearOperator([id(u1)]; regions = [1], factor = f[1])) - assign_operator!(PD, LinearOperator([id(u2)]; regions = [2], factor = f[2])) - if use_LM - assign_unknown!(PD, p) - assign_operator!(PD, BilinearOperator(interface_condition_LM!, [id(p)], [id(u1), id(u2)]; regions = [5], transposed_copy = 1, entities = ON_FACES, kwargs...)) - else - assign_operator!(PD, BilinearOperator(interface_condition!, [id(u1), id(u2)]; regions = [5], factor = τ, entities = ON_FACES, kwargs...)) - end - assign_operator!(PD, InterpolateBoundaryData(u1, boundary_conditions!; regions = 1:4)) - assign_operator!(PD, InterpolateBoundaryData(u2, boundary_conditions!; regions = 1:4)) - - sol = solve(PD, use_LM ? [FES1, FES2, FES3] : [FES1, FES2]) - - plt = plot([id(u1), id(u2), dofgrid(u1), dofgrid(u2), grid(u1)], sol; Plotter = Plotter) - - return sol, plt -end - -generateplots = ExtendableFEM.default_generateplots(Example206_CoupledSubGridProblems, "example206.png") #hide - -function jump_l2norm!(result, u, qpinfo) #hide - result[1] = (u[1] - u[2])^2 #hide - return nothing #hide -end #hide -function runtests() #hide - ## test if jump at interface vanishes for large penalty #hide - sol, plt = main(; τ = 1.0e9, nrefs = 2, order = 2) #hide - jump_integrator = ItemIntegrator(jump_l2norm!, [id(1), id(2)]; entities = ON_BFACES, regions = [5], resultdim = 1, quadorder = 4) #hide - jump_error = sqrt(sum(evaluate(jump_integrator, sol))) #hide - @info "||[u_1 - u_2]|| = $(jump_error)" #hide - @test jump_error < 1.0e-8 #hide - return nothing #hide -end #hide -end #module diff --git a/examples/Example207_AdvectionUpwindDG.jl b/examples/Example207_AdvectionUpwindDG.jl deleted file mode 100644 index e6bc499b..00000000 --- a/examples/Example207_AdvectionUpwindDG.jl +++ /dev/null @@ -1,186 +0,0 @@ -#= - -# 207 : Advection Upwind-DG -([source code](@__SOURCE_URL__)) - -This example computes the solution ``u`` of the two-dimensional advection equation -```math -\begin{aligned} -\mathrm{div} (\beta u) & = 0 \quad \text{in } \Omega -\end{aligned} -``` -with some given (divergence-free) advection field ``\beta`` and inhomogeneous -Dirichlet boundary conditions at the inflow boundary -(where ``\beta \cdot n < 0`` with ``n`` being the outer normal vector). - -In the example below the field ``\beta(x,y) = (-y, x)`` and the inflow data -```math - u(x,y) = -\begin{cases} - 1 & \text{for } x \in [0,r] \& y = 0,\\ - 0 & \text{for } x \in (r,1] \& y = 0,\\ - 0 & \text{for } x = 1 \& y \in [0,1]. -\end{cases} -``` -is employed. The expected solution is a piecewise constant function -that assumes the value one in the circle of radius ``r`` and zero elsewhere. - -Moreover, the upwind discontinuous Galerkin method for arbitrary polynomial degree -is used for the discretization of the problem, but the continuous Galerkin -method can be switched on with dg = false for comparison. For piecewise constants -the DG method satisfies the maximum principle. - -The grid (which is heavily refined along the interface of the circle) and the -computed solution looks like this: - -![](example207.png) -=# - -module Example207_AdvectionUpwindDG - -using ExtendableFEM -using ExtendableGrids -using Symbolics -using LinearAlgebra -using SimplexGridFactory -using Triangulate -using Test #hide - -## wind = advection field β -function β!(result, qpinfo) - x = qpinfo.x - result[1] = - x[2] - result[2] = x[1] - return nothing -end - -## exact solution -function exact_u!(result, qpinfo) - x = qpinfo.x - r = qpinfo.params[1] - result[1] = sqrt(x[1]^2 + x[2]^2) <= r ? 1 : 0 - return nothing -end - -## integrand of the advection bilinearform -function advection_kernel!(result, input, qpinfo) - β!(result, qpinfo) # evaluate wind β - result .*= input[1] # multiply with u_h - return nothing -end - -function outflow_kernel!(xgrid) - beta = zeros(Float64, 2) - return function closure(result, input, qpinfo) - face = qpinfo.item - β!(beta, qpinfo) - result[1] = dot(beta, qpinfo.normal) * input[1] - return nothing - end -end - -function upwind_kernel!(xgrid) - beta = zeros(Float64, 2) - return function closure(result, input, qpinfo) - face = qpinfo.item - β!(beta, qpinfo) - result[1] = dot(beta, qpinfo.normal) - return if result[1] > 0 ## wind blows this -> other - result[1] *= input[1] # upwind value = this - else ## wind blows this <- other - result[1] *= input[2] # upwind value = other - end - end -end - -## prepare error calculation -function exact_error!(result, u, qpinfo) - exact_u!(result, qpinfo) - result[1] = (result[1] - u[1])^2 - return nothing -end - -function main(; nref = 4, order = 0, r = 0.5, dg = true, Plotter = nothing, kwargs...) - - ## grid - xgrid = make_grid(nref, r) - - ## problem description - PD = ProblemDescription("advection equation") - u = Unknown("u"; name = "species") - assign_unknown!(PD, u) - - ## advection operator - assign_operator!(PD, BilinearOperator(advection_kernel!, [grad(u)], [id(u)]; factor = -1, bonus_quadorder = 1, kwargs...)) - if dg - assign_operator!(PD, BilinearOperatorDG(upwind_kernel!(xgrid), [jump(id(u))], [this(id(u)), other(id(u))]; entities = ON_IFACES, bonus_quadorder = 1, kwargs...)) - end - - ## outflow boundary (regions [3,4]) and inflow boundary (regions [5,6]) - assign_operator!(PD, BilinearOperator(outflow_kernel!(xgrid), [id(u)]; entities = ON_BFACES, regions = [3, 4])) - assign_operator!(PD, InterpolateBoundaryData(u, exact_u!; regions = [5, 6], params = [r], kwargs...)) - assign_operator!(PD, HomogeneousBoundaryData(u; regions = [1, 2], kwargs...)) - - ## solve - FES = FESpace{order == 0 ? L2P0{1} : H1Pk{1, 2, order}}(xgrid; broken = dg) - sol = solve(PD, FES; kwargs...) - - ## calculate L2 error and min/max value - ErrorIntegrator = ItemIntegrator(exact_error!, [id(u)]; quadorder = 2 * order, params = [r], kwargs...) - L2error = sqrt(sum(view(evaluate(ErrorIntegrator, sol), 1, :))) - @info "L2 error = $L2error" - @info "extrema = $(extrema(sol.entries))" - - ## plot - plt = plot([grid(u), id(u)], sol; Plotter = Plotter) - - return sol, plt -end - -## grid generator script using SimplexGridBuilder/Triangulate -function make_grid(nref = 4, radius = 0.5) - builder = SimplexGridBuilder(Generator = Triangulate) - - ## define outer boundary nodes and regions - p1 = point!(builder, 0, 0) - p12 = point!(builder, radius, 0) - p2 = point!(builder, 1, 0) - p3 = point!(builder, 1, 1) - p4 = point!(builder, 0, 1) - p41 = point!(builder, 0, radius) - - facetregion!(builder, 5) - facet!(builder, p1, p12) - facetregion!(builder, 1) - facet!(builder, p12, p2) - facetregion!(builder, 2) - facet!(builder, p2, p3) - facetregion!(builder, 3) - facet!(builder, p3, p4) - facetregion!(builder, 4) - facet!(builder, p4, p41) - facetregion!(builder, 6) - facet!(builder, p41, p1) - - ## add interior interface (quarter circle) - n = 4^(nref + 1) - points = [point!(builder, radius * sin(t), radius * cos(t)) for t in range(0, π / 2, length = n)] - facetregion!(builder, 7) - for i in 2:(n - 2) - facet!(builder, points[i], points[i + 1]) - end - facet!(builder, p41, points[1]) - facet!(builder, points[end], p12) - - ## generate - return simplexgrid(builder, maxvolume = 1) -end - -generateplots = ExtendableFEM.default_generateplots(Example207_AdvectionUpwindDG, "example207.png") #hide -function runtests() #hide - ## test if P0-DG solution stays within bounds #hide - sol, ~ = main(; order = 0, nrefs = 2) #hide - @test norm(extrema(sol.entries) .- (0, 1)) < 1.0e-12 #hide - return nothing #hide -end #hide -end # module diff --git a/examples/Example210_LshapeAdaptivePoissonProblem.jl b/examples/Example210_LshapeAdaptivePoissonProblem.jl deleted file mode 100644 index 2fc55bef..00000000 --- a/examples/Example210_LshapeAdaptivePoissonProblem.jl +++ /dev/null @@ -1,191 +0,0 @@ -#= - -# 210 : Poisson L-shape Adaptive Mesh Refinement -([source code](@__SOURCE_URL__)) - -This example computes the standard-residual error estimator for the $H^1$ error ``e = u - u_h`` of some $H^1$-conforming -approximation ``u_h`` to the solution ``u`` of some Poisson problem ``-\Delta u = f`` on an L-shaped domain, i.e. -```math -\eta^2(u_h) := \sum_{T \in \mathcal{T}} \lvert T \rvert \| f + \Delta u_h \|^2_{L^2(T)} -+ \sum_{F \in \mathcal{F}} \lvert F \rvert \| [[\nabla u_h \cdot \mathbf{n}]] \|^2_{L^2(F)} -``` -This example script showcases the evaluation of 2nd order derivatives like the Laplacian and adaptive mesh refinement. - -The resulting mesh and error convergence history for the default parameters looks like: - -![](example210.png) - -=# - -module Example210_LshapeAdaptivePoissonProblem - -using ExtendableFEM -using GridVisualize -using ExtendableGrids -using LinearAlgebra -using Test #hide - -## exact solution u for the Poisson problem -function u!(result, qpinfo) - x = qpinfo.x - r2 = x[1]^2 + x[2]^2 - φ = atan(x[2], x[1]) - if φ < 0 - φ += 2 * pi - end - result[1] = r2^(1 / 3) * sin(2 * φ / 3) - return nothing -end - -## gradient of exact solution -function ∇u!(result, qpinfo) - x = qpinfo.x - φ = atan(x[2], x[1]) - r2 = x[1]^2 + x[2]^2 - if φ < 0 - φ += 2 * pi - end - ∂r = 2 / 3 * r2^(-1 / 6) * sin(2 * φ / 3) - ∂φ = 2 / 3 * r2^(-1 / 6) * cos(2 * φ / 3) - result[1] = cos(φ) * ∂r - sin(φ) * ∂φ - result[2] = sin(φ) * ∂r + cos(φ) * ∂φ - return nothing -end - -## kernel for exact error calculation -function exact_error!(result, u, qpinfo) - u!(result, qpinfo) - ∇u!(view(result, 2:3), qpinfo) - result .-= u - result .= result .^ 2 - return nothing -end - -## kernel for face interpolation of normal jumps of gradient -function gradnormalflux!(result, ∇u, qpinfo) - result[1] = dot(∇u, qpinfo.normal) - return nothing -end - -## kernel for face refinement indicator -function η_face!(result, gradjump, qpinfo) - result .= qpinfo.volume * gradjump .^ 2 - return nothing -end - -## kernel for cell refinement indicator -function η_cell!(result, Δu, qpinfo) - result .= qpinfo.volume * Δu .^ 2 - return nothing -end - -function main(; maxdofs = 4000, θ = 0.5, μ = 1.0, nrefs = 1, order = 2, Plotter = nothing, kwargs...) - - ## problem description - PD = ProblemDescription("Poisson problem") - u = Unknown("u"; name = "u") - assign_unknown!(PD, u) - assign_operator!(PD, BilinearOperator([grad(u)]; factor = μ, kwargs...)) - assign_operator!(PD, InterpolateBoundaryData(u, u!; regions = 2:7, bonus_quadorder = 4, kwargs...)) - assign_operator!(PD, HomogeneousBoundaryData(u; regions = [1, 8])) - - ## discretize - xgrid = uniform_refine(grid_lshape(Triangle2D), nrefs) - - ## define interpolators and item integrators for error estimation and calculation - NormalJumpProjector = FaceInterpolator(gradnormalflux!, [jump(grad(u))]; resultdim = 1, order = order, only_interior = true, kwargs...) - ErrorIntegratorFace = ItemIntegrator(η_face!, [id(1)]; quadorder = 2 * order, entities = ON_FACES, kwargs...) - ErrorIntegratorCell = ItemIntegrator(η_cell!, [Δ(1)]; quadorder = 2 * (order - 2), entities = ON_CELLS, kwargs...) - ErrorIntegratorExact = ItemIntegrator(exact_error!, [id(1), grad(1)]; quadorder = 2 * order, kwargs...) - - NDofs = zeros(Int, 0) - ResultsL2 = zeros(Float64, 0) - ResultsH1 = zeros(Float64, 0) - Resultsη = zeros(Float64, 0) - sol = nothing - ndofs = 0 - level = 0 - while ndofs < maxdofs - level += 1 - - ## SOLVE : create a solution vector and solve the problem - println("------- LEVEL $level") - if ndofs < 1000 - println(stdout, unicode_gridplot(xgrid)) - end - @time begin - ## solve - FES = FESpace{H1Pk{1, 2, order}}(xgrid) - sol = ExtendableFEM.solve(PD, FES; u = [u], kwargs...) - ndofs = length(sol[1]) - push!(NDofs, ndofs) - println("\t ndof = $ndofs") - print("@time solver =") - end - - ## ESTIMATE : calculate local error estimator contributions - @time begin - ## calculate error estimator - Jumps4Faces = evaluate!(NormalJumpProjector, sol) - η_F = evaluate(ErrorIntegratorFace, Jumps4Faces) - - η_T = evaluate(ErrorIntegratorCell, sol) - facecells = xgrid[FaceCells] - for face in 1:size(facecells, 2) - η_F[face] += η_T[facecells[1, face]] - if facecells[2, face] > 0 - η_F[face] += η_T[facecells[2, face]] - end - end - - ## calculate total estimator - push!(Resultsη, sqrt(sum(η_F))) - print("@time η eval =") - end - - ## calculate exact L2 error, H1 error - @time begin - error = evaluate(ErrorIntegratorExact, sol) - push!(ResultsL2, sqrt(sum(view(error, 1, :)))) - push!(ResultsH1, sqrt(sum(view(error, 2, :)) + sum(view(error, 3, :)))) - print("@time e eval =") - end - - if ndofs >= maxdofs - break - end - - ## MARK+REFINE : mesh refinement - @time begin - if θ >= 1 ## uniform mesh refinement - xgrid = uniform_refine(xgrid) - else ## adaptive mesh refinement - ## refine by red-green-blue refinement (incl. closuring) - facemarker = bulk_mark(xgrid, view(η_F, :), θ; indicator_AT = ON_FACES) - xgrid = RGB_refine(xgrid, facemarker) - end - print("@time refine =") - end - println("\t η = $(Resultsη[level])\n\t e = $(ResultsH1[level])") - end - - ## plot - plt = GridVisualizer(; Plotter = Plotter, layout = (2, 2), clear = true, size = (1000, 1000)) - scalarplot!(plt[1, 1], id(u), sol; levels = 7, title = "u_h") - plot_convergencehistory!(plt[1, 2], NDofs, [ResultsL2 ResultsH1 Resultsη]; add_h_powers = [order, order + 1], X_to_h = X -> order * X .^ (-1 / 2), ylabels = ["|| u - u_h ||", "|| ∇(u - u_h) ||", "η"]) - gridplot!(plt[2, 1], xgrid; linewidth = 1) - gridplot!(plt[2, 2], xgrid; linewidth = 1, xlimits = [-0.0005, 0.0005], ylimits = [-0.0005, 0.0005]) - - ## print 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 -end - -generateplots = ExtendableFEM.default_generateplots(Example210_LshapeAdaptivePoissonProblem, "example210.png") #hide -function runtests() #hide - sol, plt = main(; maxdofs = 1000, order = 2) #hide - @test length(sol.entries) == 1007 #hide - return nothing #hide -end #hide -end # module diff --git a/examples/Example211_LshapeAdaptiveEQPoissonProblem.jl b/examples/Example211_LshapeAdaptiveEQPoissonProblem.jl deleted file mode 100644 index 2798d9b7..00000000 --- a/examples/Example211_LshapeAdaptiveEQPoissonProblem.jl +++ /dev/null @@ -1,372 +0,0 @@ -#= - -# 211 : Poisson L-shape Local Equilibrated Fluxes -([source code](@__SOURCE_URL__)) - -This example computes a local equilibration error estimator for the $H^1$ error of some $H^1$-conforming -approximation ``u_h`` to the solution ``u`` of some Poisson problem ``-\Delta u = f`` on an L-shaped domain, i.e. -```math -\eta^2(\sigma_h) := \| \sigma_h - \nabla u_h \|^2_{L^2(T)} -``` -where ``\sigma_h`` discretisates the exact ``\sigma`` in the dual mixed problem -```math -\sigma - \nabla u = 0 -\quad \text{and} \quad -\mathrm{div}(\sigma) + f = 0 -``` -by some local equilibration strategy, see reference below for details. - -This examples demonstrates the use of low-level structures to assemble individual problems -and a strategy to solve several small problems in parallel by use of non-overlapping node patch groups. - -!!! reference - - ''A posteriori error estimates for efficiency and error control in numerical simulations'' - Lecture Notes by M. Vohralik - [>Link<](https://who.rocq.inria.fr/Martin.Vohralik/Enseig/APost/a_posteriori.pdf) - - -The resulting mesh and error convergence history for the default parameters looks like: - -![](example211.png) - -=# - -module Example211_LshapeAdaptiveEQPoissonProblem - -using ExtendableFEM -using ExtendableFEMBase -using ExtendableGrids -using ExtendableSparse -using GridVisualize -using Test #hide - -## exact solution u for the Poisson problem -function u!(result, qpinfo) - x = qpinfo.x - r2 = x[1]^2 + x[2]^2 - φ = atan(x[2], x[1]) - if φ < 0 - φ += 2 * pi - end - result[1] = r2^(1 / 3) * sin(2 * φ / 3) - return nothing -end - -## gradient of exact solution -function ∇u!(result, qpinfo) - x = qpinfo.x - φ = atan(x[2], x[1]) - r2 = x[1]^2 + x[2]^2 - if φ < 0 - φ += 2 * pi - end - ∂r = 2 / 3 * r2^(-1 / 6) * sin(2 * φ / 3) - ∂φ = 2 / 3 * r2^(-1 / 6) * cos(2 * φ / 3) - result[1] = cos(φ) * ∂r - sin(φ) * ∂φ - result[2] = sin(φ) * ∂r + cos(φ) * ∂φ - return nothing -end - -## kernel for exact error calculation -function exact_error!(result, u, qpinfo) - u!(result, qpinfo) - ∇u!(view(result, 2:3), qpinfo) - result .-= u - result .= result .^ 2 - return nothing -end - -## kernel for equilibration error estimator -function eqestimator_kernel!(result, input, qpinfo) - σ_h, divσ_h, ∇u_h = view(input, 1:2), input[3], view(input, 4:5) - result[1] = norm(σ_h .- ∇u_h)^2 + divσ_h^2 - return nothing -end - -## unknowns for primal and dual problem -u = Unknown("u"; name = "u") -σ = Unknown("σ"; name = "equilibrated fluxes / dual stress") - -## everything is wrapped in a main function -function main(; maxdofs = 4000, μ = 1, order = 2, nlevels = 16, θ = 0.5, Plotter = nothing, kwargs...) - - ## initial grid - xgrid = grid_lshape(Triangle2D) - - ## choose some finite elements for primal and dual problem (= for equilibrated fluxes) - FEType = H1Pk{1, 2, order} - FETypeDual = HDIVRTk{2, order} - - ## setup Poisson problem - PD = ProblemDescription("Poisson problem") - assign_unknown!(PD, u) - assign_operator!(PD, BilinearOperator([grad(u)]; factor = μ, kwargs...)) - assign_operator!(PD, InterpolateBoundaryData(u, u!; regions = 2:7, bonus_quadorder = 4, kwargs...)) - assign_operator!(PD, HomogeneousBoundaryData(u; regions = [1, 8])) - - ## define error estimator : || σ_h - ∇u_h ||^2 + || div σ_h ||^2 - EQIntegrator = ItemIntegrator(eqestimator_kernel!, [id(σ), div(σ), grad(u)]; resultdim = 1, quadorder = 2 * order) - - ## setup exact error evaluations - ErrorIntegrator = ItemIntegrator(exact_error!, [id(u), grad(u)]; quadorder = 2 * order, kwargs...) - - ## refinement loop (only uniform for now) - NDofs = zeros(Int, 0) - NDofsDual = zeros(Int, 0) - ResultsL2 = zeros(Float64, 0) - ResultsH1 = zeros(Float64, 0) - Resultsη = zeros(Float64, 0) - sol = nothing - η4cell = nothing - level = 0 - while (true) - level += 1 - - ## create a solution vector and solve the problem - FES = FESpace{FEType}(xgrid) - sol = solve(PD, FES) - push!(NDofs, length(view(sol[u]))) - println("\n SOLVE LEVEL $level") - println(" ndofs = $(NDofs[end])") - - ## evaluate equilibration error estimator and append it to sol vector (for plotting etc.) - local_equilibration_estimator!(sol, FETypeDual) - η4cell = evaluate(EQIntegrator, sol) - push!(Resultsη, sqrt(sum(view(η4cell, 1, :)))) - - ## calculate L2 error, H1 error, estimator, dual L2 error and write to results - push!(NDofsDual, length(view(sol[σ]))) - error = evaluate(ErrorIntegrator, sol) - push!(ResultsL2, sqrt(sum(view(error, 1, :)))) - push!(ResultsH1, sqrt(sum(view(error, 2, :)) + sum(view(error, 3, :)))) - println(" ESTIMATE") - println(" ndofsDual = $(NDofsDual[end])") - println(" estim H1 error = $(Resultsη[end])") - println(" exact H1 error = $(ResultsH1[end])") - println(" exact L2 error = $(ResultsL2[end])") - - if NDofs[end] >= maxdofs - break - end - - ## mesh refinement - if θ >= 1 ## uniform mesh refinement - xgrid = uniform_refine(xgrid) - else ## adaptive mesh refinement - facemarker = bulk_mark(xgrid, view(η4cell, :), θ; indicator_AT = ON_CELLS) - xgrid = RGB_refine(xgrid, facemarker) - end - end - - ## plot - plt = GridVisualizer(; Plotter = Plotter, layout = (2, 2), clear = true, resolution = (1000, 1000)) - scalarplot!(plt[1, 1], id(u), sol; levels = 11, title = "u_h") - plot_convergencehistory!(plt[1, 2], NDofs, [ResultsL2 ResultsH1 Resultsη]; add_h_powers = [order, order + 1], X_to_h = X -> order * X .^ (-1 / 2), ylabels = ["|| u - u_h ||", "|| ∇(u - u_h) ||", "η"]) - gridplot!(plt[2, 1], xgrid; linewidth = 1) - gridplot!(plt[2, 2], xgrid; linewidth = 1, xlimits = [-0.0005, 0.0005], ylimits = [-0.0005, 0.0005]) - - ## 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, η4cell], plt -end - -## this function computes the local equilibrated fluxes -## by solving local problems on (disjunct groups of) node patches -function local_equilibration_estimator!(sol, FETypeDual) - ## needed grid stuff - xgrid = sol[u].FES.xgrid - xCellNodes::Array{Int32, 2} = xgrid[CellNodes] - xCellVolumes::Array{Float64, 1} = xgrid[CellVolumes] - xNodeCells::Adjacency{Int32} = atranspose(xCellNodes) - nnodes::Int = num_sources(xNodeCells) - - ## get node patch groups that can be solved in parallel - group4node = xgrid[NodePatchGroups] - - ## init equilibration space (and Lagrange multiplier space) - FESDual = FESpace{FETypeDual}(xgrid) - xItemDofs::Union{VariableTargetAdjacency{Int32}, SerialVariableTargetAdjacency{Int32}, Array{Int32, 2}} = FESDual[CellDofs] - xItemDofs_uh::Union{VariableTargetAdjacency{Int32}, SerialVariableTargetAdjacency{Int32}, Array{Int32, 2}} = sol[u].FES[CellDofs] - - ## append block in solution vector for equilibrated fluxes - append!(sol, FESDual; tag = σ) - - ## partition of unity and their gradients = P1 basis functions - POUFES = FESpace{H1P1{1}}(xgrid) - POUqf = QuadratureRule{Float64, Triangle2D}(0) - - ## quadrature formulas - qf = QuadratureRule{Float64, Triangle2D}(2 * get_polynomialorder(FETypeDual, Triangle2D)) - weights::Array{Float64, 1} = qf.w - - ## some constants - offset::Int = sol[u].offset - div_penalty::Float64 = 1.0e5 # divergence constraint is realized by penalisation - bnd_penalty::Float64 = 1.0e60 # penalty for non-involved dofs of a group - maxdofs::Int = max_num_targets_per_source(xItemDofs) - maxdofs_uh::Int = max_num_targets_per_source(xItemDofs_uh) - - ## redistribute groups for more equilibrated thread load (first groups are larger) - maxgroups = maximum(group4node) - groups = Array{Int, 1}(1:maxgroups) - for j::Int in 1:floor(maxgroups / 2) - a = groups[j] - groups[j] = groups[2 * j] - groups[2 * j] = a - end - X = Array{Array{Float64, 1}, 1}(undef, maxgroups) - - function solve_patchgroup!(group) - ## temporary variables - graduh = zeros(Float64, 2) - coeffs_uh = zeros(Float64, maxdofs_uh) - Alocal = zeros(Float64, maxdofs, maxdofs) - blocal = zeros(Float64, maxdofs) - - ## init system - A = ExtendableSparseMatrix{Float64, Int64}(FESDual.ndofs, FESDual.ndofs) - b = zeros(Float64, FESDual.ndofs) - - ## init FEBasiEvaluators - FEE_∇φ = FEEvaluator(POUFES, Gradient, POUqf) - FEE_xref = FEEvaluator(POUFES, Identity, qf) - FEE_∇u = FEEvaluator(sol[u].FES, Gradient, qf) - FEE_div = FEEvaluator(FESDual, Divergence, qf) - FEE_id = FEEvaluator(FESDual, Identity, qf) - idvals = FEE_id.cvals - divvals = FEE_div.cvals - xref_vals = FEE_xref.cvals - ∇φvals = FEE_∇φ.cvals - - ## find dofs at boundary of current node patches - ## and in interior of cells outside of current node patch group - is_noninvolveddof = zeros(Bool, FESDual.ndofs) - outside_cell::Bool = false - for cell in 1:num_cells(xgrid) - outside_cell = true - for k in 1:3 - if group4node[xCellNodes[k, cell]] == group - outside_cell = false - break - end - end - if (outside_cell) # mark interior dofs of outside cell - for j in 1:maxdofs - is_noninvolveddof[xItemDofs[j, cell]] = true - end - end - end - - for node in 1:nnodes - if group4node[node] == group - for c in 1:num_targets(xNodeCells, node) - cell = xNodeCells[c, node] - - ## find local node number of global node z - ## and evaluate (constant) gradient of nodal basis function phi_z - localnode = 1 - while xCellNodes[localnode, cell] != node - localnode += 1 - end - FEE_∇φ.citem[] = cell - update_basis!(FEE_∇φ) - - ## read coefficients for discrete flux - for j in 1:maxdofs_uh - coeffs_uh[j] = sol.entries[offset + xItemDofs_uh[j, cell]] - end - - ## update other FE evaluators - FEE_∇u.citem[] = cell - FEE_div.citem[] = cell - FEE_id.citem[] = cell - update_basis!(FEE_∇u) - update_basis!(FEE_div) - update_basis!(FEE_id) - - ## assembly on this cell - for i in eachindex(weights) - weight = weights[i] * xCellVolumes[cell] - - ## evaluate grad(u_h) and nodal basis function at quadrature point - fill!(graduh, 0) - eval_febe!(graduh, FEE_∇u, coeffs_uh, i) - - ## compute residual -f*phi_z + grad(u_h) * grad(phi_z) at quadrature point i ( f = 0 in this example !!! ) - temp2 = div_penalty * sqrt(xCellVolumes[cell]) * weight - temp = temp2 * dot(graduh, view(∇φvals, :, localnode, 1)) - for dof_i in 1:maxdofs - ## right-hand side for best-approximation (grad(u_h)*phi) - blocal[dof_i] += dot(graduh, view(idvals, :, dof_i, i)) * xref_vals[1, localnode, i] * weight - ## mass matrix Hdiv - for dof_j in dof_i:maxdofs - Alocal[dof_i, dof_j] += dot(view(idvals, :, dof_i, i), view(idvals, :, dof_j, i)) * weight - end - ## div-div matrix Hdiv * penalty (quick and dirty to avoid Lagrange multiplier) - blocal[dof_i] += temp * divvals[1, dof_i, i] - temp3 = temp2 * divvals[1, dof_i, i] - for dof_j in dof_i:maxdofs - Alocal[dof_i, dof_j] += temp3 * divvals[1, dof_j, i] - end - end - end - - ## write into global A and b - for dof_i in 1:maxdofs - dofi = xItemDofs[dof_i, cell] - b[dofi] += blocal[dof_i] - for dof_j in 1:maxdofs - dofj = xItemDofs[dof_j, cell] - if dof_j < dof_i # use that Alocal is symmetric - _addnz(A, dofi, dofj, Alocal[dof_j, dof_i], 1) - else - _addnz(A, dofi, dofj, Alocal[dof_i, dof_j], 1) - end - end - end - - ## reset local A and b - fill!(Alocal, 0) - fill!(blocal, 0) - end - end - end - - ## penalize dofs that are not involved - for j in 1:FESDual.ndofs - if is_noninvolveddof[j] - A[j, j] = bnd_penalty - b[j] = 0 - end - end - - ## solve local problem - return A \ b - end - - ## solve equilibration problems on vertex patches (in parallel) - Threads.@threads for group in groups - grouptime = @elapsed begin - @info " Starting equilibrating patch group $group on thread $(Threads.threadid())... " - X[group] = solve_patchgroup!(group) - end - @info "Finished equilibration patch group $group on thread $(Threads.threadid()) in $(grouptime)s " - end - - ## write local solutions to global vector (sequentially) - for group in 1:maxgroups - view(sol[σ]) .+= X[group] - end - return -end - -generateplots = ExtendableFEM.default_generateplots(Example211_LshapeAdaptiveEQPoissonProblem, "example211.png") #hide -function runtests() #hide - results, plt = main(; maxdofs = 21, order = 2) #hide - @test length(results[1].entries) == 96 #hide - @test results[2] ≈ [0.0005948849237161765 0.02155094069716629 0.0012248396092485094 0.01840063545660059 0.02155094068517524 0.0005948849232996949] - return nothing #hide -end #hide -end diff --git a/examples/Example212_PeriodicElasticity2D.jl b/examples/Example212_PeriodicElasticity2D.jl deleted file mode 100644 index 3b6487a5..00000000 --- a/examples/Example212_PeriodicElasticity2D.jl +++ /dev/null @@ -1,200 +0,0 @@ -#= - -# 212 : Periodic Boundary 2D -([source code](@__SOURCE_URL__)) - -This is a simple demonstration and validation of the generic periodic boundary operator. - -We construct an unstructured periodic 2D grid and solve a simple linear elastic problem -with periodic boundary coupling along the x-axis. - -![](example212.png) -=# - -module Example212_PeriodicElasticity2D - -using ExtendableFEM -using ExtendableGrids -using SimplexGridFactory -using GridVisualize -using Triangulate -using UnicodePlots -using StaticArrays -using LinearAlgebra -using Test #hide -using SparseArrays - -## enumerate the boundary regions -const reg_left = 4 -const reg_right = 2 -const reg_dirichlet = 1 -const reg_default = 3 - -## 2D reduction of the material used in Example 312 -## in Voigt notation -function material_tensor() - c11 = 396.0 - c12 = 137.0 - c44 = 116.0 - - return @SArray [ - c11 c12 0 - c12 c11 0 - 0 0 c44 - ] -end - -## generate the kernels for the linear problem -## 𝐂: Hooke tensor, 𝑓: body force -function make_kernels(𝐂, 𝑓) - - ## linear stress-strain mapping - bilinear_kernel!(σ, εv, qpinfo) = mul!(σ, 𝐂, εv) - - ## plain body force - linear_kernel!(result, qpinfo) = (result .= 𝑓) - - return bilinear_kernel!, linear_kernel! -end - - -""" - create a 2D grid with Dirichlet boundary region at the bottom center -""" -function create_grid(; h, height, width) - builder = SimplexGridBuilder(; Generator = Triangulate) - - ## bottom points - b1 = point!(builder, 0, 0) - b2 = point!(builder, 0.45 * width, 0) - b3 = point!(builder, 0.5 * width, 0) - b4 = point!(builder, 0.55 * width, 0) - b5 = point!(builder, width, 0) - - ## top points - t1 = point!(builder, 0, height) - t2 = point!(builder, width / 2, height) - t3 = point!(builder, width, height) - - ## default faces - facetregion!(builder, reg_default) - facet!(builder, b1, b2) - facet!(builder, b4, b5) - facet!(builder, t1, t2) - facet!(builder, t2, t3) - - ## left face - facetregion!(builder, reg_left) - facet!(builder, b1, t1) - - ## right face - facetregion!(builder, reg_right) - facet!(builder, b5, t3) - - ## Dirichlet face - facetregion!(builder, reg_dirichlet) - facet!(builder, b3, b4) - facet!(builder, b2, b3) - - ## divider - facetregion!(builder, reg_default) - facet!(builder, t2, b3) - - - cellregion!(builder, 1) - maxvolume!(builder, h) - regionpoint!(builder, width / 3, height / 2) - - cellregion!(builder, 2) - # much finer grid on the right half to make periodic coupling non-trivial - maxvolume!(builder, 0.1 * h) - regionpoint!(builder, 2width / 3, height / 2) - - return simplexgrid(builder) -end - -function main(; - order = 1, - periodic = true, - use_LM_restrictions = true, - Plotter = nothing, - force = 10.0, - h = 5.0e-2, - width = 6.0, - height = 1.0, - threads = 1, - kwargs... - ) - ## print options for better logs - @info "selected options" use_LM_restrictions threads - - xgrid = create_grid(; h, width, height) - - ## create finite element space and solution vector - if order == 1 - FES = FESpace{H1P1{2}}(xgrid) - elseif order == 2 - FES = FESpace{H1P2{2, 2}}(xgrid) - end - - ## problem description - PD = ProblemDescription() - u = Unknown("u"; name = "displacement") - assign_unknown!(PD, u) - - 𝐂 = material_tensor() - 𝑓 = force * [0, 1] - - bilinear_kernel!, linear_kernel! = make_kernels(𝐂, 𝑓) - assign_operator!(PD, BilinearOperator(bilinear_kernel!, [εV(u, 1.0)]; kwargs...)) - assign_operator!(PD, LinearOperator(linear_kernel!, [id(u)]; kwargs...)) - - assign_operator!(PD, HomogeneousBoundaryData(u; regions = [reg_dirichlet])) - - if periodic - function give_opposite!(y, x) - y .= x - y[1] = width - x[1] - return nothing - end - - @showtime coupling_matrix = get_periodic_coupling_matrix(FES, reg_left, reg_right, give_opposite!; parallel = threads > 1, threads) - if use_LM_restrictions - assign_restriction!(PD, CoupledDofsRestriction(coupling_matrix)) - else - assign_operator!(PD, CombineDofs(u, u, coupling_matrix; kwargs...)) - end - end - - sol, SC = solve(PD, FES, return_config = true) - residual(SC) < 1.0e-10 || error("Residual is not zero!") - - if use_LM_restrictions - @info "Lagrange residuals" SC.statistics[:restriction_residuals] - end - - plt = GridVisualizer(; Plotter, size = (1300, 800)) - - magnification = 1 - displaced_grid = deepcopy(xgrid) - displace_mesh!(displaced_grid, sol[1], magnify = magnification) - gridplot!(plt, displaced_grid, linewidth = 1, title = "displaced mesh, $(magnification)x magnified", scene3d = :LScene) - - return sol, plt -end - -generateplots = ExtendableFEM.default_generateplots(Example212_PeriodicElasticity2D, "example212.png") #hide -function runtests() #hide - sol1, _ = main(use_LM_restrictions = false, threads = 1) #hide - @test abs(maximum(view(sol1[1])) - 1.2638113215466982) < 1.0e-3 #hide - - sol2, _ = main(use_LM_restrictions = false, threads = 4) #hide - @test sol1.entries ≈ sol2.entries #hide - - sol3, _ = main(use_LM_restrictions = true, threads = 4) #hide - @test sol1.entries ≈ sol3.entries #hide - - return nothing #hide -end #hide - -end # module diff --git a/examples/Example220_ReactionConvectionDiffusion.jl b/examples/Example220_ReactionConvectionDiffusion.jl deleted file mode 100644 index 2fe84484..00000000 --- a/examples/Example220_ReactionConvectionDiffusion.jl +++ /dev/null @@ -1,151 +0,0 @@ -#= - -# 220 : Reaction-Convection-Diffusion-Problem -([source code](@__SOURCE_URL__)) - -This example computes the solution of some convection-diffusion problem -```math --\nu \Delta u + \mathbf{\beta} \cdot \nabla u + \alpha u = f \quad \text{in } \Omega -``` -with some diffusion coefficient ``\nu``, some vector-valued function ``\mathbf{\beta}``, some scalar-valued function ``\alpha`` and inhomogeneous Dirichlet boundary data. - -We prescribe an analytic solution with ``\mathbf{\beta} := (1,0)`` and ``\alpha = 0.1`` and check the L2 and H1 error convergence of the method on a series of uniformly refined meshes. -We also compare with the error of a simple nodal interpolation and plot the solution and the norm of its gradient. - -For small ``\nu``, the convection term dominates and pollutes the accuracy of the method. For demonstration some -simple gradient jump (interior penalty) stabilisation is added to improve things. - -The computed solution for the default parameters looks like this: - -![](example220.png) - -=# - -module Example220_ReactionConvectionDiffusion - -using ExtendableFEM -using ExtendableGrids -using LinearAlgebra -using Test #hide - -const α = 0.01 -const β = [1.0, 0] -const ν = 1.0e-5 - -function u!(result, qpinfo) - x = qpinfo.x - result[1] = x[1] * x[2] * (x[1] - 1) * (x[2] - 1) + x[1] - return nothing -end -function ∇u!(result, qpinfo) - x = qpinfo.x - result[1] = x[2] * (2 * x[1] - 1) * (x[2] - 1) + 1 - result[2] = x[1] * (2 * x[2] - 1) * (x[1] - 1) - return nothing -end -function Δu!(result, qpinfo) - x = qpinfo.x - result[1] = 2 * (x[2] * (x[2] - 1) + x[1] * (x[1] - 1)) - return nothing -end - -function rhs() - ∇u = zeros(Float64, 2) - Δu = zeros(Float64, 1) - u = zeros(Float64, 1) - return function closure(result, qpinfo) - ∇u!(∇u, qpinfo) - u!(u, qpinfo) - Δu!(Δu, qpinfo) - result[1] = -ν * Δu[1] + α * u[1] + dot(β, ∇u) - return nothing - end -end - -function kernel_DCR!(result, input, qpinfo) - u, ∇u = view(input, 1), view(input, 2:3) - result[1] = α * u[1] + dot(β, ∇u) - result[2] = ν * ∇u[1] - result[3] = ν * ∇u[2] - return nothing -end - -## kernel for exact error calculation -function exact_error!(result, u, qpinfo) - u!(result, qpinfo) - ∇u!(view(result, 2:3), qpinfo) - result .-= u - result .= result .^ 2 - return nothing -end - -## stab_kernel! -function stab_kernel!(result, ∇u, qpinfo) - result .= ∇u .* qpinfo.volume^2 - return nothing -end - - -## everything is wrapped in a main function -function main(; Plotter = nothing, τ = 1.0e-2, nlevels = 5, order = 2, kwargs...) - - ## create problem description - PD = ProblemDescription("reaction-convection-diffusion problem") - u = Unknown("u") - assign_unknown!(PD, u) - assign_operator!(PD, BilinearOperator(kernel_DCR!, [id(u), grad(u)]; bonus_quadorder = 1, kwargs...)) - assign_operator!(PD, LinearOperator(rhs(), [id(u)]; bonus_quadorder = 2, kwargs...)) - assign_operator!(PD, InterpolateBoundaryData(u, u!; regions = 1:4, kwargs...)) - - ## add a gradient jump (interior penalty) stabilisation for dominant convection - if τ > 0 - assign_operator!(PD, BilinearOperatorDG(stab_kernel!, [jump(grad(u))]; entities = ON_IFACES, factor = τ)) - end - - ## prepare error calculation - ErrorIntegratorExact = ItemIntegrator(exact_error!, [id(1), grad(1)]; quadorder = 2 * order, kwargs...) - Results = zeros(Float64, nlevels, 4) - NDofs = zeros(Int, nlevels) - - ## refinement loop over levels - sol = nothing - xgrid = grid_unitsquare(Triangle2D) # initial grid - for level in 1:nlevels - ## uniform mesh refinement - xgrid = uniform_refine(xgrid) - - ## generate FESpace and solve - FES = FESpace{H1Pk{1, 2, order}}(xgrid) - sol = solve(PD, FES) - - ## compute L2 and H1 errors and save data - NDofs[level] = length(sol.entries) - error = evaluate(ErrorIntegratorExact, sol) - Results[level, 1] = sqrt(sum(view(error, 1, :))) - Results[level, 3] = sqrt(sum(view(error, 2, :)) + sum(view(error, 3, :))) - - ## interpolate (just for comparison) - I = FEVector(FES) - interpolate!(I[1], u!) - error = evaluate(ErrorIntegratorExact, I) - Results[level, 2] = sqrt(sum(view(error, 1, :))) - Results[level, 4] = sqrt(sum(view(error, 2, :)) + sum(view(error, 3, :))) - end - - ## plot - plt = plot([id(u), grad(u)], sol; add = 1, ncols = 3, Plotter = Plotter) - plot_convergencehistory!(plt[1, 3], NDofs, Results; add_h_powers = [order, order + 1], X_to_h = X -> X .^ (-1 / 2), legend = :lb, ylabels = ["|| u - u_h ||", "|| u - Iu ||", "|| ∇(u - u_h) ||", "|| ∇(u - Iu) ||"], limits = (1.0e-8, 1.0e-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) ||"]) - - return Results, plt -end - -generateplots = ExtendableFEM.default_generateplots(Example220_ReactionConvectionDiffusion, "example220.png") #hide -function runtests() #hide - Results, plt = main(; nlevels = 2) #hide - @test Results[end, 1] ≈ 0.0001510021661291585 #hide - return nothing #hide -end #hide -end # module diff --git a/examples/Example225_ObstacleProblem.jl b/examples/Example225_ObstacleProblem.jl deleted file mode 100644 index 7fb67f21..00000000 --- a/examples/Example225_ObstacleProblem.jl +++ /dev/null @@ -1,80 +0,0 @@ -#= - -# 225 : Obstacle Problem -([source code](@__SOURCE_URL__)) - -This example computes the solution ``u`` of the nonlinear obstacle problem that seeks the minimiser of the energy functional -```math -\begin{aligned} - E(u) = \frac{1}{2} \int_\Omega \lvert \nabla u \rvert^2 dx - \int_\Omega f u dx -\end{aligned} -``` -with some right-hand side ``f`` within the set of admissible functions that lie above an obstacle ``\chi`` -```math -\begin{aligned} - \mathcal{K} := \lbrace u \in H^1_0(\Omega) : u \geq \chi \rbrace. -\end{aligned} -``` - -The obstacle constraint is realised via a penalty term -```math -\begin{aligned} - \frac{1}{\epsilon} \| \min(0, u - \chi) \|^2_{L^2} -\end{aligned} -``` -that is added to the energy above and is automatically differentiated for a Newton scheme. -The computed solution for the default parameters looks like this: - -![](example225.png) -=# - -module Example225_ObstacleProblem - -using ExtendableFEM -using ExtendableGrids -using Test #hide - -## define obstacle and penalty kernel -const χ! = (result, x) -> (result[1] = (cos(4 * x[1] * π) * cos(4 * x[2] * π) - 1) / 20) -function obstacle_penalty_kernel!(result, input, qpinfo) - χ!(result, qpinfo.x) # eval obstacle - result[1] = min(0, input[1] - result[1]) - return nothing -end - -function main(; Plotter = nothing, ϵ = 1.0e-4, nrefs = 6, order = 1, parallel = false, npart = 8, kwargs...) - - ## choose initial mesh - xgrid = uniform_refine(grid_unitsquare(Triangle2D), nrefs) - if parallel - xgrid = partition(xgrid, RecursiveMetisPartitioning(npart = npart)) - end - - ## problem description - PD = ProblemDescription() - u = Unknown("u"; name = "potential") - assign_unknown!(PD, u) - assign_operator!(PD, NonlinearOperator(obstacle_penalty_kernel!, [id(u)]; factor = 1 / ϵ, parallel = parallel, kwargs...)) - assign_operator!(PD, BilinearOperator([grad(u)]; store = true, parallel = parallel, kwargs...)) - assign_operator!(PD, LinearOperator([id(u)]; store = true, parallel = parallel, factor = -1, kwargs...)) - assign_operator!(PD, HomogeneousBoundaryData(u; regions = 1:4, kwargs...)) - - ## create finite element space - FES = FESpace{H1Pk{1, 2, order}}(xgrid) - - ## solve - sol = solve(PD, FES; kwargs...) - - ## plot - plt = plot([id(u), grad(u)], sol; Plotter = Plotter, ncols = 3) - - return sol, plt -end - -generateplots = ExtendableFEM.default_generateplots(Example225_ObstacleProblem, "example225.png") #hide -function runtests() #hide - sol, plt = main(; μ = 1.0, nrefs = 2, order = 2) #hide - @test maximum(sol.entries) ≈ 0.0033496680638875204 #hide - return nothing #hide -end #hide -end # module diff --git a/examples/Example226_Thermoforming.jl b/examples/Example226_Thermoforming.jl deleted file mode 100644 index 8b3dcf4c..00000000 --- a/examples/Example226_Thermoforming.jl +++ /dev/null @@ -1,149 +0,0 @@ -#= - -# 226 : Thermoforming -([source code](@__SOURCE_URL__)) - -This implements the thermoforming example taken from https://arxiv.org/abs/1802.03564 Section 6.4. -The computed solution for the default parameters looks like this: - -![](example226.png) -=# - -module Example226_Thermoforming - -using ExtendableFEM -using ExtendableGrids -using SparseArrays -using LinearAlgebra -using Test #hide - -function w(r) - if 0.1 ≤ r ≤ 0.3 - return 5.0 * r - 0.5 - elseif 0.3 < r < 0.7 - return 1.0 - elseif 0.7 <= r <= 0.9 - return 4.5 - 5.0 * r - else - return 0.0 - end -end - -## initial mould -function Φ0(x) - return w(x[1]) * w(x[2]) -end - - -function g(r, κ, s) - if r <= 0.0 - return κ - elseif r <= 0.25 * s - return κ - 8.0 * κ * r^2 / (3.0 * s^2) - elseif r <= 0.75 * s - return 7.0 / 6.0 * κ - 4.0 / 3.0 * κ * r / s - elseif r <= s - return 8.0 / 3.0 * (s - r)^2 / s^2 - else - return 0.0 - end -end - - -## The smooth bump function in [0,1] -bump(x) = (0.0 <= x <= 1.0) ? exp(-0.25 / (x - x^2)) : 0.0 - -## Bump in [0,1]^2 -bumpInUnitSquare(x) = begin - r = sqrt((x[1] - 0.5)^2 + (x[2] - 0.5)^2) - return bump(0.5 + r) -end - - -## nonlinear kernel -function nonlinear_kernel!(result, input, qpinfo) - ## results and input contain 7 variables (u,∇u,T,∇T,y) - u = view(input, 1) - ∇u = view(input, 2:3) - T = view(input, 4) - ∇T = view(input, 5:6) - y = view(input, 7) - - α = qpinfo.params[1] - k = qpinfo.params[2] - f = qpinfo.params[3] - β = qpinfo.params[4] - κ = qpinfo.params[5] - s = qpinfo.params[6] - - result[1] = α * max(0, u[1] - y[1]) - f # pattern: 1 7 - result[2:3] = ∇u # pattern: 2 / 3 - result[4] = k * T[1] - g(y[1] - u[1], κ, s) # pattern: 1 4 7 - result[5:6] = ∇T # pattern: 5 / 6 - result[7] = y[1] - Φ0(qpinfo.x) - β * bumpInUnitSquare(qpinfo.x) * T[1] # pattern: 4 7 - return nothing -end - -## custom sparsity pattern for the jacobians of the nonlinear_kernel (Symbolcs cannot handle conditional jumps) -## note: jacobians are defined row-wise -rows = [1, 1, 2, 3, 4, 4, 4, 5, 6, 7, 7] -cols = [1, 7, 2, 3, 1, 4, 7, 5, 6, 4, 7] -vals = ones(Bool, length(cols)) -sparsity_pattern = sparse(rows, cols, vals) - -function main(; - κ = 10, - s = 1, - α = 1.0e8, - k = 1, - β = 5.25e-3, - f = 100, - N = 32, - order = 1, - Plotter = nothing, - kwargs... - ) - - ## choose mesh, - h = 1 / (N + 1) - xgrid = simplexgrid(0:h:1, 0:h:1) - - ## problem description - PD = ProblemDescription() - u = Unknown("u"; name = "membrane position") - y = Unknown("y"; name = "mould") - T = Unknown("T"; name = "temperature") - assign_unknown!(PD, u) - assign_unknown!(PD, y) - assign_unknown!(PD, T) - assign_operator!(PD, NonlinearOperator(nonlinear_kernel!, [id(u), grad(u), id(T), grad(T), id(y)]; bonus_quadorder = 2, params = [α, k, f, β, κ, s], sparse_jacobians_pattern = sparsity_pattern, kwargs...)) - assign_operator!(PD, HomogeneousBoundaryData(u; regions = 1:4, kwargs...)) - assign_operator!(PD, HomogeneousBoundaryData(y; regions = 1:4, kwargs...)) - - ## create finite element space - FES = FESpace{H1Pk{1, 2, order}}(xgrid) - FESs = [FES, FES, FES] - sol = FEVector(FESs; tags = [u, y, T]) - - ## initial guess for Newton - interpolate!(sol[u], (result, qpinfo) -> (result[1] = 0.9 * Φ0(qpinfo.x))) - interpolate!(sol[T], (result, qpinfo) -> (result[1] = 0.2)) - interpolate!(sol[y], (result, qpinfo) -> (result[1] = 10.0)) - - ## solve - sol = solve(PD, FESs; init = sol, maxiterations = 420, target_residual = 1.0e-8, kwargs...) - - ## plot - plt = plot([id(u), id(T), id(y)], sol; Plotter = Plotter) - - return sol, plt -end - -generateplots = ExtendableFEM.default_generateplots(Example226_Thermoforming, "example226.png") #hide -function runtests() #hide - sol1, ~ = Example226_Thermoforming.main(; sparse_jacobians = true, N = 20) #hide - sol2, ~ = Example226_Thermoforming.main(; sparse_jacobians = false, N = 20) #hide - @test norm(sol1.entries - sol2.entries) < 1.0e-10 #hide - return nothing #hide -end #hide -end # module diff --git a/examples/Example227_ObstacleProblemLVPP.jl b/examples/Example227_ObstacleProblemLVPP.jl deleted file mode 100644 index cf169ee5..00000000 --- a/examples/Example227_ObstacleProblemLVPP.jl +++ /dev/null @@ -1,164 +0,0 @@ -#= - -# 227 : Obstacle Problem LVPP -([source code](@__SOURCE_URL__)) - -This example computes the solution ``u`` of the nonlinear obstacle problem that seeks -the minimiser of the energy functional -```math -\begin{aligned} - E(u) = \frac{1}{2} \int_\Omega \lvert \nabla u \rvert^2 dx - \int_\Omega f u dx -\end{aligned} -``` -with some right-hand side ``f`` within the set of admissible functions that lie above an obstacle ``\chi`` -```math -\begin{aligned} - \mathcal{K} := \lbrace u \in H^1_0(\Omega) : u \geq \chi \rbrace. -\end{aligned} -``` - -Opposite to Example225 the solution is computed by the latent variable proximal point (LVPP) method -that solves the problem via a series of nonlinear mixed problems that guarantee the decay of the -energy. Given ``\alpha_k`` and initial guesses ``u_0`` and ``\psi_0``, the subproblem for ``k \geq 1`` -seeks a solution ``u_{k} \in V := H^1_0(\Omega)`` and ``\psi_{k} \in W := L^\infty(\Omega)`` such that -```math -\begin{aligned} -\alpha_k (\nabla u_k, \nabla v)_{L^2} + (\psi_k, v)_{L^2} & = (\alpha_k f + \psi_{k-1},v)_{L^2} -&& \text{for all } v \in V\\ -(u_k, w)_{L^2} - (\chi + \exp(\psi_k), w)_{L^2} & = 0 -&& \text{for all } w \in W -\end{aligned} -``` -The parameter ``\alpha_k`` is initialized with ``\alpha_0 = 1`` and updated according to -``\alpha_k = \min(\max(r^(q^k) - α), 10^3)`` with ``r = q = 1.5``. The problem for each ``k`` -is solved by the Newton method. This implements Algorithm 3 in the reference below. - - -!!! reference - - ''Proximal Galerkin: A Structure-Preserving Finite Element Method for Pointwise Bound Constraints'', - Brendan Keith, Thomas M. Surowiec, Found Comput Math (2024), - [>Link<](https://doi.org/10.1007/s10208-024-09681-8) - - -![](example227.png) -=# - -module Example227_ObstacleProblemLVPP - -using ExtendableFEM -using ExtendableFEMBase -using ExtendableGrids -using LinearAlgebra -using Metis -using Test #hide - -## define obstacle -const b = 9 // 20 -const d = sqrt(1 // 4 - b^2) -function χ(x) - r = sqrt(x[1]^2 + x[2]^2) - return r <= b ? sqrt(1 // 4 - r^2) : d + b^2 / d - b * r / d -end - -## transformation of latent variable ψ to constrained variable u -function ∇R!(result, input, qpinfo) - return result[1] = χ(qpinfo.x) + exp(input[1]) -end - -## boundary data for latent variable ψ (such that ̃u := χ + ∇R(ψ) satisfies Dirichlet boundary conditions) -function bnd_ψ!(result, qpinfo) - return result[1] = log(-χ(qpinfo.x)) -end - -function main(; - nrefs = 5, - α0 = 1.0, - order = 1, - parallel = false, - npart = 8, - tol = 1.0e-12, - Plotter = nothing, - kwargs... - ) - - ## choose initial mesh - xgrid = uniform_refine(grid_unitsquare(Triangle2D; scale = (2, 2), shift = (-0.5, -0.5)), nrefs) - if parallel - xgrid = partition(xgrid, RecursiveMetisPartitioning(npart = npart)) - end - - ## create finite element space - FES = [FESpace{H1Pk{1, 2, order}}(xgrid), FESpace{H1Pk{1, 2, order}}(xgrid)] - - ## init proximal parameter - α = α0 - - ## prepare Laplacian and mass matrix - L = FEMatrix(FES[1], FES[1]) - assemble!(L, BilinearOperator([grad(1)], [grad(1)]; parallel = parallel)) - function scaled_laplacian!(A, b, args; assemble_matrix = true, assemble_rhs = true, kwargs...) - return if assemble_matrix - ## add Laplacian scaled by α - ExtendableFEMBase.add!(A, L.entries; factor = α) - end - end - M = FEMatrix(FES[1], FES[2]) - b = FEVector(FES[1]) - assemble!(M, BilinearOperator([id(1)], [id(1)]; parallel = parallel)) - - ## problem description - PD = ProblemDescription() - u = Unknown("u"; name = "solution") - ψ = Unknown("ψ"; name = "latent variable") - assign_unknown!(PD, u) - assign_unknown!(PD, ψ) - assign_operator!(PD, CallbackOperator(scaled_laplacian!, [u]; kwargs...)) - assign_operator!(PD, BilinearOperator([id(u)], [(id(ψ))]; transposed_copy = 1, store = true, parallel = parallel, kwargs...)) - assign_operator!(PD, NonlinearOperator(∇R!, [id(ψ)], [id(ψ)]; parallel = parallel, factor = -1, kwargs...)) - assign_operator!(PD, HomogeneousBoundaryData(u; regions = 1:4, kwargs...)) - assign_operator!(PD, InterpolateBoundaryData(ψ, bnd_ψ!; regions = 1:4, kwargs...)) - assign_operator!(PD, LinearOperator(b, [u]; kwargs...)) - - ## solve - sol = FEVector(FES; tags = PD.unknowns) - sol_prev = FEVector(FES; tags = PD.unknowns) - SC = nothing - r, q = 3 // 2, 3 // 2 - converged = false - k = 0 - while !converged - k += 1 - @info "Step $k: α = $(α)" - - ## save previous solution and update right-hand side - b.entries .= M.entries * view(sol[ψ]) - sol_prev.entries .= sol.entries - - ## solve nonlinear problem - sol, SC = solve(PD, FES, SC; init = sol, maxiterations = 20, verbosity = -1, timeroutputs = :hide, return_config = true, kwargs...) - niterations = length(ExtendableFEM.residuals(SC)) - - ## compute distance - dist = norm(view(sol[u]) .- view(sol_prev[u])) - @info "dist = $dist, niterations = $(niterations - 1)" - if dist < tol - converged = true - else ## increase proximal parameter - α = min(max(r^(q^k) - α), 10^3) - end - end - - ## postprocess latent variable v = ϕ + exp ψ = ∇R!(ψ) - u2 = Unknown("̃u"; name = "solution 2") - append!(sol, FES[1]; tag = u2) - lazy_interpolate!(sol[u2], sol, [id(ψ)]; postprocess = ∇R!) - - ## plot - plt = plot([id(u), id(ψ), id(u2)], sol; Plotter = Plotter, ncols = 3) - - return sol, plt -end - -generateplots = ExtendableFEM.default_generateplots(Example227_ObstacleProblemLVPP, "example227.png") #hide -end # module diff --git a/examples/Example230_NonlinearElasticity.jl b/examples/Example230_NonlinearElasticity.jl deleted file mode 100644 index 7fbf80ab..00000000 --- a/examples/Example230_NonlinearElasticity.jl +++ /dev/null @@ -1,183 +0,0 @@ -#= - -# 230 : Nonlinear Elasticity -([source code](@__SOURCE_URL__)) - -This example computes the displacement field ``u`` of the nonlinear elasticity problem -```math -\begin{aligned} --\mathrm{div}(\mathbb{C} (\epsilon(u)-\epsilon_T)) & = 0 \quad \text{in } \Omega -\end{aligned} -``` -where an isotropic stress tensor ``\mathbb{C}`` is applied to the nonlinear strain ``\epsilon(u) := \frac{1}{2}(\nabla u + (\nabla u)^T + (\nabla u)^T \nabla u)`` -and a misfit strain ``\epsilon_T := \Delta T \alpha`` due to thermal load caused by temperature(s) ``\Delta T`` and thermal expansion coefficients ``\alpha`` (that may be different) -in the two regions of the bimetal. - -This example demonstrates how to setup a (parameter- and region-dependent) nonlinear expression and how to assign it to the problem description. - -The computed solution for the default parameters looks like this: - -![](example230.png) -=# - -module Example230_NonlinearElasticity - -using ExtendableFEM -using ExtendableFEMBase -using ExtendableGrids -using GridVisualize -using UnicodePlots -using Test #hide - -## parameter-dependent nonlinear operator uses a callable struct to reduce allocations -mutable struct nonlinear_operator{T} - λ::Vector{T} - μ::Vector{T} - ϵT::Vector{T} -end - -function strain!(result, input) - result[1] = input[1] - result[2] = input[4] - result[3] = input[2] + input[3] - - ## add nonlinear part of the strain 1/2 * (grad(u)'*grad(u)) - result[1] += 1 // 2 * (input[1]^2 + input[3]^2) - result[2] += 1 // 2 * (input[2]^2 + input[4]^2) - result[3] += input[1] * input[2] + input[3] * input[4] - return nothing -end - -## kernel for nonlinear operator -(op::nonlinear_operator)(result, input, qpinfo) = ( - ## input = grad(u) written as a vector - ## compute strain and subtract thermal strain (all in Voigt notation) - region = qpinfo.region; - strain!(result, input); - result[1] -= op.ϵT[region]; - result[2] -= op.ϵT[region]; - - ## multiply with isotropic stress tensor - ## (stored in input[5:7] using Voigt notation) - a = op.λ[region] * (result[1] + result[2]) + 2 * op.μ[region] * result[1]; - b = op.λ[region] * (result[1] + result[2]) + 2 * op.μ[region] * result[2]; - c = 2 * op.μ[region] * result[3]; - - ## write strain into result - result[1] = a; - result[2] = c; - result[3] = c; - result[4] = b; - return nothing -) - -const op = nonlinear_operator([0.0, 0.0], [0.0, 0.0], [0.0, 0.0]) - -## everything is wrapped in a main function -function main(; - ν = [0.3, 0.3], # Poisson number for each region/material - E = [2.1, 1.1], # Elasticity modulus for each region/material - ΔT = [580, 580], # temperature for each region/material - α = [1.3e-5, 2.4e-4], # thermal expansion coefficients - scale = [20, 500], # scale of the bimetal, i.e. [thickness, width] - nrefs = 0, # refinement levels - order = 2, # finite element order - periodic = false, # use periodic boundary conditions? - Plotter = nothing, - kwargs... - ) - - ## compute Lame' coefficients μ and λ from ν and E - ## and thermal misfit strain and assign to operator operator - @. op.μ = E / (2 * (1 + ν)) - @. op.λ = E * ν / ((1 - 2 * ν) * (1 + ν)) - @. op.ϵT = ΔT * α - - ## generate bimetal mesh - xgrid = bimetal_strip2D(; scale = scale, n = 2 * (nrefs + 1)) - println(stdout, unicode_gridplot(xgrid)) - - ## create finite element space and solution vector - FES = FESpace{H1Pk{2, 2, order}}(xgrid) - - ## problem description - PD = ProblemDescription() - u = Unknown("u"; name = "displacement") - assign_unknown!(PD, u) - assign_operator!(PD, NonlinearOperator(op, [grad(u)]; kwargs...)) - if periodic - ## periodic boundary conditions - ## 1) couple dofs left (bregion 1) and right (bregion 3) in y-direction - function give_opposite!(y, x) - y .= x - y[1] = -x[1] - return nothing - end - coupling_matrix = get_periodic_coupling_matrix(FES, xgrid, 1, 3, give_opposite!; mask = [0, 1]) - assign_operator!(PD, CombineDofs(u, u, coupling_matrix; kwargs...)) - - ## 2) find and fix point at [0, scale[1]] - xCoordinates = xgrid[Coordinates] - closest::Int = 0 - mindist::Float64 = 1.0e30 - for j in 1:num_nodes(xgrid) - dist = xCoordinates[1, j]^2 + (xCoordinates[2, j] - scale[1])^2 - if dist < mindist - mindist = dist - closest = j - end - end - assign_operator!(PD, FixDofs(u; dofs = [closest], vals = [0])) - else - assign_operator!(PD, HomogeneousBoundaryData(u; regions = [1], mask = [1, 0], kwargs...)) - end - - ## solve - sol = solve(PD, FES; kwargs...) - - ## displace mesh and plot - plt = GridVisualizer(; Plotter = Plotter, layout = (3, 1), clear = true, size = (1000, 1500)) - grad_nodevals = nodevalues(grad(u), sol) - strain_nodevals = zeros(Float64, 3, num_nodes(xgrid)) - for j in 1:num_nodes(xgrid) - strain!(view(strain_nodevals, :, j), view(grad_nodevals, :, j)) - end - scalarplot!(plt[1, 1], xgrid, view(strain_nodevals, 1, :), levels = 3, colorbarticks = 7, xlimits = [-scale[2] / 2 - 10, scale[2] / 2 + 10], ylimits = [-30, scale[1] + 20], title = "ϵ(u)_xx + displacement") - scalarplot!(plt[2, 1], xgrid, view(strain_nodevals, 2, :), levels = 1, colorbarticks = 7, xlimits = [-scale[2] / 2 - 10, scale[2] / 2 + 10], ylimits = [-30, scale[1] + 20], title = "ϵ(u)_yy + displacement") - vectorplot!(plt[1, 1], xgrid, eval_func_bary(PointEvaluator([id(u)], sol)), rasterpoints = 20, clear = false) - vectorplot!(plt[2, 1], xgrid, eval_func_bary(PointEvaluator([id(u)], sol)), rasterpoints = 20, clear = false) - displace_mesh!(xgrid, sol[u]) - gridplot!(plt[3, 1], xgrid, linewidth = 1, title = "displaced mesh") - println(stdout, unicode_gridplot(xgrid)) - - return strain_nodevals, plt -end - -## grid -function bimetal_strip2D(; scale = [1, 1], n = 2, anisotropy_factor::Int = Int(ceil(scale[2] / (2 * scale[1])))) - X = linspace(-scale[2] / 2, 0, (n + 1) * anisotropy_factor) - X2 = linspace(0, scale[2] / 2, (n + 1) * anisotropy_factor) - append!(X, X2[2:end]) - Y = linspace(0, scale[1], 2 * n + 1) - xgrid = simplexgrid(X, Y) - cellmask!(xgrid, [-scale[2] / 2, 0.0], [scale[2] / 2, scale[1] / 2], 1) - cellmask!(xgrid, [-scale[2] / 2, scale[1] / 2], [scale[2] / 2, scale[1]], 2) - bfacemask!(xgrid, [-scale[2] / 2, 0.0], [-scale[2] / 2, scale[1] / 2], 1) - bfacemask!(xgrid, [-scale[2] / 2, scale[1] / 2], [-scale[2] / 2, scale[1]], 1) - bfacemask!(xgrid, [-scale[2] / 2, 0.0], [scale[2] / 2, 0.0], 2) - bfacemask!(xgrid, [-scale[2] / 2, scale[1]], [scale[2] / 2, scale[1]], 2) - bfacemask!(xgrid, [scale[2] / 2, 0.0], [scale[2] / 2, scale[1]], 3) - return xgrid -end - -generateplots = ExtendableFEM.default_generateplots(Example230_NonlinearElasticity, "example230.png") #hide -function runtests() #hide - strain, plt = main(periodic = false) #hide - @test maximum(strain) ≈ 0.17289633483008537 #hide - - strain, plt = main(periodic = true) #hide - @test maximum(strain) ≈ 0.17307543557837013 #hide - - return nothing #hide -end #hide -end diff --git a/examples/Example235_StokesIteratedPenalty.jl b/examples/Example235_StokesIteratedPenalty.jl deleted file mode 100644 index 860b0488..00000000 --- a/examples/Example235_StokesIteratedPenalty.jl +++ /dev/null @@ -1,125 +0,0 @@ -#= - -# 235 : Stokes iterated penalty method -([source code](@__SOURCE_URL__)) - -This example computes a velocity ``\mathbf{u}`` and pressure ``\mathbf{p}`` of the incompressible Stokes problem -```math -\begin{aligned} -- \mu \Delta \mathbf{u} + \nabla p & = \mathbf{0}\\ -\mathrm{div}(\mathbf{u}) & = 0 -\end{aligned} -``` -with some μ parameter ``\mu``. - -Here we solve the simple Hagen-Poiseuille flow on the two-dimensional unit square domain with the iterated penalty method -suggested in the reference below adapted to the Bernardi--Raugel finite element method. -Given intermediate solutions ``\mathbf{u}_h`` and ``p_h`` the next approximations are computed by the two equations - -```math -\begin{aligned} -(\nabla \mathbf{u}_h^{next}, \nabla \mathbf{v}_h) + \lambda (\mathrm{div}_h(\mathbf{u}^{next}_h) ,\mathrm{div}_h(\mathbf{v}_h)) & = (\mathbf{f},\mathbf{v}_h) + (p_h,\mathrm{div}(\mathbf{v}_h)) -&& \text{for all } \mathbf{v}_h \in \mathbf{V}_h\\ -(p^{next}_h,q_h) & = (p_h,q_h) - \lambda (\mathrm{div}(\mathbf{u}_h^{next}),q_h) && \text{for all } q_h \in Q_h -\end{aligned} -``` - -This is done consecutively until the residual of both equations is small enough. -The discrete divergence is computed via a RT0 reconstruction operator that preserves the discrete divergence. -(another way would be to compute ``B M^{-1} B^T`` where ``M`` is the mass matrix of the pressure and ``B`` is the matrix for the div-pressure block). - -!!! reference - - "An iterative penalty method for the finite element solution of the stationary Navier-Stokes equations",\ - R. Codina,\ - Computer Methods in Applied Mechanics and Engineering Volume 110, Issues 3–4 (1993),\ - [>Journal-Link<](https://doi.org/10.1016/0045-7825(93)90163-R) - - -The computed solution for the default parameters looks like this: - -![](example235.png) -=# - -module Example235_StokesIteratedPenalty - -using ExtendableFEM -using ExtendableGrids -using Test #hide - -## data for Hagen-Poiseuille flow -function p!(result, qpinfo) - x = qpinfo.x - μ = qpinfo.params[1] - result[1] = μ * (-2 * x[1] + 1.0) - return nothing -end -function u!(result, qpinfo) - x = qpinfo.x - result[1] = x[2] * (1.0 - x[2]) - result[2] = 0.0 - return nothing -end -## kernel for div projection -function div_projection!(result, input, qpinfo) - result[1] = input[1] - qpinfo.params[1] * input[2] - return nothing -end - -## everything is wrapped in a main function -function main(; Plotter = nothing, λ = 1.0e4, μ = 1.0, nrefs = 5, kwargs...) - - ## initial grid - xgrid = uniform_refine(grid_unitsquare(Triangle2D), nrefs) - - ## Bernardi--Raugel element with reconstruction operator - FETypes = (H1BR{2}, L2P0{1}) - PenaltyDivergence = Reconstruct{HDIVRT0{2}, Divergence} - - ## generate two problems - ## one for velocity, one for pressure - u = Unknown("u"; name = "velocity") - p = Unknown("p"; name = "pressure") - PDu = ProblemDescription("Stokes IPM - velocity update") - assign_unknown!(PDu, u) - assign_operator!(PDu, BilinearOperator([grad(u)]; factor = μ, store = true, kwargs...)) - assign_operator!(PDu, BilinearOperator([apply(u, PenaltyDivergence)]; store = true, factor = λ, kwargs...)) - assign_operator!(PDu, LinearOperator([div(u)], [id(p)]; factor = 1, kwargs...)) - assign_operator!(PDu, InterpolateBoundaryData(u, u!; regions = 1:4, params = [μ], bonus_quadorder = 4, kwargs...)) - - PDp = ProblemDescription("Stokes IPM - pressure update") - assign_unknown!(PDp, p) - assign_operator!(PDp, BilinearOperator([id(p)]; store = true, kwargs...)) - assign_operator!(PDp, LinearOperator(div_projection!, [id(p)], [id(p), div(u)]; params = [λ], factor = 1, kwargs...)) - - ## show and solve problem - FES = [FESpace{FETypes[1]}(xgrid), FESpace{FETypes[2]}(xgrid)] - sol = FEVector([FES[1], FES[2]]; tags = [u, p]) - SC1 = SolverConfiguration(PDu; init = sol, maxiterations = 1, target_residual = 1.0e-8, constant_matrix = true, kwargs...) - SC2 = SolverConfiguration(PDp; init = sol, maxiterations = 1, target_residual = 1.0e-8, constant_matrix = true, kwargs...) - sol, nits = iterate_until_stationarity([SC1, SC2]; init = sol, kwargs...) - @info "converged after $nits iterations" - - ## plot - plt = plot([id(u), id(p)], sol; Plotter = Plotter) - - return sol, plt -end - -generateplots = ExtendableFEM.default_generateplots(Example235_StokesIteratedPenalty, "example235.png") #hide -function exact_error!(result, u, qpinfo) #hide - u!(result, qpinfo) #hide - p!(view(result, 3), qpinfo) #hide - return result .= (result .- u) .^ 2 #hide -end #hide -function runtests(; μ = 1.0) #hide - sol, plt = main(; μ = μ) #hide - ErrorIntegratorExact = ItemIntegrator(exact_error!, [id(1), id(2)]; quadorder = 4, params = [μ]) #hide - error = evaluate(ErrorIntegratorExact, sol) #hide - error_u = sqrt(sum(view(error, 1, :)) + sum(view(error, 2, :))) #hide - error_p = sqrt(sum(view(error, 3, :))) #hide - @test error_u ≈ 3.990987355891888e-5 #hide - @test error_p ≈ 0.010437891104305222 #hide - return nothing #hide -end #hide -end diff --git a/examples/Example240_SVRTEnrichment.jl b/examples/Example240_SVRTEnrichment.jl deleted file mode 100644 index fac1a236..00000000 --- a/examples/Example240_SVRTEnrichment.jl +++ /dev/null @@ -1,465 +0,0 @@ -#= - -# 240 : Stokes ``RT`` enrichment -([source code](@__SOURCE_URL__)) - -This example computes the velocity ``\mathbf{u}`` and pressure ``\mathbf{p}`` of the incompressible Stokes problem -```math -\begin{aligned} -- \mu \Delta \mathbf{u} + \nabla p & = \mathbf{f}\\ -\mathrm{div}(\mathbf{u}) & = 0 -\end{aligned} -``` -with exterior force ``\mathbf{f}`` and some parameter ``\mu`` and inhomogeneous Dirichlet boundary data. - -The problem will be solved by a ``(P_k \oplus RTenrichment) \times P_{k-1}`` scheme, which can be seen as an inf-sup stabilized Scott-Vogelius variant -that works with general meshes, see references below. -Therein, the velocity space employs continuous ``P_{k}`` functions plus certain (only H(div)-conforming) Raviart-Thomas functions and a discontinuous ``P_{k-1}`` pressure space -leading to an exactly divergence-free discrete velocity. In a reduction step (that can be triggered with the reduce switch) all higher order pressure dofs and the -enrichment dofs can be eliminated from the system. - -!!! reference - - "A low-order divergence-free H(div)-conforming finite element method for Stokes flows",\ - X. Li, H. Rui,\ - IMA Journal of Numerical Analysis (2021),\ - [>Journal-Link<](https://doi.org/10.1093/imanum/drab080) - [>Preprint-Link<](https://arxiv.org/abs/2012.01689) - - "Inf-sup stabilized Scott--Vogelius pairs on general simplicial grids by Raviart--Thomas enrichment",\ - V. John, X. Li, C. Merdon, H. Rui,\ - [>Preprint-Link<](https://arxiv.org/abs/2206.01242) - - -The computed solution for the default parameters looks like this: - -![](example240.png) - -=# - -module Example240_SVRTEnrichment - -using ExtendableFEM -using GridVisualize -using ExtendableGrids -using ExtendableSparse -using Triangulate -using SimplexGridFactory -using Symbolics -using Test #hide - -## exact data for problem generated by Symbolics -function prepare_data(; μ = 1) - - @variables x y - - ## stream function ξ - ξ = -sin(2 * pi * x) * cos(2 * pi * y) - - ## velocity u = curl ξ - ∇ξ = Symbolics.gradient(ξ, [x, y]) - u = [-∇ξ[2], ∇ξ[1]] - - ## pressure - p = (cos(4 * pi * x) - cos(4 * pi * y)) / 4 - - ## gradient of velocity - ∇u = Symbolics.jacobian(u, [x, y]) - ∇u_reshaped = [∇u[1, 1], ∇u[1, 2], ∇u[2, 1], ∇u[2, 2]] - - ## Laplacian - Δu = [ - (Symbolics.gradient(∇u[1, 1], [x]) + Symbolics.gradient(∇u[1, 2], [y]))[1], - (Symbolics.gradient(∇u[2, 1], [x]) + Symbolics.gradient(∇u[2, 2], [y]))[1], - ] - - ## right-hand side - f = -μ * Δu + Symbolics.gradient(p, [x, y]) - - ## build functions - p_eval = build_function(p, x, y, expression = Val{false}) - u_eval = build_function(u, x, y, expression = Val{false}) - ∇u_eval = build_function(∇u_reshaped, x, y, expression = Val{false}) - f_eval = build_function(f, x, y, expression = Val{false}) - - return f_eval[2], u_eval[2], ∇u_eval[2], p_eval -end - -## grid generator function -function get_grid2D(nref; uniform = false, barycentric = false) - if uniform || barycentric - gen_ref = 0 - else - gen_ref = nref - end - grid = simplexgrid( - Triangulate; - points = [0 0; 0 1; 1 1; 1 0]', - bfaces = [1 2; 2 3; 3 4; 4 1]', - bfaceregions = [1, 2, 3, 4], - regionpoints = [0.5 0.5;]', - regionnumbers = [1], - regionvolumes = [4.0^(-gen_ref - 1)] - ) - if uniform - grid = uniform_refine(grid, nref) - end - if barycentric - grid = barycentric_refine(grid) - end - return grid -end - -## kernel for Stokes operator -function kernel_stokes_standard!(result, u_ops, qpinfo) - ∇u, p = view(u_ops, 1:4), view(u_ops, 5) - μ = qpinfo.params[1] - result[1] = μ * ∇u[1] - p[1] - result[2] = μ * ∇u[2] - result[3] = μ * ∇u[3] - result[4] = μ * ∇u[4] - p[1] - result[5] = -(∇u[1] + ∇u[4]) - return nothing -end - -function main(; nrefs = 5, μ = 1, α = 1, order = 2, Plotter = nothing, enrich = true, reduce = true, time = 0.5, bonus_quadorder = 5, kwargs...) - - ## prepare problem data - f_eval, u_eval, ∇u_eval, p_eval = prepare_data(; μ = μ) - rhs!(result, qpinfo) = (f_eval(result, qpinfo.x[1], qpinfo.x[2])) - exact_p!(result, qpinfo) = (result[1] = p_eval(qpinfo.x[1], qpinfo.x[2])) - exact_u!(result, qpinfo) = (u_eval(result, qpinfo.x[1], qpinfo.x[2])) - exact_∇u!(result, qpinfo) = (∇u_eval(result, qpinfo.x[1], qpinfo.x[2])) - - ## prepare unknowns - u = Unknown("u"; name = "velocity", dim = 2) - pfull = Unknown("p"; name = "pressure (full)", dim = 1) - pE = Unknown("p⟂"; name = "pressure (enriched)", dim = 1) - p0 = Unknown("p0"; name = "pressure (reduced)", dim = 1) # only used if enrich && reduced - uR = Unknown("uR"; name = "velocity enrichment", dim = 2) # only used if enrich == true - - ## prepare plots - plt = GridVisualizer(; Plotter = Plotter, layout = (2, 2), clear = true, size = (1000, 1000)) - - ## prepare error calculations - function exact_error!(result, u, qpinfo) - exact_u!(view(result, 1:2), qpinfo) - exact_∇u!(view(result, 3:6), qpinfo) - result .-= u - return result .= result .^ 2 - end - function exact_error_p!(result, p, qpinfo) - exact_p!(view(result, 1), qpinfo) - result .-= p - return result .= result .^ 2 - end - ErrorIntegratorExact = ItemIntegrator(exact_error!, [id(u), grad(u)]; quadorder = 2 * (order + 1), kwargs...) - ErrorIntegratorPressure = ItemIntegrator(exact_error_p!, [order == 1 ? id(p0) : id(pfull)]; quadorder = 2 * (order + 1), kwargs...) - L2NormIntegratorE = L2NormIntegrator([id(uR)]; quadorder = 2 * order) - function kernel_div!(result, u, qpinfo) - return result .= sum(u) .^ 2 - end - DivNormIntegrator = ItemIntegrator(kernel_div!, enrich ? [div(u), div(uR)] : [div(u)]; quadorder = 2 * order) - NDofs = zeros(Int, nrefs) - Results = zeros(Float64, nrefs, 5) - - for lvl in 1:nrefs - - ## grid - xgrid = get_grid2D(lvl) - - ## define and assign unknowns - PD = ProblemDescription("Stokes problem") - assign_unknown!(PD, u) - p = reduce * enrich ? p0 : pfull - assign_unknown!(PD, p) - - ################ - ### FESPACES ### - ################ - if order == 1 - FES_enrich = FESpace{HDIVRT0{2}}(xgrid) - else - FES_enrich = FESpace{HDIVRTkENRICH{2, order - 1, reduce}}(xgrid) - end - FES = Dict( - u => FESpace{H1Pk{2, 2, order}}(xgrid), - pfull => FESpace{order == 1 ? L2P0{1} : H1Pk{1, 2, order - 1}}(xgrid; broken = true), - p0 => FESpace{L2P0{1}}(xgrid; broken = true), - uR => enrich ? FES_enrich : nothing - ) - - ###################### - ### STANDARD TERMS ### - ###################### - assign_operator!(PD, LinearOperator(rhs!, [id(u)]; bonus_quadorder = bonus_quadorder, kwargs...)) - assign_operator!(PD, BilinearOperator(kernel_stokes_standard!, [grad(u), id(p)]; params = [μ], kwargs...)) - assign_operator!(PD, InterpolateBoundaryData(u, exact_u!; regions = 1:4, bonus_quadorder = bonus_quadorder)) - assign_operator!(PD, FixDofs(p; dofs = [1], vals = [0])) - - ################## - ### ENRICHMENT ### - ################## - if enrich - if reduce - if order == 1 - @info "... preparing condensation of RT0 dofs" - AR = FEMatrix(FES_enrich) - BR = FEMatrix(FES[p], FES_enrich) - bR = FEVector(FES_enrich) - assemble!(AR, BilinearOperator([div(1)]; lump = true, factor = α * μ, kwargs...)) - for bface in xgrid[BFaceFaces] - AR.entries[bface, bface] = 1.0e60 - end - assemble!(BR, BilinearOperator([id(1)], [div(1)]; factor = -1, kwargs...)) - assemble!(bR, LinearOperator(rhs!, [id(1)]; bonus_quadorder = 5, kwargs...); time = time) - ## invert AR (diagonal matrix) - AR.entries.cscmatrix.nzval .= 1 ./ AR.entries.cscmatrix.nzval - C = -BR.entries.cscmatrix * AR.entries.cscmatrix * BR.entries.cscmatrix' - c = -BR.entries.cscmatrix * AR.entries.cscmatrix * bR.entries - assign_operator!(PD, BilinearOperator(C, [p], [p]; kwargs...)) - assign_operator!(PD, LinearOperator(c, [p]; kwargs...)) - else - @info "... preparing removal of enrichment dofs" - BR = FEMatrix(FES[p], FES_enrich) - A1R = FEMatrix(FES_enrich, FES[u]) - bR = FEVector(FES_enrich) - assemble!(BR, BilinearOperator([id(1)], [div(1)]; factor = -1, kwargs...)) - assemble!(bR, LinearOperator(rhs!, [id(1)]; bonus_quadorder = 5, kwargs...); time = time) - assemble!(A1R, BilinearOperator([id(1)], [Δ(1)]; factor = -μ, kwargs...)) - F, DD_RR = div_projector(FES[u], FES_enrich) - C = F.entries.cscmatrix * A1R.entries.cscmatrix - assign_operator!(PD, BilinearOperator(C, [u], [u]; factor = 1, transposed_copy = -1, kwargs...)) - assign_operator!(PD, LinearOperator(F.entries.cscmatrix * bR.entries, [u]; kwargs...)) - end - else - assign_unknown!(PD, uR) - assign_operator!(PD, LinearOperator(rhs!, [id(uR)]; bonus_quadorder = 5, kwargs...)) - assign_operator!(PD, BilinearOperator([id(p)], [div(uR)]; transposed_copy = 1, factor = -1, kwargs...)) - if order == 1 - assign_operator!(PD, BilinearOperator([div(uR)]; lump = true, factor = μ, kwargs...)) - assign_operator!(PD, HomogeneousBoundaryData(uR; regions = 1:4)) - else - assign_operator!(PD, BilinearOperator([Δ(u)], [id(uR)]; factor = μ, transposed_copy = -1, kwargs...)) - end - end - end - - ############# - ### SOLVE ### - ############# - sol = solve(PD, FES; time = time, kwargs...) - NDofs[lvl] = length(sol.entries) - - ## move integral mean of pressure - pintegrate = ItemIntegrator([id(p)]) - pmean = sum(evaluate(pintegrate, sol)) / sum(xgrid[CellVolumes]) - view(sol[p]) .-= pmean - - ###################### - ### POSTPROCESSING ### - ###################### - if enrich && reduce - append!(sol, FES_enrich; tag = uR) - if order == 1 - ## compute enrichment part of velocity - view(sol[uR]) .= AR.entries.cscmatrix * (bR.entries - BR.entries.cscmatrix' * view(sol[p])) - else - ## compute enrichment part of velocity - view(sol[uR]) .= F.entries.cscmatrix' * view(sol[u]) - end - - ## compute higher order pressure dofs - if reduce && order > 1 - ## add blocks for higher order pressures to sol vector - VR = FES_enrich - append!(sol, VR; tag = pE) - append!(sol, FES[pfull]; tag = pfull) - sol_pE = view(sol[pE]) - sol_pfull = view(sol[pfull]) - sol_p0 = view(sol[p0]) - - res = FEVector(VR) - addblock_matmul!(res[1], A1R[1, 1], sol[u]) - celldofs_VR::VariableTargetAdjacency{Int32} = VR[CellDofs] - ndofs_VR = max_num_targets_per_source(celldofs_VR) - Ap = zeros(Float64, ndofs_VR, ndofs_VR) - bp = zeros(Float64, ndofs_VR) - xp = zeros(Float64, ndofs_VR) - for cell in 1:num_cells(xgrid) - ## solve local pressure reconstruction - ## (p_h, div VR) = - (f,VR) + a_h(u_h,VR) - for dof_j in 1:ndofs_VR - dof = celldofs_VR[dof_j, cell] - bp[dof_j] = -bR.entries[dof] + res.entries[dof] - for dof_k in 1:ndofs_VR - dof2 = celldofs_VR[dof_k, cell] - Ap[dof_j, dof_k] = DD_RR.entries[dof, dof2] - end - end - - ## solve for coefficients of div(RT1bubbles) - xp = Ap \ bp - - ## save in block id_pk - for dof_j in 1:ndofs_VR - dof = celldofs_VR[dof_j, cell] - sol_pE[dof] = xp[dof_j] - end - end - - ## interpolate into Pk basis (= same pressure basis as in full scheme) - PF = FES[pfull] - append!(sol, PF; tag = pfull) - celldofs_PF::SerialVariableTargetAdjacency{Int32} = PF[CellDofs] - ndofs_PF::Int = max_num_targets_per_source(celldofs_PF) - - ## compute local mass matrix of full pressure space - MAMA = FEMatrix(PF) - assemble!(MAMA, BilinearOperator([id(1)])) - MAMAE::ExtendableSparseMatrix{Float64, Int64} = MAMA.entries - - ## full div-pressure matrix - PFxVR = FEMatrix(PF, VR) - assemble!(PFxVR, BilinearOperator([id(1)], [div(1)])) - PFxVRE::ExtendableSparseMatrix{Float64, Int64} = PFxVR.entries - bp = zeros(Float64, ndofs_PF) - xp = zeros(Float64, ndofs_PF) - locMAMA = zeros(Float64, ndofs_PF, ndofs_PF) - for cell in 1:num_cells(xgrid) - ## solve local pressure reconstruction - fill!(bp, 0) - for dof_k in 1:ndofs_PF - dof2 = celldofs_PF[dof_k, cell] - for dof_j in 1:ndofs_VR - dof = celldofs_VR[dof_j, cell] - bp[dof_k] += PFxVRE[dof2, dof] * sol_pE[dof] - end - for dof_j in 1:ndofs_PF - dof = celldofs_PF[dof_j, cell] - locMAMA[dof_k, dof_j] = MAMAE[dof2, dof] - end - end - - ## solve for coefficients of div(RT1bubbles) - xp = locMAMA \ bp - for dof_j in 1:ndofs_PF - dof = celldofs_PF[dof_j, cell] - sol_pfull[dof] = sol_p0[cell] + xp[dof_j] - end - end - elseif reduce && order == 1 - pfull = p0 - end - end - - ######################## - ### ERROR EVALUATION ### - ######################## - error = evaluate(ErrorIntegratorExact, sol) - L2errorU = sqrt(sum(view(error, 1, :)) + sum(view(error, 2, :))) - H1errorU = sqrt(sum(view(error, 3, :)) + sum(view(error, 4, :)) + sum(view(error, 5, :)) + sum(view(error, 6, :))) - @info "L2error(u) = $L2errorU" - @info "L2error(∇u) = $H1errorU" - evaluate!(error, ErrorIntegratorPressure, sol) - L2errorP = sqrt(sum(view(error, 1, :))) - @info "L2error(p) = $L2errorP" - Results[lvl, 4] = L2errorP - if enrich - fill!(error, 0) - evaluate!(error, L2NormIntegratorE, sol) - L2normUR = sqrt(sum(view(error, 1, :)) + sum(view(error, 2, :))) - @info "L2norm(uR) = $L2normUR" - end - fill!(error, 0) - evaluate!(error, DivNormIntegrator, sol) - L2normDiv = sqrt(sum(view(error, 1, :))) - @info "L2norm(div(u+uR)) = $L2normDiv" - - Results[lvl, 1] = L2errorU - Results[lvl, 2] = H1errorU - Results[lvl, 3] = L2normUR - Results[lvl, 5] = L2normDiv - - ############# - ### PLOTS ### - ############# - scalarplot!(plt[1, 1], id(u), sol; abs = true) - scalarplot!(plt[1, 2], id(pfull), sol) - if order == 1 && enrich - scalarplot!(plt[2, 2], id(uR), sol) - end - end - plot_convergencehistory!( - plt[2, 1], - NDofs, - Results[:, 1:4]; - add_h_powers = [order, order + 1], - X_to_h = X -> 8 * X .^ (-1 / 2), - legend = :best, - ylabels = ["|| u - u_h ||", "|| ∇(u - u_h) ||", "|| uR ||", "|| p - p_h ||", "|| div(u + uR) ||"], - ) - - print_convergencehistory(NDofs, Results; X_to_h = X -> X .^ (-1 / 2), ylabels = ["|| u - u_h ||", "|| ∇(u - u_h) ||", "|| uR ||", "|| p - p_h ||", "|| div(u + uR) ||"], xlabel = "ndof") - - return Results, plt -end - -function div_projector(V1, VR) - - ## setup interpolation matrix - celldofs_V1 = V1[CellDofs] - celldofs_VR = VR[CellDofs] - ndofs_V1 = max_num_targets_per_source(celldofs_V1) - ndofs_VR = max_num_targets_per_source(celldofs_VR) - - DD_RR = FEMatrix(VR) - assemble!(DD_RR, BilinearOperator([div(1)])) - DD_RRE = DD_RR.entries - DD_1R = FEMatrix(V1, VR) - assemble!(DD_1R, BilinearOperator([div(1)])) - DD_1RE = DD_1R.entries - Ap = zeros(Float64, ndofs_VR, ndofs_VR) - bp = zeros(Float64, ndofs_VR) - xp = zeros(Float64, ndofs_VR) - ncells = num_sources(celldofs_V1) - F = FEMatrix(V1, VR) - FE = F.entries - for cell in 1:ncells - - ## solve local pressure reconstruction for RTk part - for dof_j in 1:ndofs_VR - dof = celldofs_VR[dof_j, cell] - for dof_k in 1:ndofs_VR - dof2 = celldofs_VR[dof_k, cell] - Ap[dof_j, dof_k] = DD_RRE[dof, dof2] - end - end - - for dof_j in 1:ndofs_V1 - dof = celldofs_V1[dof_j, cell] - for dof_k in 1:ndofs_VR - dof2 = celldofs_VR[dof_k, cell] - bp[dof_k] = -DD_1RE[dof, dof2] - end - - xp = Ap \ bp - - for dof_k in 1:ndofs_VR - dof2 = celldofs_VR[dof_k, cell] - FE[dof, dof2] = xp[dof_k] - end - end - end - flush!(FE) - return F, DD_RR -end - -generateplots = ExtendableFEM.default_generateplots(Example240_SVRTEnrichment, "example240.png") #hide -function runtests() #hide - Results, plt = main(; nrefs = 2) #hide - @test Results[end, 1] ≈ 0.09600693353585522 #hide - @test Results[end, 5] < 1.0e-9 #hide - return nothing #hide -end #hide -end # module diff --git a/examples/Example245_NSEFlowAroundCylinder.jl b/examples/Example245_NSEFlowAroundCylinder.jl index fd4534e3..0eddf1a3 100644 --- a/examples/Example245_NSEFlowAroundCylinder.jl +++ b/examples/Example245_NSEFlowAroundCylinder.jl @@ -30,6 +30,7 @@ using SimplexGridFactory using ExtendableGrids using GridVisualize using LinearAlgebra +using UnicodePlots; import Term using Test #hide ## inlet data for Karman vortex street example @@ -86,7 +87,7 @@ end ## everything is wrapped in a main function -function main(; Plotter = nothing, μ = 1.0e-3, maxvol = 1.0e-3, reconstruct = true, parallel = false, npart = 8, kwargs...) +function main(; Plotter = UnicodePlots, μ = 1.0e-3, maxvol = 1.0e-3, reconstruct = true, parallel = false, npart = 8, kwargs...) ## load grid (see function below) xgrid = make_grid(W, H; n = Int(ceil(sqrt(1 / maxvol))), maxvol = maxvol) @@ -123,11 +124,13 @@ function main(; Plotter = nothing, μ = 1.0e-3, maxvol = 1.0e-3, reconstruct = t ## plots via GridVisualize plt = GridVisualizer(; Plotter = Plotter, layout = (4, 1), clear = true, size = (800, 1200)) - gridplot!(plt[1, 1], xgrid, cellcoloring = :partitions, linewidth = 1) - gridplot!(plt[2, 1], xgrid, cellcoloring = :partitions, linewidth = 1, xlimits = [0, 0.3], ylimits = [0.1, 0.3]) - scalarplot!(plt[3, 1], xgrid, nodevalues(sol[u]; abs = true)[1, :]) - vectorplot!(plt[3, 1], xgrid, eval_func_bary(PointEvaluator([id(u)], sol)), rasterpoints = 20, clear = false) - scalarplot!(plt[4, 1], xgrid, view(nodevalues(sol[p]), 1, :), levels = 11, title = "p_h") + # gridplot!(plt[1, 1], xgrid) + gridplot!(plt[1, 1], xgrid, title = "grid") + # gridplot!(plt[2, 1], xgrid, cellcoloring = :partitions, linewidth = 1, xlimits = [0, 0.3], ylimits = [0.1, 0.3], title = "grid around obstacle") + # scalarplot!(plt[3, 1], xgrid, nodevalues(sol[u]; abs = true)[1, :]) + # vectorplot!(plt[3, 1], xgrid, eval_func_bary(PointEvaluator([id(u)], sol)), rasterpoints = 20, clear = false, title = "|u_h| + quiver") + # scalarplot!(plt[4, 1], xgrid, view(nodevalues(sol[p]), 1, :), levels = 11, title = "p_h") + reveal(plt) return [draglift[1], draglift[2], pdiff[1]], plt end diff --git a/examples/Example250_NSELidDrivenCavity.jl b/examples/Example250_NSELidDrivenCavity.jl deleted file mode 100644 index 7fb0c061..00000000 --- a/examples/Example250_NSELidDrivenCavity.jl +++ /dev/null @@ -1,159 +0,0 @@ -#= - -# 250 : Navier--Stokes Lid-driven cavity -([source code](@__SOURCE_URL__)) - -This example computes the velocity ``\mathbf{u}`` and pressure ``\mathbf{p}`` of the incompressible Navier--Stokes problem -```math -\begin{aligned} -- \mu \Delta \mathbf{u} + \left(\mathbf{u} \cdot \nabla\right) \mathbf{u}+ \nabla p & = \mathbf{f}\\ -\mathrm{div}(\mathbf{u}) & = 0 -\end{aligned} -``` -in a lid driven cavity example over a cone and plots the solution and the formed eddies. - -The computed solution for the default parameters looks like this: - -![](example250.png) - -=# - -module Example250_NSELidDrivenCavity - -using ExtendableFEM -using Triangulate -using ExtendableGrids -using SimplexGridFactory -using LinearAlgebra -using GridVisualize -using Test #hide - - -function create_cone(h) - builder = SimplexGridBuilder(; Generator = Triangulate) - - ## points - p1 = point!(builder, -1, 0) - p2 = point!(builder, 1, 0) - p3 = point!(builder, 0, -2) - - ## top face - facetregion!(builder, 1) - facet!(builder, p1, p2) - - ## other faces - facetregion!(builder, 2) - facet!(builder, p2, p3) - facet!(builder, p3, p1) - - cellregion!(builder, 1) - maxvolume!(builder, h) - regionpoint!(builder, 0, -0.5) - - return simplexgrid(builder) -end - -const 𝕀 = [1 0; 0 1] - -function NSE_kernel!(result, u_ops, qpinfo) - - u = tensor_view(u_ops, 1, TDVector(2)) - v = tensor_view(result, 1, TDVector(2)) - ∇u = tensor_view(u_ops, 3, TDMatrix(2)) - ∇v = tensor_view(result, 3, TDMatrix(2)) - p = tensor_view(u_ops, 7, TDScalar()) - q = tensor_view(result, 7, TDScalar()) - μ = qpinfo.params[1] - - tmul!(v, ∇u, u) - ∇v .= μ .* ∇u .- p[1] .* 𝕀 - q[1] = -dot(∇u, 𝕀) - - return nothing -end - - -## boundary function -function u_boundary!(result, qpinfo) - result[1] = 1 - result[2] = 0 - return nothing -end - - -function main(; - μ_final = 0.0005, # flow parameter - order = 2, # FE order of the flow field (pressure order is order-1) - h = 1.0e-3, # grid cell volume - nrefs = 1, # additional grid refinements - Plotter = nothing, - kwargs... - ) - - ## prepare parameter field - extra_params = Array{Float64, 1}([max(μ_final, 0.05)]) - - ## problem description - PD = ProblemDescription() - u = Unknown("u"; name = "velocity") - p = Unknown("p"; name = "pressure") - - assign_unknown!(PD, u) - assign_unknown!(PD, p) - assign_operator!(PD, NonlinearOperator(NSE_kernel!, [id(u), grad(u), id(p)]; params = extra_params, kwargs...)) - - ## boundary data - assign_operator!(PD, InterpolateBoundaryData(u, u_boundary!; regions = 1)) - assign_operator!(PD, HomogeneousBoundaryData(u; regions = [2])) - - ## grid - xgrid = uniform_refine(create_cone(h), nrefs) - - ## prepare FESpace - FES = [ - FESpace{H1Pk{2, 2, order}}(xgrid), - FESpace{H1Pk{1, 2, order - 1}}(xgrid), - ] - - ## solve by μ embedding - step = 0 - sol = nothing - SC = nothing - PE = PointEvaluator([id(1)]) - while (true) - step += 1 - @info "Step $step : solving for μ=$(extra_params[1])" - sol, SC = ExtendableFEM.solve( - PD, - FES, - SC; - return_config = true, - target_residual = 1.0e-10, - maxiterations = 20, - kwargs... - ) - - if step == 1 - initialize!(PE, sol) - end - - if extra_params[1] <= μ_final - break - else - extra_params[1] = max(μ_final, extra_params[1] / 2) - end - end - - plt = plot([id(u), streamlines(u)], sol; Plotter = Plotter, rasterpoints = 30, width = 800, height = 400, title_add = " (μ = $(μ_final))") - - return sol, plt -end - -generateplots = ExtendableFEM.default_generateplots(Example250_NSELidDrivenCavity, "example250.png") #hide -function runtests() #hide - sol, plt = main(; μ_final = 0.005) #hide - sum(view(sol[1])) ≈ 237.24628017878518 #hide - return nothing #hide -end #hide - -end # module diff --git a/examples/Example252_NSEPlanarLatticeFlow.jl b/examples/Example252_NSEPlanarLatticeFlow.jl deleted file mode 100644 index e9946da8..00000000 --- a/examples/Example252_NSEPlanarLatticeFlow.jl +++ /dev/null @@ -1,152 +0,0 @@ -#= - -# 252 : Navier--Stokes Planar Lattice Flow -([source code](@__SOURCE_URL__)) - -This example computes an approximation to the planar lattice flow test problem of the Navier-Stokes equations -```math -\begin{aligned} -- \nu \Delta \mathbf{u} + (\mathbf{u} \cdot \nabla) \mathbf{u} + \nabla p & = \mathbf{f}\\ -\mathrm{div}(\mathbf{u}) & = 0 -\end{aligned} -``` -with an exterior force ``\mathbf{f}`` and some viscosity parameter ``\nu`` and Dirichlet boundary data for ``\mathbf{u}``. - -Here the exact data for the planar lattice flow -```math -\begin{aligned} -\mathbf{u}(x,y,t) & := \exp(-8 \pi^2 \nu t) \begin{pmatrix} \sin(2 \pi x) \sin(2 \pi y) \\ \cos(2 \pi x) \cos(2 \pi y) \end{pmatrix}\\ -p(x,y,t) & := \exp(-8 \pi^2 \nu t) ( \cos(4 \pi x) - \cos(4 \pi y)) / 4 -\end{aligned} -``` -is prescribed at fixed time ``t = 0`` with ``\mathbf{f} = - \nu \Delta \mathbf{u}``. - -In this example the Navier-Stokes equations are solved with a pressure-robust variant of the Bernardi--Raugel finite element method -and the nonlinear convection term (that involves reconstruction operators) is automatically differentiated for a Newton iteration. - -The computed solution for the default parameters looks like this: - -![](example252.png) - -=# - -module Example252_NSEPlanarLatticeFlow - -using ExtendableFEM -using ExtendableGrids -using LinearAlgebra -using Test #hide - -## exact velocity (and Dirichlet data) -function u!(result, qpinfo) - x = qpinfo.x - result[1] = sin(2 * pi * x[1]) * sin(2 * pi * x[2]) - result[2] = cos(2 * pi * x[1]) * cos(2 * pi * x[2]) - return nothing -end - -## right-hand side f := -μ Δu + (u⋅∇)u + ∇p -function f!(μ) - α = [0, 0] - return function closure(result, qpinfo) - x = qpinfo.x - result[1] = (μ * 8 * pi^2 + α[1]) * sin(2 * pi * x[1]) * sin(2 * pi * x[2]) - result[2] = (μ * 8 * pi^2 + α[2]) * cos(2 * pi * x[1]) * cos(2 * pi * x[2]) - return nothing - end -end - -## exact pressure -function p!(result, qpinfo) - x = qpinfo.x - result[1] = (cos(4 * pi * x[1]) - cos(4 * pi * x[2])) / 4 - return nothing -end - -function kernel_nonlinear!(result, u_ops, qpinfo) - u, ∇u, p = view(u_ops, 1:2), view(u_ops, 3:6), view(u_ops, 7) - μ = qpinfo.params[1] - result[1] = dot(u, view(∇u, 1:2)) - result[2] = dot(u, view(∇u, 3:4)) - result[3] = μ * ∇u[1] - p[1] - result[4] = μ * ∇u[2] - result[5] = μ * ∇u[3] - result[6] = μ * ∇u[4] - p[1] - result[7] = -(∇u[1] + ∇u[4]) - return nothing -end - -function exact_error!(result, u, qpinfo) - u!(result, qpinfo) - p!(view(result, 3), qpinfo) - result .-= u - return result .= result .^ 2 -end - -function main(; - μ = 0.001, - nrefs = 5, - reconstruct = true, - Plotter = nothing, - use_LM_restrictions = true, - kwargs... - ) - - ## problem description - PD = ProblemDescription() - u = Unknown("u"; name = "velocity") - p = Unknown("p"; name = "pressure") - id_u = reconstruct ? apply(u, Reconstruct{HDIVBDM1{2}, Identity}) : id(u) - - assign_unknown!(PD, u) - assign_unknown!(PD, p) - assign_operator!(PD, NonlinearOperator(kernel_nonlinear!, [id_u, grad(u), id(p)]; params = [μ], kwargs...)) - assign_operator!(PD, LinearOperator(f!(μ), [id_u]; kwargs...)) - assign_operator!(PD, InterpolateBoundaryData(u, u!; regions = 1:4)) - if use_LM_restrictions - assign_restriction!(PD, ZeroMeanValueRestriction(p)) - else - assign_operator!(PD, FixDofs(p; dofs = [1])) - end - - ## grid - xgrid = uniform_refine(grid_unitsquare(Triangle2D), nrefs) - - ## prepare FESpace - FES = [FESpace{H1BR{2}}(xgrid), FESpace{L2P0{1}}(xgrid)] - - ## solve - sol = solve(PD, FES; kwargs...) - - pintegrate = ItemIntegrator([id(p)]) - pmean = sum(evaluate(pintegrate, sol)) / sum(xgrid[CellVolumes]) - if use_LM_restrictions - @show pmean - else - ## move integral mean of pressure - view(sol[p]) .-= pmean - end - - ## error calculation - ErrorIntegratorExact = ItemIntegrator(exact_error!, [id(u), id(p)]; quadorder = 4, params = [μ], kwargs...) - error = evaluate(ErrorIntegratorExact, sol) - L2errorU = sqrt(sum(view(error, 1, :)) + sum(view(error, 2, :))) - L2errorP = sqrt(sum(view(error, 3, :))) - @info "L2error(u) = $L2errorU" - @info "L2error(p) = $L2errorP" - - ## plot - plt = plot([id(u), id(p)], sol; Plotter = Plotter) - - return L2errorU, plt -end - -generateplots = ExtendableFEM.default_generateplots(Example252_NSEPlanarLatticeFlow, "example252.png") #hide -function runtests() #hide - L2errorU, _ = main(; use_LM_restrictions = false, nrefs = 4) #hide - @test L2errorU ≈ 0.010616923333947861 #hide - L2errorU, _ = main(; use_LM_restrictions = true, nrefs = 4) #hide - @test L2errorU ≈ 0.010616923333947861 #hide - return nothing #hide -end #hide -end # module diff --git a/examples/Example260_AxisymmetricNavierStokesProblem.jl b/examples/Example260_AxisymmetricNavierStokesProblem.jl deleted file mode 100644 index 8f20d4c6..00000000 --- a/examples/Example260_AxisymmetricNavierStokesProblem.jl +++ /dev/null @@ -1,160 +0,0 @@ -#= - -# 260 : Axisymmetric Stokes -([source code](@__SOURCE_URL__)) - -This example solves the 3D stagnation point flow via -the 2.5D axisymmetric formulation of the Navier--Stokes problem that seeks a velocity ``\mathbf{u} = (u_z, u_r)`` -and pressure ``p`` such that -```math -\begin{aligned} -- \mu\left(\partial^2_r + r^{-1} \partial_r + \partial^2_z - r^{-2} \right) u_r -+ (u_r \partial_r + u_z \partial_z) u_r + \partial_r p & = \mathbf{f}_r\\ -- \mu\left(\partial^2_r + r^{-1} \partial_r + \partial^2_z \right) u_z -+ (u_r \partial_r + u_z \partial_z) u_z + \partial_z p & = \mathbf{f}_z\\ -(\partial_r + r^{-1})u_r + \partial_z u_z & = 0 -\end{aligned} -``` -with exterior force ``\mathbf{f}`` and some viscosity parameter ``\mu``. - -The axisymmetric formulation assumes that the velocity in some -3D-domain, that is obtained by rotation of a 2D domain ``\Omega``, -only depends on the distance ``r`` to the rotation axis and the -``z``-coordinate tangential to the x-axis, but not on the angular coordinate -of the cylindric coordinates. -The implementation employs ``r``-dependent bilinear forms and a -Cartesian grid for the 2D ``(z,r)`` domain that is assumed to be rotated -around the ``r=0``-axis. - -This leads to the weak formulation -```math -\begin{aligned} -a(u,v) + b(p,v) & = (f,v) \\ - b(q,u) & = 0 -\end{aligned} -``` -with the bilinear forms -```math -\begin{aligned} -a(u,v) := \int_{\Omega} \left( \nabla u : \nabla v + r^{-2} u_r v_r \right) r dr dz\\ -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. The factor ``2\pi`` from the -integral over the rotation angle drops out on both sides. - -The computed solution for the default parameters looks like this: - -![](example260.png) - -=# - -module Example260_AxisymmetricNavierStokesProblem - -using ExtendableFEM -using ExtendableGrids -using SimplexGridFactory -using Triangulate -using Test #hide - - -function kernel_convection!(result, input, qpinfo) - u, ∇u = view(input, 1:2), view(input, 3:6) - r = qpinfo.x[1] - result[1] = r * (∇u[1] * u[1] + ∇u[2] * u[2]) - result[2] = r * (∇u[3] * u[1] + ∇u[4] * u[2]) - return nothing -end - -function kernel_stokes_axisymmetric!(result, u_ops, qpinfo) - u, ∇u, p = view(u_ops, 1:2), view(u_ops, 3:6), view(u_ops, 7) - r = qpinfo.x[1] - μ = qpinfo.params[1] - ## add Laplacian - result[1] = μ / r * u[1] - p[1] - result[2] = 0 - result[3] = μ * r * ∇u[1] - r * p[1] - result[4] = μ * r * ∇u[2] - result[5] = μ * r * ∇u[3] - result[6] = μ * r * ∇u[4] - r * p[1] - result[7] = -(r * (∇u[1] + ∇u[4]) + u[1]) - return nothing -end - -function u!(result, qpinfo) - x = qpinfo.x - result[1] = x[1] - result[2] = -2 * x[2] - return nothing -end - -function kernel_l2div(result, u_ops, qpinfo) - u, divu = view(u_ops, 1:2), view(u_ops, 3) - result[1] = (qpinfo.x[1] * divu[1] + u[1])^2 - return nothing -end - - -function main(; μ = 0.1, nrefs = 4, nonlinear = false, uniform = false, Plotter = nothing, kwargs...) - - ## problem description - PD = ProblemDescription() - u = Unknown("u"; name = "velocity") - p = Unknown("p"; name = "pressure") - assign_unknown!(PD, u) - assign_unknown!(PD, p) - assign_operator!(PD, BilinearOperator(kernel_stokes_axisymmetric!, [id(u), grad(u), id(p)]; params = [μ], kwargs...)) #; jacobian = kernel_jacobian!)) - if nonlinear - assign_operator!(PD, NonlinearOperator(kernel_convection!, [id(u)], [id(u), grad(u)]; bonus_quadorder = 1, kwargs...)) #; jacobian = kernel_jacobian!)) - end - assign_operator!(PD, InterpolateBoundaryData(u, u!; regions = [3])) - assign_operator!(PD, HomogeneousBoundaryData(u; regions = [4], mask = (1, 0, 1))) - assign_operator!(PD, HomogeneousBoundaryData(u; regions = [1], mask = (0, 1, 1))) - - ## grid - if uniform - xgrid = uniform_refine(grid_unitsquare(Triangle2D), nrefs) - else - xgrid = simplexgrid( - Triangulate; - points = [0 0 ; 5 0 ; 5 1 ; 0 1]', - bfaces = [1 2 ; 2 3 ; 3 4 ; 4 1 ]', - bfaceregions = [1, 2, 3, 4], - regionpoints = [0.5 0.5;]', - regionnumbers = [1], - regionvolumes = [4.0^(-nrefs - 1)] - ) - end - - ## solve - FES = [FESpace{H1P2B{2, 2}}(xgrid), FESpace{L2P1{1}}(xgrid)] - sol = ExtendableFEM.solve(PD, FES; kwargs...) - - ## compute divergence in cylindrical coordinates by volume integrals - DivIntegrator = ItemIntegrator(kernel_l2div, [id(u), div(u)]; quadorder = 4, resultdim = 1) - div_error = sqrt(sum(evaluate(DivIntegrator, sol))) - @info "||div(u_h)|| = $div_error" - - ## compute L2error - function kernel_l2error(result, u_ops, qpinfo) - u!(result, qpinfo) - return result .= (result - u_ops) .^ 2 - end - ErrorIntegratorExact = ItemIntegrator(kernel_l2error, [id(1)]; entities = ON_BFACES, regions = [3], quadorder = 4, kwargs...) - error = evaluate(ErrorIntegratorExact, sol) - L2error = sqrt(sum(view(error, 1, :)) + sum(view(error, 2, :))) - @info "||u - u_h|| = $L2error" - - ## plot - plt = plot([id(u)], sol; Plotter = Plotter) - - return [div_error, L2error], plt -end - -generateplots = ExtendableFEM.default_generateplots(Example260_AxisymmetricNavierStokesProblem, "example260.png") #hide -function runtests() #hide - errors, plt = main(; nrefs = 1) #hide - @test all(errors .<= 1.0e-12) #hide - return nothing #hide -end #hide -end # module diff --git a/examples/Example264_StokesDarcy.jl b/examples/Example264_StokesDarcy.jl deleted file mode 100644 index 3e34fe76..00000000 --- a/examples/Example264_StokesDarcy.jl +++ /dev/null @@ -1,304 +0,0 @@ -#= - -# 264 : Stokes+Darcy -([source code](@__SOURCE_URL__)) - -This example solves the coupled Stokes-Darcy problem in a domain -``\Omega := \Omega_\text{FF} \cup \Omega_\text{PM}`` that has a free flow region -``\Omega_\text{FF}`` and a porous media regions ``\Omega_\text{PM}``. -In the free flow region a Stokes problem is solved that seeks a velocity -``\mathbf{u}_\text{FF}`` and a pressure ``\mathbf{p}_\text{FF}`` -such that -```math -\begin{aligned} -- 2\mu \mathrm{div}(\epsilon(\mathbf{u}_\text{FF}) - p_\text{FF}I) & = \mathbf{f}_\text{FF}\\ -\mathrm{div}(\mathbf{u}_\text{FF}) & = 0. -\end{aligned} -``` -In the porous media region the Darcy problem is solved that seeks a velocity ``\mathbf{u}_\text{PM}`` -and a pressure ``\mathbf{p}_\text{PM}`` such that -```math -\begin{aligned} -\mathbf{u}_\text{PM} + k \nabla p_\text{PM} & = 0\\ -\mathrm{div}(\mathbf{u}_\text{PM}) & = f_\text{PM}. -\end{aligned} -``` -On the interface ``\Gamma := \partial \Omega_\text{FF} \cap \partial \Omega_\text{PM}`` -the two velocities are coupled via several conditions, i.e., the conservation of mass -```math -\mathbf{u}_\text{FF} \cdot \mathbf{n} = \mathbf{u}_\text{PM} \cdot \mathbf{n} \quad \text{on } \Gamma, -``` -the balance of normal forces -```math -p_\text{FF} - 2\mu \epsilon(\mathbf{u}_\text{FF}) \mathbf{n}_\text{FF} \cdot \mathbf{n}_\text{FF} = p_\text{PM}, -``` -and the Beavers-Joseph-Saffman condition -```math -\mathbf{u}_\text{FF} \cdot \mathbf{\tau} -= -\frac{\sqrt{\mu k}}{\mu \alpha} 2 \epsilon(\mathbf{u}_\text{FF}) \mathbf{n}_\text{FF} \cdot \mathbf{\tau}. -``` -The interface condition for the normal fluxes is realized weakly via a Lagrange multiplier ``\lambda`` -that only lives on the interface. - -The weak formulation leads to the problem: -seek ``(\mathbf{u}_\text{FF}, \mathbf{u}_\text{PM}, p_\text{FF}, p_\text{PM}, \lambda)`` -such that, for all ``(\mathbf{v}_\text{FF}, \mathbf{v}_\text{PM}, q_\text{FF}, q_\text{PM}, \chi)``, -```math -\begin{aligned} -a_1(\mathbf{u}_\text{FF},\mathbf{v}_\text{FF}) + a_2(\mathbf{u}_\text{PM}, \mathbf{v}_\text{PM}) -+ b_1(p_\text{FF},\mathbf{v}_\text{FF}) + b_2(p_\text{PM},\mathbf{v}_\text{PM}) -+ b_{\Gamma}(\mathbf{v}_\text{FF} - \mathbf{v}_\text{PM}, \lambda) -& = (\mathbf{f}_\text{FF}, \mathbf{v}_\text{FF})_{L^2(\Omega_\text{FF})}\\ -b_1(q_\text{FF},\mathbf{u}_\text{FF}) + b_2(q_\text{PM},\mathbf{u}_\text{PM}) & = (f_\text{PM}, q_\text{PM})_{L^2(\Omega_\text{PM})}\\ -b_{\Gamma}(\mathbf{u}_\text{FF} - \mathbf{u}_\text{PM}, \chi) & = 0. -\end{aligned} -``` - -The bilinearforms read -```math -\begin{aligned} -a_1(\mathbf{u}_\text{FF},\mathbf{v}_\text{FF}) & := 2\mu (\epsilon(\mathbf{u}_\text{FF}), \epsilon(\mathbf{v}_\text{FF}))_{L^2(\Omega_\text{FF})} -+ \frac{αμ}{\sqrt{μk}} (\mathbf{u}_\text{FF} \cdot \mathbf{\tau},\mathbf{v}_\text{FF} \cdot \mathbf{\tau})_{L^2(\Gamma)}\\ -a_2(\mathbf{u}_\text{PM}, \mathbf{v}_\text{PM}) & := (\mathbf{u}_\text{PM}, \mathbf{v}_\text{PM})_{L^2(\Omega_\text{PM})}\\ -b_1(q_\text{FF},\mathbf{v}_\text{FF}) & := -(\mathrm{div} \mathbf{v}_\text{FF}, q_\text{FF})_{L^2(\Omega_\text{FF})}\\ -b_2(q_\text{PM},\mathbf{v}_\text{PM}) & := -(\mathrm{div} \mathbf{v}_\text{PM}, q_\text{PM})_{L^2(\Omega_\text{PM})}\\ -b_{\Gamma}(\mathbf{v}, \lambda) & := (\mathbf{v} \cdot \mathbf{n}, \lambda)_{L^2(\Gamma)} -\end{aligned} -``` - -Details on the model can be found e.g. in the reference below. - -!!! reference - - ''Coupling Fluid Flow with Porous Media Flow'' - SIAM Journal on Numerical Analysis 2002 40:6, 2195-2218 - [>Link<](https://doi.org/10.1137/S0036142901392766) - -In this example an analytic benchmark problem is solved with the Taylor--Hood FEM in the free flow -domain and the Raviart--Thomas FEM in the porous media domain. - -The computed solution for the default parameters looks like this: - -![](example264.png) - -=# - -module Example264_StokesDarcy - -using ExtendableFEM -using ExtendableFEMBase -using ExtendableGrids -using GridVisualize -using SimplexGridFactory -using Triangulate -using Metis -using Test #hide - -## exact solution and data functions -function u!(result, qpinfo) - x = qpinfo.x - result[1] = -cos(pi * x[1]) * sin(pi * x[2]) - result[2] = sin(pi * x[1]) * cos(pi * x[2]) - return nothing -end -function f_FF!(result, qpinfo) - x = qpinfo.x - result[1] = -(4 * pi^2 + 1) * cos(pi * x[1]) * sin(pi * x[2]) - return result[2] = - sin(pi * x[1]) * cos(pi * x[2]) -end -function p!(result, qpinfo) - x = qpinfo.x - result[1] = -sin(pi * x[1])sin(pi * x[2]) * (1 / pi + 2 * pi) + 4 / pi - return nothing -end -function q!(result, qpinfo) - k = qpinfo.params[1] - x = qpinfo.x - result[1] = -sin(pi * x[1])sin(pi * x[2]) / pi + 4 / pi - return nothing -end -function v!(result, qpinfo) # = -K∇q - k = qpinfo.params[1] - x = qpinfo.x - result[1] = cos(pi * x[1]) * sin(pi * x[2]) - result[2] = sin(pi * x[1]) * cos(pi * x[2]) - return nothing -end -function f_PM!(result, qpinfo) # = div(v) - x = qpinfo.x - return result[1] = -2 * pi * sin(pi * x[1]) * sin(pi * x[2]) -end - -## kernel functions for operators -function kernel_stokes_standard!(result, u_ops, qpinfo) - ∇u, p = view(u_ops, 1:4), view(u_ops, 5) - μ = 2 * qpinfo.params[1] - result[1] = μ * ∇u[1] - p[1] - result[2] = μ * (∇u[2] + ∇u[3]) - result[3] = μ * (∇u[2] + ∇u[3]) - result[4] = μ * ∇u[4] - p[1] - result[5] = -(∇u[1] + ∇u[4]) - return nothing -end - -function coupling_normal!(result, input, qpinfo) - return result[1] = input[2] - input[1] -end - -function exact_error!(f!::Function) - return function closure(result, u, qpinfo) - f!(result, qpinfo) - result .-= u - return result .= result .^ 2 - end -end - -## everything is wrapped in a main function -function main(; - μ = 1, # viscosity - k = 1, # permeability - α = 1, # parameter in interface condition - coupled = true, # solve coupled problem with interface conditions or decoupled problems on subgrids ? - order_u = 2, # polynomial order for free flow velocity - order_p = order_u - 1, # polynomial order for free flow pressure - order_v = order_u - 1, # polynomial order for porous media velocity - nrefs = 4, # number of mesh refinements - Plotter = nothing, # backend for Plotting (e.g. GLMakie) - parallel = false, # do parallel assembly? - npart = 8, # number of partitions for grid coloring (if parallel = true) - kwargs... - ) - - ## region numbers - rFF = 1 # free flow region - rPM = 2 # porous media region - - ## load mesh and refine - xgrid = - uniform_refine( - simplexgrid( - Triangulate; - points = [0 0; 0 -1; 1 -1; 1 0; 1 1; 0 1]', - bfaces = [1 2; 2 3; 3 4; 4 5; 5 6; 6 1; 1 4]', - bfaceregions = [1, 2, 3, 4, 5, 6, 7], - regionpoints = [0.5 0.5; 0.5 -0.5]', - regionnumbers = [1, 2], - regionvolumes = [1.0, 1.0] - ), nrefs - ) - - if parallel - xgrid = partition(xgrid, RecursiveMetisPartitioning(npart = npart)) - end - - ## define unknowns - u = Unknown("u"; name = "velocity FF", dim = 2) - p = Unknown("p"; name = "pressure FF", dim = 1) - v = Unknown("q"; name = "velocity PM", dim = 2) - q = Unknown("q"; name = "pressure PM", dim = 2) - λ = Unknown("λ"; name = "interface LM", dim = 1) - - ## problem description - PD = ProblemDescription("Stokes-Darcy problem") - assign_unknown!(PD, u) - assign_unknown!(PD, p) - assign_unknown!(PD, v) - assign_unknown!(PD, q) - assign_unknown!(PD, λ) - if coupled - ## when coupled no boundary conditions on boundary region 7 - bnd_stokes = [4, 5, 6] - bnd_darcy = [1, 2, 3] - else - bnd_stokes = [4, 5, 6, 7] - bnd_darcy = [1, 2, 3, 7] - end - - ## define free flow problem: Stokes equations to solve for velocity u and pressure p - assign_operator!(PD, BilinearOperator(kernel_stokes_standard!, [grad(u), id(p)]; params = [μ], parallel = parallel, regions = [rFF], kwargs...)) - assign_operator!(PD, LinearOperator(f_FF!, [id(u)]; bonus_quadorder = 4, parallel = parallel, regions = [rFF], kwargs...)) - assign_operator!(PD, InterpolateBoundaryData(u, u!; bonus_quadorder = 4, regions = bnd_stokes, kwargs...)) - - ## define porous media flow: Darcy problem to solve for pressure q - assign_operator!(PD, BilinearOperator([id(v)]; parallel = parallel, regions = [rPM], kwargs...)) - assign_operator!(PD, BilinearOperator([div(v)], [id(q)]; factor = -1, transposed_copy = -1, parallel = parallel, regions = [rPM], kwargs...)) - assign_operator!(PD, LinearOperator(f_PM!, [id(q)]; bonus_quadorder = 4, parallel = parallel, regions = [rPM], kwargs...)) - - ## coupling conditions on interface (bregion 7) - if coupled - ## interface condition for tangential part - assign_operator!(PD, BilinearOperator([tangentialflux(u)]; factor = α * μ / sqrt(μ * k), entities = ON_FACES, regions = [7])) - ## interface condition for normal part via Lagrange multiplier - assign_operator!(PD, BilinearOperator(coupling_normal!, [id(λ)], [normalflux(u), normalflux(v)]; transposed_copy = 1, entities = ON_FACES, regions = [7])) - else - ## dummy problem for Lagrange multiplier - assign_operator!(PD, BilinearOperator([id(λ)]; entities = ON_FACES, parallel = parallel, regions = [7])) - ## Dirichlet condition for Darcy problem with different factor (due to different direction of normal) - assign_operator!(PD, LinearOperator(q!, [normalflux(v)]; bonus_quadorder = 4, factor = 1, entities = ON_FACES, params = [k], regions = [7], kwargs...)) - end - assign_operator!(PD, LinearOperator(q!, [normalflux(v)]; bonus_quadorder = 4, factor = -1, entities = ON_BFACES, params = [k], regions = setdiff(bnd_darcy, 7), kwargs...)) - - ## generate FESpaces and a solution vector for all 3 unknowns - FEType_λ = order_v == 0 ? L2P0{1} : H1Pk{1, 1, min(order_v, order_u)} - FES = [ - FESpace{H1Pk{2, 2, order_u}}(xgrid; regions = [rFF]), - FESpace{H1Pk{1, 2, order_p}}(xgrid; regions = [rFF]), - FESpace{HDIVRTk{2, order_v}}(xgrid; regions = [rPM]), - FESpace{H1Pk{1, 2, order_v}}(xgrid; regions = [rPM], broken = true), - FESpace{FEType_λ, ON_FACES}(xgrid; regions = [7], broken = true), - ] - sol = FEVector(FES; tags = PD.unknowns) - - ## solve the coupled problem - sol = solve(PD, [FES[1], FES[2], FES[3], FES[4], FES[5]]; kwargs...) - - ## correct integral mean of pressures - pmean_exact = integrate(xgrid, ON_CELLS, p!, 1; quadorder = 10, regions = [rFF]) - qmean_exact = integrate(xgrid, ON_CELLS, q!, 1; params = [k], quadorder = 10, regions = [rPM]) - pintegrate = ItemIntegrator([id(p)]; regions = [rFF]) - qintegrate = ItemIntegrator([id(q)]; regions = [rPM]) - pmean = sum(evaluate(pintegrate, sol)) - qmean = sum(evaluate(qintegrate, sol)) - if coupled - view(sol[p]) .-= (pmean + qmean) - view(sol[q]) .-= (pmean + qmean) - else - view(sol[p]) .-= pmean - pmean_exact - view(sol[q]) .-= qmean - qmean_exact - end - - ## error calculation - ErrorIntegratorExactU = ItemIntegrator(exact_error!(u!), [id(u)]; bonus_quadorder = 4, regions = [rFF], kwargs...) - ErrorIntegratorExactV = ItemIntegrator(exact_error!(v!), [id(v)]; params = [k], bonus_quadorder = 4, regions = [rPM], kwargs...) - ErrorIntegratorExactP = ItemIntegrator(exact_error!(p!), [id(p)]; bonus_quadorder = 4, regions = [rFF], kwargs...) - ErrorIntegratorExactQ = ItemIntegrator(exact_error!(q!), [id(q)]; params = [k], bonus_quadorder = 4, regions = [rPM], kwargs...) - errorU = evaluate(ErrorIntegratorExactU, sol) - errorV = evaluate(ErrorIntegratorExactV, sol) - errorP = evaluate(ErrorIntegratorExactP, sol) - errorQ = evaluate(ErrorIntegratorExactQ, sol) - L2errorU = sqrt(sum(view(errorU, 1, :)) + sum(view(errorU, 2, :))) - L2errorV = sqrt(sum(view(errorV, 1, :)) + sum(view(errorV, 2, :))) - L2errorP = sqrt(sum(view(errorP, 1, :))) - L2errorQ = sqrt(sum(view(errorQ, 1, :))) - @info "L2error(u) = $L2errorU" - @info "L2error(v) = $L2errorV" - @info "L2error(p) = $L2errorP" - @info "L2error(q) = $L2errorQ" - - ## plot - plt = plot([id(u), id(v), id(p), id(q), grid(u)], sol; Plotter = Plotter, ncols = 2, rasterpoints = 12, width = 800, height = 1200) - - return [L2errorU, L2errorV, L2errorP, L2errorQ], plt -end - -generateplots = ExtendableFEM.default_generateplots(Example264_StokesDarcy, "example264.png") #hide -function runtests() #hide - errors, plt = main(; μ = 1, k = 1, σ = 1, nrefs = 3, coupled = true) #hide - @info errors - @test errors[1] <= 0.000804 #hide - @test errors[2] <= 0.004468 #hide - @test errors[3] <= 0.043904 #hide - @test errors[4] <= 0.002966 #hide - return nothing #hide -end #hide -end # module diff --git a/examples/Example265_FlowTransport.jl b/examples/Example265_FlowTransport.jl deleted file mode 100644 index 201250d4..00000000 --- a/examples/Example265_FlowTransport.jl +++ /dev/null @@ -1,235 +0,0 @@ -#= - -# 265 : Flow + Transport -([source code](@__SOURCE_URL__)) - -This example solve the Stokes problem in an Omega-shaped pipe and then uses the velocity in a transport equation for a species with a certain inlet concentration. -Altogether, we are looking for a velocity ``\mathbf{u}``, a pressure ``\mathbf{p}`` and a stationary species concentration ``\mathbf{c}`` such that -```math -\begin{aligned} -- \mu \Delta \mathbf{u} + \nabla p & = 0\\ -\mathrm{div}(\mathbf{u}) & = 0\\ -\mathbf{c}_t - \kappa \Delta \mathbf{c} + \mathbf{u} \cdot \nabla \mathbf{c} & = 0 -\end{aligned} -``` -with some viscosity parameter and diffusion parameter ``\kappa``. - -The diffusion coefficient for the species is chosen (almost) zero such that the isolines of the concentration should stay parallel from inlet to outlet. -For the discretisation of the convection term in the transport equation three possibilities can be chosen: - -1. Classical Bernardi--Raugel stationary finite element discretisations ``\mathbf{u}_h \cdot \nabla \mathbf{c}_h`` - [set FVtransport = false, reconstruct = false] -2. As in 1. but with divergence-free reconstruction operator in convection term ``\Pi_\text{reconst} \mathbf{u}_h \cdot \nabla \mathbf{c}_h`` - [set FVtransport = false, reconstruct = true] -3. Time-dependent upwind finite volume discretisation for ``\kappa = 0`` based on normal fluxes along the faces [set FVtransport = true] - -Observe that the divergence-free postprocessing helps a lot for mass conservation, but is still not perfect. The finite volume -upwind discretisation ensures mass conservation. - -Note, that the transport equation is very convection-dominated and no stabilisation in the finite element discretisations was used here (but instead a nonzero ``\kappa``). -Also note, that only the finite volume discretisation perfectly obeys the maximum principle for the concentration but the isolines do no stay -parallel until the outlet is reached, possibly due to artificial diffusion. - -The computed solution for the default parameters looks like this: - -![](example265.png) - -=# - -module Example265_FlowTransport - -using ExtendableFEM -using ExtendableFEMBase -using ExtendableGrids -using SimplexGridFactory -using Triangulate -using Test #hide - -## boundary data -function u_inlet!(result, qpinfo) - x = qpinfo.x - result[1] = 4 * x[2] * (1 - x[2]) - result[2] = 0 - return nothing -end -function c_inlet!(result, qpinfo) - result[1] = (1 - qpinfo.x[2]) * qpinfo.x[2] - return nothing -end - -function kernel_stokes_standard!(result, u_ops, qpinfo) - ∇u, p = view(u_ops, 1:4), view(u_ops, 5) - μ = qpinfo.params[1] - result[1] = μ * ∇u[1] - p[1] - result[2] = μ * ∇u[2] - result[3] = μ * ∇u[3] - result[4] = μ * ∇u[4] - p[1] - result[5] = -(∇u[1] + ∇u[4]) - return nothing -end - -function kernel_convection!(result, ∇T, u, qpinfo) - result[1] = ∇T[1] * u[1] + ∇T[2] * u[2] - return nothing -end - -function kernel_inlet!(result, input, qpinfo) - c_inlet!(result, qpinfo) - result[1] *= -input[1] - return nothing -end - - -## everything is wrapped in a main function -function main(; nrefs = 4, Plotter = nothing, reconstruct = true, FVtransport = true, parallel = false, npart = 8, μ = 1, kwargs...) - - ## load mesh and refine - xgrid = uniform_refine( - simplexgrid( - Triangulate; - points = [0 0; 3 0; 3 -3; 7 -3; 7 0; 10 0; 10 1; 6 1; 6 -2; 4 -2; 4 1; 0 1]', - bfaces = [1 2; 2 3; 3 4; 4 5; 5 6; 6 7; 7 8; 8 9; 9 10; 10 11; 11 12; 12 1]', - bfaceregions = [1; 1; 1; 1; 1; 2; 3; 3; 3; 3; 3; 4], - regionpoints = [0.5 0.5;]', - regionnumbers = [1], - regionvolumes = [1.0] - ), nrefs - ) - - if parallel - xgrid = partition(xgrid, RecursiveMetisPartitioning(npart = npart)) - end - - ## define unknowns - u = Unknown("u"; name = "velocity", dim = 2) - p = Unknown("p"; name = "pressure", dim = 1) - T = Unknown("T"; name = "temperature", dim = 1) - - id_u = reconstruct ? apply(u, Reconstruct{HDIVBDM1{2}, Identity}) : id(u) - - ## define first sub-problem: Stokes equations to solve for velocity u - PD = ProblemDescription("Stokes problem") - assign_unknown!(PD, u) - assign_unknown!(PD, p) - assign_operator!(PD, BilinearOperator(kernel_stokes_standard!, [grad(u), id(p)]; params = [μ], parallel = parallel, kwargs...)) - assign_operator!(PD, InterpolateBoundaryData(u, u_inlet!; regions = 4, parallel = parallel, kwargs...)) - assign_operator!(PD, HomogeneousBoundaryData(u; regions = [1, 3], kwargs...)) - - ## add transport equation of species - PDT = ProblemDescription("transport problem") - assign_unknown!(PDT, T) - if FVtransport ## FVM discretisation of transport equation (pure upwind convection) - τ = 1.0e3 - assign_operator!(PDT, CallbackOperator(assemble_fv_operator!(), [u]; kwargs...)) - assign_operator!(PDT, BilinearOperator([id(T)]; store = true, factor = 1 / τ, parallel = parallel, kwargs...)) - assign_operator!(PDT, LinearOperator([id(T)], [id(T)]; factor = 1 / τ, parallel = parallel, kwargs...)) - else ## FEM discretisation of transport equation (with small diffusion term) - assign_operator!(PDT, BilinearOperator([grad(T)]; factor = 1.0e-6, parallel = parallel, kwargs...)) - assign_operator!(PDT, BilinearOperator(kernel_convection!, [id(T)], [grad(T)], [id_u]; parallel = parallel, kwargs...)) - assign_operator!(PDT, InterpolateBoundaryData(T, c_inlet!; regions = [4], kwargs...)) - end - - ## generate FESpaces and a solution vector for all 3 unknowns - FETypes = [H1BR{2}, L2P0{1}, FVtransport ? L2P0{1} : H1P1{1}] - FES = [FESpace{FETypes[j]}(xgrid) for j in 1:3] - sol = FEVector(FES; tags = [u, p, T]) - - ## solve the two problems separately - sol = solve(PD; init = sol, kwargs...) - sol = solve(PDT; init = sol, maxiterations = 20, target_residual = 1.0e-12, constant_matrix = true, kwargs...) - - ## print minimal and maximal concentration to check max principle (should be in [0,1]) - println("\n[min(c),max(c)] = [$(minimum(view(sol[T]))),$(maximum(view(sol[T])))]") - - ## plot - plt = plot([id(u), id(T)], sol; Plotter = Plotter, ncols = 1, rasterpoints = 40, width = 800, height = 800) - - return sol, plt -end - -## pure convection finite volume operator for transport -function assemble_fv_operator!() - - BndFluxIntegrator = ItemIntegrator(kernel_inflow!, [normalflux(1)]; entities = ON_BFACES) - FluxIntegrator = ItemIntegrator([normalflux(1)]; entities = ON_FACES) - fluxes::Matrix{Float64} = zeros(Float64, 1, 0) - - return function closure(A, b, args; assemble_matrix = true, assemble_rhs = true, kwargs...) - - ## prepare grid and stash - xgrid = args[1].FES.xgrid - nfaces = size(xgrid[FaceCells], 2) - if size(fluxes, 2) < nfaces - fluxes = zeros(Float64, 1, nfaces) - end - - ## right-hand side = boundary inflow fluxes if velocity points inward - if assemble_rhs - fill!(fluxes, 0) - evaluate!(fluxes, BndFluxIntegrator, [args[1]]) - facecells = xgrid[FaceCells] - bface2face = xgrid[BFaceFaces] - for bface in 1:lastindex(bface2face) - b[facecells[1, bface2face[bface]]] -= fluxes[bface] - end - end - - ## assemble upwind finite volume fluxes over cell faces into matrix - if assemble_matrix - ## integrate normalfux of velocity - fill!(fluxes, 0) - evaluate!(fluxes, FluxIntegrator, [args[1]]) - - cellfaces = xgrid[CellFaces] - cellfacesigns = xgrid[CellFaceSigns] - for cell in 1:num_cells(xgrid) - nfaces4cell = num_targets(cellfaces, cell) - for cf in 1:nfaces4cell - face = cellfaces[cf, cell] - other_cell = facecells[1, face] - if other_cell == cell - other_cell = facecells[2, face] - end - flux = fluxes[face] * cellfacesigns[cf, cell] - if (other_cell > 0) - flux *= 1 // 2 # because it will be accumulated on two cells - end - if flux > 0 # flow from cell to other_cell or out of domain - _addnz(A, cell, cell, flux, 1) - if other_cell > 0 - _addnz(A, other_cell, cell, -flux, 1) - ## otherwise flow goes out of domain - end - else # flow from other_cell into cell or into domain - _addnz(A, cell, cell, 1.0e-16, 1) # add zero to keep pattern for LU - if other_cell > 0 # flow comes from neighbour cell - _addnz(A, other_cell, other_cell, -flux, 1) - _addnz(A, cell, other_cell, flux, 1) - end - ## otherwise flow comes from outside into domain, handled in rhs side loop above - end - end - end - end - return nothing - end -end - - -function kernel_inflow!(result, input, qpinfo) - return if input[1] < 0 # if velocity points into domain - c_inlet!(result, qpinfo) - result[1] *= input[1] - else - result[1] = 0 - end -end - -generateplots = ExtendableFEM.default_generateplots(Example265_FlowTransport, "example265.png") #hide -function runtests() #hide - sol, plt = main() #hide - @test minimum(view(sol[3])) >= 0 #hide - @test maximum(view(sol[3])) <= 0.25 #hide - return nothing #hide -end #hide -end # module diff --git a/examples/Example270_NaturalConvectionProblem.jl b/examples/Example270_NaturalConvectionProblem.jl deleted file mode 100644 index 7f07d1f0..00000000 --- a/examples/Example270_NaturalConvectionProblem.jl +++ /dev/null @@ -1,160 +0,0 @@ -#= - -# 270 : Natural convection -([source code](@__SOURCE_URL__)) - -Seek velocity ``\mathbf{u}``, pressure ``p`` and temperature ``\theta`` such that -```math -\begin{aligned} - - \mu \Delta u + (\mathbf{u} \cdot \nabla) \mathbf{u} + \nabla p & = Ra \, \theta \, g \\ - - \Delta \theta + \mathbf{u} \cdot \nabla \theta & = 0 -\end{aligned} -``` -on a given domain ``\Omega`` (here a triangle) and boundary conditions -```math -\begin{aligned} - \mathbf{u} & = 0 && \quad \text{along } \partial \Omega\\ - T & = T_\text{bottom} &&\quad \text{along } y = 0\\ - T & = 0 &&\quad \text{along } x = 0 -\end{aligned} -``` - -The weak formulation seeks ``(\mathbf{u},p,\theta) \in V \times Q \times X \subseteq H^1_0(\Omega)^2 \times L^2_0(\Omega) \times H^1_D(\Omega)`` such that -```math -\begin{aligned} - (\mu \nabla \mathbf{u}, \nabla \mathbf{v}) + ((\mathbf{u} \cdot \nabla) \mathbf{u}, \mathbf{v}) - (\mathrm{div} \mathbf{v}, p) & = (\mathbf{v}, Ra g \, \theta) && \quad \text{for all } \mathbf{v} \in V,\\ -(\mathrm{div} \mathbf{u}, q) & = 0 && \quad \text{for all } q \in Q,\\ - (\nabla \theta, \nabla \varphi) + (u \cdot \nabla \theta, \varphi) & = 0 - && \quad \text{for all } \varphi \in X. -\end{aligned} -``` - -To render the discrete method pressure-robust, a reconstruction operator is applied to all identity evaluations of ``\mathbf{u}`` and ``\mathbf{v}`` -(when the switch reconstruct is set to true). -Further explanations and discussion on this example can be found in the reference below. - -!!! reference - - "On the divergence constraint in mixed finite element methods for incompressible flows",\ - V. John, A. Linke, C. Merdon, M. Neilan and L. Rebholz,\ - SIAM Review 59(3) (2017),\ - [>Journal-Link<](https://doi.org/10.1137/15M1047696) - [>Preprint-Link<](http://www.wias-berlin.de/publications/wias-publ/run.jsp?template=abstract&type=Preprint&year=2015&number=2177) - -The computed solution for the default parameters looks like this: - -![](example270.png) - -=# - - -module Example270_NaturalConvectionProblem - -using ExtendableFEM -using ExtendableGrids -using GridVisualize -using LinearAlgebra -using Test #hide - -function kernel_nonlinear!(result, u_ops, qpinfo) - u, ∇u, p, ∇T, T = view(u_ops, 1:2), view(u_ops, 3:6), view(u_ops, 7), view(u_ops, 8:9), view(u_ops, 10) - Ra, μ, ϵ = qpinfo.params[1], qpinfo.params[2], qpinfo.params[3] - result[1] = dot(u, view(∇u, 1:2)) - result[2] = dot(u, view(∇u, 3:4)) - Ra * T[1] - result[3] = μ * ∇u[1] - p[1] - result[4] = μ * ∇u[2] - result[5] = μ * ∇u[3] - result[6] = μ * ∇u[4] - p[1] - result[7] = -(∇u[1] + ∇u[4]) - result[8] = ϵ * ∇T[1] - result[9] = ϵ * ∇T[2] - result[10] = dot(u, ∇T) - return nothing -end - -function T_bottom!(result, qpinfo) - x = qpinfo.x - result[1] = 2 * (1 - cos(2 * π * x[1])) - return nothing -end - -function main(; - nrefs = 5, - μ = 1.0, - ϵ = 1.0, - Ra_final = 1.0e6, - reconstruct = true, - Plotter = nothing, - kwargs... - ) - - ## problem description - PD = ProblemDescription() - u = Unknown("u"; name = "velocity") - p = Unknown("p"; name = "pressure") - T = Unknown("T"; name = "temperature") - assign_unknown!(PD, u) - assign_unknown!(PD, p) - assign_unknown!(PD, T) - id_u = reconstruct ? apply(u, Reconstruct{HDIVBDM1{2}, Identity}) : id(u) - assign_operator!(PD, NonlinearOperator(kernel_nonlinear!, [id_u, grad(u), id(p), grad(T), id(T)]; kwargs...)) - assign_operator!(PD, HomogeneousBoundaryData(u; regions = 1:3)) - assign_operator!(PD, FixDofs(p; dofs = [1], vals = [0])) - assign_operator!(PD, HomogeneousBoundaryData(T; regions = 3)) - assign_operator!(PD, InterpolateBoundaryData(T, T_bottom!; regions = 1)) - - ## grid - xgrid = uniform_refine(reference_domain(Triangle2D), nrefs) - - ## FESpaces - FES = Dict( - u => FESpace{H1BR{2}}(xgrid), - p => FESpace{L2P0{1}}(xgrid), - T => FESpace{H1P1{1}}(xgrid) - ) - - ## prepare plots - plt = GridVisualizer(; Plotter = Plotter, layout = (1, 3), clear = true, size = (1200, 400)) - - ## solve by Ra embedding - params = Array{Float64, 1}([min(Ra_final, 4000), μ, ϵ]) - sol = nothing - SC = nothing - step = 0 - while (true) - ## solve (params are given to all operators) - sol, SC = ExtendableFEM.solve(PD, FES, SC; init = sol, return_config = true, target_residual = 1.0e-6, params = params, kwargs...) - - ## plot - scalarplot!(plt[1, 1], id(u), sol; levels = 0, colorbarticks = 7, abs = true) - vectorplot!(plt[1, 1], id(u), sol; clear = false, title = "|u| + quiver (Ra = $(params[1]))") - scalarplot!(plt[1, 2], id(T), sol; title = "T (Ra = $(params[1]))") - scalarplot!(plt[1, 3], id(p), sol; title = "p (Ra = $(params[1]))") - - ## stop if Ra_final is reached - if params[1] >= Ra_final - break - end - - ## increase Ra - params[1] = min(Ra_final, params[1] * 3) - step += 1 - @info "Step $step : solving for Ra=$(params[1])" - end - - ## compute Nusselt number along bottom (= boundary region 1) - ∇T_faces = FaceInterpolator([jump(grad(T))]; order = 0, kwargs...) - NuIntegrator = ItemIntegrator((result, input, qpinfo) -> (result[1] = -input[2]), [id(1)]; entities = ON_FACES, regions = [1]) - Nu = sum(evaluate(NuIntegrator, evaluate!(∇T_faces, sol))) - @info "Nu = $Nu" - - return Nu, plt -end - -generateplots = ExtendableFEM.default_generateplots(Example270_NaturalConvectionProblem, "example270.png") #hide -function runtests() #hide - Nu, plt = main(; nrefs = 4) #hide - @test Nu ≈ 17.641450080135293 #hide - return nothing #hide -end #hide -end # module diff --git a/examples/Example275_OptimalControlStokes.jl b/examples/Example275_OptimalControlStokes.jl deleted file mode 100644 index ff20340a..00000000 --- a/examples/Example275_OptimalControlStokes.jl +++ /dev/null @@ -1,134 +0,0 @@ -#= - -# 275 : Optimal Control Stokes -([source code](@__SOURCE_URL__)) - -This example studies the optimal control problem for the Stokes operator -with divergence-free velocity space ``\mathbf{V}_0 \subset \mathbf{H}^1_0``, i.e., -for given data ``\mathbf{u}^d`` minimize the functional -```math -\begin{aligned} -\min_{(\mathbf{u},\mathbf{q}) \in \mathbf{V}_0 \times \mathbf{L}^2} \| \mathbf{u} - \mathbf{u}^d \|^2 + \frac{\alpha}{2} \| \mathbf{q} \|^2 -\quad \text{s.t. } (\mu \nabla \mathbf{u}, \nabla \mathbf{v}) = (\mathbf{q}, \mathbf{v}) \quad \text{for all } \mathbf{v} \in \mathbf{V}_0 -\end{aligned} -``` - -This results in the set of variational equations that seeks -``(\mathbf{u}, \mathbf{z}, p , \lambda)`` such that -```math -\begin{aligned} -(\mu \nabla \mathbf{u}, \nabla \mathbf{v}) + (p, \mathrm{div} \mathbf{v}) & = - \alpha^{-1/2} (\mathbf{z}, \mathbf{v})\\ - (q, \mathrm{div} \mathbf{u}) & = 0\\ -(\mu \nabla \mathbf{z}, \nabla \mathbf{w}) + (λ, \mathrm{div} \mathbf{w}) & = \alpha^{-1/2} (\mathbf{u} - \mathbf{u}^d, \mathbf{w})\\ - (φ, \mathrm{div} \mathbf{z}) & = 0. -\end{aligned} -``` -for all test functions ``(\mathbf{v}, \mathbf{w}, q , \varphi)``. - -Here, we study pressure-robustness with the given data -```math -\mathbf{u}^d := \mathrm{curl} \left(x^4y^4(x-1)^4(y-1)^4\right) + \epsilon \nabla(\cos(x)\sin(y)) -``` -with a gradient field distortion that can be steered by ``ϵ \geq 0`` which was -an example in the reference below. - -!!! reference - - "Pressure-robustness in the context of optimal control",\ - C. Merdon and W. Wollner,\ - SIAM Journal on Control and Optimization 61:1, 342-360 (2023),\ - [>Journal-Link<](https://epubs.siam.org/doi/10.1137/22M1482603) - [>Preprint-Link<](https://arxiv.org/abs/2203.02224) - -The computed solution for the default parameters looks like this: - -![](example275.png) - -=# - -module Example275_OptimalControlStokes - -using ExtendableFEM -using ExtendableGrids -using Symbolics - -function prepare_data!(; ϵ = 0) - @variables x y - - ## stream function ξ - ξ = x^4 * y^4 * (x - 1)^4 * (y - 1)^4 - ∇ξ = Symbolics.gradient(ξ, [x, y]) - - ## irrotational perturbation (to study pressure-robustness) - ϕ = cos(x) * sin(y) - ∇ϕ = Symbolics.gradient(ϕ, [x, y]) - - ## final data = curl ξ + ϵ ∇ϕ - d = [-∇ξ[2], ∇ξ[1]] + ϵ * ∇ϕ - d_eval = build_function(d, x, y, expression = Val{false}) - - return d_eval[2] -end - -## standard Stokes kernel -function kernel_stokes_standard!(result, u_ops, qpinfo) - ∇u, p = view(u_ops, 1:4), view(u_ops, 5) - μ = qpinfo.params[1] - result[1] = μ * ∇u[1] + p[1] - result[2] = μ * ∇u[2] - result[3] = μ * ∇u[3] - result[4] = μ * ∇u[4] + p[1] - result[5] = (∇u[1] + ∇u[4]) - return nothing -end - -## everything is wrapped in a main function -function main(; nrefs = 4, Plotter = nothing, reconstruct = true, μ = 1, α = 1.0e-6, ϵ = 0, kwargs...) - - ## prepare target data - d_eval = prepare_data!(; ϵ = ϵ) - data!(result, qpinfo) = (d_eval(result, qpinfo.x[1], qpinfo.x[2]);) - - ## load mesh and refine - xgrid = uniform_refine(grid_unitsquare(Triangle2D), nrefs) - - ## define unknowns - u = Unknown("u"; name = "velocity", dim = 2) - z = Unknown("z"; name = "control", dim = 2) - p = Unknown("p"; name = "pressure", dim = 1) - λ = Unknown("λ"; name = "control pressure", dim = 1) - - ## prepare reconstruction operator (if reconstruct = true) - idR(u) = reconstruct ? apply(u, Reconstruct{HDIVBDM1{2}, Identity}) : id(u) - - ## define optimal control problem - PD = ProblemDescription("Stokes optimal control problem") - assign_unknown!(PD, u) - assign_unknown!(PD, z) - assign_unknown!(PD, p) - assign_unknown!(PD, λ) - assign_operator!(PD, BilinearOperator(kernel_stokes_standard!, [grad(u), id(p)]; params = [μ], kwargs...)) - assign_operator!(PD, BilinearOperator(kernel_stokes_standard!, [grad(z), id(λ)]; params = [μ], kwargs...)) - assign_operator!(PD, BilinearOperator([idR(z)], [idR(u)]; factor = -1 / sqrt(α), transposed_copy = -1, kwargs...)) - assign_operator!(PD, LinearOperator(data!, [idR(z)]; factor = -1 / sqrt(α), bonus_quadorder = 5, kwargs...)) - assign_operator!(PD, HomogeneousBoundaryData(u; regions = 1:4, kwargs...)) - assign_operator!(PD, HomogeneousBoundaryData(z; regions = 1:4, kwargs...)) - - ## solve with Bernardi--Raugel method - FETypes = [H1BR{2}, L2P0{1}] - FES = [FESpace{FETypes[j]}(xgrid) for j in 1:2] - sol = solve(PD, [FES[1], FES[1], FES[2], FES[2]]; kwargs...) - - ## plot solution - plt = plot([id(u), id(p), id(z), id(λ)], sol; add = 1, Plotter = Plotter) - - ## plot target data - I = FEVector(FES[1]; name = "u^d") - interpolate!(I[1], data!) - plot!(plt, [id(1)], I; keep = 1:4) - - return sol, plt -end - -generateplots = ExtendableFEM.default_generateplots(Example275_OptimalControlStokes, "example275.png") #hide -end diff --git a/examples/Example280_CompressibleStokes.jl b/examples/Example280_CompressibleStokes.jl deleted file mode 100644 index 8bf45a38..00000000 --- a/examples/Example280_CompressibleStokes.jl +++ /dev/null @@ -1,343 +0,0 @@ -#= - -# 280 : Compressible Stokes -([source code](@__SOURCE_URL__)) - -This example solves the two-dimensional -compressible Stokes equations where one seeks a (vector-valued) velocity -``\mathbf{u}``, a density ``\varrho`` and a pressure ``p`` such that -```math -\begin{aligned} -- \mu \Delta \mathbf{u} + \lambda \nabla(\mathrm{div}(\mathbf{u})) + \nabla p & = \mathbf{f} + \varrho \mathbf{g}\\ -\mathrm{div}(\varrho \mathbf{u}) & = 0\\ - p & = eos(\varrho)\\ - \int_\Omega \varrho \, dx & = M\\ - \varrho & \geq 0. -\end{aligned} -``` -Here, eos ``eos`` is some equation of state function that describes the dependence of the pressure on the density -(and further physical quantities like temperature in a more general setting). -Moreover, ``\mu`` and ``\lambda`` are Lame parameters and ``\mathbf{f}`` and ``\mathbf{g}`` are given right-hand side data. - -There are two testcases. The first testcase solves an analytical toy problem with the prescribed solution -```math -\begin{aligned} -\mathbf{u}(\mathbf{x}) & =0\\ -\varrho(\mathbf{x}) & = \exp(-y/c) \\ -p &= eos(\varrho) := c \varrho^\gamma -\end{aligned} -``` -such that ``\mathbf{f} = 0`` and ``\mathbf{g}`` nonzero to match the prescribed solution. The second testcase -tests an analytical nonzero velocity benchmark problem with the same density. - -This example is designed to study the well-balanced property of a discretisation. The gradient-robust discretisation -approximates the well-balanced state much better, i.e. has a much smaller L2 velocity error. For larger c -(= smaller Mach number) the problem gets more incompressible which reduces the error further as then the -right-hand side is a perfect gradient also when evaluated with the (now closer to a constant) discrete density. -See reference below for more details. - -!!! reference - - "A gradient-robust well-balanced scheme for the compressible isothermal Stokes problem",\ - M. Akbas, T. Gallouet, A. Gassmann, A. Linke and C. Merdon,\ - Computer Methods in Applied Mechanics and Engineering 367 (2020),\ - [>Journal-Link<](https://doi.org/10.1016/j.cma.2020.113069) - [>Preprint-Link<](https://arxiv.org/abs/1911.01295) - -The computed solution for the default parameters looks like this: - -![](example280.png) - -=# - -module Example280_CompressibleStokes - -using ExtendableFEM -using ExtendableGrids -using Triangulate -using SimplexGridFactory -using GridVisualize -using Symbolics -using LinearAlgebra -using Test #hide - -## everything is wrapped in a main function -## testcase = 1 : well-balanced test (stratified no-flow over mountain) -## testcase = 2 : vortex example (ϱu is div-free p7 vortex) -function main(; - testcase = 1, - nrefs = 4, - M = 1, - c = 1, - ufac = 100, - pressure_stab = 0, - laplacian_in_rhs = false, # for data in example 2 - maxsteps = 5000, - target_residual = 1.0e-11, - Plotter = nothing, - reconstruct = true, - μ = 1, - order = 1, - kwargs... - ) - - ## load data for testcase - grid_builder, kernel_gravity!, kernel_rhs!, u!, ∇u!, ϱ!, τfac = load_testcase_data(testcase; laplacian_in_rhs = laplacian_in_rhs, M = M, c = c, μ = μ, ufac = ufac) - xgrid = grid_builder(nrefs) - - ## define unknowns - u = Unknown("u"; name = "velocity", dim = 2) - ϱ = Unknown("ϱ"; name = "density", dim = 1) - p = Unknown("p"; name = "pressure", dim = 1) - - ## define reconstruction operator - if order == 1 - FETypes = [H1BR{2}, L2P0{1}, L2P0{1}] - id_u = reconstruct ? apply(u, Reconstruct{HDIVRT0{2}, Identity}) : id(u) - elseif order == 2 - FETypes = [H1P2B{2, 2}, L2P1{1}, L2P1{1}] - id_u = reconstruct ? apply(u, Reconstruct{HDIVRT1{2}, Identity}) : id(u) - end - - ## define first sub-problem: Stokes equations to solve for velocity u - PD = ProblemDescription("Stokes problem") - assign_unknown!(PD, u) - assign_operator!(PD, BilinearOperator([grad(u)]; factor = μ, store = true, kwargs...)) - assign_operator!(PD, LinearOperator([div(u)], [id(ϱ)]; factor = c, kwargs...)) - assign_operator!(PD, HomogeneousBoundaryData(u; regions = 1:4, kwargs...)) - if kernel_rhs! !== nothing - assign_operator!(PD, LinearOperator(kernel_rhs!, [id_u]; factor = 1, store = true, bonus_quadorder = 3 * order, kwargs...)) - end - assign_operator!(PD, LinearOperator(kernel_gravity!, [id_u], [id(ϱ)]; factor = 1, bonus_quadorder = 3 * order, kwargs...)) - - ## FVM for continuity equation - τ = μ / (order^2 * M * sqrt(τfac)) # time step for pseudo timestepping - @info "timestep = $τ" - PDT = ProblemDescription("continuity equation") - assign_unknown!(PDT, ϱ) - if order > 1 - assign_operator!(PDT, BilinearOperator(kernel_continuity!, [grad(ϱ)], [id(ϱ)], [id(u)]; quadorder = 2 * order, factor = -1, kwargs...)) - end - if pressure_stab > 0 - psf = pressure_stab #* xgrid[CellVolumes][1] - assign_operator!(PDT, BilinearOperator(stab_kernel!, [jump(id(ϱ))], [jump(id(ϱ))], [id(u)]; entities = ON_IFACES, factor = psf, kwargs...)) - end - assign_operator!(PDT, BilinearOperator([id(ϱ)]; quadorder = 2 * (order - 1), factor = 1 / τ, store = true, kwargs...)) - assign_operator!(PDT, LinearOperator([id(ϱ)], [id(ϱ)]; quadorder = 2 * (order - 1), factor = 1 / τ, kwargs...)) - assign_operator!(PDT, BilinearOperatorDG(kernel_upwind!, [jump(id(ϱ))], [this(id(ϱ)), other(id(ϱ))], [id(u)]; quadorder = order + 1, entities = ON_IFACES, kwargs...)) - - ## prepare error calculation - EnergyIntegrator = ItemIntegrator(energy_kernel!, [id(u)]; resultdim = 1, quadorder = 2 * (order + 1), kwargs...) - ErrorIntegratorExact = ItemIntegrator(exact_error!(u!, ∇u!, ϱ!), [id(u), grad(u), id(ϱ)]; resultdim = 9, quadorder = 2 * (order + 1), kwargs...) - NDofs = zeros(Int, nrefs) - Results = zeros(Float64, nrefs, 7) - - sol = nothing - xgrid = nothing - op_upwind = 0 - for lvl in 1:nrefs - xgrid = grid_builder(lvl) - @show xgrid - FES = [FESpace{FETypes[j]}(xgrid) for j in 1:3] - sol = FEVector(FES; tags = [u, ϱ, p]) - - ## initial guess - fill!(sol[ϱ], M) - interpolate!(sol[u], u!) - interpolate!(sol[ϱ], ϱ!) - NDofs[lvl] = length(sol.entries) - - ## solve the two problems iteratively [1] >> [2] >> [1] >> [2] ... - SC1 = SolverConfiguration(PD; init = sol, maxiterations = 1, target_residual = target_residual, constant_matrix = true, kwargs...) - SC2 = SolverConfiguration(PDT; init = sol, maxiterations = 1, target_residual = target_residual, kwargs...) - sol, nits = iterate_until_stationarity([SC1, SC2]; energy_integrator = EnergyIntegrator, maxsteps = maxsteps, init = sol, kwargs...) - - residual_momentum = residual(SC1) - residual_continuity = residual(SC2) - @info "final residual momentum = $(residual_momentum)" - @info "final residual continuity = $(residual_continuity)" - - ## calculate error - error = evaluate(ErrorIntegratorExact, sol) - Results[lvl, 1] = sqrt(sum(view(error, 1, :)) + sum(view(error, 2, :))) - Results[lvl, 2] = sqrt(sum(view(error, 3, :)) + sum(view(error, 4, :)) + sum(view(error, 5, :)) + sum(view(error, 6, :))) - Results[lvl, 3] = sqrt(sum(view(error, 7, :))) - Results[lvl, 4] = sqrt(sum(view(error, 8, :)) + sum(view(error, 9, :))) - Results[lvl, 5] = nits - Results[lvl, 6] = residual_momentum - Results[lvl, 7] = residual_continuity - - ## print results - print_convergencehistory(NDofs[1:lvl], Results[1:lvl, :]; X_to_h = X -> X .^ (-1 / 2), ylabels = ["|| u - u_h ||", "|| ∇(u - u_h) ||", "|| ϱ - ϱ_h ||", "|| ϱu - ϱu_h ||", "#its"], xlabel = "ndof") - end - - - ## plot - plt = GridVisualizer(; Plotter = Plotter, layout = (2, 2), clear = true, size = (1000, 1000)) - scalarplot!(plt[1, 1], xgrid, view(nodevalues(sol[u]; abs = true), 1, :), levels = 0, colorbarticks = 7) - vectorplot!(plt[1, 1], xgrid, eval_func_bary(PointEvaluator([id(u)], sol)), rasterpoints = 10, clear = false, title = "u_h (abs + quiver)") - scalarplot!(plt[2, 1], xgrid, view(nodevalues(sol[ϱ]), 1, :), levels = 11, title = "ϱ_h") - plot_convergencehistory!(plt[1, 2], NDofs, Results[:, 1:4]; add_h_powers = [order, order + 1], X_to_h = X -> 0.2 * X .^ (-1 / 2), legend = :best, ylabels = ["|| u - u_h ||", "|| ∇(u - u_h) ||", "|| ϱ - ϱ_h ||", "|| ϱu - ϱu_h ||", "#its"]) - gridplot!(plt[2, 2], xgrid) - - return Results, plt -end - -function stab_kernel!(result, p, u, qpinfo) - result[1] = p[1] #*abs(u[1] + u[2]) - return nothing -end - -## kernel for (uϱ, ∇λ) ON_CELLS in continuity equation -function kernel_continuity!(result, ϱ, u, qpinfo) - result[1] = ϱ[1] * u[1] - result[2] = ϱ[1] * u[2] - return nothing -end - -## kernel for (u⋅n ϱ^upw, λ) ON_IFACES in continuity equation -function kernel_upwind!(result, input, u, qpinfo) - flux = dot(u, qpinfo.normal) # u * n - return if flux > 0 - result[1] = input[1] * flux # rho_left * flux - else - result[1] = input[2] * flux # rho_righ * flux - end -end - -## kernel for exact error calculation -function exact_error!(u!, ∇u!, ϱ!) - return function closure(result, u, qpinfo) - u!(view(result, 1:2), qpinfo) - ∇u!(view(result, 3:6), qpinfo) - ϱ!(view(result, 7), qpinfo) - result[8] = result[1] * result[7] - result[9] = result[2] * result[7] - view(result, 1:7) .-= u - result[8] -= u[1] * u[7] - result[9] -= u[2] * u[7] - return result .= result .^ 2 - end -end - -## kernel for gravity term in testcase 1 -function standard_gravity!(result, ϱ, qpinfo) - result[1] = 0 - result[2] = -ϱ[1] - return nothing -end - -function energy_kernel!(result, u, qpinfo) - result[1] = dot(u, u) / 2 - return nothing -end - -function load_testcase_data(testcase::Int = 1; laplacian_in_rhs = true, M = 1, c = 1, μ = 1, ufac = 100) - if testcase == 1 - grid_builder = (nref) -> simplexgrid( - Triangulate; - points = [0 0; 0.2 0; 0.3 0.2; 0.45 0.05; 0.55 0.35; 0.65 0.2; 0.7 0.3; 0.8 0; 1 0; 1 1 ; 0 1]', - bfaces = [1 2; 2 3; 3 4; 4 5; 5 6; 6 7; 7 8; 8 9; 9 10; 10 11; 11 1]', - bfaceregions = ones(Int, 11), - regionpoints = [0.5 0.5;]', - regionnumbers = [1], - regionvolumes = [4.0^-(nref) / 2] - ) - xgrid = grid_builder(3) - u1!(result, qpinfo) = (fill!(result, 0);) - ∇u1!(result, qpinfo) = (fill!(result, 0);) - M_exact = integrate(xgrid, ON_CELLS, (result, qpinfo) -> (result[1] = exp(-qpinfo.x[2] / c) / M;), 1; quadorder = 20) - area = sum(xgrid[CellVolumes]) - ϱ1!(result, qpinfo) = (result[1] = exp(-qpinfo.x[2] / c) / (M_exact / area);) - return grid_builder, standard_gravity!, nothing, u1!, ∇u1!, ϱ1!, 1 - elseif testcase == 2 - grid_builder = (nref) -> simplexgrid( - Triangulate; - points = [0 0; 1 0; 1 1 ; 0 1]', - bfaces = [1 2; 2 3; 3 4; 4 1]', - bfaceregions = ones(Int, 4), - regionpoints = [0.5 0.5;]', - regionnumbers = [1], - regionvolumes = [4.0^-(nref)] - ) - - xgrid = grid_builder(3) - M_exact = integrate(xgrid, ON_CELLS, (result, qpinfo) -> (result[1] = exp(-qpinfo.x[1]^3 / (3 * c));), 1; quadorder = 20) - ϱ_eval, g_eval, f_eval, u_eval, ∇u_eval = prepare_data!(; laplacian_in_rhs = laplacian_in_rhs, M = M_exact, c = c, μ = μ, ufac = ufac) - ϱ2!(result, qpinfo) = (result[1] = ϱ_eval(qpinfo.x[1], qpinfo.x[2]);) - - M_exact = integrate(xgrid, ON_CELLS, ϱ2!, 1) - area = sum(xgrid[CellVolumes]) - - function kernel_gravity!(result, input, qpinfo) - g_eval(result, qpinfo.x[1], qpinfo.x[2]) - return result .*= input[1] - end - - function kernel_rhs!(result, qpinfo) - return f_eval(result, qpinfo.x[1], qpinfo.x[2]) - end - - u2!(result, qpinfo) = (u_eval(result, qpinfo.x[1], qpinfo.x[2]);) - ∇u2!(result, qpinfo) = (∇u_eval(result, qpinfo.x[1], qpinfo.x[2]);) - return grid_builder, kernel_gravity!, f_eval === nothing ? nothing : kernel_rhs!, u2!, ∇u2!, ϱ2!, ufac - end -end - -## exact data for testcase 2 computed by Symbolics -function prepare_data!(; M = 1, c = 1, μ = 1, ufac = 100, laplacian_in_rhs = true) - - @variables x y - - ## density - ϱ = exp(-x^3 / (3 * c)) / M - - ## stream function ξ - ## sucht that ϱu = curl ξ - ξ = x^2 * y^2 * (x - 1)^2 * (y - 1)^2 * ufac - - ∇ξ = Symbolics.gradient(ξ, [x, y]) - - ## velocity u = curl ξ / ϱ - u = [-∇ξ[2], ∇ξ[1]] ./ ϱ - - ## gradient of velocity - ∇u = Symbolics.jacobian(u, [x, y]) - ∇u_reshaped = [∇u[1, 1], ∇u[1, 2], ∇u[2, 1], ∇u[2, 2]] - - ## Laplacian - Δu = [ - (Symbolics.gradient(∇u[1, 1], [x]) + Symbolics.gradient(∇u[1, 2], [y]))[1], - (Symbolics.gradient(∇u[2, 1], [x]) + Symbolics.gradient(∇u[2, 2], [y]))[1], - ] - - ## gravity ϱg = - Δu + ϱ∇log(ϱ) - - if laplacian_in_rhs - f = - μ * Δu - g = c * Symbolics.gradient(log(ϱ), [x, y]) - else - g = - μ * Δu / ϱ + c * Symbolics.gradient(log(ϱ), [x, y]) - f = 0 - end - - #Δu = Symbolics.derivative(∇u[1,1], [x]) + Symbolics.derivative(∇u[2,2], [y]) - - ϱ_eval = build_function(ϱ, x, y, expression = Val{false}) - u_eval = build_function(u, x, y, expression = Val{false}) - ∇u_eval = build_function(∇u_reshaped, x, y, expression = Val{false}) - g_eval = build_function(g, x, y, expression = Val{false}) - f_eval = build_function(f, x, y, expression = Val{false}) - - return ϱ_eval, g_eval[2], f == 0 ? nothing : f_eval[2], u_eval[2], ∇u_eval[2] -end - -generateplots = ExtendableFEM.default_generateplots(Example280_CompressibleStokes, "example280.png") #hide -function runtests() #hide - Results, plt = main(; nrefs = 2) #hide - @test Results[end, 6] <= 1.0e-11 #hide - @test Results[end, 7] <= 1.0e-11 #hide - @test Results[end, 1] ≈ 6.732891488265023e-7 #hide - return nothing #hide -end #hide -end diff --git a/examples/Example282_IncompressibleMHD.jl b/examples/Example282_IncompressibleMHD.jl deleted file mode 100644 index 134b0881..00000000 --- a/examples/Example282_IncompressibleMHD.jl +++ /dev/null @@ -1,117 +0,0 @@ -#= - -# 282 : Incompressible MHD -([source code](@__SOURCE_URL__)) - -This example yields a prototype for the stationary incompressible viscous MHD equations that seek -a velocity field ``\mathbf{u}``, a pressure field ``p`` and a divergence-free magnetic field ``\mathbf{B}`` such that -```math -\begin{aligned} -- \mu \Delta \mathbf{u} + \nabla \cdot (\mathbf{u} \otimes \mathbf{u} - \mathbf{B} \otimes \mathbf{B}) + \nabla (p + \frac{1}{2} \mathbf{B} \cdot \mathbf{B}) -& = 0\\ -\mathrm{div}(\mathbf{u}) & = 0\\ -- \eta \Delta \mathbf{B} + \nabla \cdot (\mathbf{u} \otimes \mathbf{B} - \mathbf{B} \otimes \mathbf{u}) & = 0\\ -\mathrm{div}(\mathbf{B}) & = 0\\ -\end{aligned} -``` -on a rectangular 2D domain. Here, ``\mu`` and ``\eta`` are the viscosity and resistivity of the fluid and the magnetic field, respectively. - -![](example282.png) - -=# - -module Example282_IncompressibleMHD - -using ExtendableFEM -using ExtendableGrids -using LinearAlgebra - -function f!(result, qpinfo) - return result .= 0 -end - -function g!(result, qpinfo) - x = qpinfo.x - result[1] = sin(2 * pi * x[2]) * cos(pi * x[1]) - result[2] = 0 - return nothing -end - -function kernel_nonlinear!(result, u_ops, qpinfo) - u, B, ∇u, ∇B, p, q = view(u_ops, 1:2), view(u_ops, 3:4), view(u_ops, 5:8), view(u_ops, 9:12), view(u_ops, 13), view(u_ops, 14) - μ = qpinfo.params[1] - η = qpinfo.params[2] - - ## viscous terms and pressures - result[5] = μ * ∇u[1] - p[1] - result[6] = μ * ∇u[2] - result[7] = μ * ∇u[3] - result[8] = μ * ∇u[4] - p[1] - result[9] = η * ∇B[1] - q[1] - result[10] = η * ∇B[2] - result[11] = η * ∇B[3] - result[12] = η * ∇B[4] - q[1] - - ## Lorentz force - result[1] = - dot(B, view(∇B, 1:2)) - result[2] = - dot(B, view(∇B, 3:4)) - BdotB = (B[1] * B[1] + B[2] * B[2]) / 2 - result[5] -= BdotB - result[8] -= BdotB - - ## convection term for u and B - result[1] += dot(u, view(∇u, 1:2)) - result[2] += dot(u, view(∇u, 3:4)) - result[3] = dot(u, view(∇B, 1:2)) - dot(B, view(∇u, 1:2)) - result[4] = dot(u, view(∇B, 3:4)) - dot(B, view(∇u, 3:4)) - - ## divergence constraint - result[13] = -(∇u[1] + ∇u[4]) - result[14] = -(∇B[1] + ∇B[4]) - return nothing -end - - -## everything is wrapped in a main function -function main(; Plotter = nothing, μ = 1.0e-3, η = 1.0e-1, nrefs = 5, kwargs...) - - ## load grid (see function below) - xgrid = uniform_refine(grid_unitsquare(Triangle2D), nrefs) - - ## problem description - PD = ProblemDescription() - u = Unknown("u"; name = "velocity") - B = Unknown("B"; name = "magnetic field") - p = Unknown("p"; name = "pressure") - q = Unknown("q"; name = "magnetic pressure") - - assign_unknown!(PD, u) - assign_unknown!(PD, B) - assign_unknown!(PD, p) - assign_unknown!(PD, q) - - assign_operator!(PD, NonlinearOperator(kernel_nonlinear!, [id(u), id(B), grad(u), grad(B), id(p), id(q)]; bonus_quadorder = 2, params = [μ, η], kwargs...)) - assign_operator!(PD, LinearOperator(f!, [id(u)])) - assign_operator!(PD, LinearOperator(g!, [id(B)])) - assign_operator!(PD, HomogeneousBoundaryData(u; regions = 1:4)) - assign_operator!(PD, HomogeneousBoundaryData(B; regions = [1])) - assign_operator!(PD, FixDofs(p; dofs = [1])) - assign_operator!(PD, FixDofs(q; dofs = [1])) - - ## P2-bubble finite element method - FETypes = [H1P2{2, 2}, H1P2{2, 2}, H1P1{1}, H1P1{1}] - - ## generate FESpaces and Solution vector - FES = [FESpace{FETypes[j]}(xgrid) for j in 1:4] - - ## solve - sol = ExtendableFEM.solve(PD, FES; target_residual = 1.0e-8, time = 0, kwargs...) - - ## plot - plt = plot([id(u), id(B), id(p), id(q)], sol; Plotter = Plotter) - - return sol, plt -end - -generateplots = ExtendableFEM.default_generateplots(Example282_IncompressibleMHD, "example282.png") #hide -end diff --git a/examples/Example284_LevelSetMethod.jl b/examples/Example284_LevelSetMethod.jl deleted file mode 100644 index 3079b1df..00000000 --- a/examples/Example284_LevelSetMethod.jl +++ /dev/null @@ -1,118 +0,0 @@ -#= - -# 284 : Level Set Method -([source code](@__SOURCE_URL__)) - -This example studies the level-set method of some level function ``\mathbf{\phi}`` convected in time via the equation -```math -\begin{aligned} -\phi_t + \mathbf{u} \cdot \nabla \phi & = 0. -\end{aligned} -``` -Here this is tested with the (conservative) initial level set function ``\phi(x) = 0.5 \tanh((\lvert x - (0.5,0.75) \rvert - 0.15)/(2ϵ) + 1)`` -such that the level ``\phi \equiv 0.5`` forms a circle which is then convected by the velocity -``\mathbf{u} = \mathrm{curl} \pi^{-1} \sin^2(\pi x) \sin^2(\pi y)``. -No reinitialisation step is performed. - -The initial condition and the solution at ``T = 1`` for the default parameters looks like this: - -![](example284.png) - -=# - -module Example284_LevelSetMethod - -using ExtendableFEM -using ExtendableGrids -using GridVisualize -using LinearAlgebra -using OrdinaryDiffEqSDIRK - -function ϕ_init!(result, qpinfo) - x = qpinfo.x - ϵ = qpinfo.params[1] - result[1] = 1 / 2 * (tanh((sqrt((x[1] - 0.5)^2 + (x[2] - 0.75)^2) - 0.15) / (2 * ϵ)) + 1) - return nothing -end - -function velocity!(result, qpinfo) - result[1] = 0.5 - result[2] = 1.0 - result[1] = -2 * cos(π * qpinfo.x[2]) * sin(π * qpinfo.x[2]) * sin(π * qpinfo.x[1])^2 - result[2] = 2 * cos(π * qpinfo.x[1]) * sin(π * qpinfo.x[1]) * sin(π * qpinfo.x[2])^2 - return nothing -end - - -function kernel_convection!() - u = zeros(Float64, 2) - return function closure(result, input, qpinfo) - velocity!(u, qpinfo) - result[1] = dot(u, input) - return nothing - end -end - -## everything is wrapped in a main function -function main(; - Plotter = nothing, ϵ = 0.05, τ = 1.0e-2, T = 1.0, order = 2, nref = 6, use_diffeq = false, - solver = ImplicitEuler(), verbosity = -1, kwargs... - ) - - ## initial grid and final time - xgrid = uniform_refine(grid_unitsquare(Triangle2D), nref) - - ## define main level set problem - PD = ProblemDescription("level set problem") - ϕ = Unknown("ϕ"; name = "level set function") - assign_unknown!(PD, ϕ) - assign_operator!(PD, BilinearOperator(kernel_convection!(), [id(ϕ)], [grad(ϕ)]; kwargs...)) - assign_operator!(PD, HomogeneousBoundaryData(ϕ; value = 1, regions = 1:4, kwargs...)) - - ## generate FESpace and solution vector and interpolate initial state - FES = FESpace{H1Pk{1, 2, order}}(xgrid) - sol = FEVector(FES; tags = PD.unknowns) - interpolate!(sol[ϕ], ϕ_init!; params = [ϵ]) - - ## prepare plot and plot init solution - plt = GridVisualizer(; Plotter = Plotter, layout = (1, 2), clear = true, resolution = (800, 400)) - scalarplot!(plt[1, 1], id(ϕ), sol; levels = [0.5], flimits = [-0.05, 1.05], colorbarticks = [0, 0.25, 0.5, 0.75, 1], title = "ϕ (t = 0)") - - if (use_diffeq) - ## generate DifferentialEquations.ODEProblem - prob = generate_ODEProblem(PD, FES, (0.0, T); init = sol, constant_matrix = true) - - ## solve ODE problem - de_sol = DifferentialEquations.solve(prob, solver, abstol = 1.0e-6, reltol = 1.0e-4, dt = τ, dtmin = 1.0e-8, adaptive = true) - @info "#tsteps = $(length(de_sol.u))" - - ## get final solution - sol.entries .= de_sol.u[end] - else - ## add backward Euler time derivative - M = FEMatrix(FES) - assemble!(M, BilinearOperator([id(1)])) - assign_operator!(PD, BilinearOperator(M, [ϕ]; factor = 1 / τ, kwargs...)) - assign_operator!(PD, LinearOperator(M, [ϕ], [ϕ]; factor = 1 / τ, kwargs...)) - - ## generate solver configuration - SC = SolverConfiguration(PD, FES; init = sol, maxiterations = 1, constant_matrix = true, verbosity = verbosity, kwargs...) - - ## iterate tspan - t = 0 - for it in 1:Int(floor(T / τ)) - 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 - - ## plot final state - 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)") - - return sol, plt -end - -generateplots = ExtendableFEM.default_generateplots(Example284_LevelSetMethod, "example284.png") #hide -end diff --git a/examples/Example285_CahnHilliard.jl b/examples/Example285_CahnHilliard.jl deleted file mode 100644 index 411cc8e2..00000000 --- a/examples/Example285_CahnHilliard.jl +++ /dev/null @@ -1,139 +0,0 @@ -#= - -# 285 : Cahn-Hilliard Equations -([source code](@__SOURCE_URL__)) - -This example studies the mixed form of the Cahn-Hilliard equations that seeks -``(c,\mu)`` such that -```math -\begin{aligned} -c_t - \mathbf{div} (M \nabla \mu) & = 0\\ -\mu - \partial f / \partial c + \lambda \nabla^2c & = 0. -\end{aligned} -``` - -with ``f(c) = 100c^2(1-c)^2``, constant parameters ``M`` and ``\lambda`` and (random) -initial concentration as defined in the code below. - -The computed solution at different timesteps -for the default parameters and a randomized initial state look like this: - -![](example285.png) - -=# - -module Example285_CahnHilliard - -using ExtendableFEM -using ExtendableGrids -using GridVisualize -using ForwardDiff -using Random -Random.seed!(135791113) - -## parameters and initial condition -const f = (c) -> 100 * c^2 * (1 - c)^2 -const dfdc = (c) -> ForwardDiff.derivative(f, c) - -function c0!(result, qpinfo) - result[1] = 0.63 + 0.02 * (0.5 - rand()) - return nothing -end - -## everything is wrapped in a main function -function main(; - order = 2, # finite element order for c and μ - nref = 4, # refinement level - periodic = true, - M = 1.0, - λ = 1.0e-2, - iterations_until_next_plot = 20, - τ = 5 / 1000000, # time step (for main evolution phase) - τ_increase = 1.1, # increase factor for τ after each plot - Plotter = nothing, # Plotter (e.g. PyPlot) - kwargs..., - ) - - ## initial grid and final time - xgrid = uniform_refine(grid_unitsquare(Triangle2D), nref) - - ## define unknowns - c = Unknown("c"; name = "concentration", dim = 1) - μ = Unknown("μ"; name = "chemical potential", dim = 1) - - ## define main level set problem - PD = ProblemDescription("Cahn-Hilliard equation") - assign_unknown!(PD, c) - assign_unknown!(PD, μ) - assign_operator!(PD, BilinearOperator([grad(c)], [grad(μ)]; factor = M, store = true)) - assign_operator!(PD, BilinearOperator([id(μ)]; store = true)) - assign_operator!(PD, BilinearOperator([grad(μ)], [grad(c)]; factor = -λ, store = true)) - - if periodic - assign_restriction!(PD, CoupledDofsRestriction(c, 1, 3)) - assign_restriction!(PD, CoupledDofsRestriction(c, 2, 4)) - assign_restriction!(PD, CoupledDofsRestriction(μ, 1, 3)) - assign_restriction!(PD, CoupledDofsRestriction(μ, 2, 4)) - end - - ## add nonlinear reaction part (= -df/dc times test function) - function kernel_dfdc!(result, input, qpinfo) - return result[1] = -dfdc(input[1]) - end - assign_operator!(PD, NonlinearOperator(kernel_dfdc!, [id(μ)], [id(c)]; bonus_quadorder = 1)) - - ## generate FESpace and solution vector and interpolate initial state - FES = FESpace{H1Pk{1, 2, order}}(xgrid) - sol = FEVector([FES, FES]; tags = PD.unknowns) - interpolate!(sol[c], c0!) - - ## init plot (if order > 1, solution is upscaled to finer grid for plotting) - plt = GridVisualizer(; Plotter = Plotter, layout = (4, 3), clear = true, resolution = (900, 1200)) - if order > 1 - xgrid_upscale = uniform_refine(xgrid, order - 1) - SolutionUpscaled = FEVector(FESpace{H1P1{1}}(xgrid_upscale)) - lazy_interpolate!(SolutionUpscaled[1], sol) - else - xgrid_upscale = xgrid - SolutionUpscaled = sol - end - nodevals = nodevalues_view(SolutionUpscaled[1]) - scalarplot!(plt[1, 1], xgrid_upscale, nodevals[1]; limits = (0.61, 0.65), xlabel = "", ylabel = "", levels = 1, title = "c (t = 0)") - - ## prepare backward Euler time derivative - M = FEMatrix(FES) - b = FEVector(FES) - assemble!(M, BilinearOperator([id(1)]; factor = 1.0 / τ)) - assign_operator!(PD, BilinearOperator(M, [c]; kwargs...)) - assign_operator!(PD, LinearOperator(b, [c]; kwargs...)) - - ## generate solver configuration - SC = SolverConfiguration(PD, [FES, FES]; init = sol, maxiterations = 50, target_residual = 1.0e-6, kwargs...) - - ## advance in time, plot from time to time - t = 0 - for j in 1:11 - ## do some timesteps until next plot - for it in 1:iterations_until_next_plot - t += τ - ## update time derivative - b.entries .= M.entries * view(sol[c]) - ExtendableFEM.solve(PD, [FES, FES], SC; time = t) - end - - ## enlarge time step a little bit - τ *= τ_increase - M.entries.cscmatrix.nzval ./= τ_increase - - ## plot at current time - if order > 1 - lazy_interpolate!(SolutionUpscaled[1], sol) - end - scalarplot!(plt[1 + Int(floor((j) / 3)), 1 + (j) % 3], xgrid_upscale, nodevals[1]; xlabel = "", ylabel = "", limits = (-0.1, 1.1), levels = 1, title = "c (t = $(Float32(t)))") - end - - return sol, plt -end - -generateplots = ExtendableFEM.default_generateplots(Example285_CahnHilliard, "example285.png") #hide -end diff --git a/examples/Example290_PoroElasticity.jl b/examples/Example290_PoroElasticity.jl deleted file mode 100644 index f7727a81..00000000 --- a/examples/Example290_PoroElasticity.jl +++ /dev/null @@ -1,221 +0,0 @@ -#= - -# 290 : Poro-Elasticity -([source code](@__SOURCE_URL__)) - -This example concerns the three-field solution ``(\mathbf{w},\mathbf{u},p)`` of Biot's consolidation -model in poroelasticity given by -```math -\begin{aligned} --(\lambda + \mu) \nabla (\mathrm{div} \mathbf{v}) - \mu \Delta \mathbf{v} + \alpha \nabla p & = f \quad \text{in } \Omega \times [0,T]\\ -\partial_t (c_o + α \mathrm{div}(\mathbf{v})) + \mathrm{div}(w\mathbf{w}) & = g \quad \text{in } \Omega \times [0,T]\\ -K^{-1} \mathbf{w} + \nabla p & = 0 \quad \text{in } \Omega \times [0,T] -\end{aligned} -``` -and suitable boundary conditions and given initial state. - -The discretisation involves an Hdiv-conforming reconstruction operator to avoid -Poisson locking which results in a scheme similar to the one suggested in the reference below. -As a test problem the first benchmark problem from the same reference is used. - -!!! reference - - "A study of two modes of locking in poroelasticity",\ - S.-Y. Yi,\ - SIAM J. Num. Anal. 55(4) (2017),\ - [>Journal-Link<](https://epubs.siam.org/doi/10.1137/16M1056109) - -The computed solution for the default parameters looks like this: - -![](example290.png) - -=# - -module Example290_PoroElasticity - -using ExtendableFEM -using ExtendableGrids -using GridVisualize -using Symbolics -using Test #hide - - -## exact data for testcase 2 computed by Symbolics -function prepare_data!(; μ = 1, λ = 1, K = 1, c0 = 1, α = 1) - - @variables x y t - - ## displacement and pressure - u = [ - exp(-t) * (sin(2 * pi * y) * (-1 + cos(2 * pi * x)) + sin(pi * x) * sin(pi * y) / (μ + λ)) - exp(-t) * (sin(2 * pi * x) * (1 - cos(2 * pi * y)) + sin(pi * x) * sin(pi * y) / (μ + λ)) - ] - p = exp(-t) * sin(pi * x) * sin(pi * y) - - ## gradient of displacement - ∇u = Symbolics.jacobian(u, [x, y]) - ∇u_reshaped = [∇u[1, 1], ∇u[1, 2], ∇u[2, 1], ∇u[2, 2]] - - ## gradient of pressure - ∇p = [Symbolics.gradient(p, [x])[1], Symbolics.gradient(p, [y])[1]] - - ## Laplacian - Δu = [ - (Symbolics.gradient(∇u[1, 1], [x]) + Symbolics.gradient(∇u[1, 2], [y]))[1], - (Symbolics.gradient(∇u[2, 1], [x]) + Symbolics.gradient(∇u[2, 2], [y]))[1], - ] - Δp = Symbolics.gradient(∇p[1], [x]) + Symbolics.gradient(∇p[2], [y]) - divu = ∇u[1, 1] + ∇u[2, 2] - ∇divu = [Symbolics.gradient(divu, [x])[1], Symbolics.gradient(divu, [y])[1]] - divu_dt = Symbolics.gradient(divu, [t]) - - f = -μ * Δu .+ α * ∇p .- (μ + λ) * ∇divu - g = c0 * Symbolics.gradient(p, [t]) - K * Δp + α * divu_dt - - u_eval = build_function(u, x, y, t, expression = Val{false}) - ∇u_eval = build_function(∇u_reshaped, x, y, t, expression = Val{false}) - g_eval = build_function(g, x, y, t, expression = Val{false}) - f_eval = build_function(f, x, y, t, expression = Val{false}) - p_eval = build_function(p, x, y, t, expression = Val{false}) - ∇p_eval = build_function(∇p, x, y, t, expression = Val{false}) - - return f_eval[2], g_eval[2], u_eval[2], ∇u_eval[2], p_eval, ∇p_eval[2] -end - -function linear_kernel!(result, input, qpinfo) - ∇u, divu, p, w, divw = view(input, 1:4), view(input, 5), view(input, 6), view(input, 7:8), view(input, 9) - μ, λ, α, K = qpinfo.params[1], qpinfo.params[2], qpinfo.params[3], qpinfo.params[4] - result[1] = μ * ∇u[1] + (λ + μ) * divu[1] - p[1] - result[2] = μ * ∇u[2] - result[3] = μ * ∇u[3] - result[4] = μ * ∇u[4] + (λ + μ) * divu[1] - p[1] - result[5] = divu[1] - result[6] = divw[1] - result[7] = w[1] / K - result[8] = w[2] / K - result[9] = -p[1] - return nothing -end - -## kernel for exact error calculation -function exact_error!(u!, ∇u!, p!) - return function closure(result, u, qpinfo) - u!(view(result, 1:2), qpinfo) - ∇u!(view(result, 3:6), qpinfo) - p!(view(result, 7), qpinfo) - view(result, 1:7) .-= u - return result .= result .^ 2 - end -end - -function main(; α = 0.93, E = 1.0e5, ν = 0.4, K = 1.0e-7, nrefs = 6, T = 0.5, τ = 1.0e-2, c0 = 1, order = 1, reconstruct = true, Plotter = nothing, kwargs...) - - ## calculate Lame' parameter - μ = E / (2 * (1 + ν)) - λ = E * ν / ((1 - 2 * ν) * (1 + ν)) - - ## initial and exact state for u and p at time t0 - f_eval, g_eval, u_eval, ∇u_eval, p_eval, ∇p_eval = prepare_data!(; μ = μ, λ = λ, K = K, c0 = c0, α = α) - f!(result, qpinfo) = (f_eval(result, qpinfo.x[1], qpinfo.x[2], qpinfo.time)) - g!(result, qpinfo) = (g_eval(result, qpinfo.x[1], qpinfo.x[2], qpinfo.time)) - exact_p!(result, qpinfo) = (result[1] = p_eval(qpinfo.x[1], qpinfo.x[2], qpinfo.time)) - exact_∇p!(result, qpinfo) = (∇p_eval(result, qpinfo.x[1], qpinfo.x[2], qpinfo.time)) - exact_u!(result, qpinfo) = (u_eval(result, qpinfo.x[1], qpinfo.x[2], qpinfo.time)) - exact_∇u!(result, qpinfo) = (∇u_eval(result, qpinfo.x[1], qpinfo.x[2], qpinfo.time)) - - ## problem description - PD = ProblemDescription("Heat Equation") - u = Unknown("u"; name = "displacement") - p = Unknown("p"; name = "pressure") - w = Unknown("w"; name = "Darcy velocity") - assign_unknown!(PD, u) - assign_unknown!(PD, p) - assign_unknown!(PD, w) - - ## prepare reconstruction operator - if reconstruct - FES_Reconst = order == 1 ? HDIVBDM1{2} : HDIVBDM2{2} - divu = apply(u, Reconstruct{FES_Reconst, Divergence}) - idu = apply(u, Reconstruct{FES_Reconst, Identity}) - else - divu = div(u) - idu = id(u) - end - - ## linear operator - assign_operator!(PD, BilinearOperator(linear_kernel!, [grad(u), divu, id(p), id(w), div(w)]; params = [μ, λ, α, K], store = true, kwargs...)) - - ## right-hand side data - assign_operator!(PD, LinearOperator(f!, [idu]; kwargs...)) - assign_operator!(PD, LinearOperator(g!, [id(p)]; kwargs...)) - - ## boundary conditions - assign_operator!(PD, InterpolateBoundaryData(u, exact_u!; regions = 1:4)) - assign_operator!(PD, InterpolateBoundaryData(p, exact_p!; regions = 1:4)) - - ## grid - xgrid = uniform_refine(grid_unitsquare(Triangle2D; scale = [4, 4], shift = [-0.5, -0.5]), nrefs) - - ## prepare solution vector - if order == 1 - FES = [FESpace{H1BR{2}}(xgrid), FESpace{L2P0{1}}(xgrid; broken = true), FESpace{HDIVRT0{2}}(xgrid)] - elseif order == 2 - FES = [FESpace{H1P2B{2, 2}}(xgrid), FESpace{H1P1{1}}(xgrid; broken = true), FESpace{HDIVRT1{2}}(xgrid)] - end - sol = FEVector(FES; tags = PD.unknowns) - - ## initial data - interpolate!(sol[u], exact_u!; bonus_quadorder = 5, time = 0) - interpolate!(sol[p], exact_p!; bonus_quadorder = 5, time = 0) - - ## init plotter and plot initial data and grid - plt = GridVisualizer(; Plotter = Plotter, layout = (3, 2), clear = true, size = (800, 1200)) - scalarplot!(plt[1, 1], id(u), sol; abs = true, title = "u_h (t = 0)") - scalarplot!(plt[2, 1], id(p), sol; title = "p_h (t = 0)") - gridplot!(plt[3, 1], xgrid; linewidth = 1) - - ## compute mass matrix - M = FEMatrix(FES) - assemble!(M, BilinearOperator([id(2)]; factor = c0)) - assemble!(M, BilinearOperator([id(2)], [div(1)]; factor = -α)) - - ## add backward Euler time derivative - assign_operator!(PD, BilinearOperator(M, [u, p, w]; factor = 1 / τ, kwargs...)) - assign_operator!(PD, LinearOperator(M, [u, p, w], [u, p, w]; factor = 1 / τ, kwargs...)) - - ## generate solver configuration - SC = SolverConfiguration(PD, FES; init = sol, maxiterations = 1, verbosity = -1, constant_matrix = true, kwargs...) - - ## iterate tspan - t = 0 - for it in 1:Int(floor(T / τ)) - t += τ - @info "t = $t" - ExtendableFEM.solve(PD, FES, SC; time = t, timeroutputs = :hide) - end - - ## error calculation - ErrorIntegrator = ItemIntegrator(exact_error!(exact_u!, exact_∇u!, exact_p!), [id(u), grad(u), id(p)]; quadorder = 2 * (order + 1), kwargs...) - error = evaluate(ErrorIntegrator, sol; time = T) - L2errorU = sqrt(sum(view(error, 1, :)) + sum(view(error, 2, :))) - H1errorU = sqrt(sum(view(error, 3, :)) + sum(view(error, 4, :)) + sum(view(error, 5, :)) + sum(view(error, 6, :))) - L2errorP = sqrt(sum(view(error, 7, :))) - @info "|| u - u_h || = $L2errorU - || ∇(u - u_h) || = $H1errorU - || p - p_h || = $L2errorP" - - ## plot final state - scalarplot!(plt[1, 2], id(u), sol; abs = true, title = "u_h (t = $T)") - scalarplot!(plt[2, 2], id(p), sol; title = "p_h (t = $T)") - scalarplot!(plt[3, 2], id(w), sol; abs = true, title = "|w_h| (t = $T)") - - return L2errorU, plt -end - -generateplots = ExtendableFEM.default_generateplots(Example290_PoroElasticity, "example290.png") #hide -function runtests() #hide - L2errorU, plt = main(; nrefs = 4) #hide - @test L2errorU ≈ 0.18232484430836826 #hide - return nothing #hide -end #hide -end # module diff --git a/examples/Example295_SlidingDroplet.jl b/examples/Example295_SlidingDroplet.jl deleted file mode 100644 index 9e36ec0e..00000000 --- a/examples/Example295_SlidingDroplet.jl +++ /dev/null @@ -1,239 +0,0 @@ -#= - -# 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 Navier-slip and capillary boundary conditions, -see the reference below for details. - -The ALE method computes a postprocessed displacement -``\mathbf{w}`` from ``\mathbf{u}`` by Laplace smoothing to prevent mesh folding. The smoothing -preserves the normal fluxes via a Lagrange multiplier. - -The weak formulation characterizes ``(\mathbf{u},p,\mathbf{w},\lambda) \in V \times Q \times W \times \Lambda`` by -```math -\begin{aligned} - (2\mu \varepsilon(\mathbf{u}), \varepsilon(\mathbf{v})) - (\mathrm{div} \mathbf{v}, p)\\ - + \beta (\mathbf{u}, \mathbf{v})_{\Gamma_{solid}} - + \delta (\mathbf{u}, \mathbf{v})_{\Gamma_{contact}} - & = \mathrm{Bo} (\mathbf{e}_x, \mathbf{v}) - (\nabla_s \mathbf{x}, \nabla_s \mathbf{v})_{\Gamma_{air}} - + \cos(\theta) (\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}) + (\mathbf{z} \cdot \mathbf{n}, \lambda)_{\partial\Omega} - & = 0 - && \quad \text{for all } \mathbf{z} \in W,\\ - (\mathbf{w} \cdot \mathbf{n}, \psi)_{\partial\Omega} - & = (\mathbf{u} \cdot \mathbf{n}, \psi)_{\partial\Omega} - && \quad \text{for all } \mathbf{\psi} \in \Lambda -\end{aligned} -``` -Here, ``Bo``, ``\beta`` and ``\delta`` are dimensionless parameters -and ``\theta`` is the equilibration contact angle, see the reference below for details. - -When the droplet speed, i.e. the mean velocity in x-direction, -becomes constant in time, a travelling wave 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, ## polynomial FEM order - β = 1.0, ## friction at liquid-solid interface - δ = 1.0, ## friction at contact line - Bo = 0.5, ## Bond number - θ = π / 2, ## equilibrium contact angle - nsteps = 800, ## number of ALE steps - τ = 0.005, ## ALE stepsize - nrefs = 4, ## mesh refinement level - stationarity_target = 1.0e-2, - Plotter = nothing, kwargs... - ) - - @info "δ = $δ, β = $(β)" - - ## 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, BilinearOperator([id(p)], [div(u)]; transposed_copy = 1, factor = -1, kwargs...)) - assign_operator!(PD, BilinearOperator([id(u)]; entities = ON_BFACES, regions = [2], factor = β, kwargs...)) - if abs(δ) > 0 - assign_operator!(PD, CallbackOperator(triple_junction_kernel!(triple_nodes, δ), [u]; kwargs...)) - end - - ## RHS - assign_operator!(PD, LinearOperator(g!, [id(u)]; factor = Bo, kwargs...)) - assign_operator!(PD, LinearOperatorDG(surface_tension!, [grad(u)], [grad(x)]; entities = ON_BFACES, quadorder = 2, factor = -1, regions = [1], kwargs...)) - if abs(cos(θ)) > 1.0e-12 - assign_operator!(PD, LinearOperatorDG(surface_tension!, [grad(u)], [grad(x)]; entities = ON_BFACES, quadorder = 2, factor = cos(θ), regions = [2], kwargs...)) - end - - ## 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 mean velocity integration - vintegrate = ItemIntegrator([id(u, 1)]; piecewise = false) - - ## time loop - time = 0.0 - mass = sum(xgrid[CellVolumes]) - v0, vprev = 0.0, 0.0 - for step in 1:nsteps - time += τ - - ## redefine x for computation of the tangential identity grad(x) - view(sol[x]) .= view(xgrid[Coordinates]', :) - - ## solve Stokes problem - solve(PD, FES[[1, 2]], SC; kwargs...) - - ## solve ALE problem - solve(PDALE, FES[[1, 3]], SCALE; kwargs...) - - ## displace mesh - displace_mesh!(xgrid, sol[w]; magnify = τ) - - ## calculate droplet speed - vprev = v0 - v0 = evaluate(vintegrate, sol)[1] / mass - - @info "STEP $step ------------- - time = $(Float16(time)) - droplet_speed = $(Float32(v0)) - stationarity = $(step > 1 ? Float32((vprev - v0) / τ) : String("init"))" - - if step > 1 && (vprev - v0 < τ * stationarity_target) - @info "detected stationarity" - break - elseif step == nsteps - @info "maximum number of ALE steps reached" - end - 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 = 21) - - 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 diff --git a/examples/Example301_PoissonProblem.jl b/examples/Example301_PoissonProblem.jl deleted file mode 100644 index 78cdde64..00000000 --- a/examples/Example301_PoissonProblem.jl +++ /dev/null @@ -1,113 +0,0 @@ -#= - -# 301 : Poisson-Problem -([source code](@__SOURCE_URL__)) - -This example computes the solution ``u`` of the two-dimensional Poisson problem -```math -\begin{aligned} --\Delta u & = f \quad \text{in } \Omega -\end{aligned} -``` -with right-hand side ``f(x,y) \equiv xy`` and homogeneous Dirichlet boundary conditions -on the unit cube domain ``\Omega`` on a given grid. The computed solution for the default -parameters looks like this: - -![](example301.png) - -This examples uses an iterative solver with an IncompleteLU preconditioner as the default solver. -It can be changed via the arguments 'method_linear' and 'precon_linear', see the runtests function -for more examples. - -=# - -module Example301_PoissonProblem - -using ExtendableFEM -using ExtendableGrids -using LinearSolve -using IncompleteLU -using LinearAlgebra -using Test - -function f!(result, qpinfo) - result[1] = qpinfo.params[1] * (1.7^2 + 3.9^2) * sin(1.7 * qpinfo.x[1]) * cos(3.9 * qpinfo.x[2]) - return nothing -end - -function u!(result, qpinfo) - result[1] = sin(1.7 * qpinfo.x[1]) * cos(3.9 * qpinfo.x[2]) - return nothing -end - -function main(; - μ = 1.0, - nrefs = 4, - method_linear = KrylovJL_GMRES(), - precon_linear = method_linear == KrylovJL_GMRES() ? IncompleteLU.ilu : nothing, - Plotter = nothing, - kwargs... - ) - - ## problem description - PD = ProblemDescription() - u = Unknown("u"; name = "potential") - assign_unknown!(PD, u) - assign_operator!(PD, BilinearOperator([grad(u)]; factor = μ)) - assign_operator!(PD, LinearOperator(f!, [id(u)]; params = [μ])) - assign_operator!(PD, InterpolateBoundaryData(u, u!; regions = 1:6)) - - ## discretize - xgrid = uniform_refine(grid_unitcube(Tetrahedron3D), nrefs) - FES = FESpace{H1P2{1, 3}}(xgrid) - - ## solve - sol = solve(PD, FES; method_linear, precon_linear, kwargs...) - - ## compute error - function exact_error!(result, u, qpinfo) - u!(result, qpinfo) - result .-= u - result .= result .^ 2 - return nothing - end - ErrorIntegratorExact = ItemIntegrator(exact_error!, [id(u)]; quadorder = 8) - - ## calculate error - error = evaluate(ErrorIntegratorExact, sol) - L2error = sqrt(sum(view(error, 1, :))) - @info "L2 error = $L2error" - - ## plot - plt = plot([id(u)], sol; Plotter = Plotter) - - return L2error, plt -end - -generateplots = ExtendableFEM.default_generateplots(Example301_PoissonProblem, "example301.png") -function runtests() - expected_error = 8.56e-5 - - ## test direct solver - L2error, plt = main(; nrefs = 4) - @test L2error <= expected_error - - ## test iterative solver with IncompleteLU (fastest) - method_linear = KrylovJL_GMRES() - precon_linear = IncompleteLU.ilu - L2error, plt = main(; method_linear, precon_linear, nrefs = 4) - @test L2error <= expected_error - - ## test other iterative solvers - method_linear = KrylovJL_GMRES(precs = (A, p) -> (Diagonal(A), I)) - precon_linear = nothing - L2error, plt = main(; method_linear, precon_linear, nrefs = 4) - @test L2error <= expected_error - - ## also working: - ## method_linear = KrylovJL_GMRES(precs = (A, p) -> (AMGCLWrap.AMGPrecon(ruge_stuben(A)), I)) - ## method_linear = KrylovJL_GMRES(precs = (A, p) -> (AlgebraicMultigrid.aspreconditioner(ruge_stuben(A)), I)) - - return nothing -end -end # module diff --git a/examples/Example310_DivFreeBasis.jl b/examples/Example310_DivFreeBasis.jl deleted file mode 100644 index 63fa942f..00000000 --- a/examples/Example310_DivFreeBasis.jl +++ /dev/null @@ -1,204 +0,0 @@ -#= - -# 310 : Div-free RT0 basis -([source code](@__SOURCE_URL__)) - -This example computes the best-approximation ``\mathbf{\psi}_h`` of a divergence-free velocity -``\mathbf{u} = \mathrm{curl} \mathbf{\psi}`` by solving -for a curl-potential ``\mathbf{\phi}_h \in N_0`` with -```math -\begin{aligned} -(\mathrm{curl} \mathbf{\phi}_h, \mathrm{curl} \mathbf{\theta}_h) & = (\mathbf{u}, \mathrm{curl} \mathbf{\theta}_h) \quad \text{for all } \mathbf{\theta} \in N_0 -\end{aligned} -``` -Here, ``N_0`` denotes the lowest-order Nedelec space which renders the problem ill-posed unless one selects -a linear independent basis. This is done with the algorithm suggested in the reference below. - -!!! reference - - "Decoupling three-dimensional mixed problems using divergence-free finite elements",\ - R. Scheichl,\ - SIAM J. Sci. Comput. 23(5) (2002),\ - [>Journal-Link<](http://www.siam.org/journals/sisc/23-5/37588.html) - -The computed solution for the default parameters looks like this: - -![](example310.png) - -=# - -module Example310_DivFreeBasis - -using ExtendableFEM -using GridVisualize -using ExtendableGrids -using ExtendableSparse -using LinearAlgebra -using Symbolics -using Test #hide - -## exact data for problem generated by symbolics -function prepare_data() - - @variables x y z - - ## stream function ξ - ξ = [x * y * z, x * y * z, x * y * z] - - ## velocity u = curl ξ - ∇ξ = Symbolics.jacobian(ξ, [x, y, z]) - u = [∇ξ[3, 2] - ∇ξ[2, 3], ∇ξ[1, 3] - ∇ξ[3, 1], ∇ξ[2, 1] - ∇ξ[1, 2]] - - ## build function - u_eval = build_function(u, x, y, z, expression = Val{false}) - - return u_eval[2] -end - -function main(; - nrefs = 4, ## number of refinement levels - bonus_quadorder = 2, ## additional quadrature order for data evaluations - divfree_basis = true, ## if true uses curl(N0), if false uses mixed FEM RT0xP0 - Plotter = nothing, ## Plotter (e.g. PyPlot) - kwargs... - ) - - ## prepare problem data - u_eval = prepare_data() - exact_u!(result, qpinfo) = (u_eval(result, qpinfo.x[1], qpinfo.x[2], qpinfo.x[3])) - - ## prepare plots - plt = GridVisualizer(; Plotter = Plotter, layout = (2, 2), clear = true, size = (800, 800)) - - ## prepare error calculation - function exact_error!(result, u, qpinfo) - exact_u!(view(result, 1:3), qpinfo) - result .-= u - return result .= result .^ 2 - end - ErrorIntegratorExact = ItemIntegrator(exact_error!, [divfree_basis ? curl3(1) : id(1)]; bonus_quadorder = 2 + bonus_quadorder, kwargs...) - NDofs = zeros(Int, nrefs) - L2error = zeros(Float64, nrefs) - - sol = nothing - for lvl in 1:nrefs - ## grid - xgrid = uniform_refine(grid_unitcube(Tetrahedron3D), lvl) - - if divfree_basis - - ## use Nedelec FESpace and determine linear independent basis - FES = FESpace{HCURLN0{3}}(xgrid) - - @time begin - ## get subset of edges, spanning the node graph - spanning_tree = get_spanning_edge_subset(xgrid) - - ## get all other edges = linear independent degrees of freedom - subset = setdiff(1:num_edges(xgrid), spanning_tree) - end - NDofs[lvl] = length(subset) - - ## assemble full Nedelec curl-curl problem... - u = Unknown("u"; name = "curl potential of velocity") - PD = ProblemDescription("curl-curl formulation") - assign_unknown!(PD, u) - assign_operator!(PD, BilinearOperator([curl3(u)])) - assign_operator!(PD, LinearOperator(exact_u!, [curl3(u)]; bonus_quadorder = bonus_quadorder)) - - ## ...and solve with subset - sol = solve(PD, FES; restrict_dofs = [subset[:]]) - else - ## use RT0 functions + side constraint for divergence - FES = [FESpace{HDIVRT0{3}}(xgrid), FESpace{L2P0{1}}(xgrid)] - NDofs[lvl] = FES[1].ndofs + FES[2].ndofs - - u = Unknown("u"; name = "velocity") - p = Unknown("u"; name = "pressure") - PD = ProblemDescription("mixed formulation") - assign_unknown!(PD, u) - assign_unknown!(PD, p) - assign_operator!(PD, BilinearOperator([id(u)])) - assign_operator!(PD, BilinearOperator([div(u)], [id(p)]; transposed_copy = 1)) - assign_operator!(PD, LinearOperator(exact_u!, [id(u)]; bonus_quadorder = bonus_quadorder)) - sol = solve(PD, FES) - end - - ## evaluate error - error = evaluate(ErrorIntegratorExact, sol) - L2error[lvl] = sqrt(sum(view(error, 1, :)) + sum(view(error, 2, :))) - if divfree_basis - @info "|| u - curl(ϕ_h) || = $(L2error[lvl])" - else - @info "|| u - u_h || = $(L2error[lvl])" - end - end - - ## plot - if divfree_basis - scalarplot!(plt[1, 1], curl3(1), sol; abs = true) - else - scalarplot!(plt[1, 1], id(1), sol; abs = true) - end - - ## print convergence history as table - print_convergencehistory(NDofs, L2error; X_to_h = X -> X .^ (-1 / 3), ylabels = ["|| u - u_h ||"], xlabel = "ndof") - - return L2error, plt -end - - -## finds a minimal subset (of dimension #nodes - 1) of edges, such that all nodes are connected -function get_spanning_edge_subset(xgrid) - nnodes = num_nodes(xgrid) - edgenodes = xgrid[EdgeNodes] - bedgenodes = xgrid[BEdgeNodes] - bedgeedges = xgrid[BEdgeEdges] - - ## boolean arrays to memorize which nodes are visited - ## and which edges belong to the spanning tree - visited = zeros(Bool, nnodes) - markededges = zeros(Bool, num_edges(xgrid)) - - function find_spanning_tree(edgenodes, remap) - nodeedges = atranspose(edgenodes) - function recursive(node) - visited[node] = true - nneighbors = num_targets(nodeedges, node) - for e in 1:nneighbors - edge = nodeedges[e, node] - for k in 1:2 - node2 = edgenodes[k, edge] - if !visited[node2] - ## mark edge - markededges[remap[edge]] = true - recursive(node2) - end - end - end - return nothing - end - return recursive(edgenodes[1]) - end - - ## find spanning tree for Neumann boundary - ## local bedges >> global edge numbers - find_spanning_tree(bedgenodes, bedgeedges) - - ## find spanning tree for remaining part - other_nodes = setdiff(1:nnodes, unique(view(bedgenodes, :))) - if length(other_nodes) > 0 - find_spanning_tree(edgenodes, 1:num_edges(xgrid)) - end - - ## return all marked edges - return findall(==(true), markededges) -end - -generateplots = ExtendableFEM.default_generateplots(Example310_DivFreeBasis, "example310.png") #hide -function runtests() #hide - L2error, plt = main(; nrefs = 2) #hide - @test L2error[2] ≈ 0.06821145277709957 #hide - return nothing #hide -end #hide -end # module diff --git a/examples/Example312_PeriodicElasticity3D.jl b/examples/Example312_PeriodicElasticity3D.jl deleted file mode 100644 index 78cb2228..00000000 --- a/examples/Example312_PeriodicElasticity3D.jl +++ /dev/null @@ -1,218 +0,0 @@ -#= - -# 312 : Periodic Boundary 3D -([source code](@__SOURCE_URL__)) - -This is a simple demonstration and validation of the generic periodic boundary operator. - -We construct an unstructured periodic 3D grid and solve a simple linear elastic problem -with periodic coupling along the x-axis. - -![](example312.png) -=# - -module Example312_PeriodicElasticity3D - -using ExtendableFEM -using ExtendableGrids -using SimplexGridFactory -using GridVisualize -using TetGen -using UnicodePlots -using StaticArrays -using LinearAlgebra -using Test #hide - -const reg_left = 1 -const reg_right = 2 -const reg_dirichlet = 3 -const reg_default = 4 - - -# define the Hooke tensor for AlN (from DOI 10.1063/1.1368156) -function material_tensor() - c11 = 396.0 - c12 = 137.0 - c13 = 108.0 - c33 = 373.0 - c44 = 116.0 - - return @SArray [ - c11 c12 c13 0 0 0 - c12 c11 c13 0 0 0 - c13 c13 c33 0 0 0 - 0 0 0 c44 0 0 - 0 0 0 0 c44 0 - 0 0 0 0 0 c44 - ] -end - -## generate the kernels for the linear problem -## 𝐂: Hooke tensor, 𝑓: body force -function make_kernels(𝐂, 𝑓) - - ## linear stress-strain mapping - bilinear_kernel!(σ, εv, qpinfo) = mul!(σ, 𝐂, εv) - - ## body force - linear_kernel!(result, qpinfo) = (result .= 𝑓) - - return bilinear_kernel!, linear_kernel! -end - - -""" - create 3D grid with Dirichlet boundary region at the bottom center -""" -function create_grid(; h, height, width, depth) - builder = SimplexGridBuilder(; Generator = TetGen) - - ## bottom points - b01 = point!(builder, 0, 0, 0) - b02 = point!(builder, 0.45 * width, 0, 0) - b03 = point!(builder, 0.55 * width, 0, 0) - b04 = point!(builder, width, 0, 0) - - b11 = point!(builder, 0, depth, 0) - b12 = point!(builder, 0.45 * width, depth, 0) - b13 = point!(builder, 0.55 * width, depth, 0) - b14 = point!(builder, width, depth, 0) - - ## top points - t01 = point!(builder, 0, 0, height) - t02 = point!(builder, width, 0, height) - - t11 = point!(builder, 0, depth, height) - t12 = point!(builder, width, depth, height) - - ## center points - c01 = point!(builder, 0.5 * width, 0, 0) - c02 = point!(builder, 0.5 * width, 0, height) - c11 = point!(builder, 0.5 * width, depth, 0) - c12 = point!(builder, 0.5 * width, depth, height) - - ## default faces - facetregion!(builder, reg_default) - facet!(builder, b01, b02, b12, b11) - facet!(builder, b03, b04, b14, b13) - facet!(builder, [t01, c02, c12, t11]) - facet!(builder, [c02, t02, t12, c12]) - facet!(builder, [b01, b02, c01, c02, t01]) - facet!(builder, [c01, b03, b04, t02, c02]) - facet!(builder, [b11, b12, c11, c12, t11]) - facet!(builder, [c11, b13, b14, t12, c12]) - facet!(builder, c01, c02, c12, c11) - - ## left face - facetregion!(builder, reg_left) - facet!(builder, b01, t01, t11, b11) - - ## right face - facetregion!(builder, reg_right) - facet!(builder, b04, t02, t12, b14) - - ## Dirichlet face - facetregion!(builder, reg_dirichlet) - facet!(builder, [b02, c01, c11, b12]) - facet!(builder, [c01, b03, b13, c11]) - - cellregion!(builder, 1) - maxvolume!(builder, h) - regionpoint!(builder, width / 3, depth / 2, height / 2) - cellregion!(builder, 2) - ## finer grid on the right half to make the periodic coupling non-trivial - maxvolume!(builder, 0.3 * h) - regionpoint!(builder, 2 * width / 3, depth / 2, height / 2) - - return simplexgrid(builder) -end - -function main(; - order = 1, - periodic_coupling = :high_level_restriction, # :restriction, :operator, :high_level_restriction - Plotter = nothing, - force = 1.0, - h = 1.0e-4, - width = 6.0, - height = 0.2, - depth = 1, - threads = 1, - kwargs... - ) - ## print options for better logs - @info "selected options" periodic_coupling threads - - xgrid = create_grid(; h, width, height, depth) - - ## create finite element space and solution vector - if order == 1 - FES = FESpace{H1P1{3}}(xgrid) - elseif order == 2 - FES = FESpace{H1P2{3, 3}}(xgrid) - end - - ## problem description - PD = ProblemDescription() - u = Unknown("u"; name = "displacement") - assign_unknown!(PD, u) - - 𝐂 = material_tensor() - 𝑓 = force * [0, 0, 1] - - bilinear_kernel!, linear_kernel! = make_kernels(𝐂, 𝑓) - assign_operator!(PD, BilinearOperator(bilinear_kernel!, [εV(u, 1.0)]; kwargs...)) - assign_operator!(PD, LinearOperator(linear_kernel!, [id(u)]; kwargs...)) - - assign_operator!(PD, HomogeneousBoundaryData(u; regions = [reg_dirichlet])) - - if periodic_coupling == :high_level_restriction - - # new high-level call - assign_restriction!(PD, CoupledDofsRestriction(u, reg_left, reg_right; parallel = threads > 1, threads)) - - elseif periodic_coupling != :none - - function give_opposite!(y, x) - y .= x - y[1] = width - x[1] - return nothing - end - @showtime coupling_matrix = get_periodic_coupling_matrix(FES, reg_left, reg_right, give_opposite!; parallel = threads > 1, threads) - if periodic_coupling == :restriction - assign_restriction!(PD, CoupledDofsRestriction(coupling_matrix)) - else # :operator - assign_operator!(PD, CombineDofs(u, u, coupling_matrix; kwargs...)) - end - end - - ## solve - sol, SC = solve(PD, FES; return_config = true, kwargs...) - residual(SC) < 1.0e-10 || error("Residual is not zero!") - - @info "Lagrange residuals" SC.statistics[:restriction_residuals] - - displace_mesh!(xgrid, sol[u]) - plt = plot([grid(u)], sol; Plotter, do_vector_plots = false, width = 1200, height = 800, title = "displaced mesh", scene3d = :LScene) - - return sol, plt - -end - -generateplots = ExtendableFEM.default_generateplots(Example312_PeriodicElasticity3D, "example312.png") #hide -function runtests() #hide - sol1, _ = main(periodic_coupling = :operator, threads = 1) #hide - @test abs(maximum(view(sol1[1])) - 1.8004602502175202) < 2.0e-3 #hide - - sol2, _ = main(periodic_coupling = :operator, threads = 4) #hide - @test sol1.entries ≈ sol2.entries #hide - - sol3, _ = main(periodic_coupling = :restriction, threads = 4) #hide - @test sol1.entries ≈ sol3.entries #hide - - sol4, _ = main(periodic_coupling = :high_level_restriction, threads = 4) #hide - @test sol1.entries ≈ sol4.entries #hide - - return nothing #hide -end #hide - -end # module diff --git a/examples/Example313_PeriodicPoisson.jl b/examples/Example313_PeriodicPoisson.jl deleted file mode 100644 index a4f8e761..00000000 --- a/examples/Example313_PeriodicPoisson.jl +++ /dev/null @@ -1,148 +0,0 @@ -#= - -# 312 : Periodic Poisson 3D -([source code](@__SOURCE_URL__)) - -This is a simple demonstration and validation of the new restriction based periodic boundary operator. - -An unstructured cube grid is coupled along two axes periodically and restricted to non-zero Dirichlet values at the other to faces. - -The result is a linear function and the error is measured directly. - -Note that the result is independent of the periodic coupling! Therefore the correctness of the -new QR based column compression is covered by this example. - -![](example313.png) -=# - -module Example313_PeriodicPoisson - -using ExtendableFEM -using ExtendableGrids -using SimplexGridFactory -using GridVisualize -using TetGen -using UnicodePlots -using StaticArrays -using LinearAlgebra -using Test #hide - -const reg_left = 1 -const reg_right = 2 -const reg_front = 3 -const reg_back = 4 -const reg_bottom = 5 -const reg_top = 6 - - -function create_grid(; h) - builder = SimplexGridBuilder(; Generator = TetGen) - - ## bottom points - b00 = point!(builder, 0, 0, 0) - b01 = point!(builder, 0, 1, 0) - b10 = point!(builder, 1, 0, 0) - b11 = point!(builder, 1, 1, 0) - - ## top points - t00 = point!(builder, 0, 0, 1) - t01 = point!(builder, 0, 1, 1) - t10 = point!(builder, 1, 0, 1) - t11 = point!(builder, 1, 1, 1) - - ## left face - facetregion!(builder, reg_left) - facet!(builder, b00, b01, t01, t00) - - ## right face - facetregion!(builder, reg_right) - facet!(builder, b10, b11, t11, t10) - - ## front face - facetregion!(builder, reg_front) - facet!(builder, b00, b10, t10, t00) - - ## back face - facetregion!(builder, reg_back) - facet!(builder, b01, b11, t11, t01) - - ## top face - facetregion!(builder, reg_top) - facet!(builder, t00, t10, t11, t01) - - ## bottom face - facetregion!(builder, reg_bottom) - facet!(builder, b00, b10, b11, b01) - - - cellregion!(builder, 1) - maxvolume!(builder, h) - regionpoint!(builder, 0.5, 0.5, 0.5) - - return simplexgrid(builder) -end - -function main(; - Plotter = nothing, - h = 1.0e-3, - periodic = true, - kwargs... - ) - ## print options for better logs - @info "selected options" periodic - - xgrid = create_grid(; h) - - FES = FESpace{H1P1{1}}(xgrid) - - ## problem description - PD = ProblemDescription("Periodic Poisson Problem") - u = Unknown("u"; name = "temperature") - assign_unknown!(PD, u) - - assign_operator!(PD, BilinearOperator([grad(u)]; kwargs...)) - - assign_restriction!(PD, BoundaryDataRestriction(u; value = -1, regions = [reg_bottom])) - assign_restriction!(PD, BoundaryDataRestriction(u; value = +1, regions = [reg_top])) - - if periodic - assign_restriction!(PD, CoupledDofsRestriction(u, reg_left, reg_right)) - assign_restriction!(PD, CoupledDofsRestriction(u, reg_front, reg_back)) - end - - ## solve - sol, SC = solve(PD, FES; return_config = true, kwargs...) - residual(SC) < 1.0e-10 || error("Residual is not zero!") - - @info "Lagrange residuals" SC.statistics[:restriction_residuals] - - function exact_error!(out, u, qpinfo) - # exact solution here is u(x,y,z) = 2z - 1 - val = qpinfo.x[3] * 2 - 1 - out[1] = (val - u[1])^2 - return nothing - end - - plt = plot([grid(u), id(u)], sol; Plotter, width = 1200, height = 800, scene3d = :LScene, slice = :y => 0.5) - - ErrorIntegrator = ItemIntegrator(exact_error!, [id(u)]; quadorder = 2) - L2error = sqrt(sum(evaluate(ErrorIntegrator, sol))) - - @show L2error - - return L2error, plt - -end - -generateplots = ExtendableFEM.default_generateplots(Example313_PeriodicPoisson, "example313.png") #hide -function runtests() #hide - error1, _ = main(periodic = true) #hide - @test error1 < 1.0e-12 #hide - - error2, _ = main(periodic = false) #hide - @test error2 < 1.0e-12 #hide - - return nothing #hide -end #hide - -end # module diff --git a/examples/Example330_HyperElasticity.jl b/examples/Example330_HyperElasticity.jl deleted file mode 100644 index 160ad8c1..00000000 --- a/examples/Example330_HyperElasticity.jl +++ /dev/null @@ -1,153 +0,0 @@ -#= - -# 330 : Hyperelasticity -([source code](@__SOURCE_URL__)) - -This examples computes the solution of a nonlinear elasticity problem for hyperelastic media -via minimisation of the (neo-Hookian) energy functional -```math -\begin{aligned} - W(u, F(\mathbf{u})) := \int_\Omega \frac{\mu}{2} (F:F - 3 - 2\log(\mathrm{det}(F))) + \frac{\lambda}{2} \log(\mathrm{det}(F))^2 - B \cdot \mathbf{u} \textit{dx} - \int_{\partial \Omega} T \cdot \mathbf{u} \textit{ds} -\end{aligned} -``` -where ``F(\mathbf{u}) := I + \nabla u`` is the deformation gradient and ``\mu`` and ``\lambda`` are the Lame parameters. -The energy is differentiated twice by automatic differentiation to setup a Newton scheme for a -Lagrangian finite element approximation of ``\mathbf{u}``, once in the code below to define the kernel for the NonlinearOperator -and this kernel is differentiated again in the assembly of the Newton scheme for the nonlinear operator. - -The deformed unit cube and the displacement for the default parameters and inhomogeneous boundary conditions as defined in the code -looks like this: - -![](example330.png) - -=# - -module Example330_HyperElasticity - -using ExtendableFEM -using DifferentiationInterface -using ForwardDiff -using LinearAlgebra -using SimplexGridFactory -using TetGen - -## inhomogeneous boundary conditions for bregion 1 -function bnd_1!(result, qpinfo) - x, y, z = qpinfo.x[1], qpinfo.x[2], qpinfo.x[3] - angle = pi / 3 - result[1] = 0.0 - result[2] = (0.5 + (y - 0.5) * cos(angle) - (z - 0.5) * sin(angle) - y) / 2.0 - return result[3] = (0.5 + (y - 0.5) * sin(angle) + (z - 0.5) * cos(angle) - x) / 2.0 -end - -## kernel for body and traction forces -function apply_force!(result, qpinfo) - return result .= qpinfo.params[1] -end - -## energy functional (only nonlinear part, without exterior forces) -function W!(result, F, qpinfo) - F[1] += 1 - F[5] += 1 - F[9] += 1 - μ, λ = qpinfo.params[1], qpinfo.params[2] - detF = -(F[3] * (F[5] * F[7] - F[4] * F[8]) + F[2] * ((-F[6]) * F[7] + F[4] * F[9]) + F[1] * (F[6] * F[8] - F[5] * F[9])) - result[1] = μ / 2 * (dot(F, F) - 3 - 2 * log(detF)) + λ / 2 * (log(detF))^2 - return nothing -end - -## derivative of energy functional (by ForwardDiff) -function nonlinkernel_DW!() - ## A dictionary with cached data for each input type - jac_prep = Dict{DataType, Any}() - scalar_dummy = Dict{DataType, Any}() - AD_backend = AutoForwardDiff() - - ## generate input -> output function for a given quadrature point - make_W(qpinfo) = (out, in) -> W!(out, in, qpinfo) - - return function closure(result, input::Vector{T}, qpinfo) where {T} - W = make_W(qpinfo) - if !haskey(jac_prep, T) - ## first initialization of cached data when type of input = F is known - scalar_dummy[T] = zeros(eltype(input), 1) - jac_prep[T] = prepare_jacobian(W, scalar_dummy[T], AD_backend, input) - end - jacobian!(W, scalar_dummy[T], result, jac_prep[T], AD_backend, input) - return nothing - end -end - - -function main(; - maxvolume = 0.001, # parameter for grid generator - E = 10, # Young modulus - ν = 0.3, # Poisson ratio - order = 3, # finite element order - B = [0, -0.5, 0], # body force - T = [0.1, 0, 0], # traction force - Plotter = nothing, - kwargs... - ) - - ## compute Lame parameters - μ = E / (2 * (1 + ν)) - λ = E * ν / ((1 + ν) * (1 - 2 * ν)) - - ## define unknowns - u = Unknown("u"; name = "displacement") - - ## define problem - PD = ProblemDescription("Hyperelasticity problem") - assign_unknown!(PD, u) - assign_operator!(PD, NonlinearOperator(nonlinkernel_DW!(), [grad(u)]; sparse_jacobians = false, params = [μ, λ])) - assign_operator!(PD, LinearOperator(apply_force!, [id(u)]; params = [B])) - assign_operator!(PD, LinearOperator(apply_force!, [id(u)]; entities = ON_BFACES, store = true, regions = 1:6, params = [T])) - assign_operator!(PD, HomogeneousBoundaryData(u; regions = [2])) - assign_operator!(PD, InterpolateBoundaryData(u, bnd_1!; regions = [1], quadorder = 10)) - - ## grid - xgrid = tetrahedralization_of_cube(maxvolume = maxvolume) - - ## solve - FES = FESpace{H1P1{3}}(xgrid) - sol = solve(PD, FES; maxiterations = 20) - - ## displace mesh and plot final result - displace_mesh!(xgrid, sol[u]) - plt = plot([grid(u), id(u)], sol; Plotter = Plotter, do_vector_plots = false) - - return sol, plt -end - -generateplots = ExtendableFEM.default_generateplots(Example330_HyperElasticity, "example330.png") #hide - -function tetrahedralization_of_cube(; maxvolume = 0.1) - builder = SimplexGridBuilder(; Generator = TetGen) - - p1 = point!(builder, 0, 0, 0) - p2 = point!(builder, 1, 0, 0) - p3 = point!(builder, 1, 1, 0) - p4 = point!(builder, 0, 1, 0) - p5 = point!(builder, 0, 0, 1) - p6 = point!(builder, 1, 0, 1) - p7 = point!(builder, 1, 1, 1) - p8 = point!(builder, 0, 1, 1) - - facetregion!(builder, 1) - facet!(builder, p1, p2, p3, p4) - facetregion!(builder, 2) - facet!(builder, p5, p6, p7, p8) - facetregion!(builder, 3) - facet!(builder, p1, p2, p6, p5) - facetregion!(builder, 4) - facet!(builder, p2, p3, p7, p6) - facetregion!(builder, 5) - facet!(builder, p3, p4, p8, p7) - facetregion!(builder, 6) - facet!(builder, p4, p1, p5, p8) - - return simplexgrid(builder; maxvolume = maxvolume) -end - -end # module diff --git a/src/ExtendableFEM.jl b/src/ExtendableFEM.jl index 0febdb08..3350f1f8 100644 --- a/src/ExtendableFEM.jl +++ b/src/ExtendableFEM.jl @@ -36,7 +36,7 @@ using ExtendableFEMBase: ExtendableFEMBase, AbstractFiniteElement, initialize!, integrate, integrate!, integrate_segment!, lazy_interpolate!, nodevalues, nodevalues!, nodevalues_subset!, nodevalues_view, - norms, unicode_gridplot, unicode_scalarplot, + norms, broken_scalarplot!, update_basis!, SymmetricGradient using ExtendableGrids: ExtendableGrids, AT_NODES, AbstractElementGeometry, Adjacency, AssemblyType, BEdgeNodes, BFaceFaces, BFaceNormals, @@ -95,14 +95,12 @@ export interpolate!, lazy_interpolate! export PointEvaluator, evaluate, evaluate!, evaluate_bary!, eval_func, eval_func_bary export SegmentIntegrator, integrate_segment!, initialize! export integrate!, integrate, QuadratureRule -export unicode_gridplot, unicode_scalarplot export CellDofs, BFaceDofs, FaceDofs, EdgeDofs, BEdgeDofs # export get_polynomialorder export displace_mesh, displace_mesh! export Reconstruct, Identity, Divergence, Gradient # export _addnz export addblock!, addblock_matmul! -export unicode_gridplot, unicode_scalarplot ## reexport stuff from GridVisualize export reveal, save diff --git a/src/plots.jl b/src/plots.jl index 4ae3fc52..af702604 100644 --- a/src/plots.jl +++ b/src/plots.jl @@ -34,14 +34,9 @@ function plot!( p::GridVisualizer, ops, sol; - rasterpoints = 10, - linewidth = 1, keep = [], ncols = size(p.subplots, 2), - do_abs = true, - do_vector_plots = true, - title_add = "", - average_broken_plots = false, + reveal = true, kwargs... ) col, row, id = 0, 1, 0 @@ -58,81 +53,79 @@ function plot!( col, row = 1, row + 1 end end - if op[2] == "grid" - 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...) + if typeof(op[1]) <: Unknown + uname = op[2] == Identity ? String(op[1].identifier) : "$(op[2])(" * String(op[1].identifier) * ")" else - ncomponents = get_ncomponents(sol[op[1]]) - edim = size(sol[op[1]].FES.xgrid[Coordinates], 1) - resultdim = Length4Operator(op[2], edim, ncomponents) - if typeof(op[1]) <: Unknown - title = op[2] == Identity ? String(op[1].identifier) : "$(op[2])(" * String(op[1].identifier) * ")" + uname = op[2] == Identity ? "$(sol[op[1]].name)" : "$(op[2])($(sol[op[1]].name))" + end + plot!(p[row, col], op, sol; reveal = false, kwargs...) + end + if reveal + GridVisualize.reveal(p) + end + return p +end + + +function plot!( + psub, + op, + sol; + rasterpoints = 10, + linewidth = 1, + do_abs = true, + do_vector_plots = true, + title_add = "", + reveal = true, + average_broken_plots = false, + kwargs... + ) + if typeof(op[1]) <: Unknown + uname = op[2] == Identity ? String(op[1].identifier) : "$(op[2])(" * String(op[1].identifier) * ")" + else + uname = op[2] == Identity ? "$(sol[op[1]].name)" : "$(op[2])($(sol[op[1]].name))" + end + if op[2] == "grid" + gridplot!(psub, sol[op[1]].FES.xgrid; linewidth = linewidth, title = uname * title_add, kwargs...) + elseif op[2] == "dofgrid" + gridplot!(psub, sol[op[1]].FES.dofgrid; linewidth = linewidth, title = uname * title_add, kwargs...) + elseif op[2] == "streamlines" + PE = PointEvaluator([apply(op[1], Identity)], sol) + streamplot!(psub, sol[op[1]].FES.dofgrid, eval_func_bary(PE); rasterpoints = rasterpoints, title = uname * title_add, kwargs...) + else + ncomponents = get_ncomponents(sol[op[1]]) + edim = size(sol[op[1]].FES.xgrid[Coordinates], 1) + resultdim = Length4Operator(op[2], edim, ncomponents) + if resultdim == 1 + if !average_broken_plots && ExtendableFEMBase.broken(sol[op[1]].FES) + ExtendableFEMBase.broken_scalarplot!(psub, sol[op[1]], op[2]; title = uname * title_add, kwargs...) else - title = op[2] == Identity ? "$(sol[op[1]].name)" : "$(op[2])($(sol[op[1]].name))" + GridVisualize.scalarplot!(psub, sol[op[1]].FES.dofgrid, view(nodevalues(sol[op[1]], op[2]; abs = false), 1, :), title = uname * title_add; kwargs...) end - if resultdim == 1 - if !average_broken_plots && ExtendableFEMBase.broken(sol[op[1]].FES) - broken_scalarplot!(p[row, col], sol[op[1]], op[2]; title = title * title_add, kwargs...) - else - GridVisualize.scalarplot!(p[row, col], sol[op[1]].FES.dofgrid, view(nodevalues(sol[op[1]], op[2]; abs = false), 1, :), title = title * title_add; kwargs...) - end - elseif do_abs == true - GridVisualize.scalarplot!(p[row, col], sol[op[1]].FES.dofgrid, view(nodevalues(sol[op[1]], op[2]; abs = true), 1, :), title = "|" * title * "|" * title_add; kwargs...) - else - nv = nodevalues(sol[op[1]], op[2]; abs = false) - for k in 1:resultdim - if k > 1 - col += 1 - if col == ncols + 1 - col, row = 1, row + 1 - end + elseif do_abs == true + GridVisualize.scalarplot!(psub, sol[op[1]].FES.dofgrid, view(nodevalues(sol[op[1]], op[2]; abs = true), 1, :), title = "|" * uname * "|" * title_add; kwargs...) + else + nv = nodevalues(sol[op[1]], op[2]; abs = false) + for k in 1:resultdim + if k > 1 + col += 1 + if col == ncols + 1 + col, row = 1, row + 1 end - GridVisualize.scalarplot!(p[row, col], sol[op[1]].FES.dofgrid, view(nv, k, :), title = title * " (component $k)" * title_add, kwargs...) end + GridVisualize.scalarplot!(psub, sol[op[1]].FES.dofgrid, view(nv, k, :), title = uname * " (component $k)" * title_add, kwargs...) end - if resultdim > 1 && do_vector_plots && do_abs == true && edim > 1 - GridVisualize.vectorplot!(p[row, col], sol[op[1]].FES.dofgrid, eval_func_bary(PointEvaluator([op], sol)); rasterpoints = rasterpoints, title = "|" * title * "|" * " + quiver" * title_add, clear = false, kwargs...) - end + end + if resultdim > 1 && do_vector_plots && do_abs == true && edim > 1 + GridVisualize.vectorplot!(psub, sol[op[1]].FES.dofgrid, eval_func_bary(PointEvaluator([op], sol)); title = "|" * uname * "|" * " + quiver" * title_add, clear = false, kwargs...) end end - return p -end - - -""" - broken_scalarplot!(vis, feVectorBlock::FEVectorBlock, operator = Identity; kwargs...) - -A "broken" scalarplot of a broken finite element vector. -Instead of averaging the discontinuous values on the grid nodes, each grid cell is plotted -independently. Thus, a discontinuous plot is generated. - -All kwargs of the calling method are transferred to the scalarplot in this method. -""" -function broken_scalarplot!(vis, feVectorBlock::FEVectorBlock, operator = Identity; kwargs...) - - dofgrid = feVectorBlock.FES.dofgrid - cell_nodes = dofgrid[CellNodes] - coords = dofgrid[Coordinates] - - all_values = nodevalues(feVectorBlock, operator; cellwise = true) # cellwise evaluation of the FE - all_coords = @views coords[:, cell_nodes[:]] - all_cells = reshape(1:length(all_values), size(all_values)) - - GridVisualize.scalarplot!(vis, simplexgrid(all_coords, all_cells, dofgrid[CellRegions]), view(all_values, :); kwargs...) - - return nothing + if reveal + GridVisualize.reveal(psub) + end + return psub end - """ ```` function plot!(p::GridVisualizer, ops, sol; Plotter = nothing, kwargs...) @@ -161,24 +154,9 @@ function plot(ops, sol; add = 0, Plotter = nothing, ncols = min(2, length(ops) + return plot!(p, ops, sol; do_abs = do_abs, kwargs...) end -""" -```` -function plot_unicode(sol; kwargs...) -```` - -Plots all blocks of sol into stdout -(via plot_scalarplot from the UnicodePlots extension of ExtendableFEMBase) - -""" -function plot_unicode(sol; kwargs...) - for u in 1:length(sol) - println(stdout, unicode_scalarplot(sol[u]; title = sol[u].name, kwargs...)) - end - return -end function GridVisualize.vectorplot!(p, xgrid, op::Tuple{Union{Unknown, Int}, DataType}, sol; title = sol[op[1]].name, kwargs...) - return GridVisualize.vectorplot!(p, xgrid, eval_func(PointEvaluator([op], sol)); title = title, kwargs...) + return GridVisualize.vectorplot!(p, xgrid, eval_func_bary(PointEvaluator([op], sol)); title = title, kwargs...) end @@ -256,7 +234,7 @@ function plot_convergencehistory!( end for p in add_h_powers label = "h^$p" - 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...) + scalarplot!(target, simplexgrid(X), X_to_h(X) .^ p; linestyle = :dot, xlabel = xlabel, ylabel = ylabel, color = :black, clear = false, markershape = :none, xscale = :log, yscale = :log, label = label, legend = legend, title = title, args...) end return end @@ -288,3 +266,36 @@ function default_generateplots(example_module, filename; kwargs...) return GridVisualize.save(joinpath(dir, filename), scene; Plotter = Plotter) end end + + +### deprecated +""" +```` +function plot_unicode(sol; kwargs...) +```` + +Plots all blocks of sol into stdout +(via plot_scalarplot from the UnicodePlots extension of ExtendableFEMBase) + +""" +function plot_unicode(sol; kwargs...) + Base.depwarn( + "plot_unicode(sol) will be deprecated in a future release. Use plot(sol; Plotter = UnicodePlots) instead.", + :plot_unicode; force = true + ) + plot(sol; Plotter = UnicodePlots, kwargs...) + return +end + +""" +```` +function plot(sol::FEVector; kwargs...) +```` + +Plots all blocks of sol. + +""" +function plot(sol::FEVector; kwargs...) + plot([id(j) for j in 1:length(sol)], sol; kwargs...) + return nothing +end diff --git a/src/solvers.jl b/src/solvers.jl index b225cfe3..3eb1ffb7 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -787,9 +787,7 @@ function CommonSolve.solve(PD::ProblemDescription, FES::Union{<:FESpace, Vector{ end if SC.parameters[:plot] - for u in unknowns - println(stdout, unicode_scalarplot(sol[u]; title = u.name, kwargs...)) - end + plot([id(u) for u in unknowns], sol; Plotter = UnicodePlots, ncols = length(unknowns)) end # Print the timings in the default way diff --git a/test/runtests.jl b/test/runtests.jl index 79f3bf91..ae1952ba 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,22 +1,5 @@ using Test -using ExtendableGrids -using ExtendableFEMBase -using ExtendableFEM using ExampleJuggler -using ExplicitImports -using Metis -using Aqua -using Triangulate -using SimplexGridFactory -using Symbolics - - -include("test_dgblf.jl") -include("test_boundary_operator.jl") -include("test_itemintegrator.jl") -include("test_timedependence.jl") -include("test_nonlinear_operator.jl") -include("test_helper_functions.jl") function run_examples() ExampleJuggler.verbose!(true) @@ -73,5 +56,5 @@ function run_all_tests() return nothing end -run_all_tests() +# run_all_tests() run_examples()