Linear Reaction-Diffusion Equations

• Linear Reaction-Diffusion Equations
• • Mathematical Details
  • Implementation
  • Using the Provided Template
  • Just the code

Next, we write a specialised solver for solving linear reaction-diffusion equations. What we produce in this section can also be accessed in FiniteVolumeMethod.LinearReactionDiffusionEquation.

Mathematical Details

To start, let's give the mathematical details. The problems we will be solving take the form

\[\pdv{u}{t} = \div\left[D(\vb x)\grad u\right] + f(\vb x)u.\]

We want to turn this into an equation of the form $\mathrm d\vb u/\mathrm dt = \vb A\vb u + \vb b$ as usual. This takes the same form as our diffusion equation example, except with the extra $f(\vb x)u$ term, which just adds an extra $f(\vb x)$ term to the diagonal of $\vb A$. See the previois sections for further mathematical details.

Implementation

Let us now implement the solver. For constructing $\vb A$, we can use FiniteVolumeMethod.triangle_contributions! as in the previous sections, but we will need an extra function to add $f(\vb x)$ to the appropriate diagonals. We can also reuse apply_dirichlet_conditions!, apply_dudt_conditions, and boundary_edge_contributions! from the diffusion equation example. Here is our implementation.

using FiniteVolumeMethod, SparseArrays, OrdinaryDiffEq, LinearAlgebra, SciMLOperators
const FVM = FiniteVolumeMethod
function linear_source_contributions!(
        A, mesh, conditions, source_function, source_parameters)
    for i in each_solid_vertex(mesh.triangulation)
        if !FVM.has_condition(conditions, i)
            x, y = get_point(mesh, i)
            A[i, i] += source_function(x, y, source_parameters)
        end
    end
end
function linear_reaction_diffusion_equation(mesh::FVMGeometry,
        BCs::BoundaryConditions,
        ICs::InternalConditions = InternalConditions();
        diffusion_function,
        diffusion_parameters = nothing,
        source_function,
        source_parameters = nothing,
        initial_condition,
        initial_time = 0.0,
        final_time)
    conditions = Conditions(mesh, BCs, ICs)
    n = DelaunayTriangulation.num_solid_vertices(mesh.triangulation)
    Afull = zeros(n + 1, n + 1)
    A = @views Afull[begin:(end - 1), begin:(end - 1)]
    b = @views Afull[begin:(end - 1), end]
    _ic = vcat(initial_condition, 1)
    FVM.triangle_contributions!(
        A, mesh, conditions, diffusion_function, diffusion_parameters)
    FVM.boundary_edge_contributions!(
        A, b, mesh, conditions, diffusion_function, diffusion_parameters)
    linear_source_contributions!(A, mesh, conditions, source_function, source_parameters)
    FVM.apply_dudt_conditions!(b, mesh, conditions)
    FVM.apply_dirichlet_conditions!(_ic, mesh, conditions)
    FVM.fix_missing_vertices!(A, b, mesh)
    Af = sparse(Afull)
    prob = ODEProblem(MatrixOperator(Af), _ic, (initial_time, final_time))
    return prob
end

linear_reaction_diffusion_equation (generic function with 2 methods)

If you go and look back at the diffusion_equation function from the diffusion equation example, you will see that this is essentially the same function except we now have linear_source_contributions! and source_function and source_parameters arguments.

Let's now test this function. We consider the problem

\[\pdv{T}{t} = \div\left[10^{-3}x^2y\grad T\right] + (x-1)(y-1)T, \quad \vb x \in [0,1]^2,\]

with $\grad T \vdot\vu n = 1$ on the boundary.

using DelaunayTriangulation
tri = triangulate_rectangle(0, 1, 0, 1, 150, 150, single_boundary = true)
mesh = FVMGeometry(tri)
BCs = BoundaryConditions(mesh, (x, y, t, u, p) -> one(x), Neumann)
diffusion_function = (x, y, p) -> p.D * x^2 * y
diffusion_parameters = (D = 1e-3,)
source_function = (x, y, p) -> (x - 1) * (y - 1)
initial_condition = [x^2 + y^2 for (x, y) in DelaunayTriangulation.each_point(tri)]
final_time = 8.0
prob = linear_reaction_diffusion_equation(mesh, BCs;
    diffusion_function, diffusion_parameters,
    source_function, initial_condition, final_time)

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
Non-trivial mass matrix: false
timespan: (0.0, 8.0)
u0: 22501-element Vector{Float64}:
 0.0
 4.504301608035674e-5
 0.00018017206432142696
 0.00040538714472321063
 0.0007206882572857079
 ⋮
 1.9601369307688843
 1.9733345344804287
 1.986622224224134
 2.0
 1.0

sol = solve(prob, Tsit5(); saveat = 2)

retcode: Success
Interpolation: 1st order linear
t: 5-element Vector{Float64}:
 0.0
 2.0
 4.0
 6.0
 8.0
u: 5-element Vector{Vector{Float64}}:
 [0.0, 4.504301608035674e-5, 0.00018017206432142696, 0.00040538714472321063, 0.0007206882572857079, 0.0011260754020089186, 0.0016215485788928425, 0.0022071077879374803, 0.0028827530291428314, 0.0036484843025088956  …  1.8955002026935723, 1.9082473762443133, 1.921084635827215, 1.9340119814422772, 1.9470294130895003, 1.9601369307688843, 1.9733345344804287, 1.986622224224134, 2.0, 1.0]
 [5.392972715241145e-10, 0.0003283969313676686, 0.0012960759990657787, 0.002877290944573625, 0.005046982038098889, 0.007780761845897447, 0.011054900287088448, 0.014846309987847453, 0.01913253192748715, 0.023891721371031464  …  1.857601258988074, 1.8669187520861212, 1.8756974896164815, 1.8847638874373243, 1.8917905879259949, 1.9028816755063145, 1.9023588034293817, 1.9299909658674983, 1.8825037936504037, 1.0]
 [7.916735417539208e-9, 0.0023942552773780603, 0.009323370071034549, 0.020421920082540512, 0.035343899135842696, 0.05376187635010648, 0.07536617230272585, 0.0998640679093906, 0.12697904481386701, 0.1564500561223992  …  1.8509136105071593, 1.8585916000623355, 1.866489589989137, 1.8731272794777627, 1.8825835384209562, 1.884116496741115, 1.9053114153277322, 1.8762173442192143, 1.9784909526952361, 1.0]
 [8.71629135976034e-8, 0.01745587228824827, 0.06706791816304722, 0.1449467115730236, 0.24751160642335715, 0.3714705521119446, 0.5138009456452567, 0.6717316202228336, 0.8427259076096044, 1.0244657148949712  …  1.8560527056615181, 1.8629729344642032, 1.8696298851466968, 1.8768276952756353, 1.882405375120809, 1.8921110896830655, 1.8912408569476789, 1.91713526139898, 1.8728687859752957, 1.0]
 [8.530436664276854e-7, 0.12726600200660892, 0.48245430288397806, 1.0287720968254592, 1.7333054836766981, 2.5666820729791886, 3.5027585194801283, 4.51833268831312, 5.592878646428223, 6.708302802132751  …  1.8700526286528925, 1.8761847240187435, 1.8825065062351811, 1.8884430769360525, 1.8956204930445022, 1.8997921589345754, 1.911893062612213, 1.9042189119218382, 1.9487418307408346, 1.0]

using CairoMakie
fig = Figure(fontsize = 38)
for j in eachindex(sol.u)
    ax = Axis(fig[1, j], width = 600, height = 600,
        xlabel = "x", ylabel = "y",
        title = "t = $(sol.t[j])")
    tricontourf!(ax, tri, sol.u[j], levels = 0:0.1:1, extendlow = :auto,
        extendhigh = :auto, colormap = :turbo)
    tightlimits!(ax)
end
resize_to_layout!(fig)
fig

Here is how we could convert this into an FVMProblem. Note that the Neumann boundary conditions are expressed as $\grad T\vdot\vu n = 1$ above, but for FVMProblem we need them in the form $\vb q\vdot\vu n = \ldots$. For this problem, $\vb q=-D\grad T$, which gives $\vb q\vdot\vu n = -D$.

_BCs = BoundaryConditions(mesh, (x, y, t, u, p) -> -p.D(x, y, p.Dp), Neumann;
    parameters = (D = diffusion_function, Dp = diffusion_parameters))
fvm_prob = FVMProblem(
    mesh,
    _BCs;
    diffusion_function = let D=diffusion_function
        (x, y, t, u, p) -> D(x, y, p)
    end,
    diffusion_parameters = diffusion_parameters,
    source_function = let S=source_function
        (x, y, t, u, p) -> S(x, y, p) * u
    end,
    final_time = final_time,
    initial_condition = initial_condition
)
fvm_sol = solve(fvm_prob, Tsit5(), saveat = 2.0)

retcode: Success
Interpolation: 1st order linear
t: 5-element Vector{Float64}:
 0.0
 2.0
 4.0
 6.0
 8.0
u: 5-element Vector{Vector{Float64}}:
 [0.0, 4.504301608035674e-5, 0.00018017206432142696, 0.00040538714472321063, 0.0007206882572857079, 0.0011260754020089186, 0.0016215485788928425, 0.0022071077879374803, 0.0028827530291428314, 0.0036484843025088956  …  1.882843115174992, 1.8955002026935723, 1.9082473762443133, 1.921084635827215, 1.9340119814422772, 1.9470294130895003, 1.9601369307688843, 1.9733345344804287, 1.986622224224134, 2.0]
 [5.392972715241148e-10, 0.0003283969313676689, 0.0012960759990657787, 0.0028772909445736266, 0.0050469820380988906, 0.007780761845897442, 0.011054900287088448, 0.014846309987847429, 0.019132531927487158, 0.02389172137103148  …  1.8480635912090662, 1.8576012582625903, 1.8669187542720505, 1.8756974833133424, 1.884763904846072, 1.8917905417452556, 1.902881793954017, 1.9023585053824077, 1.9299917234623891, 1.8825017475754235]
 [7.91673541753924e-9, 0.00239425527737807, 0.009323370071034566, 0.020421920082540564, 0.03534389913584279, 0.05376187635010652, 0.07536617230272606, 0.09986406790939072, 0.12697904481386738, 0.15645005612239948  …  1.8429322427231234, 1.8509136099818762, 1.8585916016450499, 1.8664895854253758, 1.8731272920825033, 1.882583504983962, 1.884116582502736, 1.9053111995279193, 1.876217892753008, 1.9784894712421874]
 [8.716291359760363e-8, 0.017455872288248315, 0.06706791816304726, 0.14494671157302377, 0.24751160642335737, 0.3714705521119439, 0.5138009456452576, 0.671731620222835, 0.8427259076096059, 1.024465714894972  …  1.8491360538216244, 1.856052707046341, 1.8629729302916802, 1.869629897178227, 1.8768276620455662, 1.8824054632712972, 1.892110863588329, 1.8912414258641936, 1.9171338152905495, 1.872872691553759]
 [8.5304366642769e-7, 0.12726600200660945, 0.4824543028839793, 1.0287720968254628, 1.7333054836767012, 2.566682072979191, 3.5027585194801394, 4.518332688313142, 5.592878646428247, 6.708302802132766  …  1.863890335744982, 1.8700526280569478, 1.8761847258143511, 1.8825065010575321, 1.8884430912362864, 1.8956204551098068, 1.899792256232221, 1.9118928177845884, 1.9042195342403732, 1.9487401500135686]

Using the Provided Template

The above code is implemented in LinearReactionDiffusionEquation in FiniteVolumeMethod.jl.

prob = LinearReactionDiffusionEquation(mesh, BCs;
    diffusion_function, diffusion_parameters,
    source_function, initial_condition, final_time)
sol = solve(prob, Tsit5(); saveat = 2)

retcode: Success
Interpolation: 1st order linear
t: 5-element Vector{Float64}:
 0.0
 2.0
 4.0
 6.0
 8.0
u: 5-element Vector{Vector{Float64}}:
 [0.0, 4.504301608035674e-5, 0.00018017206432142696, 0.00040538714472321063, 0.0007206882572857079, 0.0011260754020089186, 0.0016215485788928425, 0.0022071077879374803, 0.0028827530291428314, 0.0036484843025088956  …  1.8955002026935723, 1.9082473762443133, 1.921084635827215, 1.9340119814422772, 1.9470294130895003, 1.9601369307688843, 1.9733345344804287, 1.986622224224134, 2.0, 1.0]
 [5.392972715241145e-10, 0.0003283969313676686, 0.0012960759990657787, 0.002877290944573625, 0.005046982038098889, 0.007780761845897447, 0.011054900287088448, 0.014846309987847453, 0.01913253192748715, 0.023891721371031464  …  1.857601258988074, 1.8669187520861212, 1.8756974896164815, 1.8847638874373243, 1.8917905879259949, 1.9028816755063145, 1.9023588034293817, 1.9299909658674983, 1.8825037936504037, 1.0]
 [7.916735417539208e-9, 0.0023942552773780603, 0.009323370071034549, 0.020421920082540512, 0.035343899135842696, 0.05376187635010648, 0.07536617230272585, 0.0998640679093906, 0.12697904481386701, 0.1564500561223992  …  1.8509136105071593, 1.8585916000623355, 1.866489589989137, 1.8731272794777627, 1.8825835384209562, 1.884116496741115, 1.9053114153277322, 1.8762173442192143, 1.9784909526952361, 1.0]
 [8.71629135976034e-8, 0.01745587228824827, 0.06706791816304722, 0.1449467115730236, 0.24751160642335715, 0.3714705521119446, 0.5138009456452567, 0.6717316202228336, 0.8427259076096044, 1.0244657148949712  …  1.8560527056615181, 1.8629729344642032, 1.8696298851466968, 1.8768276952756353, 1.882405375120809, 1.8921110896830655, 1.8912408569476789, 1.91713526139898, 1.8728687859752957, 1.0]
 [8.530436664276854e-7, 0.12726600200660892, 0.48245430288397806, 1.0287720968254592, 1.7333054836766981, 2.5666820729791886, 3.5027585194801283, 4.51833268831312, 5.592878646428223, 6.708302802132751  …  1.8700526286528925, 1.8761847240187435, 1.8825065062351811, 1.8884430769360525, 1.8956204930445022, 1.8997921589345754, 1.911893062612213, 1.9042189119218382, 1.9487418307408346, 1.0]

Here is a benchmark comparison of LinearReactionDiffusionEquation versus FVMProblem.

using BenchmarkTools
using Sundials
@btime solve($prob, $CVODE_BDF(linear_solver = :GMRES); saveat = $2);

  48.360 ms (1087 allocations: 1.58 MiB)

@btime solve($fvm_prob, $CVODE_BDF(linear_solver = :GMRES); saveat = $2);

  163.686 ms (83267 allocations: 90.84 MiB)

Just the code

An uncommented version of this example is given below. You can view the source code for this file here.

using FiniteVolumeMethod, SparseArrays, OrdinaryDiffEq, LinearAlgebra, SciMLOperators
const FVM = FiniteVolumeMethod
function linear_source_contributions!(
        A, mesh, conditions, source_function, source_parameters)
    for i in each_solid_vertex(mesh.triangulation)
        if !FVM.has_condition(conditions, i)
            x, y = get_point(mesh, i)
            A[i, i] += source_function(x, y, source_parameters)
        end
    end
end
function linear_reaction_diffusion_equation(mesh::FVMGeometry,
        BCs::BoundaryConditions,
        ICs::InternalConditions = InternalConditions();
        diffusion_function,
        diffusion_parameters = nothing,
        source_function,
        source_parameters = nothing,
        initial_condition,
        initial_time = 0.0,
        final_time)
    conditions = Conditions(mesh, BCs, ICs)
    n = DelaunayTriangulation.num_solid_vertices(mesh.triangulation)
    Afull = zeros(n + 1, n + 1)
    A = @views Afull[begin:(end - 1), begin:(end - 1)]
    b = @views Afull[begin:(end - 1), end]
    _ic = vcat(initial_condition, 1)
    FVM.triangle_contributions!(
        A, mesh, conditions, diffusion_function, diffusion_parameters)
    FVM.boundary_edge_contributions!(
        A, b, mesh, conditions, diffusion_function, diffusion_parameters)
    linear_source_contributions!(A, mesh, conditions, source_function, source_parameters)
    FVM.apply_dudt_conditions!(b, mesh, conditions)
    FVM.apply_dirichlet_conditions!(_ic, mesh, conditions)
    FVM.fix_missing_vertices!(A, b, mesh)
    Af = sparse(Afull)
    prob = ODEProblem(MatrixOperator(Af), _ic, (initial_time, final_time))
    return prob
end

using DelaunayTriangulation
tri = triangulate_rectangle(0, 1, 0, 1, 150, 150, single_boundary = true)
mesh = FVMGeometry(tri)
BCs = BoundaryConditions(mesh, (x, y, t, u, p) -> one(x), Neumann)
diffusion_function = (x, y, p) -> p.D * x^2 * y
diffusion_parameters = (D = 1e-3,)
source_function = (x, y, p) -> (x - 1) * (y - 1)
initial_condition = [x^2 + y^2 for (x, y) in DelaunayTriangulation.each_point(tri)]
final_time = 8.0
prob = linear_reaction_diffusion_equation(mesh, BCs;
    diffusion_function, diffusion_parameters,
    source_function, initial_condition, final_time)

sol = solve(prob, Tsit5(); saveat = 2)

using CairoMakie
fig = Figure(fontsize = 38)
for j in eachindex(sol.u)
    ax = Axis(fig[1, j], width = 600, height = 600,
        xlabel = "x", ylabel = "y",
        title = "t = $(sol.t[j])")
    tricontourf!(ax, tri, sol.u[j], levels = 0:0.1:1, extendlow = :auto,
        extendhigh = :auto, colormap = :turbo)
    tightlimits!(ax)
end
resize_to_layout!(fig)
fig

_BCs = BoundaryConditions(mesh, (x, y, t, u, p) -> -p.D(x, y, p.Dp), Neumann;
    parameters = (D = diffusion_function, Dp = diffusion_parameters))
fvm_prob = FVMProblem(
    mesh,
    _BCs;
    diffusion_function = let D=diffusion_function
        (x, y, t, u, p) -> D(x, y, p)
    end,
    diffusion_parameters = diffusion_parameters,
    source_function = let S=source_function
        (x, y, t, u, p) -> S(x, y, p) * u
    end,
    final_time = final_time,
    initial_condition = initial_condition
)
fvm_sol = solve(fvm_prob, Tsit5(), saveat = 2.0)

prob = LinearReactionDiffusionEquation(mesh, BCs;
    diffusion_function, diffusion_parameters,
    source_function, initial_condition, final_time)
sol = solve(prob, Tsit5(); saveat = 2)

This page was generated using Literate.jl.