# [ANN] MovingBoundaryProblems1D.jl

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)


3 Likes