[ANN] MovingBoundaryProblems1D.jl

[DanielVandH]

Hi all,

I’ve just released my package MovingBoundaryProblems1D.jl, and it’s now registered. This is a package for solving one-dimensional (single phase) moving boundary problems of the form:

\begin{align*} \begin{array}{rcll} \dfrac{\partial u(x, t)}{\partial t} & = & \dfrac{\partial}{\partial x}\left(D\left(u(x, t), x, t\right)\dfrac{\partial u(x, t)}{\partial x}\right) + R\left(u(x, t), x, t\right) & 0 < x < L(t),\, t>0, \\[9pt] \dfrac{\partial u(0, t)}{\partial x} & = & a_0\left(u(0, t), t\right) & x = 0,\,t>0, \\[9pt] \dfrac{\partial u(L(t), t)}{\partial x} & = & a_1\left(u(L(t), t), t\right) & x=L(t),\,t>0, \\[9pt] \dfrac{\mathrm dL}{\mathrm dt} & = & a_2\left(u(L(t), t), t\right) + b_2\left(u(L(t), t), t\right)\dfrac{\partial u(L(t), t)}{\partial x} & x = L(t),\,t>0, \\[9pt] u(x, 0) & = & u_0(x) & 0 \leq x \leq L(0), \\[9pt] L(0) & = & L_0. \end{array} \end{align*}

Dirichlet boundary conditions are also supported with Dirichlet. The finite volume method together with the transform \xi = x/L(t) is used to solve these problems. The package is built on top of the SciMLBase interface.

There are plenty of examples in the docs, where I also talk about computing steady states and solving Stefan problems, but here is a demonstration (this is Example V in the docs - see there for more detail). The following moving boundary problem comes up when considering the dynamics of epithelial sheets (other epithelial examples are also shown here for my other work EpithelialDynamics1D.jl if you are interested):

\begin{align*} \begin{array}{rcll} \dfrac{\partial q}{\partial t} & = & \dfrac{\partial}{\partial x}\left(D(q)\dfrac{\partial q}{\partial x}\right) + qG\left(\dfrac{1}{q}\right) & 0 < x < L(t),\,t>0,\\[9pt] \dfrac{\partial q}{\partial x} & = & 0 & x = 0,\,t>0,\\[9pt] \dfrac{1}{\eta}F\left(\dfrac{1}{q}\right) + \dfrac{D(q)}{2q}\dfrac{\partial q}{\partial x} & = & 0 & x = L(t),\,t>0,\\[9pt] \dfrac{\mathrm dL}{\mathrm dt} & = & -\dfrac{D(q)}{q}\dfrac{\partial q}{\partial x} & x = L(t),\,t>0, \\[9pt] q(x, 0) & = & q_0(x) & 0 \leq x \leq L(0). \end{array} \end{align*}

Where D(q) = k/(\eta q^2), F(q) = k(s-q), G(q) = \beta, and q_0(x) is a Gaussian density. We are interested in (1) the evolution of the density over time, (2) the number of cells N(t) = \int_0^{L(t)} q(x, t)\,\mathrm dx, and (3) the position of the leading edge L(t). The following code does all these computations, noting that (1/\eta)F(1/q) + [D(q)/(2q)]\partial q/\partial x = 0 had to be rearranged into \partial q(L(t), t)/\partial x = -2qF(1/q)/[\eta D(q)]:

using MovingBoundaryProblems1D, SpecialFunctions 

## Define the parameters 
k, s, η, β = 10.0, 1.0, 1.0, 0.00577
F = (q, p) -> p.k * (p.s - q)
D = (q, p) -> p.k / (p.η * q^2)
G = (q, p) -> p.β
L₀ = 10.0
N₀ = 40.0
σ = sqrt(3)

## Define the initial condition 
mesh_points = LinRange(0, L₀, 1000)
q₀ = x -> N₀ * exp(-(2x - L₀)^2 / (8σ^2)) / (erf(L₀ * sqrt(2) / (4σ)) * sqrt(2π * σ^2))
initial_condition = q₀.(mesh_points)

## Define the PDE 
diffusion_function = (q, x, t, p) -> p.D(q, p)
diffusion_parameters = (D=D, k=k, η=η)
reaction_function = (q, x, t, p) -> q * p.G(inv(q), p)
reaction_parameters = (G=G, β=β)

## Define the boundary conditions 
lhs = Neumann(0.0)
rhs_f = (q, t, p) -> -2q * p.F(inv(q), p) / (p.η * p.D(q, p))
rhs_p = (F=F, η=η, D=D, s=s, k=k)
rhs = Neumann(rhs_f, rhs_p)
moving_boundary_f = (q, t, p) -> (zero(q), -p.D(q, p) / q)
moving_boundary_p = (D=D, k=k, η=η)
moving_boundary = Robin(moving_boundary_f, moving_boundary_p)

## Define the problem 
prob = MBProblem(mesh_points, lhs, rhs, moving_boundary;
    diffusion_function,
    diffusion_parameters,
    reaction_function,
    reaction_parameters,
    initial_condition,
    initial_endpoint=L₀,
    final_time=400.0)

## Solve the problem 
using OrdinaryDiffEq, LinearSolve
sol = solve(prob, TRBDF2(linsolve=KLUFactorization()))

## Computing the number of cells 
using DataInterpolations
function integrate_solution(prob, sol) # could also consider simple quadrature rules
    N = zeros(length(sol))
    x = prob.geometry.mesh_points
    for i in eachindex(sol)
        q = @views sol.u[i][begin:(end-1)]
        L = sol.u[i][end]
        interp = LinearInterpolation(q, x) 
        N[i] = DataInterpolations.integral(interp, 0.0, 1.0) * L # ∫₀ᴸ q(x) dx = L∫₀¹ q(ξ) dξ 
    end
    return N
end
Nt = integrate_solution(prob, sol)

## Plot 
using CairoMakie 

fig = Figure(fontsize=33)
colors = (:red, :blue, :black, :magenta, :darkgreen)

ax1 = Axis(fig[1, 1], xlabel=L"x", ylabel=L"q(x, t)",
    title=L"(a): $q(x, t)$", titlealign=:left,
    width=600, height=300)
ax2 = Axis(fig[1, 2], xlabel=L"t", ylabel=L"N(t)",
    title=L"(b): $N(t)$", titlealign=:left,
    width=600, height=300)
ax3 = Axis(fig[1, 3], xlabel=L"t", ylabel=L"L(t)",
    title=L"(c): $L(t)$", titlealign=:left,
    width=600, height=300)

t = [0.0, 100.0, 200.0, 300.0, 400.0]
qL = sol.(t)
q = [qL[begin:(end-1)] for qL in qL]
L = [qL[end] for qL in qL]
ξ_grid = prob.geometry.mesh_points
[lines!(ax1, ξ_grid .* L[i], q[i], color=colors[i]) for i in eachindex(q)]
lines!(ax2, sol.t, Nt, color=:black, linewidth=3)
lines!(ax3, sol.t, sol[end, :], color=:black, linewidth=3)
resize_to_layout!(fig)