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), "8f95e526-3241-11f0-08cc-6110dc25d8f1")
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), "90550f00-3241-11f0-334a-eb8e6593354a")
begin
## PDE data
const μ = 1.0e-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
endp! (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
endExtendableGrids.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
end2-element Vector{FESpace{Float64, Int32, FEType, ON_CELLS} where FEType<:AbstractFiniteElement}:
FESpace information
===================
name = velocity space
FEType = H1Pk{2,2,2} (ON_CELLS)
FEClass = ExtendableFEMBase.AbstractH1FiniteElement
ndofs = 8450 (coffset = 4225)
xgrid = ExtendableGrids.ExtendableGrid{Float64, Int32}(dim=2, nnodes=1089, ncells=2048, nbfaces=128)
dofgrid = xgrid
DofMaps
=======
FESpace information
===================
name = pressure space
FEType = H1Pk{1,2,1} (ON_CELLS)
FEClass = ExtendableFEMBase.AbstractH1FiniteElement
ndofs = 1089
xgrid = ExtendableGrids.ExtendableGrid{Float64, Int32}(dim=2, nnodes=1089, ncells=2048, nbfaces=128)
dofgrid = xgrid
DofMaps
=======
begin
## call low level solver
sol, u_init = solve_stokes_lowlevel(FES, μ, f!)
sol
endFEVector information ==================== block | starts | ends | length | min / max | FEType (name/tag) [ 1] | 1 | 8450 | 8450 | -1.00e+00/1.00e+00 | velocity space (#1) [ 2] | 8451 | 9539 | 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 in 1:nbfaces
for j in 1:num_targets(BFaceDofs, 1)
dof_j = BFaceDofs[j, bface]
push!(fixed_dofs, dof_j)
end
end
end
for it in 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] = 1.0e60
b.entries[dof] = 1.0e60 * 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(1.0e-12, 20 * linres)
break
end
end
return Solution, u_init
endsolve_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]
return 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 in 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 in 1:ndofs4cell_u, k in 1:ndofs4cell_u
temp = 0
for qp in 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 in 1:ndofs4cell_u, k in 1:ndofs4cell_p
temp = 0
for qp in 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 in 1:ndofs4cell_u
## right-hand side
temp = 0
for qp in 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 in 1:nweights
fill!(input, 0)
for j in 1:ndofs4cell_u
dof_j = CellDofs_u[j, cell]
for d in 1:2
input[d] += Solution[dof_j] * iduvals[d, j, qp]
end
for d in 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 in 1:ndofs4cell_u
# multiply ansatz function with local jacobian
fill!(tempV, 0)
for d in 1:2
tempV[1] += jac[1, d] * iduvals[d, j, qp]
tempV[2] += jac[2, d] * iduvals[d, j, qp]
end
for d in 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 in 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 in 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 in 1:ndofs4cell_u
dof_j = CellDofs_u[j, cell]
for k in 1:ndofs4cell_u
dof_k = CellDofs_u[k, cell]
rawupdateindex!(A, +, Aloc[j, k], dof_j, dof_k)
end
if (linear)
for k in 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
return
end
function update_system!(linear::Bool, nonlinear::Bool)
barrier(EG, L2G, linear, nonlinear)
return flush!(A)
end
return update_system!
endprepare_assembly! (generic function with 2 methods)
Built with Julia 1.11.5 and
DiffResults 1.1.0ExtendableFEMBase 1.1.0
ExtendableGrids 1.13.1
ExtendableSparse 1.7.1
ForwardDiff 1.0.1
GridVisualize 1.12.0
LinearAlgebra 1.11.0
PlutoVista 1.2.1