# \[ANN\] MovingBoundaryProblems1D.jl

**URL:** https://discourse.julialang.org/t/ann-movingboundaryproblems1d-jl/100094
**Category:** Package Announcements
**Tags:** package, announcement, pde, differentialequation
**Created:** [June 9, 2023, 10:39am UTC](https://discourse.julialang.org/t/ann-movingboundaryproblems1d-jl/100094 "2023-06-09T10:39:48Z")
**Posts on this page:** 1
**Page:** 1

<div class="post-metadata">

### Author: ![DanielVandH](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielvandh/32/31134_2.png) [@DanielVandH](https://discourse.julialang.org/u/DanielVandH)
#### Post date: [June 9, 2023, 10:39am UTC](https://discourse.julialang.org/t/ann-movingboundaryproblems1d-jl/100094/1 "2023-06-09T10:39:48Z")

</div>

Hi all,

I’ve just released my package [MovingBoundaryProblems1D.jl](https://github.com/DanielVandH/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](https://danielvandh.github.io/MovingBoundaryProblems1D.jl/dev/), 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](https://danielvandh.github.io/EpithelialDynamics1D.jl/dev/examples/) for my other work [EpithelialDynamics1D.jl](https://github.com/DanielVandH/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)]:

```julia
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)

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/4/1/41ccd8539a7649c8251389d673ad73e99e0047b3.png)