Diffusion Equation in a Wedge with Mixed Boundary Conditions
In this example, we consider a diffusion equation on a wedge with angle $\alpha$ and mixed boundary conditions:
\[\begin{equation*} \begin{aligned} \pdv{u(r, \theta, t)}{t} &= \grad^2u(r,\theta,t), & 0<r<1,\,0<\theta<\alpha,\,t>0,\\[6pt] \pdv{u(r, 0, t)}{\theta} & = 0 & 0<r<1,\,t>0,\\[6pt] u(1, \theta, t) &= 0 & 0<\theta<\alpha,\,t>0,\\[6pt] \pdv{u(r,\alpha,t)}{\theta} & = 0 & 0<\theta<\alpha,\,t>0,\\[6pt] u(r, \theta, 0) &= f(r,\theta) & 0<r<1,\,0<\theta<\alpha, \end{aligned} \end{equation*}\]
where we take $f(r,\theta) = 1-r$ and $\alpha=\pi/4$.
Note that the PDE is provided in polar form, but Cartesian coordinates are assumed for the operators in our code. The conversion is easy, noting that the two Neumann conditions are just equations of the form $\grad u \vdot \vu n = 0$. Moreover, although the right-hand side of the PDE is given as a Laplacian, recall that $\grad^2 = \div\grad$, so we can write the PDE as $\partial u/\partial t + \div \vb q = 0$, where $\vb q = -\grad u$.
Let us now setup the problem. To define the geometry, we need to be careful that the Triangulation recognises that we need to split the boundary into three parts, one part for each boundary condition. This is accomplished by providing a single vector for each part of the boundary as follows (and as described in DelaunayTriangulation.jl's documentation), where we also refine! the mesh to get a better mesh. For the arc, we use the CircularArc so that the mesh knows that it is triangulating a certain arc in that area.
using DelaunayTriangulation, FiniteVolumeMethod, ElasticArrays
α = π / 4
points = [(0.0, 0.0), (1.0, 0.0), (cos(α), sin(α))]
bottom_edge = [1, 2]
arc = CircularArc((1.0, 0.0), (cos(α), sin(α)), (0.0, 0.0))
upper_edge = [3, 1]
boundary_nodes = [bottom_edge, [arc], upper_edge]
tri = triangulate(points; boundary_nodes)
A = get_area(tri)
refine!(tri; max_area=1e-4A)
mesh = FVMGeometry(tri)FVMGeometry with 8176 control volumes, 15985 triangles, and 24160 edgesThis is the mesh we've constructed.
using CairoMakie
fig, ax, sc = triplot(tri)
fig
To confirm that the boundary is now in three parts, see:
get_boundary_nodes(tri)3-element Vector{Vector{Int64}}:
[1, 6612, 5358, 6615, 3442, 6037, 4731, 6829, 439, 4576 … 5461, 136, 6959, 963, 6646, 145, 5702, 4190, 5398, 2]
[2, 5394, 148, 6669, 4703, 7724, 140, 4696, 7168, 174 … 6540, 2350, 4906, 99, 4889, 751, 6243, 4786, 5480, 3]
[3, 5301, 4782, 7099, 710, 6023, 4771, 4790, 85, 4795 … 6521, 438, 6915, 4488, 5360, 3441, 6446, 5357, 6611, 1]We now need to define the boundary conditions. For this, we need to provide Tuples, where the ith element of the Tuples refers to the ith part of the boundary. The boundary conditions are thus:
lower_bc = arc_bc = upper_bc = (x, y, t, u, p) -> zero(u)
types = (Neumann, Dirichlet, Neumann)
BCs = BoundaryConditions(mesh, (lower_bc, arc_bc, upper_bc), types)BoundaryConditions with 3 boundary conditions with types (Neumann, Dirichlet, Neumann)Now we can define the PDE. We use the reaction-diffusion formulation, specifying the diffusion function as a constant.
f = (x, y) -> 1 - sqrt(x^2 + y^2)
D = (x, y, t, u, p) -> one(u)
initial_condition = [f(x, y) for (x, y) in DelaunayTriangulation.each_point(tri)]
final_time = 0.1
prob = FVMProblem(mesh, BCs; diffusion_function=D, initial_condition, final_time)FVMProblem with 8176 nodes and time span (0.0, 0.1)If you did want to use the flux formulation, you would need to provide
flux = (x, y, t, α, β, γ, p) -> (-α, -β)#9 (generic function with 1 method)which replaces u with αx + βy + γ so that we approximate $\grad u$ by $(\alpha,\beta)^{\mkern-1.5mu\mathsf{T}}$, and the negative is needed since $\vb q = -\grad u$.
We now solve the problem. We provide the solver for this problem. In my experience, I've found that TRBDF2(linsolve=KLUFactorization()) typically has the best performance for these problems.
using OrdinaryDiffEq, LinearSolve
sol = solve(prob, TRBDF2(linsolve=KLUFactorization()), saveat=0.01, parallel=Val(false))retcode: Success
Interpolation: 1st order linear
t: 11-element Vector{Float64}:
0.0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
0.1
u: 11-element Vector{Vector{Float64}}:
[1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.5, 0.2742608203909226, 0.25 … 0.21102390215660594, 0.23064488154843021, 0.24518994207244638, 0.4565742808527762, 0.03802231065392825, 0.7635611349574392, 0.07869015399181345, 0.10952685520892969, 0.695094259650801, 0.16752299089587674]
[0.8216930251744623, 0.0, 0.0, 0.0, 0.0, 0.0, 0.47952443673617756, 0.4795259124903454, 0.2742608203909226, 0.23680452099197022 … 0.19875982421928057, 0.21790088551594727, 0.23210246444904253, 0.43781800582489017, 0.034085456087435115, 0.7165630572135004, 0.07168229728060796, 0.10078878739700568, 0.6598376275417672, 0.15648696940996745]
[0.7467921384789642, 0.0, 0.0, 0.0, 0.0, 0.0, 0.4579250873463284, 0.45792607330593005, 0.2742608203909226, 0.2246926380170625 … 0.18812974244430883, 0.20650631285573123, 0.2201659885530376, 0.4182853361760805, 0.03199034959545989, 0.668729551137776, 0.06735718447359969, 0.09483197439827178, 0.6218118509806899, 0.14770948145351825]
[0.6920196642268366, 0.0, 0.0, 0.0, 0.0, 0.0, 0.43544968415761237, 0.43545023344757033, 0.2742608203909226, 0.21363285315565658 … 0.17875120989325263, 0.19627674251557997, 0.20931179563101934, 0.39816095112088806, 0.030336363240709464, 0.6251460853096077, 0.0638746886797474, 0.08995863913947615, 0.5843728861782704, 0.14024560125515995]
[0.6446033689709099, 0.0, 0.0, 0.0, 0.0, 0.0, 0.41292311941321946, 0.41292339074035683, 0.2742608203909226, 0.20314304981203704 … 0.169951838009694, 0.18662831135459362, 0.19903160030056585, 0.3780004068592615, 0.028824366260305982, 0.5858359691170152, 0.060698189751909666, 0.08549386154690733, 0.5494426521332185, 0.13331662958551746]
[0.602981220262972, 0.0, 0.0, 0.0, 0.0, 0.0, 0.3907729024794227, 0.390773015185898, 0.2742608203909226, 0.1930702855682513 … 0.16156000765112, 0.17739567308806442, 0.189168748950216, 0.3581202800551191, 0.02740864677995577, 0.5501463996340893, 0.05771817745685726, 0.08129446805467515, 0.5169864425450424, 0.12675311640807538]
[0.5650750604385473, 0.0, 0.0, 0.0, 0.0, 0.0, 0.3694938456366688, 0.36949385055152484, 0.2742608203909226, 0.18331143101428285 … 0.153440722173054, 0.1684564066101971, 0.17961467011792365, 0.33894397681069993, 0.026046631464241864, 0.5170260015223731, 0.054847473456861744, 0.07724568383444007, 0.4865888495292707, 0.12041346657432943]
[0.5302188873233703, 0.0, 0.0, 0.0, 0.0, 0.0, 0.34917622722049185, 0.3491761580217371, 0.2742608203909226, 0.17386217446663574 … 0.145578067861775, 0.15979976884089436, 0.17036334063497477, 0.32056186715614937, 0.024729518219448916, 0.48617296844133323, 0.05207019388627076, 0.07332752955850709, 0.458098238797051, 0.11427554220841263]
[0.498180431695116, 0.0, 0.0, 0.0, 0.0, 0.0, 0.3298217347852667, 0.3298216197406723, 0.2742608203909226, 0.16473448493905865 … 0.13797975900503506, 0.15143557497388135, 0.16142619242965217, 0.3029860915918949, 0.023457047083576335, 0.4574699876059712, 0.049386444569434584, 0.0695408796931405, 0.43143658199766594, 0.10834317739562538]
[0.4687274243314588, 0.0, 0.0, 0.0, 0.0, 0.0, 0.3114320558853681, 0.3114319178524672, 0.2742608203909226, 0.15594033144530023 … 0.13065350989058622, 0.1433736402028831, 0.15281465743023853, 0.28622879061836365, 0.022228958094623507, 0.4307997462312883, 0.04679633133070348, 0.0658866087046046, 0.4065258507803975, 0.10262020622126805]
[0.44162759601007323, 0.0, 0.0, 0.0, 0.0, 0.0, 0.2940088780751707, 0.29400873450125914, 0.2742608203909226, 0.14749168299910925 … 0.1236070348061804, 0.1356237797216247, 0.14454016756501661, 0.27030210473598304, 0.02104499129058984, 0.4060449315322861, 0.04429995999442773, 0.06236559105916366, 0.3832880167945281, 0.097110462770641]using CairoMakie
fig = Figure(fontsize=38)
for (i, j) in zip(1:3, (1, 6, 11))
local ax
ax = Axis(fig[1, i], width=600, height=600,
xlabel="x", ylabel="y",
title="t = $(sol.t[j])",
titlealign=:left)
tricontourf!(ax, tri, sol.u[j], levels=0:0.01:1, colormap=:matter)
tightlimits!(ax)
end
resize_to_layout!(fig)
fig
Just the code
An uncommented version of this example is given below. You can view the source code for this file here.
using DelaunayTriangulation, FiniteVolumeMethod, ElasticArrays
α = π / 4
points = [(0.0, 0.0), (1.0, 0.0), (cos(α), sin(α))]
bottom_edge = [1, 2]
arc = CircularArc((1.0, 0.0), (cos(α), sin(α)), (0.0, 0.0))
upper_edge = [3, 1]
boundary_nodes = [bottom_edge, [arc], upper_edge]
tri = triangulate(points; boundary_nodes)
A = get_area(tri)
refine!(tri; max_area=1e-4A)
mesh = FVMGeometry(tri)
using CairoMakie
fig, ax, sc = triplot(tri)
fig
get_boundary_nodes(tri)
lower_bc = arc_bc = upper_bc = (x, y, t, u, p) -> zero(u)
types = (Neumann, Dirichlet, Neumann)
BCs = BoundaryConditions(mesh, (lower_bc, arc_bc, upper_bc), types)
f = (x, y) -> 1 - sqrt(x^2 + y^2)
D = (x, y, t, u, p) -> one(u)
initial_condition = [f(x, y) for (x, y) in DelaunayTriangulation.each_point(tri)]
final_time = 0.1
prob = FVMProblem(mesh, BCs; diffusion_function=D, initial_condition, final_time)
flux = (x, y, t, α, β, γ, p) -> (-α, -β)
using OrdinaryDiffEq, LinearSolve
sol = solve(prob, TRBDF2(linsolve=KLUFactorization()), saveat=0.01, parallel=Val(false))
using CairoMakie
fig = Figure(fontsize=38)
for (i, j) in zip(1:3, (1, 6, 11))
local ax
ax = Axis(fig[1, i], width=600, height=600,
xlabel="x", ylabel="y",
title="t = $(sol.t[j])",
titlealign=:left)
tricontourf!(ax, tri, sol.u[j], levels=0:0.01:1, colormap=:matter)
tightlimits!(ax)
end
resize_to_layout!(fig)
figThis page was generated using Literate.jl.