Reaction-Diffusion Equation with a Time-dependent Dirichlet Boundary Condition on a Disk
In this tutorial, we consider a reaction-diffusion equation on a disk with a boundary condition of the form $\mathrm du/\mathrm dt = u$:
\[\begin{equation*} \begin{aligned} \pdv{u(r, \theta, t)}{t} &= \div[u\grad u] + u(1-u) & 0<r<1,\,0<\theta<2\pi,\\[6pt] \dv{u(1, \theta, t)}{t} &= u(1,\theta,t) & 0<\theta<2\pi,\,t>0,\\[6pt] u(r,\theta,0) &= \sqrt{I_0(\sqrt{2}r)} & 0<r<1,\,0<\theta<2\pi, \end{aligned} \end{equation*}\]
where $I_0$ is the modified Bessel function of the first kind of order zero. For this problem the diffusion function is $D(\vb x, t, u) = u$ and the source function is $R(\vb x, t, u) = u(1-u)$, or equivalently the force function is
\[\vb q(\vb x, t, \alpha,\beta,\gamma) = \left(-\alpha(\alpha x + \beta y + \gamma), -\beta(\alpha x + \beta y + \gamma)\right)^{\mkern-1.5mu\mathsf{T}}.\]
As usual, we start by generating the mesh.
using FiniteVolumeMethod, DelaunayTriangulation, ElasticArrays
r = 1.0
circle = CircularArc((0.0, r), (0.0, r), (0.0, 0.0))
points = NTuple{2, Float64}[]
boundary_nodes = [circle]
tri = triangulate(points; boundary_nodes)
A = get_area(tri)
refine!(tri; max_area=1e-4A)
mesh = FVMGeometry(tri)FVMGeometry with 8888 control volumes, 17518 triangles, and 26405 edgesusing CairoMakie
triplot(tri)
Now we define the boundary conditions and the PDE.
using Bessels
BCs = BoundaryConditions(mesh, (x, y, t, u, p) -> u, Dudt)BoundaryConditions with 1 boundary condition with type Dudtf = (x, y) -> sqrt(besseli(0.0, sqrt(2) * sqrt(x^2 + y^2)))
D = (x, y, t, u, p) -> u
R = (x, y, t, u, p) -> u * (1 - u)
initial_condition = [f(x, y) for (x, y) in DelaunayTriangulation.each_point(tri)]
final_time = 0.10
prob = FVMProblem(mesh, BCs;
diffusion_function=D,
source_function=R,
final_time,
initial_condition)FVMProblem with 8888 nodes and time span (0.0, 0.1)We can now solve.
using OrdinaryDiffEq, LinearSolve
alg = FBDF(linsolve=UMFPACKFactorization(), autodiff=false)
sol = solve(prob, alg, saveat=0.01)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.2514323512504983, 1.2514323512504983, 1.2514323512504983, 1.2514323512504983, 1.2514323512504983, 1.2514323512504983, 1.2514323512504983, 1.2514323512504983, 1.0, 1.0732643239064652 … 1.0314926207281716, 1.0303458702182084, 1.134584129951588, 1.1452063869613691, 1.150274117652823, 1.044973178267007, 1.2189095906775154, 1.120129034752343, 1.1154275787943209, 1.1776687251690439]
[1.26401009909278, 1.26401009909278, 1.26401009909278, 1.26401009909278, 1.26401009909278, 1.26401009909278, 1.26401009909278, 1.26401009909278, 1.010079320110259, 1.0840136612418794 … 1.0418458314876649, 1.040663282007412, 1.1459741255773648, 1.1567263403183645, 1.161843867498476, 1.0554640129088042, 1.2311801207173887, 1.1313820339349707, 1.1266676160413327, 1.1895217319370281]
[1.2767133779263626, 1.2767133779263626, 1.2767133779263626, 1.2767133779263626, 1.2767133779263626, 1.2767133779263626, 1.2767133779263626, 1.2767133779263626, 1.020196398075027, 1.0949559342528987 … 1.0523330184839967, 1.0511698891367502, 1.1575128260999403, 1.1683465423733979, 1.1735160262342053, 1.0660849132049408, 1.2435369061958719, 1.1427627078687648, 1.1379604244506434, 1.2014595033529263]
[1.289547019854722, 1.289547019854722, 1.289547019854722, 1.289547019854722, 1.289547019854722, 1.289547019854722, 1.289547019854722, 1.289547019854722, 1.0304826679285801, 1.1059576759755656 … 1.0628796443877715, 1.0617085289813408, 1.1691255903711848, 1.1800956156267774, 1.1853184161207906, 1.0767936790911095, 1.2560230711660976, 1.1542346782756807, 1.1493899500162217, 1.2135459019988146]
[1.3025106827520256, 1.3025106827520256, 1.3025106827520256, 1.3025106827520256, 1.3025106827520256, 1.3025106827520256, 1.3025106827520256, 1.3025106827520256, 1.0408137785634688, 1.1170945907596395 … 1.073590871292195, 1.0724149953671376, 1.180900109375207, 1.1919566205236611, 1.1972308709963917, 1.0876269489114563, 1.2686551082539554, 1.165852036599947, 1.1609447072973786, 1.2257362285702487]
[1.315601713292154, 1.315601713292154, 1.315601713292154, 1.315601713292154, 1.315601713292154, 1.315601713292154, 1.315601713292154, 1.315601713292154, 1.051266300801148, 1.1282682838171225 … 1.0844153341284695, 1.083177941337161, 1.1927921131480979, 1.2039416611689822, 1.2092678976338478, 1.098564971077452, 1.2814494933789509, 1.1775878833542042, 1.1726694261675212, 1.2380690459531225]
[1.328819827612282, 1.328819827612282, 1.328819827612282, 1.328819827612282, 1.328819827612282, 1.328819827612282, 1.328819827612282, 1.328819827612282, 1.0618438302980409, 1.1392193053314696 … 1.095454487192969, 1.0938800647394191, 1.2048653026857834, 1.2160757903351667, 1.2214496513820128, 1.1096151097494944, 1.2945938378616926, 1.189498351272095, 1.1848057280029078, 1.2506103835619344]
[1.3421681270338466, 1.3421681270338466, 1.3421681270338466, 1.3421681270338466, 1.3421681270338466, 1.3421681270338466, 1.3421681270338466, 1.3421681270338466, 1.072529852614277, 1.1501593081345876 … 1.1066463592791635, 1.1046318898632317, 1.2170842987380335, 1.2283408790531452, 1.2337606586706193, 1.120776831076538, 1.3079522098534944, 1.201550457438544, 1.1971711664673919, 1.2633058443877871]
[1.3556499728606162, 1.3556499728606162, 1.3556499728606162, 1.3556499728606162, 1.3556499728606162, 1.3556499728606162, 1.3556499728606162, 1.3556499728606162, 1.083315345853519, 1.1614114307468806 … 1.1178754296660967, 1.1155868191830278, 1.2293789057857074, 1.2407088577583074, 1.2461780165923515, 1.1320439860123472, 1.3213028248692755, 1.2136820644964874, 1.2094749746720208, 1.2760746352464698]
[1.36926872639636, 1.36926872639636, 1.36926872639636, 1.36926872639636, 1.36926872639636, 1.36926872639636, 1.36926872639636, 1.36926872639636, 1.0941912881194278, 1.1732988116887526 … 1.129026177632813, 1.1268982551732385, 1.2416789283096654, 1.253151656886044, 1.2586788222398944, 1.1434104255106858, 1.3344238984239547, 1.2258310350888597, 1.221426385727842, 1.288835962953773]
[1.383027748944847, 1.383027748944847, 1.383027748944847, 1.383027748944847, 1.383027748944847, 1.383027748944847, 1.383027748944847, 1.383027748944847, 1.1051486575156662, 1.1861445894806077 … 1.139983082458356, 1.1387196003082936, 1.2539141707907666, 1.2656412068717444, 1.2712401727059321, 1.1548700005253179, 1.347093646032451, 1.2379352318585968, 1.232734632745903, 1.301509034325486]fig = Figure(fontsize=38)
for (i, j) in zip(1:3, (1, 6, 11))
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=1:0.01:1.4, colormap=:matter)
tightlimits!(ax)
end
resize_to_layout!(fig)
fig
!DelaunayTriangulation.has_vertex(tri, i) && continueJust the code
An uncommented version of this example is given below. You can view the source code for this file here.
using FiniteVolumeMethod, DelaunayTriangulation, ElasticArrays
r = 1.0
circle = CircularArc((0.0, r), (0.0, r), (0.0, 0.0))
points = NTuple{2, Float64}[]
boundary_nodes = [circle]
tri = triangulate(points; boundary_nodes)
A = get_area(tri)
refine!(tri; max_area=1e-4A)
mesh = FVMGeometry(tri)
using CairoMakie
triplot(tri)
using Bessels
BCs = BoundaryConditions(mesh, (x, y, t, u, p) -> u, Dudt)
f = (x, y) -> sqrt(besseli(0.0, sqrt(2) * sqrt(x^2 + y^2)))
D = (x, y, t, u, p) -> u
R = (x, y, t, u, p) -> u * (1 - u)
initial_condition = [f(x, y) for (x, y) in DelaunayTriangulation.each_point(tri)]
final_time = 0.10
prob = FVMProblem(mesh, BCs;
diffusion_function=D,
source_function=R,
final_time,
initial_condition)
using OrdinaryDiffEq, LinearSolve
alg = FBDF(linsolve=UMFPACKFactorization(), autodiff=false)
sol = solve(prob, alg, saveat=0.01)
fig = Figure(fontsize=38)
for (i, j) in zip(1:3, (1, 6, 11))
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=1:0.01:1.4, colormap=:matter)
tightlimits!(ax)
end
resize_to_layout!(fig)
fig
!DelaunayTriangulation.has_vertex(tri, i) && continueThis page was generated using Literate.jl.