Evolving finite elements for advection diffusion with an evolving interface

Technical Summary: Evolving Finite Elements for Advection Diffusion with an Evolving Interface

Problem Statement
This paper addresses the numerical approximation of parabolic equations defined on a stationary domain $\Omega$ containing a moving interface $\Gamma(t)$. The interface divides the domain into two subdomains, $\Omega_1(t)$ and $\Omega_2(t)$. The governing equations model the transport and diffusion of scalar fields $u_1$ and $u_2$ within these subdomains, subject to:

• Advection-diffusion-reaction dynamics within each subdomain with time-dependent coefficients ($A_i, B_i, C_i$).
• A Dirichlet boundary condition on the outer boundary $\partial\Omega$.
• Continuity of the solution across the interface ($u_1 = u_2$).
• A prescribed flux jump across the interface ($A_1 \partial_{\nu} u_1 - A_2 \partial_{\nu} u_2 = g$).

The primary challenges identified are the evolution of the subdomains over time and the presence of the flux jump across the moving interface. While applications exist in fluid dynamics, materials science, and cell biology, the paper focuses on the mathematical and numerical formulation of the problem.

Methodology
The authors propose an Arbitrary Lagrangian-Eulerian (ALE) approach utilizing Evolving Finite Element Spaces (ESFEM) based on the abstract framework established in [21].

1. Functional Analytic Setting: The continuous problem is formulated within the context of evolving Hilbert spaces. The authors define compatible pairs of spaces $(X(t), \phi_t)$ where $\phi_t$ is a flow map transporting functions from a reference domain to the evolving domain. This allows for the definition of weak material derivatives and the derivation of a well-posed weak formulation (Theorem 2.16).
2. Mesh Construction:
   • Initial Mesh: The domain is initially triangulated using piecewise linear simplices. Crucially, the vertices of interface facets are constrained to lie exactly on the initial interface $\Gamma(0)$ (Assumption M5).
   • Isoparametric Mapping: To achieve higher-order accuracy, the linear mesh is mapped to an isoparametric mesh using a diffeomorphism $\Psi_h$. This map pushes the interface facets onto the exact curved interface $\Gamma(0)$ while preserving the continuity of the mesh.
   • Time Evolution: The mesh evolves by moving the Lagrangian nodes according to a known, smooth global velocity field $w$. The discrete flow $\Phi^h$ is constructed as a polynomial interpolation of the continuous flow, ensuring the mesh remains a conforming admissible subdivision at all times.
3. Lift Mapping: A lift operator $\Lambda_h$ is introduced to map functions from the discrete computational domain (which approximates the curved interface) back to the exact continuous domain. This is essential for comparing the discrete solution to the exact solution in error analysis.
4. Discrete Scheme: The method employs isoparametric finite elements of order $k$. The discrete variational problem is formulated using the discrete material derivative and the lifted coefficients. The scheme is solved using a time-stepping method (specifically Backward Difference Formulas, BDF, in the numerical experiments).

Key Contributions

• Well-Posedness: The paper provides a rigorous functional analytic setting proving the existence and uniqueness of the solution to the continuous problem (Theorem 2.16), including additional regularity results under specific assumptions on the data.
• Isoparametric ESFEM Framework: The authors extend the abstract evolving finite element framework to handle two-sided interface problems with flux jumps. They detail the construction of isoparametric elements that accurately approximate the curved boundary, enabling higher-order convergence.
• Optimal Error Bounds: The main theoretical result (Theorem 5.1) establishes optimal order error bounds. Specifically, the error in the $L^2$ norm is bounded by $O(h^{k+1})$, and the error in the $H^1$ semi-norm is bounded by $O(h^k)$, where $h$ is the mesh size and $k$ is the polynomial degree. These bounds hold up to a constant and match the interpolation error of a known smooth solution.
• Implementation and Verification: The paper includes a robust implementation using the Firedrake package and GMSH. Numerical experiments in 2D and 3D verify the theoretical convergence rates for various polynomial orders.

Results
The numerical experiments demonstrate the method's effectiveness on a model problem with a known exact solution that is continuous but not differentiable across the interface.

• Convergence: For spatial discretization of order $k$, the numerical results confirm convergence of order $k+1$ in the $L^2$ norm and order $k$ in the $H^1$ semi-norm.
• Robustness: The results hold for both 2D and 3D domains, even when the meshes are generated independently (rather than via uniform refinement), indicating the method's stability.
• Code Availability: The authors provide the simulation code to allow for reproducibility and further use by the community.

Significance and Scope
The paper claims to provide a complete theoretical and practical framework for solving advection-diffusion problems with moving interfaces using high-order isoparametric finite elements. The significance lies in bridging the gap between abstract evolving space theory and practical implementation for interface problems with flux jumps.

The authors maintain a modest scope regarding the velocity field: they assume a global, smooth velocity field is given a priori. They explicitly state that the work does not address how to compute such a velocity field (e.g., via re-parametrization to prevent mesh tangling), noting that such methods exist in the literature [8, 18, 19]. Furthermore, while the method could theoretically apply to coupled PDE systems where domain evolution is part of the solution, the current analysis does not cover such coupled problems. The focus remains strictly on the approximation of the parabolic transmission problem given a known domain evolution.