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 8946 control volumes, 17634 triangles, and 26579 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 8946 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.0195288169650498, 1.049951240619799, 1.051222108035346, 1.2034832383083245, 1.1659603065940944, 1.1713660564109187, 1.048821102568457, 1.0908021706611652, 1.0927000706710919, 1.0526110220030425]
[1.2640100443849882, 1.2640100443849882, 1.2640100443849882, 1.2640100443849882, 1.2640100443849882, 1.2640100443849882, 1.2640100443849882, 1.2640100443849882, 1.0100508315497618, 1.0841896054826454 … 1.0297242815471341, 1.0604099668415994, 1.061777867555827, 1.2155776618855119, 1.1776591227650801, 1.1831181204608627, 1.0593085831863873, 1.1017516723701355, 1.1036340875976107, 1.0631736356403894]
[1.276713393923123, 1.276713393923123, 1.276713393923123, 1.276713393923123, 1.276713393923123, 1.276713393923123, 1.276713393923123, 1.276713393923123, 1.0202000216233924, 1.0946977288286597 … 1.040195555561176, 1.0713237281692425, 1.0724768172554613, 1.2278090678815954, 1.1895451718270802, 1.195099760749267, 1.0700931305020651, 1.1128633856660775, 1.114868366128756, 1.0739127160855582]
[1.2895471779387442, 1.2895471779387442, 1.2895471779387442, 1.2895471779387442, 1.2895471779387442, 1.2895471779387442, 1.2895471779387442, 1.2895471779387442, 1.0304556333880142, 1.105004483466503 … 1.0508738325790115, 1.0825624647849255, 1.0833147687632716, 1.2401735887594356, 1.2015906224361652, 1.2072775551579549, 1.081106408201479, 1.124124806219305, 1.1263444330796104, 1.0848157635039926]
[1.3025106254398593, 1.3025106254398593, 1.3025106254398593, 1.3025106254398593, 1.3025106254398593, 1.3025106254398593, 1.3025106254398593, 1.3025106254398593, 1.0408323496279388, 1.1165456688876156 … 1.0613383027214451, 1.0932113676696646, 1.094210617425147, 1.2526153415928174, 1.213622473318858, 1.219300184132052, 1.0918687211989409, 1.1354000537042388, 1.137550578849258, 1.0957014768739086]
[1.3156015580932814, 1.3156015580932814, 1.3156015580932814, 1.3156015580932814, 1.3156015580932814, 1.3156015580932814, 1.3156015580932814, 1.3156015580932814, 1.0512650710432427, 1.1287545730983366 … 1.071631800441959, 1.10344212251076, 1.105060829481845, 1.2651977955595999, 1.2256901350919183, 1.2313493381547955, 1.1023937301743354, 1.1466750023633563, 1.1485227044736595, 1.1064844001833058]
[1.328822551519961, 1.328822551519961, 1.328822551519961, 1.328822551519961, 1.328822551519961, 1.328822551519961, 1.328822551519961, 1.328822551519961, 1.061822177710086, 1.1403435723173516 … 1.0823099026963456, 1.1143403369763303, 1.1161264569240705, 1.2779114865999546, 1.237977578124719, 1.2436704234326137, 1.1133592061672737, 1.1581637430720322, 1.159945459193002, 1.1175421884620396]
[1.342176117979284, 1.342176117979284, 1.342176117979284, 1.342176117979284, 1.342176117979284, 1.342176117979284, 1.342176117979284, 1.342176117979284, 1.0725059478847305, 1.1515162039113929 … 1.0933077468673893, 1.1257733876145677, 1.1273948560601306, 1.2907525185408903, 1.2504598369478839, 1.2562223168561963, 1.1246927576527224, 1.1698477529520201, 1.1717467559403292, 1.1288490568914136]
[1.355664352105631, 1.355664352105631, 1.355664352105631, 1.355664352105631, 1.355664352105631, 1.355664352105631, 1.355664352105631, 1.355664352105631, 1.0832986095293549, 1.162731072722484 … 1.1044360164827367, 1.1373722824311348, 1.1387865633184018, 1.3037239970163141, 1.2630733177215387, 1.2689213263987897, 1.136167714737864, 1.1816587504341278, 1.1837007015679175, 1.1402871244823456]
[1.3692874036955431, 1.3692874036955431, 1.3692874036955431, 1.3692874036955431, 1.3692874036955431, 1.3692874036955431, 1.3692874036955431, 1.3692874036955431, 1.0941884391143673, 1.1741425238185619 … 1.1156239500344547, 1.148995537941922, 1.1502611670379599, 1.316830786960187, 1.275797674279118, 1.2817383941559402, 1.1476931603960185, 1.193566637670295, 1.195712547349068, 1.1517992924411675]
[1.3830463861860862, 1.3830463861860862, 1.3830463861860862, 1.3830463861860862, 1.3830463861860862, 1.3830463861860862, 1.3830463861860862, 1.3830463861860862, 1.1051814830726971, 1.1861591692243643 … 1.1267616062052712, 1.1603869112547631, 1.1617960836542283, 1.3300589115963801, 1.2885906908740354, 1.2945684069204564, 1.1591388069790087, 1.2055339761079609, 1.2076355901704512, 1.1633324744120304]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
Just 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)
figThis page was generated using Literate.jl.