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

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.

## 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.$

## 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.

• 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 .

• 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)}.$

where

\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 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.

• 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:

• integration by parts;
• a nice boundary duality argument .

# 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).

# Summary

• 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 .

• Numerical results demonstrate the convergence.

• Questions remain over the quality of the moving mesh.

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.