Modify NSE example to include iterative GMRES and ILU preconditioning

Author: pjaap

CHANGELOG.md

@@ -4,6 +4,11 @@
 
 ## v1.0.0 April 7, 2025
 
+### Added
+
+  - example `Example250` now also demonstrates a nonlinear problem solved by an iterative linear solver
+    with preconditioning and an initial solution.
+
 ### Changed
 
   - `solve` uses now the residual equation for the linear systems.

Project.toml

@@ -33,6 +33,7 @@ ExtendableGrids = "1.10.3"
 ExtendableSparse = "1.5.3"
 ForwardDiff = "0.10.35,1"
 GridVisualize = "1.8.1"
+IncompleteLU = "0.2.1"
 LinearAlgebra = "1.9"
 LinearSolve = "2, 3"
 Metis = "1.5.0"
@@ -55,6 +56,7 @@ julia = "1.9"
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 ExampleJuggler = "3bbe58f8-ed81-4c4e-a134-03e85fcf4a1a"
 ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
+IncompleteLU = "40713840-3770-5561-ab4c-a76e7d0d7895"
 Metis = "2679e427-3c69-5b7f-982b-ece356f1e94b"
 OrdinaryDiffEqRosenbrock = "43230ef6-c299-4910-a778-202eb28ce4ce"
 OrdinaryDiffEqSDIRK = "2d112036-d095-4a1e-ab9a-08536f3ecdbf"
@@ -65,4 +67,4 @@ TetGen = "c5d3f3f7-f850-59f6-8a2e-ffc6dc1317ea"
 Triangulate = "f7e6ffb2-c36d-4f8f-a77e-16e897189344"
 
 [targets]
-test = ["Aqua", "ExampleJuggler", "ExplicitImports", "Metis", "OrdinaryDiffEqRosenbrock", "OrdinaryDiffEqSDIRK", "SimplexGridFactory", "StaticArrays", "Test", "TetGen", "Triangulate"]
+test = ["Aqua", "ExampleJuggler", "ExplicitImports", "IncompleteLU", "Metis", "OrdinaryDiffEqRosenbrock", "OrdinaryDiffEqSDIRK", "SimplexGridFactory", "StaticArrays", "Test", "TetGen", "Triangulate"]

docs/Project.toml

@@ -8,6 +8,7 @@ ExtendableGrids = "cfc395e8-590f-11e8-1f13-43a2532b2fa8"
 ExtendableSparse = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 GridVisualize = "5eed8a63-0fb0-45eb-886d-8d5a387d12b8"
+IncompleteLU = "40713840-3770-5561-ab4c-a76e7d0d7895"
 KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"

examples/Example250_NSELidDrivenCavity.jl

@@ -21,25 +21,62 @@ The computed solution for the default parameters looks like this:
 module Example250_NSELidDrivenCavity
 
 using ExtendableFEM
-using GridVisualize
+using Triangulate
 using ExtendableGrids
+using SimplexGridFactory
 using LinearAlgebra
+using GridVisualize
+using LinearSolve
+using IncompleteLU
 using Test #hide
 
-function kernel_nonlinear!(result, u_ops, qpinfo)
-    u, ∇u, p = view(u_ops, 1:2), view(u_ops, 3:6), view(u_ops, 7)
+
+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]
-    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])
+
+    tmul!(v, ∇u, u)
+    ∇v .= μ .* ∇u .- p[1] .* 𝕀
+    q[1] = -dot(∇u, 𝕀)
+
     return nothing
 end
 
-function boundarydata!(result, qpinfo)
+
+## boundary function
+function u_boundary!(result, qpinfo)
     result[1] = 1
     result[2] = 0
     return nothing
@@ -58,7 +95,16 @@ function initialgrid_cone()
     return xgrid
 end
 
-function main(; μ_final = 0.0005, order = 2, nrefs = 5, Plotter = nothing, kwargs...)
+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
+        method_linear = nothing, # linear solver ("nothing" invokes the default solver)
+        precon_linear = nothing, # preconditioner
+        Plotter = nothing,
+        kwargs...
+    )
 
     ## prepare parameter field
     extra_params = Array{Float64, 1}([max(μ_final, 0.05)])
@@ -70,30 +116,52 @@ function main(; μ_final = 0.0005, order = 2, nrefs = 5, Plotter = nothing, kwar
 
     assign_unknown!(PD, u)
     assign_unknown!(PD, p)
-    assign_operator!(PD, NonlinearOperator(kernel_nonlinear!, [id(u), grad(u), id(p)]; params = extra_params, kwargs...))
-    assign_operator!(PD, InterpolateBoundaryData(u, boundarydata!; regions = 3))
-    assign_operator!(PD, HomogeneousBoundaryData(u; regions = [1, 2]))
-    assign_operator!(PD, FixDofs(p; dofs = [1]))
+    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(initialgrid_cone(), nrefs)
+    xgrid = uniform_refine(create_cone(h), nrefs)
 
     ## prepare FESpace
-    FES = [FESpace{H1Pk{2, 2, order}}(xgrid), FESpace{H1Pk{1, 2, order - 1}}(xgrid)]
+    FES = [
+        FESpace{H1Pk{2, 2, order}}(xgrid),
+        FESpace{H1Pk{1, 2, order - 1}}(xgrid),
+    ]
 
     ## prepare plots
     plt = GridVisualizer(; Plotter = Plotter, layout = (1, 2), clear = true, size = (1600, 800))
 
+
+    ## prepare an initial solution matching the boundary data
+    sol = FEVector(FES; tags = PD.unknowns)
+    interpolate!(sol[u], ON_BFACES, u_boundary!; regions = [1])
+
+
     ## 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...)
+        sol, SC = ExtendableFEM.solve(
+            PD,
+            FES,
+            SC;
+            init = sol,
+            method_linear,
+            # use new preconditioner API: https://docs.sciml.ai/LinearSolve/stable/basics/Preconditioners/#Specifying-Preconditioners
+            # this does currently now work with ILU, see https://github.com/WIAS-PDELib/ExtendableFEM.jl/pull/47#issuecomment-2796849931
+            precon_linear,
+            return_config = true,
+            target_residual = 1.0e-10,
+            maxiterations = 20,
+            kwargs...
+        )
+
         if step == 1
             initialize!(PE, sol)
         end
@@ -116,9 +184,18 @@ function main(; μ_final = 0.0005, order = 2, nrefs = 5, Plotter = nothing, kwar
 end
 
 generateplots = ExtendableFEM.default_generateplots(Example250_NSELidDrivenCavity, "example250.png") #hide
-function runtests() #hide
-    sol, plt = main(; nrefs = 3, μ_final = 0.005) #hide
-    @test sum(view(sol[1])) ≈ 9.501630403050289 #hide
-    return nothing #hide
-end #hide
+function runtests()                                                                                  #hide
+    sol, plt = main(; μ_final = 0.005)                                                               #hide
+    sum(view(sol[1])) ≈ 237.24628017878518                                                           #hide
+
+    method_linear = KrylovJL_GMRES()                                                                 #hide
+    precon_linear = IncompleteLU.ilu                                                                 #hide
+
+    sol, plt = main(; μ_final = 0.005, method_linear, precon_linear)                                 #hide
+    @test sum(view(sol[1])) ≈ 237.24628017878518                                                     #hide
+
+    return nothing                                                                                   #hide
+end                                                                                                  #hide
+
+
 end # module