\[ \renewcommand{\vec}[1]{\boldsymbol{#1}} \renewcommand{\tilde}[1]{\widetilde{#1}} \newcommand{\N}{\mathcal{N}} \newcommand{\R}{\mathbb{R}} \newcommand{\T}{\mathcal{T}} \newcommand{\dd}{\,\mathrm{d}} \newcommand{\dt}{\frac{\mathrm{d}}{\mathrm{d}t}} \newcommand{\md}{\partial^\bullet} \newcommand{\abs}[1]{\left|#1\right|} \newcommand{\norm}[1]{\|#1\|} \newcommand{\diam}{\mathop{\rm diam}\nolimits} \]

Numerical analysis of a coupled bulk-surface problem in an evolving domain

Tom Ranner

joint work with Charles M. Elliott and Pierre Stepanov

Slides available at tomranner.org/fbp2021

The problem

Today: error analysis

Given initial condition \(u_0\), (possibly discontinuous) diffusion coefficient \(A\), and data \(f_1, f_2, g\), find \(u \colon \Omega \to \R\) such that

\[ \begin{aligned} u_t - \nabla \cdot( A|_{\Omega_i(t)} \nabla u ) & = f_i && \mbox{ in } \Omega_i(t) \\ [A \partial_\nu u ] = g \quad \mbox{ and } \quad [u] & = 0 && \mbox{ on } \Gamma(t) \\ u & = 0 && \mbox{ on } \partial \Omega \\ u(\cdot, 0) & = u_0 && \mbox{ in } \Omega. \end{aligned} \]

\([\eta]\) denotes jump in \(\eta\) across \(\Gamma(t)\).

Formulation of problem

  • We assume the (nice, smooth) flow of the domain given by the flow \(\Phi^t \colon \Omega \to \Omega\).

  • We assume the flow preserves partition of the domain so that \(\Phi^t(\Gamma_0) = \Gamma(t)\).

  • We define a global velocity field \(\vec{w}\) by \[ \dt \Phi_t(x) = w(\Phi_t(x), t) \qquad \mbox{ for all } x \in \Omega. \]

Moving spaces

  • We identify a pair \(v = (v_1, v_2)\) for a function \(v\) defined on \(\Omega\) by \(v_i = v|_{\Omega_i}\).

  • We define the spaces: \[ \begin{aligned} H(t) & = L^2(\Omega_1(t)) \times L^2(\Omega_2(t)) \\ V(t) & = \{ v \in H^1(\Omega_1(t)) \times H^1(\Omega_2(t)) \,|\, u_1|_{\Gamma(t)} = u_2|_{\Gamma(t)} \mbox{ and } u_2|_{\partial \Omega} = 0 \}. \end{aligned} \]

  • \(V(t) \subset H(t) \subset V^*(t)\) form an evolving Hilbert triple.

(Alphonse et al., 2015)

Moving spaces (ii)

  • Each spaces forms a compatible pair with the push-forward (or pull-back) map given by \[ \begin{aligned} \phi_t v &= (v_1(\Phi_{-t}), v_2(\Phi_{-t})) && \mbox{ for } v \in H(0) \\ \phi_{-t} v &= (v_1(\Phi_{t}), v_2(\Phi_{t})) && \mbox{ for } v \in H(t). \end{aligned} \]

  • This allows us to define moving Sobolev spaces (for \(X = V, H\)) \[ L^2_X := \{ v \colon [0,T] \to \bigcup_{t \in [0,T]} X(t) \times \{ t \} \,|\, t \mapsto (\hat{v}(t), t), \phi_{-t} \hat{v}(t) \in L^2(0, T; X_0) \} \]

  • and a material derivative for smooth function \(v\) by \[ \md_t v := \phi_t \left( \dt \phi_{-t} v \right) = \partial_t v + \vec{w} \cdot \nabla \eta. \]

(Alphonse et al., 2015)

The variational form

Let \(f \in L^2_{V^*}\) and \(g \in L^2_{H^{-1/2}(\Gamma)}\), find \(u \in \{ v \in L^2_V : \md_t v \in L^2_{V^*}\}\) such that \[ \begin{aligned} \int_0^T \langle \md_t u, v \rangle_{V(t)} \dd t + \underbrace{\sum_i \int_0^T \int_{\Omega_i(t)} A_i \nabla u_i \cdot \nabla v_i - w \cdot \nabla u_i v_i \dd x}_{ =: a(t; u, v)} \dd t \\ = \int_0^T \underbrace{\langle f, v \rangle_{V(t)} + \langle g, v \rangle_{H^{1/2}(\Gamma(t))}}_{ =: l(t; v)} \dd t, \end{aligned} \] with the initial condition \(u(0) = u_0 \in H_0\).

Well posedness

The variational form has a unique solution \(u\) that satisfies \[ \norm{u}_{L^2_V}^2 + \norm{\md_t u}_{L^2_{V^*}}^2 \le C \left( \norm{f}_{L^2_{V^*}}^2 + \norm{g}_{L^2_{H^{-1/2}(\Gamma)}}^2 + \norm{u_0}_{H_0}^2 \right). \]

Suppose further that \(f \in L^2_H\), \(g \in L^2_{H^{1/2}(\Gamma)}\) and \(u_0 \in V(0)\) then the solution \(u \in \{v \in L^2_V, \md_t v \in L^2_H\}\) and satisfies \[ \norm{u}_{L^2_V}^2 + \norm{\md_t u}_{L^2_{H}}^2 \le C \left( \norm{f}_{L^2_{H}}^2 + \norm{g}_{L^2_{H^{1/2}(\Gamma)}}^2 + \norm{u_0}_{V_0}^2 \right). \]

The finite element scheme

We construct a moving isoparametric finite element method of order \(k\).

Discretisation of the domain (i)

  • We assume we are given an initial (nice) piecewise linear triangulation \(\tilde{\T}_h\) of \(\Omega\) which approximates the domains at time \(t = 0\).
  • We denote by \(\tilde{\Gamma}_h\) the initial approximation of \(\Gamma_0\).

  • We construct a bijection \(\Psi_h\colon \Omega \to \Omega\) such that

    • \(\Psi_h|_{\tilde{\Gamma}_h} = p|_{\tilde{\Gamma}_h}\).
    • \(\Psi_h|_{K}\) is smooth (on the closed element).

\[ \Psi_h(x) = x + (\mu^*(x))^{\textcolor{green}{k+2}} (p(y(x), 0) - y(x)) \quad \mbox{ if } \mu^*(x) \neq 0 \]

Discretisation of the domain (ii)

The initial mesh \(F_{K}(\hat{x}, 0)\) is defined by the interpolation of \(\Psi_h\) of order \(k\)

Time dependent domain

  • We move each Lagrange node by the Lagrange interpolation (order \(k\)) of the smooth velocity \(w\).

  • We assume that this gives a good mesh.

  • On going work - how to generalise/relax this condition.

Lifted element

  • We construct a bijection between the computational domain and the exact domain. We call this the lift denoted by \(\Lambda_h \colon \Omega \to \Omega\).

  • We can also lift functions defined on the triangulation: \(v_h \mapsto v_h^\ell\).

  • The process also defines a lifted material derivative. (Dziuk & Elliott, 2012)

  • The procedure is similar to construction of the initial domain.

  • Important: The lift is the identity mapping for all elements away from the interface.

Finite element spaces

\[ \begin{aligned} S_h^i(t) := \Big\{ & \chi_h = (\chi_K)_{K \in \T_h^i(t)} \in \prod_{K \in \T_h^i(t)} \{ \hat\chi \circ F^{-1}_{K} : \hat\chi \in P_k(\hat{K}) \} : \\ & \chi_K(a(t)) = \chi_{K'}(a(t)) \mbox{ for all } K, K' \in \T^i(a(t)), \\ & \qquad \mbox{ for all } a(t) \in \N_h^i(t) \Big\} \end{aligned} \]

\[ \begin{aligned} S_h(t) := \big\{ & \chi_h = (\chi^1_h, \chi^2_h) \in S_h^1(t) \times S_h^2(t) : \\ & \chi_h^1(a(t)) = \chi_h^2(a(t)) \mbox{ for all } a(t) \in \N_h^\Gamma(t), \\ & \chi_h^2(a) = 0 \mbox{ for all } a \in \N_h^{\partial \Omega} \big\} \end{aligned} \]

Discrete problem

Given \(f_h\), \(g_h\), diffusion coefficients \(A_h^i\) and initial data \(U_{h,0}\), find \(U_h \in \{ v_h \in L^2_{S_h} : \md_h v_h \in L^2_{S_h}\}\) such that \[ \begin{aligned} \dt \int_{\Omega(t)} U_h v_h \dd x + \underbrace{\sum_i \int_{\Omega_{h}^i(t)} A_{h}^i \nabla U_{h}^i \cdot \nabla v_h - W_h \cdot \nabla U_h^i v_h \dd x}_{ =: a_h(t; U_h, v_h)} \\ = \underbrace{\int_{\Omega_h(t)} f_h v_h \dd x + \int_{\Gamma_h(t)} g_h v_h \dd \sigma_h}_{ =: l_h(t; v_h)} + \int_{\Omega_h(t)} U_h \md_h v_h \end{aligned} \] with the initial condition \(U_h(0) = U_{h,0}\).

Analysis of scheme

The Discrete problem has a unique solution \(U_h \in C^1_{S_h}\) and \[ \sup_{t \in [0,T]} \norm{U_h}_{H_h(t)}^2 + \int_0^T \norm{U_h}_{V_h(t)}^2 \dd t \le C(T) \norm{U_{h,0}}_{H_h(0)}^2. \]

\[ \begin{aligned} & \sup_{t \in [0,T]} \norm{u - U_h^\ell}_{H(t)}^2 + h^2 \int_0^T \norm{u - U_h^\ell}_{V(t)}^2 \\ & \qquad \le \norm{u_0 - U_{h,0}^\ell}_{H(t)}^2 + c h^{2k+2} C(u). \end{aligned} \]

Ideas of proof

  • We apply the ideas from (Elliott & TR, 2020).

  • We use the Ritz projection \(\Pi_h \colon V(t) \to S^h(t)\) given by \[ a_h(t; \Pi_h z, v_h) = a(t; z, v_h^\ell) \qquad \mbox{ for all } v_h \in S^h(t). \]

  • Split errors: \[ u - U_h^\ell = \rho + \theta = (u - (\Pi_h u)^\ell) + (\Phi_h u - U_h)^\ell. \]

  • Estimates for \(\rho\) and \(\theta\) both require geometric estimates, interpolation estimates, duality arguments.

  • The geometric errors compare a lifted problem (with lifted material derivative) to the continuous problem.

Geometric estimates (example) i

We want to show: \[ \tag{P1} \abs{m(t; \eta_h^\ell, v_h^\ell) - m_h(t; \eta_h, v_h)} \le c h^{k+1} \norm{\eta_h^\ell}_{V(t)} \norm{v_h^\ell}_{V(t)}. \]


\[ \begin{aligned} m(t; \eta, v) &= \sum_i \int_{\Omega^i} \eta^i v^i \dd x \\ m_h(t; \eta_h, v_h) &= \sum_i \int_{\Omega_h^i} \eta_h^i v_h^i \dd x. \end{aligned} \]

Geometric estimates (example) ii

\[ \tag{P1} \abs{m(t; \eta_h^\ell, v_h^\ell) - m_h(t; \eta_h, v_h)} \le c h^{k+1} \norm{\eta_h^\ell}_{V(t)} \norm{v_h^\ell}_{V(t)}. \]

  • Let \(J_h\) be the Jacobian \(\sqrt{\det( \nabla \Lambda_h^T \nabla \Lambda h)}\) then \[ \sup_{t \in [0, T]} \norm{J_h - 1}_{L^\infty(\Omega_{h,i}(t))} \le c h^k. \]

  • \[ \begin{aligned} & \abs{m(t; \eta_h^\ell, v_h^\ell) - m_h(t; \eta_h, v_h)} = \abs{\int_{\Omega} \eta_h^\ell v_h^\ell (J_h - 1)} \\ & \qquad \le c h^k \norm{ \eta_h }_{L^2(\Omega)} \norm{ v_h }_{L^2(\Omega)}. \end{aligned} \]

Geometric estimates (example) iii

  • Denote by \([J_j \neq 1]\) be the elements where \(J_h \neq 1\) then \[ \begin{aligned} & \abs{m(t; \eta_h^\ell, v_h^\ell) - m_h(t; \eta_h, v_h)} = \abs{\int_{\Omega} \eta_h^\ell v_h^\ell (J_h - 1)} \\ & \le c h^k \left(\int_{[J_h \neq 1]} (\eta_h^\ell)^2 \right)^{1/2} \left(\int_{[J_h \neq 1]} (v_h^\ell)^2 \right)^{1/2}. \end{aligned} \]

  • Applying the narrow band trace inequality(Elliott & TR, 2012) we get the result: For any \(\eta \in H^1(\Omega)\) \[ \norm{ \eta }_{L^2([J_h \neq 1])} \le c h^{1/2} \norm{ \eta }_{H^1(\Omega)}. \]

  • \[ \abs{m(t; \eta_h^\ell, v_h^\ell) - m_h(t; \eta_h, v_h)} \le c h^{k+1} \norm{ \eta_h }_{V(t)} \norm{ v_h }_{V(t)}. \]

Geometric estimates - what about \(a\) form?

\[ \tag{P2} \abs{a(t; \eta_h^\ell, v_h^\ell) - a_h(t; \eta_h, v_h)} \le c h^{k+1} \norm{\eta_h^\ell}_{Z_0(t)} \norm{v_h^\ell}_{Z_0(t)}. \]

  • We can apply the same ideas as for \(m\) bilinear forms to get order \(k+1\) estimate.

  • But this needs careful application to ensure that we can apply this for smooth functions.

  • We do this by using our Ritz projection so that (P2) is only applied to smooth functions.

Additional duality argument

  • Define the bilinear form \(b \colon V(t) \times V(t) \to \R\) by \[ \dt a(t; \eta, v) = a(t; \md_t \eta, v) + a(t; \eta, \md_t v) + b(t; \eta, v). \]

  • Let \(z \in Z_k(t)\) and \(\eta = z - (\Pi_h z)^\ell\). We need to show \[ \abs{b(t; \eta, \zeta)} \le c \left( \norm{\eta}_{H(t)} + h \norm{\eta}_{V(t)} + h^{k+1} \norm{z}_{Z_k(t)} \right) \norm{\zeta}_{Z_0(t)}. \tag{B3} \]

  • To show this we use:

Numerical results

  • Implemented using the high-level finite element package firedrake.

  • Allows quick implementation using ufl but with efficient, parallel implementation using PETSc.

  • Use appropriate order BDF scheme for time stepping (including mesh movement).

Convergence result


  • We have presented a numerical method and accompanying analysis for a moving coupled bulk-surface problem.

  • The method and analysis is based on the fundamental and abstract frameworks from (Elliott & TR, 2020).

  • Numerical results demonstrate the convergence.

  • Questions remain over the quality of the moving mesh.

Thank you for your attention!

C.M. Elliott & TR. A unified theory for continuous-in-time evolving finite element space approximations to partial differential equations in evolving domains, IMA Journal of Numerical Analysis, Volume 41, Issue 3, July 2021, Pages 1696–1845, https://doi.org/10.1093/imanum/draa062


Alphonse, A., Elliott, C., & Stinner, B. (2015) An abstract framework for parabolic PDEs on evolving spaces. Port. Math., 72, 1–46, https://doi.org/10.4171/pm/1955.
Douglas, J. & Dupont, T. (1973) Galerkin methods for parabolic equations with nonlinear boundary conditions. Numer. Math., 20, 213–237, https://doi.org/10.1007/bf01436565.
Dziuk, G. & Elliott, C.M. (2012) \(L^2\)-estimates for the evolving surface finite element method. Math. Comp., 82, 1–24, https://doi.org/10.1090/s0025-5718-2012-02601-9.
Elliott, C.M. & TR (2012) Finite element analysis for a coupled bulk-surface partial differential equation. IMA J. Numer. Anal., 33, 377–402, https://doi.org/10.1093/imanum/drs022.
Elliott, C.M. & TR (2020) A unified theory for continuous-in-time evolving finite element space approximations to partial differential equations in evolving domains. IMA J. Numer. Anal., 41, 1696–1845, https://doi.org/10.1093/imanum/draa062.