An unfitted discontinuous Galerkin scheme for conservation laws on evolving surfaces††thanks: The work was partially supported by the German Research Foundation (DFG), grant EN 1042/4-1 and by the UK Engineering and Physical Sciences Research Council, grant EPSRC EP/J004057/1.
Motivated by considering partial differential equations arising from conservation laws posed on evolving surfaces, a new numerical method for an advection problem is developed and simple numerical tests are performed. The method is based on an unfitted discontinuous Galerkin approach where the surface is not explicitly tracked by the mesh which means the method is extremely flexible with respect to geometry. Furthermore, the discontinuous Galerkin approach is well-suited to capture the advection driven by the evolution of the surface without the need for a space-time formulation, back-tracking trajectories or streamline diffusion. The method is illustrated by a one-dimensional example and numerical results are presented that show good convergence properties for a simple test problem.
Key words. Discontinuous Garlerkin; unfitted finite elements; surface partial differential equations; conservation laws.
AMS subject classifications. 35L02, 35L65, 35Q90, 35R01, 65N12, 65N30
Interest in the study of partial differential equations on evolving surfaces has grown in recent years with applications in many areas including materials science (e.g. diffusion of species along grain boundaries ), fluid dynamics (e.g. surface active agents along the interface between two fluids ) and cell biology (e.g. cell motility involving the processes on the cell membrane ), for example. See the review of  for a more detailed list. In this work, we derive a new numerical scheme for an essential conservation law on evolving hypersurfaces and present first numerical results.
As a model problem, we consider the evolution of a conserved quantity on an evolving curve or surface. Fix and let be a time-dependent, connected, compact, smooth -dimensional hypersurface embedded in for , denote a field of unit normal vectors to , and write . We consider two different descriptions of the surface. First, we say is defined by a diffeomorphic parametrization so that for . The map defines a field which describes the velocity of , precisely by . Second, we describe as the zero level set of a smooth function so that . We will use the diffeomorphic parametrization to define the equations we solve and the level set function to define the geometry in our computational scheme.
Consider a control volume which is the image of a control volume under the flow , i.e. . Let denote a time-dependent field on which satisfies a conservation law of the form
where is a tangential flux on and is the co-normal vector to ( is normal to and ). Applying a transport formula to the left hand side of LABEL:eq:balance and the divergence theorem for hypersurfaces to the right hand side yields
where is the material derivative of and denotes the tangential divergence operator. See LABEL:sec:prelim for details. We describe equation LABEL:eq:local-cons as a local conservation law for or a global conservation law in the case .
Since LABEL:eq:local-cons holds for any , we arrive at a pointwise conservation law. We wish to find a time-dependent surface field with
For more details on the notation, see .
In this work, we aim to derive a new numerical method using an implicit representation of the surface for which we can show a discrete analogue to the global conservation law. To solve problems of this form, we are required to consider two separate problems. We have to discretize the surface partial differential equation at each given time step and we have to evaluate the transport from the old surface to the new surface. Either these two components are computed in a splitting approach and discretized separately, or computed directly in a coupled manner. In either case, we must understand the effects of the flux and the advection of the moving surface separately. Previously (see the following section), attention has been paid to the stationary case and the treatment of the flux . In this work, we are interested in the advection driven by the moving surface. We thus restrict to the case and a single time step, which simplifies considerations. Extensions to the truly time-dependent case and will be considered in future work.
1.2 Previous approaches
Previous computational approaches lie in two categories. The first is based on a moving mesh approach and the second is based on a static mesh.
For surface PDEs, the basic ideas of the moving mesh approach date back to the ideas of the finite element method for elliptic PDEs on stationary surfaces which was presented in . A smooth surface is approximated by a union of elements (usually simplices but also possibly quadradrilaterals) whose vertices move according to a globally defined smooth velocity. It was applied by , using a finite element discretization for an advection-diffusion problem, and by , using a finite volume method for conservation laws. Moving mesh schemes have proven to be useful in many practical situations and can be shown analytically to be stable and accurate, or preserve various conservative properties. Their main disadvantage is that the moving meshes may degenerate, leading to frequent remeshing . Extensions have been presented recently by [15, 16] that reduce the need for remeshing by allowing mesh points to move with a different non-physical velocity. However, in many situations preserving a good mesh remains very challenging.
In the static mesh approach, an arbitrary fixed background mesh is used and the evolution of the surface is defined implicitly. One approach is to use a level set function together with an unfitted finite element method where the computational domain consists of partial cut-cell elements [6, 20, 7, 22]. The authors of  apply a space-time formulation to produce a stable, accurate method, but the method is only approximately globally conservative and requires the use of complicated four-dimensional space-time elements which are not available for geometries other than simplices. An alternative method proposed by  uses a semi-Lagrangian formulation that employs a non-local right hand side. The resulting scheme recovers a globally mass conservation property and performs well in practice for a test problem on a curve, but lacks analysis and the non-local term is expensive to compute, especially in three space dimensions. A phase field representation of the interface was used by [14, 23]. The authors of  use a narrow band formulation and recover a discrete analogue to the continuous level global conservation law but must use stream line diffusion to stabilize the scheme.
A key difficulty to overcome with the static mesh approach is to derive a stable scheme that adequately treats the advection driven by the evolution of the surface. The moving mesh approach implicitly deals with this problem by using moving basis functions. The advective flux usually has a component orthogonal to the surface. Therefore, any other fluxes (e.g. diffusive surface fluxes) can not stabilize the method. Furthermore, the previous ideas based on finite element methods are not well-suited to advection-dominated problems. Extra complications such as space-time formulations, semi-Lagrangian terms or streamline diffusion are required. In this work, we use a more suitable discontinuous Galerkin (DG) method which can naturally handle advection or advection-dominated problems.
2 Computational approach
Our approach is based on reformulating problem LABEL:eq:problem as a sequence of bulk advection problems with singular source and sink terms and applying the unfitted discontinuous Galerkin (UDG) method  on a static bulk mesh. This mesh can be chosen independent of the evolving surface, the surface is not explicitly represented by the mesh. Our formulation guarantees a globally conservative scheme.
Let be a polygonal domain in which contains for all times . We suppose that for each time , is defined as the zero level set of a smooth function and assume for , . Let be the field of unit normal vectors to oriented by . The normal component of the velocity field is given by .
The tangential gradient of a smooth surface field can be defined as
where is a spatially differentiable extension of to a neighborhood of and is its gradient with respect to the ambient Cartesian coordinates. The tangential gradient has components . For a smooth surface vector field , the tangential divergence is given by
The material derivative of a smooth, time-dependent surface field is defined as
where again and are Cartesian derivatives of a differentiable extension of to a space-time neighborhood of . The material derivative describes the variation of with respect to the evolution of the surface. It also has a key role in the following transport relation for evolving material surfaces, see e.g. [5, 10]:
Proceeding formally, we multiply problem LABEL:eq:problem (with ) by a smooth function , integrate over for and apply transport relation LABEL:eq:transport to see
Since does not depend on time, we have .
We fix a time and small enough that . Integrating in time over the time interval yields
We will approximate the space-time integral in the third term by a scaled bulk integral over the spatial-only domain
which corresponds to the projection of the space-time domain to spatial-only coordinates. An example of the geometric setup is shown in LABEL:fig:geometry.
Denoting an extension of a function by , we have
where we employ the right-hand rectangle method, an approximation of the resulting surface integral using the level sets of and apply the coarea formula, while assuming that the time step size as well as with and are small. Note that the assumption on is reasonable for a small time step size. If is constant and is a signed distance function, i.e. , corresponds to the travel distance of an individual point on the surface.
In order to be able to discretize using flux-based numerical schemes, such as unfitted DG schemes, we use integration by parts to reformulate this integral via the divergence of an advective flux tested with . As an approximate reformulation of problem LABEL:eq:problem over the time interval , we therefore obtain the following stationary problem: seek with
for all smooth functions . Here and in the following, we use the notation , and , and denotes the outward pointing unit normal vector field to .
In the special case (e.g. a geometrical setup as in LABEL:fig:domain-ex), problem LABEL:eq:problem_reformulated is consistent with the following bulk advection problem in strong form:
2.3 Unfitted discrete geometry
In order to discretize problem LABEL:eq:problem_reformulated, we start by discretizing the geometry. Let be a shape regular decomposition of into closed elements, either tetrahedra or hexahedra for , triangles or quadrilaterals for , and denote by the maximum element size. Let be the space of piecewise linear (for simplices), bilinear (for quadrilaterals) or trilinear (for hexahedra) continuous functions over and denote by interpolation of functions in into . We will write for a discrete level set function and set
Note that since is a closed set. We also use the notation
An example is shown in LABEL:fig:domain-ex.
It will be useful to consider the restriction of the decomposition to the computational domain :
We note that the are arbitrary shaped elements and call those elements cut cells. In general, they are not either shape regular or even convex. We will also consider the internal faces of the cut cell mesh , denoted by . The set is often called the skeleton of the mesh. To each internal face , which is the intersection of two elements , we assign a unit normal vector field and arbitrarily choose .
Our solution will be defined over the union of elements in , but we are only interested in the values of our solution variable over the sharp interface approximation . The integrals in the method will be computed over either the sharp interface approximation or the unfitted cut cell mesh . This is a similar approach to , but different to , thus avoiding difficulties in defining our discrete spaces and constructing a natural basis.
We introduce the discrete spaces given by
where denotes the space of piecewise polynomials of degree over the volumetric domain . Note that functions in are discontinuous and do not take a unique value along the faces , in general. On each face with adjacent elements as defined above, we define the jump of a function as and its average as .
We discretize problem LABEL:eq:problem_reformulated by approximating by , the test functions by , by and by a discrete approximation (which is given in more detail in the sequel). Integrals over and are approximated by integrals over and , respectively, and integrals over are approximated by integrals over . Finally, we split the bulk integral into a sum of integrals over the cut cells and integrate by parts on each . Using classical upwind stabilization (as in ), we obtain the following unfitted DG scheme:
Given , find , such that
where denotes the upwind solution of on which is given by
Note that if a section of or is part of a face , we choose as and as , respectively.
We note that, since this construction implies that is an admissible test function, a global mass conservation law is recovered:
There are many different options for the choice of . In this work, we will simply use the continuous normal velocity defined by the level set function and assume there is no tangential component to the velocity field. It is possible to use a backward difference to extract a normal velocity from the two discrete level set functions and by
3 Understanding the method in one dimension
In this section, we want to give a better interpretation of the method and its limitations, using a simple 1D example. We consider a point moving along a one-dimensional axis with positive speed , i.e. , . Furthermore, we consider a time step with .
Let be the interval and be its decomposition into sub-intervals of width , with . Let and , such that and . Discretizing using Scheme 1 and a piecewise constant discrete space , i.e. polynomial degree , we obtain a finite volume type scheme and equation LABEL:eq:scheme simplifies to
By fixing the basis of which consists of characteristic functions , , we obtain a discrete system. Denoting the vector of unkowns by and supposing that is a given scalar at , we wish to find with
|This system is uniquely solvable, yielding a piecewise constant solution with|
For implementation reasons one might want to compute on a domain larger than , say a computational domain . For the domain in LABEL:fig:domain-ex, for example, is not easily reconstructed without using knowledge of all intermediate curves. As we will illustrate now, the construction of the scheme unfortunately does not immediately carry over to this case. Unique solvability requires that resolves as a domain boundary. Furthermore, if does not meet this requirement, we cannot guarantee mass conservation any more.
We take a slightly larger domain with the same mesh size and extend the domain by — to the left and to the right. This means that and , i.e. and lie in the inner part of and , respectively. The system now changes to
|which is underdetermined, due to the pure Neumann boundary conditions. Analogous to pseudo time stepping, we regularize the system with an additional mass term of order , the limit yields the solution|
For , the solution is not well defined and diverges in the –limit. For , we observe that mass conservation is not fulfilled anymore and for the special case we obtain . Thus it is necessary to resolve with the computational domain .
4 Numerical results
The presented scheme has been implemented in the DUNE framework [1, 2] using the dune-UDG library . Besides providing unfitted DG spaces like the discrete space from LABEL:sec:method, dune-UDG enables the evaluation of integrals over cut cells and their faces , which is needed for the assembly of system matrices in unfitted DG schemes.
The fundamental mesh (see LABEL:sec:discrete_geometry) used in our code is a structured, Cartesian, quadradrilateral mesh over a freely choosable domain . In the following numerical experiments, we use the two-dimensional domain , decomposed into an equal number of quadrilaterals in - and -direction. Furthermore, we employ spaces of polynomial degree . With respect to an a priori known analytical solution on , we compute errors , and .
4.1 A circle shrinking with constant normal speed
We consider a circle of initial radius , centered at the origin, which is shrinking with constant normal velocity . Using the description as a level set function, can be described as the zero level set of , , , . Note that is a signed distance function, i.e. . Considering a time step with , the spatial-only domain from LABEL:sec:motivation takes the form . Its discrete reconstruction for different values of can be seen in LABEL:fig:annulus_solution_h_refinement.
Results for -refinement using , on and a corresponding analytical solution on are shown in the top row of LABEL:fig:annulus_solution_h_refinement and in LABEL:tab:annulus_solution_h_refinement. Note that #ncells is the number of cells in the fundamental mesh . It does not refer to the number of cut cells . Furthermore, since the scheme is globally mass conservative up to machine precision, the mass written in LABEL:tab:annulus_solution_h_refinement is both the mass of on and the mass of on . Note that it is an indicator for the geometrical error coming from the reconstruction of , which determines the amount of mass entering the discrete system. For the continuous problem and , the total mass in the system equals . We observe convergence of order in the -norm and the -norm, limited by geometrical errors for coarse meshes. For the -norm, the convergence rate is not that clear but also seems to be approaching order for small values of .
4.2 Non-constant initial concentration
To illustrate the effect of a non-constant , we compute the same setup as in LABEL:subsec:annulus_solution_h_refinement but use a with a binary distribution. As the worst case scenario, we consider value for and else. In this scenario, the numerical diffusion is expected to have the largest impact and for small values of the jump sharpens, see bottom row of LABEL:fig:annulus_solution_h_refinement. The numerical diffusion can be further reduced by using higher-order methods with flux-limiters.
We presented a new approach to solve the advection problem driven by the evolution of an evolving surface. The method shows attractive properties which justify further investigation. It is mass conservative and relatively easy to implement. We avoid constructing space-time meshes or following characteristics by reformulating the original problem as a classical transport problem on an unfitted domain.
In future work, we plan to extend the method to allow computations on and high-order shape functions. Furthermore, we will apply the method to truly time-dependent problems and will also consider more general equations with .
-  P. Bastian, M. Blatt, A. Dedner, C. Engwer, R. Klöfkorn , M. Ohlberger, and O. Sander. A generic grid interface for parallel and adaptive scientific computing. Part I: Abstract framework. Computing, 82(2–3):103–119, 7 2008.
-  P. Bastian, M. Blatt, A. Dedner, C. Engwer, R. Klöfkorn , R. Kornhuber, M. Ohlberger, and O. Sander. A generic grid interface for parallel and adaptive scientific computing. Part II: Implementation and tests in DUNE. Computing, 82(2–3):121–138, 7 2008.
-  P. Bastian and C. Engwer. An unfitted finite element method using discontinuous galerkin. International Journal for Numerical Methods in Engineering, 79(12):1557–1576, 2009.
-  Peter Bastian, Christian Engwer, Jorrit Fahlke, and Olaf Ippisch. An unfitted discontinuous galerkin method for pore-scale simulations of solute transport. Mathematics and Computers in Simulation, 81(10):2051 – 2061, 2011. MAMERN 2009: 3rd International Conference on Approximation Methods and Numerical Modeling in Environment and Natural Resources.
-  Paolo Cermelli, Eliot Fried, and Morton E. Gurtin. Transport relations for surface integrals arising in the formulation of balance laws for evolving fluid interfaces. Journal of Fluid Mechanics, 544:339–351, 12 2005.
-  Klaus Deckelnick, Gerhard Dziuk, Charles M. Elliott, and Claus-Justus Heine. An -narrow band finite-element method for elliptic equations on implicit surfaces. IMA Journal of Numerical Analysis, 30(2):351–376, 2010.
-  Klaus Deckelnick, Charles M. Elliott, and Thomas Ranner. Unfitted finite element methods using bulk meshes for surface partial differential equations. SIAM Journal on Numerical Analysis, 52(4):2137–2162, 2014.
-  Klaus Deckelnick, Charles M Elliott, and Vanessa Styles. Numerical diffusion-induced grain boundary motion. Interfaces and Free Boundaries, 3(4):393–414, 2001.
-  Gerhard Dziuk. Finite elements for the Beltrami operator on arbitrary surfaces. In Stefan Hildebrandt and Rolf Leis, editors, Partial Differential Equations and Calculus of Variations, volume 1357 of Lecture Notes in Mathematics, pages 142–155. 1988.
-  Gerhard Dziuk and Charles M. Elliott. Finite elements on evolving surfaces. IMA Journal of Numerical Analysis, 27(2):262–292, 2007.
-  Gerhard Dziuk and Charles M. Elliott. Finite element methods for surface PDEs. Acta Numerica, 22:289–396, 2013.
-  Gerhard Dziuk, Dietmar Kröner, and Thomas Müller. Scalar conservation laws on moving hypersurfaces. Interfaces and Free Boundaries, 15(2):203–236, 2013.
-  Carston Eilks and Charles M. Elliott. Numerical simulation of dealloying by surface dissolution via the evolving surface finite element method. Journal of Computational Physics, 227(23):9727–9741, 2008.
-  Charles M. Elliott, Björn Stinner, Vanessa Styles, and Richard Welford. Numerical computation of advection and diffusion on evolving diffuse interfaces. IMA Journal of Numerical Analysis, 31(3):786–812, 2011.
-  Charles M. Elliott and Vanessa Styles. An ALE ESFEM for solving PDEs on evolving surfaces. Milan Journal of Mathematics, 80(2):469–501, 2012.
-  Charles M Elliott and Chandrasekhar Venkataraman. Error analysis for an ale evolving surface finite element method. Numerical Methods for Partial Differential Equations, 31(2):459–499, 2015.
-  C. Engwer and F. Heimann. Dune-udg: A cut-cell framework for unfitted discontinuous galerkin methods. In Proceedings of the DUNE User Meeting 2010, 2011.
-  Ashley J. James and John Lowengrub. A surfactant-conserving volume-of-fluid method for interfacial flows with insoluble surfactant. Journal of Computational Physics, 201(2):685–722, December 2004.
-  Matthew P. Neilson, John A. Mackenzie, Steven D. Webb, and Robert H. Insall. Modeling Cell Movement and Chemotaxis Using Pseduopod-Based Feedback. SIAM Journal on Scientific Computing, 33:1035–1057, 2011.
-  Maxim A. Olshanskii, Arnold Reusken, and Jörg Grande. A finite element method for elliptic equations on surfaces. SIAM Journal on Numerical Analysis, 47(5):3339–3358, 2009.
-  Maxim A. Olshanskii, Arnold Reusken, and Xianmin Xu. An eulerian space-time finite element method for diffusion problems on evolving surfaces. SIAM Journal on Numerical Analysis, 52(3):1354–1377, 2014.
-  Maxim A. Olshanskii, Arnold Reusken, and Xianmin Xu. A stabilized finite element method for advection–diffusion equations on surfaces. IMA Journal of Numerical Analysis, 34(2):732–758, 2014.
-  Knut Erik Teigen, Xiangrong Li, John Lowengrub, Fan Wang, and Axel Voigt. A diffuse-interface approach for modeling transport, diffusion and adsorption/desorption of material quantities on a deformable interface. Communications in mathematical sciences, 4(7):1009, 2009.