GitHub

Tutorial notebook: Navier–Stokes problem

Consider the Navier-Stokes problem that seeks \(u\) and \(p\) such that

$$\begin{aligned} - \mu \Delta u + (u \cdot \nabla) u + \nabla p &= f\\ \mathrm{div}(u) & = 0. \end{aligned}$$

The weak formulation seeks \(u \in V := H^1_0(\Omega)\) and \(p \in Q := L^2_0(\Omega)\) such that

$$\begin{aligned} \mu (\nabla u, \nabla v) + ((u \cdot \nabla) u, v) - (p, \mathrm{div}(v)) & = (f, v) & \text{for all } v \in V\\ (q, \mathrm{div}(u)) & = 0 & \text{for all } q \in Q\\ \end{aligned}$$

This tutorial notebook compute a planar lattice flow with inhomogeneous Dirichlet boundary conditions (which requires some modification above). Newton's method with automatic differentation is used to handle the nonlinear convection term.

This is a plot of the computed velocity and pressure:

2-element Vector{PlutoVTKPlot}:
 PlutoVTKPlot(Dict{String, Any}("2axisfontsize" => 10, "1" => "tricontour", "2ylabel" => "y", "2" => "axis", "cbar_levels" => [6.123233995736766e-17, 0.16667347690578882, 0.3333469538115776, 0.5000204307173663, 0.6666939076231551, 0.833367384528944, 1.0000408614347327], "cbar" => 1, "2xlabel" => "x", "2zlabel" => "z", "cbar_stops" => [0.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09  …  0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99, 1.0], "2zoom" => 1.0…), 300.0, 300.0, false, (title = "", titlefontsize = 12, axisfontsize = 10, tickfontsize = 10, xlabel = "x", ylabel = "y", zlabel = "z", aspect = 1.0, zoom = 1.0, legendfontsize = 10, colorbarticks = :default, clear = false, levels = 5), "2dd28466-a722-11ef-39c3-e1ffc42915d8")
 PlutoVTKPlot(Dict{String, Any}("2axisfontsize" => 10, "1" => "tricontour", "2ylabel" => "y", "2" => "axis", "cbar_levels" => [-0.5064542228750328, -0.33763730567358835, -0.1688203884721439, -3.4712706993289544e-6, 0.16881344593074513, 0.3376303631321897, 0.5064472803336342], "cbar" => 1, "2xlabel" => "x", "2zlabel" => "z", "cbar_stops" => [0.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09  …  0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99, 1.0], "2zoom" => 1.0…), 300.0, 300.0, false, (title = "", titlefontsize = 12, axisfontsize = 10, tickfontsize = 10, xlabel = "x", ylabel = "y", zlabel = "z", aspect = 1.0, zoom = 1.0, legendfontsize = 10, colorbarticks = :default, clear = false, levels = 5), "2e858a34-a722-11ef-3434-d348d9ce5eca")
begin
    ## PDE data
    const μ = 1e-2
    function f!(fval, x, t) # right-hand side
        fval[1] = 8.0*π*π*μ * exp(-8.0*π*π*μ*t) * sin(2.0*π*x[1])*sin(2.0*π*x[2])
        fval[2] = 8.0*π*π*μ * exp(-8.0*π*π*μ*t) * cos(2.0*π*x[1])*cos(2.0*π*x[2])
        return nothing
    end

    # exact velocity (for boundary data and error calculation)
    function u!(uval, qpinfo) 
        x = qpinfo.x
        t = qpinfo.time
        uval[1] = exp(-8.0*π*π*μ*t) * sin(2.0*π*x[1])*sin(2.0*π*x[2])
        uval[2] = exp(-8.0*π*π*μ*t) * cos(2.0*π*x[1])*cos(2.0*π*x[2])
        return nothing
    end
    
    ## discretization parameters
    const nref = 5
    const teval = 0
    const order = 2
    
    ## prepare error calculation
    function p!(pval,x,t) # exact pressure (for error calculation)
        pval[1] = exp(-16*pi*pi*μ*t)*(cos(4*pi*x[1])-cos(4*pi*x[2]))/4
        return nothing
    end
end
p! (generic function with 1 method)
begin
    ## create grid
    X = LinRange(0,1,2^nref+1)
    Y = LinRange(0,1,2^nref+1)
    println("Creating grid...")
    @time xgrid = simplexgrid(X,Y)
    println("Preparing FaceNodes...")
    @time xgrid[FaceNodes]
    println("Preparing CellVolumes...")
    @time xgrid[CellVolumes]
    xgrid
end
ExtendableGrids.ExtendableGrid{Float64, Int32}
      dim =       2
   nnodes =    1089
   ncells =    2048
  nbfaces =     128
begin
    ## create finite element space (Taylor--Hood)
    FETypes = [H1Pk{2,2,order}, H1Pk{1,2,order-1}]

    ## prepare finite element space and dofmaps
    println("Creating FESpace...")
    @time FES = [FESpace{FETypes[1]}(xgrid; name = "velocity space"),                  	             FESpace{FETypes[2]}(xgrid; name = "pressure space")]
    FES
end
2-element Vector{FESpace{Float64, Int32, FEType, ON_CELLS} where FEType<:AbstractFiniteElement}:
 
FESpace information
===================
     name = velocity space
   FEType = H1Pk{2,2,2}
  FEClass = ExtendableFEMBase.AbstractH1FiniteElement
    ndofs = 8450


DofMaps
==========

 
FESpace information
===================
     name = pressure space
   FEType = H1Pk{1,2,1}
  FEClass = ExtendableFEMBase.AbstractH1FiniteElement
    ndofs = 1089


DofMaps
==========

begin
    ## call low level solver
    sol, u_init = solve_stokes_lowlevel(FES, μ, f!)
    sol
end
FEVector information
====================
   block  |  ndofs 	|     min  /  max    	| FEType 		 (name/tag)
 [    1]  |    8450	| -1.00e+00/1.00e+00  	| velocity space  	 (#1)
 [    2]  |    1089	| -5.06e-01/5.06e-01  	| pressure space  	 (#2)
function solve_stokes_lowlevel(FES, μ, f!)
    
    println("Initializing system...")
    Solution = FEVector(FES)
    A = FEMatrix(FES)
    b = FEVector(FES)
    @time update_system! = prepare_assembly!(A, b, FES[1], FES[2], Solution, f!, μ)
    @time update_system!(true, false)
    Alin = deepcopy(A) # = keep linear part of system matrix
    blin = deepcopy(b) # = keep linear part of right-hand side
    
    println("Prepare boundary conditions...")
    @time begin
        u_init = FEVector(FES)
        interpolate!(u_init[1], u!; time = teval)
        
        fixed_dofs = [size(A.entries,1)] # fix one pressure dof = last dof
        BFaceDofs::Adjacency{Int32} = FES[1][ExtendableFEMBase.BFaceDofs]
        nbfaces::Int = num_sources(BFaceDofs)
        AM::ExtendableSparseMatrix{Float64,Int64} = A.entries
        dof_j::Int = 0
        for bface = 1 : nbfaces
            for j = 1 : num_targets(BFaceDofs,1)
                dof_j = BFaceDofs[j, bface]
                push!(fixed_dofs, dof_j)
            end
        end
    end

    
    for it = 1 : 20
        ## solve
        println("\nITERATION $it\n=============")
        println("Solving linear system...")
        @time copyto!(Solution.entries, A.entries \ b.entries)
        res = A.entries.cscmatrix * Solution.entries .- b.entries
        for dof in fixed_dofs
            res[dof] = 0
        end
        linres = norm(res)
        println("linear residual = $linres")
        
        fill!(A.entries.cscmatrix.nzval,0)
        fill!(b.entries,0)
        println("Updating linear system...")
        @time begin
            update_system!(false,true)
            A.entries.cscmatrix += Alin.entries.cscmatrix
            b.entries .+= blin.entries
        end
        
        ## fix boundary dofs
        for dof in fixed_dofs
            AM[dof,dof] = 1e60
            b.entries[dof] = 1e60 * u_init.entries[dof]
        end
        ExtendableSparse.flush!(A.entries)

        ## calculate nonlinear residual
        res = A.entries.cscmatrix * Solution.entries .- b.entries
        for dof in fixed_dofs
            res[dof] = 0
        end
        nlres = norm(res)
        println("nonlinear residual = $nlres")
        if nlres < max(1e-12, 20*linres)
            break
        end
    end

    return Solution, u_init
end
solve_stokes_lowlevel (generic function with 1 method)
function prepare_assembly!(A, b, FESu, FESp, Solution, f, μ = 1)

    A = A.entries
    b = b.entries
    Solution = Solution.entries
    xgrid = FESu.xgrid
    EG = xgrid[UniqueCellGeometries][1]
    FEType_u = eltype(FESu)
    FEType_p = eltype(FESp)
    L2G = L2GTransformer(EG, xgrid, ON_CELLS)
    cellvolumes = xgrid[CellVolumes]
    ncells::Int = num_cells(xgrid)

    ## dofmap
    CellDofs_u = FESu[ExtendableFEMBase.CellDofs]
    CellDofs_p = FESp[ExtendableFEMBase.CellDofs]
    offset_p = FESu.ndofs
    
    ## quadrature formula
    qf = QuadratureRule{Float64, EG}(3*get_polynomialorder(FEType_u, EG)-1)
    weights::Vector{Float64} = qf.w
    xref::Vector{Vector{Float64}} = qf.xref
    nweights::Int = length(weights)
    
    ## FE basis evaluator
    FEBasis_∇u = FEEvaluator(FESu, Gradient, qf)
    ∇uvals = FEBasis_∇u.cvals
    FEBasis_idu = FEEvaluator(FESu, Identity, qf)
    iduvals = FEBasis_idu.cvals
    FEBasis_idp = FEEvaluator(FESp, Identity, qf)
    idpvals = FEBasis_idp.cvals

    ## prepare automatic differentation of convection operator
    function operator!(result, input)
        # result = (u ⋅ ∇)u
        result[1] = input[1]*input[3]+input[2]*input[4]
        result[2] = input[1]*input[5]+input[2]*input[6]
    end
    result = Vector{Float64}(undef,2)
    input = Vector{Float64}(undef,6)
    tempV = zeros(Float64, 2)
    Dresult = DiffResults.JacobianResult(result, input)
    cfg = ForwardDiff.JacobianConfig(operator!, result, input, ForwardDiff.Chunk{6}())
    jac = DiffResults.jacobian(Dresult)
    value = DiffResults.value(Dresult)
    
   
    ## ASSEMBLY LOOP
    function barrier(EG, L2G::L2GTransformer, linear::Bool, nonlinear::Bool)
        ## barrier function to avoid allocations caused by L2G
        
    	ndofs4cell_u::Int = get_ndofs(ON_CELLS, FEType_u, EG)
    	ndofs4cell_p::Int = get_ndofs(ON_CELLS, FEType_p, EG)
    	Aloc = zeros(Float64, ndofs4cell_u, ndofs4cell_u)
    	Bloc = zeros(Float64, ndofs4cell_u, ndofs4cell_p)
    	dof_j::Int, dof_k::Int = 0, 0
        fval::Vector{Float64} = zeros(Float64,2)
        x::Vector{Float64} = zeros(Float64, 2)
        
        for cell = 1 : ncells
            ## update FE basis evaluators
            update_basis!(FEBasis_∇u, cell)
            update_basis!(FEBasis_idu, cell)
            update_basis!(FEBasis_idp, cell) 
    
            ## assemble local stiffness matrix (symmetric)
            if (linear)
                for j = 1 : ndofs4cell_u, k = 1 : ndofs4cell_u
                    temp = 0
                    for qp = 1 : nweights
                        temp += weights[qp] * dot(view(∇uvals,:,j,qp), view(∇uvals,:,k,qp))
                    end
                    Aloc[k,j] = μ * temp
                end

                ## assemble div-pressure coupling
                for j = 1 : ndofs4cell_u, k = 1 : ndofs4cell_p
                    temp = 0
                    for qp = 1 : nweights
                        temp -= weights[qp] * (∇uvals[1,j,qp] + ∇uvals[4,j,qp]) * 
                        idpvals[1,k,qp]
                    end
                    Bloc[j,k] = temp
                end
                Bloc .*= cellvolumes[cell]
                
                ## assemble right-hand side
                update_trafo!(L2G, cell)
                for j = 1 : ndofs4cell_u
                    ## right-hand side
                    temp = 0
                    for qp = 1 : nweights
                        ## get global x for quadrature point
                        eval_trafo!(x, L2G, xref[qp])
                        ## evaluate (f(x), v_j(x))
                        f!(fval, x, teval)
                        temp += weights[qp] * dot(view(iduvals,: , j, qp), fval)
                    end
                    ## write into global vector
                    dof_j = CellDofs_u[j, cell]
                    b[dof_j] += temp * cellvolumes[cell]
                end
            end

            ## assemble nonlinear term
            if (nonlinear)
                for qp = 1 : nweights
                    fill!(input,0)
                    for j = 1 : ndofs4cell_u
                        dof_j = CellDofs_u[j, cell]
                        for d = 1 : 2
                            input[d] += Solution[dof_j] * iduvals[d,j,qp]
                        end
                        for d = 1 : 4
                            input[2+d] += Solution[dof_j] * ∇uvals[d,j,qp]
                        end
                    end
                    
                	## evaluate jacobian
                    ForwardDiff.chunk_mode_jacobian!(Dresult, operator!, result, input, cfg)
                    
                    # update matrix
                    for j = 1 : ndofs4cell_u
                        # multiply ansatz function with local jacobian
                        fill!(tempV,0)
                        for d = 1 : 2
                            tempV[1] += jac[1,d] * iduvals[d,j,qp]
                            tempV[2] += jac[2,d] * iduvals[d,j,qp]
                        end
                        for d = 1 : 4
                            tempV[1] += jac[1,2+d] * ∇uvals[d,j,qp]
                            tempV[2] += jac[2,2+d] * ∇uvals[d,j,qp]
                        end
    
                        # multiply test function operator evaluation
                        for k = 1 : ndofs4cell_u
                            Aloc[k,j] += dot(tempV,view(iduvals,:,k,qp)) * weights[qp]
                        end
                    end 
    
                    # update rhs
                    mul!(tempV, jac, input)
                    tempV .-= value
                    for j = 1 : ndofs4cell_u
                		dof_j = CellDofs_u[j, cell]
                		b[dof_j] += dot(tempV, view(iduvals,:,j,qp)) * weights[qp] * cellvolumes[cell]
                    end
                end
            end
            
            ## add local matrices to global matrix
            Aloc .*= cellvolumes[cell]
            for j = 1 : ndofs4cell_u
                dof_j = CellDofs_u[j, cell]
                for k = 1 : ndofs4cell_u
                    dof_k = CellDofs_u[k, cell]
                    rawupdateindex!(A, +, Aloc[j,k], dof_j, dof_k)
                end
                if (linear)
                    for k = 1 : ndofs4cell_p
                        dof_k = CellDofs_p[k, cell] + offset_p
                        rawupdateindex!(A, +, Bloc[j,k], dof_j, dof_k) 
                        rawupdateindex!(A, +, Bloc[j,k], dof_k, dof_j)
                    end
                end
            end
            fill!(Aloc, 0)
            fill!(Bloc, 0)
        end
    end

    function update_system!(linear::Bool, nonlinear::Bool)
        barrier(EG, L2G, linear, nonlinear)
        flush!(A)
    end
    update_system!
end
prepare_assembly! (generic function with 2 methods)

Built with Julia 1.11.1 and

DiffResults 1.1.0
ExtendableFEMBase 0.8.1
ExtendableGrids 1.11.0
ExtendableSparse 1.6.0
ForwardDiff 0.10.38
GridVisualize 1.8.2
LinearAlgebra 1.11.0
PlutoVista 1.0.1