Energy stable discretization of Allen-Cahn type problems modeling the motion of phase boundaries
We study the systematic numerical approximation of a class of Allen-Cahn type problems modeling the motion of phase interfaces. The common feature of these models is an underlying gradient flow structure which gives rise to a decay of an associated energy functional along solution trajectories. We first study the discretization in space by a conforming Galerkin approximation of a variational principle which characterizes smooth solutions of the problem. Well-posedness of the resulting semi-discretization is established and the energy decay along discrete solution trajectories is proven. A problem adapted implicit time-stepping scheme is then proposed and we establish its well-posed and decay of the free energy for the fully discrete scheme. Some details about the numerical realization by finite elements are discussed, in particular the iterative solution of the nonlinear problems arising in every time-step. The theoretical results are illustrated by numerical tests which also provide further evidence for asymptotic expansions of the interface velocities derived by Alber et al.
Keywords: Allen-Cahn equation, phase-field models, gradient systems, mean curvature flow, finite elements, implicit time stepping
AMS-classification (2000): 65M60,74N20,35J93,53C44
We consider the systematic numerical approximation of a class of Allen-Cahn type equations describing, for instance, the motion of anti-phase boundaries in crystalline solids , the evolution of interfaces in more general phase field models [6, 23, 36], or the geometric evolution by mean curvature [12, 24]. For motivation of our considerations, let us briefly consider two particular examples which have been studied in detail in [1, 2]. As a first model problem, we consider the generalized Allen-Cahn equation
Here is an order parameter which takes values close to zero in one phase of the medium and close to one in the other, is the double well potential, is a source term due to external loading, and , are model parameters. The coefficient plays the role of a regularization parameter that allows to widen the transition zone between the phases and thus to alleviate the numerical simulation. For the choice , , , and , one obtains the standard Allen-Cahn equation
which has been studied intensively in the literature; see e.g. [18, 26, 33] and the references given there. It is well-known  that (1.2), and therefore also (1.1), can be understood as a gradient flow for an associated free energy functional, which for the problem (1.1) has the form
The multiplicative term on the right hand side of (1.1) yields a scaling of time that leads to a non-trivial behavior for , which is called the sharp interface limit. It is shown in [1, 2] that the velocity of the interface can be expanded as
where can be expressed expicitly in terms of the problem data. Here denotes the jump of the Eshelby tensor [16, 17] which is related to the external forces and is the normal vector while is the mean curvature of the interface . In the limit , the velocity of the interface motion can thus be described by the well-known relation [26, 33]
Choosing the regularization parameter leads to an incorrect speed of the interface motion, and according to (1.4), one should choose to keep the consistency error small. For parameter the interface is widened to a transition zone of width . In order to resolve the smoothed interface in simulations, the mesh size should therefore be chosen at least as small as . More restrictive assumptions on the mesh size are required to really prove a good approximation of the interface evolution; see [19, 30] for details. If is small, then a very small mesh size is thus required to obtain a good approximation.
This equation, which has been termed hybrid model in [1, 2], describes the motion of the levelsets of the order parameter with normal speed proportional to the term in parenthesis. Following the arguments of [9, 13, 31], the problem (1.6) can again be interpreted as a gradient flow for an associated energy functional, which here reads
The gradient here has to be defined with respect to a metric induced by the solution itselve; we refer to [9, 13] and to Section 2 for details. Using asymptotic analysis, the propagation speed of the interface has been shown in [1, 2] to behave like
where the parameter can again be computed from the model data.
A comparison of the two formulas (1.4) and (1.8) reveals that for a proper choice of the model and scaling parameters, the hybrid model leads to the same propagation speed for the interface as in the Allen-Cahn model up to higher order terms. To obtain a good agreement of the results, one should choose and in this case width of the smoothed interface in the hybrid model (1.6) can be shown to behave like ; see [1, 2] for details. When is small, the smoothed interface of the hybrid model is thus substantially larger then that of the corresponding Allen-Cahn model and a mesh size of should be sufficient to resolve the interface in computations, instead of required for the corresponding Allen-Cahn model. The hybrid model should therefore allow to compute the interface evolution on a much coarser mesh, which was the main motivation for the proposal of this method.
Due to the close connection to mean curvature flow, there has been much interest in the numerical approximation for the Allen-Cahn equation (1.2). Various discretization schemes have been proposed and analyzed, see e.g. [7, 35] and the references given there. Particular aspects of the efficient numerical approximation have been addressed in [20, 21, 34] and a-posteriori error estimates have been derived in [5, 22, 28]. The convergence of numerical approximations of the Allen-Cahn equation in the sharp interface limit towards the generalized motion by mean curvature has been established in [19, 30]; see also  for further numerical tests. The direct numerical approximation of mean curvature flow has also been investigated intensively; see for instance [11, 37] and refer to [12, 25] for a review about available results and further references.
The focus of the current manuscript is somewhat different to these previous works: Here we want to study the systematic numerical approximation of a large class of Allen-Cahn type models which includes (1.1), (1.2), and (1.6) as special cases. The common feature of all these models is the underlying gradient flow structure with respect to an associated energy functional, and we aim to strictly preserve this structure in the approximation process on the semi-discrete and the fully discrete level. With this in mind, we automatically obtain discretization schemes that are uniformly energy stable and thus consistent with the second law of thermodynamics.
In this paper we establish well-posedness and qualitative properties of a discretization strategy appropriate for a large class of Allen-Cahn type problems. The models (1.1) and (1.6) introduced above will be considered as particular examples. As a by-product of our investigations, we also obtain further numerical evidence for the asymptotic expansions (1.4) and (1.8) of the interface velocities. Our main arguments can be extended to elastic Allen-Cahn type models which describe the evolution of phase interfaces in binary alloys [2, 23, 24]. The required coupling to the additional elasticity system will be investigated in a forth-comming publication.
The remainder of the manuscript is organized as follows: In Section 2, we first introduce a unified formulation of the problems discussed in the introduction and then present a variational characterization of smooth solutions. In addition, we verify that the respective energies decrease along solution trajectories which is due to the underlying gradient flow structure. In Section 3, we then discuss the Galerkin approximation of the underlying variational principles. We establish the well-posedness of the resulting semi-discrete schemes and verify the energy decay on the discrete level. In Section 4, we then discuss the subsequent discretization in time by a problem adapted implicit time stepping scheme. We establish the well-posedness of the fully discrete problem and again prove the decay in energy. A particular space discretization by finite elements is discussed in Seciton 5, and some further aspects of the implementation are presented. For illustration of our theoretical results, we report in Section 6 about some numerical tests,in which we also verify the asymptotic expansions for the interface velocities given in (1.4) and (1.8). We close with a short discussion of our results and open problems for future research.
2. A unified variational framework
For the rest of the manuscript, we consider the following general model problem:
Here , is a bounded Lipschitz domain, are model parameters, and the functions , are given a-priori. For our analysis later on, we will assume that
with , , and for all .
Some of these technical assumptions are made only for convenience and could further be relaxed without problems. This will become clear from our analysis below and we will make some remarks in this direction. To complete the model description, we additionally assume that
Other boundary conditions can also be treated without much difficulty. The above problem is of quasilinear parabolic type and existence of solutions for appropriate initial values can be proven with standard arguments; see e.g. [4, 29, 32].
Here denotes the space of continuous functions that are continuously differentiable with respect to and twice continuously differentiable with respect to . Note that acording to our definition, smooth solutions and their spatial derivatives are uniformly bounded.
The Allen-Cahn model (1.1) can be phrased in the form (2.1) by setting , , , and defining the potential as . Note that the second derivative of the potential is unbounded here. Using maximum principles, the solution can however be shown to remain bounded and therefore could be modified for large in order to satisfy (A2) without changing the solution; see  for other appropriate conditions on the potential. The hybrid model (1.6) can also be cast into the form (2.1) with , , , and . To avoid problems for vanishing gradient, one can set for some instead, and thus obtain a regularized version of (1.6). The limit may be studied in the framework of -convergence [10, 27].
for all and . The last term is the negative directional derivative of the energy functional in direction . These formal arguments motivate the following variational principle.
Lemma 2.3 (Variational characterization).
As a direct consequence of the gradient flow structure, we obtain
3. Galerkin approximation in space
Motivated by the results of the previous section, we now consider the Galerkin approximation of the variational principle (2.6) to construct semi-discretizations in space. Let us choose some finite dimensional subspace and consider the following discrete variational problem.
Problem 3.1 (Galerkin semi-discretization).
Let and be given. Find such that and
for all test functions and for all .
Since we employ a conforming Galerkin approximation, the gradient flow structure and thus also the energy decay estimates are inherited automatically by the semi-discrete problem.
Let denote a solution of Problem 3.1. Then
The assertion follows with literally the same arguments as in Lemma 2.4. ∎
As a next step, let us establish the well-posedness of the semi-discrete problem.
Lemma 3.3 (Well-posedness).
For any and , Problem 3.1 has a unique solution.
After choosing a basis for , we can rewrite the discrete problem as initial value problem
where here denotes the coordinate vector when expanding in the chosen basis. Due to assumptions (A1)–(A2) on the coefficients and since is finite dimensional, the matrix functions and can be seen to be Lipschitz continuous. From the lower bounds on we further obtain that is symmetric positive definite. Existence of a unique local solution then follows from the Picard-Lindelöf theorem. From (2.3) and (3.2), we can see that
Here we used that due to assumption (A2). This yields and by the growth condition and the positivity assumption for , we further obtain
By the Poincaré inequality and positivity of , we thus get with constant only depending on the domain, on the bounds in the assumptions, and on the energy of the initial value. This implies that remains uniformly bounded and hence the solution can be extended uniquely to arbitrarily large . ∎
Under assumptions (A1)–(A2), we thus obtain global existence of a unique semi-discrete solution . The constant in the a-priori bound only depends on the initial energy , on the values of and , and on the lower bound for the growth rate of the potential . Let us note that the energy decay estimate also provides a uniform bound for the time derivative of the solution which might be useful for a full convergence analysis for the Galerkin approximation; we refer to [7, 19] for results in this direction.
4. Time discretization
As a second step in the approximation process, let us now discuss the discretization in time. Given some time step , we define for , and we denote by
the backward difference quotient at time . As will become clear from our analysis, non-uniform time steps could be considered as well. For the time discretization of the Galerkin approximation defined in Problem 3.1, we then consider the following implicit time-stepping scheme.
Problem 4.1 (Fully discrete scheme).
Set with initial value as in Problem 3.1. For any , find such that
for all . We set wherever .
Let us note that a similar treatment of the nonlinear term has been proposed in the context of mean-curvature flow in . The particular of the discrete variational problem (4.1) allows us to establish a decay estimate for the energy of the fully discrete approximation with very similar arguments as used for the analysis on the continuous level.
Let denote a solution of Problem 4.1. Then
Using definition (2.3), we can decompose the discrete energy as
We now utilize the fully discrete variational principle with to replace the first and last term together by which directly yields the result. ∎
The second term on the right hand side of the above energy estimate stems from numerical dissipation of the implicit time-stepping scheme. The particular time discretization thus enhances the stability of the numerical solution. As further theoretical backup, we next establish the well-posedness of the fully discrete scheme.
Let be given. Then (4.1) has at least one solution , and the solution is unique for all time steps with some sufficiently small.
Existence of a solution follows from the energy estimate given in the previous section and the Brouwer fixed point theorem, and we know that . Now let and be two solutions of (4.1). Then testing with and using assumption (A1) leads to
For ease of notation, we utilized the symbol
here to abbreviate the difference quotient. By the fundamental theorem of calculus and using the bounds of assumption (A2) for the coefficient , we further get
A combination of the two estimates then directly leads to
For small enough, the right hand side can be absorbed by the terms on the left side and we thus also obtain the uniqueness of the solution. ∎
Before we proceed, let us briefly indicate at this point how the assumption (A1)–(A2) could be possibly relaxed, which can be deduced by inspection of the previous proof.
The bounds for are not required for all here, but only for values that are actually attained during the evolution. If a bound is available, which may be proved by a discrete comparison principles , then one could replace the constant in the above proof by . If the space admits an inverse inequality
then the a-priori estimate implies with . The constant can be estimated explicitly for many discretization schemes, e.g., for finite element methods which discussed in the following section. In one space dimension, the above inverse inequality holds true for any choice of due to the continuous embedding of into . More general conditions on the potential have also been considered in .
5. Details on the implementation
Before we turn to numerical tests, let us briefly discuss some details of the implementation that will be used in our computations later on. This particularly involves the choice of the approximation spaces and the solution of the nonlinear systems (4.1) in every time step.
5.1. A finite element method
as the space of continuous piecewise polynomials of degree . In our numerical tests we will actually only utilize the lowest order case , but the methods can in principle be extended to higher order without difficulty. Similar approximations in space have been considered in [7, 14].
5.2. Nonlinear potential
Without loss of generality, one may assume that the potential is a polynomial function. Otherwise, can be replaced by a polynomial and due to the Weierstrass theorem, the error of this approximation can be made arbitrarily small on bounded intervals. In our numerical tests, we will utilize a function of the form
The first part is the usual double-well potential and the linear term amounts to external forces. Note that the potential does not strictly satisfy the assumptions (A2), in particular, the upper bounds on the second derivative are violated. As noted in Remark 4.4, this assumption was made for convenience and is not required as long as stays uniformly bounded. For the above choice of the potential , the difference quotient arising in (4.1) can be expressed as
The nonlinear problem (4.1) that has to be solved in every time step of the fully discrete scheme thus has a polynomial nonlinearity and is non-degenerate. Standard iterative methods can therefore be utilized for an efficient solution.
5.3. Iterative solver
The following fixed-point iteration will be used in our numerical tests for the solution of the nonlinear problem that arises in every time step of the fully discrete scheme.
Problem 5.1 (Fixed-point iteration).
Let and , be given and define . Further set and for , find such that
for all . If convergence is reached after iteration , set .
Choosing allows to improve the convergence behavior also in case ; see  for a similar argument. As an immediate consequence of the particular construction, we obtain
For any , the above fixed-point iteration is well-defined. If is sufficiently small, then the iteration converges to the unique solution of Problem 4.1.
First note that the problem that has to be solved in every step of the fixed-point iteration of Problem 5.1 is linear. Moreover, the left hand side in the above iteration defines a positive definite quadratic form. This implies the existence of a unique solution in every step of the iteration. For sufficiently small, one can show that the map is contractive, and convergence follows by the Banach fixed-point theorem. ∎
For time step sufficiently small, the contraction constant in the above iteration becomes small and convergence will typically be observed within a few iterations. To speed up convergence, one may alternatively also utilize Newton-type iterations.
6. Numerical illustration
For our computations, we consider the Allen-Cahn equation (1.1) and the hybrid model (1.6). The model parameters are set to in all tests, and the double-well potential is defined as . As potential for the external forces, we choose with constant to be specified below. The constant in the Allen-Cahn equation is chosen as . For this setting, the interface velocities of Allen-Cahn and the hybrid model given in (1.4) and (1.8) agree up to higher order terms; see [1, 2] for details.
For the space discretization, we utilize the finite element method described in Section 5 with polynomial degree . The time integration is performed by the implicit time-stepping scheme (4.1) with constant time step , and the nonlinear systems arising in every time step are solved with the fixed-point iteration outlined above.
6.1. A quasi one-dimensional problem
As a first test case, we consider a quasi one-dimensional geometric setting. We choose and set
where is the characteristic function of the set . The interface between the two phases at time is a straight line here and thue the curvature is zero. This can be shown to remain true for all by a symmetry argument. The evolution of the interface is therefore only driven by external forces in this example. For our computations, we choose with which yields a constant driving force acting on the whole domain.
As mentioned above, the solution can be shown to be independent of for all . We thus obtain with defined by
Here is the one-dimensional cross-section at arbitrary and . The initial condition for the one-dimensional problem is . In Figure 6.1, we display a few snapshots of the solution obtained with the Allen-Cahn equation (1.1) in two dimensions and the corresponding solutions and for the Allen-Cahn problem and the hybrid model in one space dimension.
The results clearly illustrate the relation between the solutions of the one- and the two-dimensional problem. Moreover, the propagation of interface for the one-dimensional Allen-Cahn and the corresponding hybrid model take place at the same speed, which is in perfect agreement with formulas (1.4) and (1.8). In Figure 6.2, we display the area surrounded by the interface and we depict the decay of the energy functional for the two models illustrating the underlying gradient flow structure.
6.2. Shrinking of a circle by mean curvature
As a second test case, we consider the evolution of a circular interface . We again choose and define the initial values as
Here we set which yields and thus the motion is only driven by the curvature of the interface. The solution of the sharp interface limit can be shown to be radially symmetric and according to (1.5), the propagation of the interface is governed by
here denotes the radial position of the interface at time . This allows to compute the area surrounded by the interface analytically, i.e.,
A brief visual inspection shows that the motion of the interface is rather similar for both models. One can further deduce from the plots that the width of the transition zone in the hybrid model is somewhat larger than for the Allen-Cahn model; compare with our discussion in the introduction. A coarser mesh might therefore be used for the simulation of the interface evolution by the hybrid model. For a better comparison of the two models, we again display in Figure 6.4 the evolution of the area and of the energy representing one phase of the system.
6.3. A non-smooth and non-convex geometry
As a last test case, we now consider the evolution of an interface which initially is non-smooth and encloses a non-convex region at time . Some snapshots of the evolution are depicted in Figure 6.5.
Like in the previous test, we set the force potental to which implies , and hence the motion is again driven only by the curvature. Note that in the vicinity of convex corners, the interface moves towards the interior of the enclosed area while at the non-convex corners, the interface propagates in the other direction, which is caused by the change of sign in the curvature. In Figure 6.6, we again depict the area surrounded by the interface and the uniform decay of the energy for the numerical solutions obtained with the two models.
Like in the previous examples, the velocity of the interface motion in the Allen-Cahn and the hybrid model are very similar; compare with the formulas (1.4) and (1.8). Due to the gradient flow structure of the fully discrete evolution, a strict decrease in the energy can again be observed. The results are very similar, and the slight discrepancies can be explained by approximation errors.
In this paper, we investigated the systematic numerical approximation of a general class of Allen-Cahn type equations. The common feature of these models was a gradient flow structure with respect to an associated energy. We proposed and analyzed semi-discrete and fully discrete numerical schemes which strictly preserve the underlying gradient flow structure and which therefore automatically yield uniformly energy stable discrete approximations. Well-posedness of the numerical schemes and energy decay was established theoretically and illustrated by numerical tests. Our computational results also provide further numerical evidence for asymptotic expansions of interface motion obtained [1, 2].
In our analysis, we used some rather strong structural assumptions on the coefficients arising in the equations, which should allow to conduct a full a-priori convergence analysis of the semi- and fully discrete scheme by following the arguments of [7, 19]. As indicated in remarks, a detailed inspection of the proofs allows to further relax these assumptions. Our arguments and methods of proof are rather general and seem to be applicable also to phase-field models for the dynamics of solid-solid phase transition which are governed by coupled Allen-Cahn and elasticity equations. These topics and further extensions are left for future research.
This work was supported by the German Research Foundation (DFG) via grants IRTG 1529, TRR 154, and Eg-331/1-1 and by the German Excellence Initiative via grant GSC 233.
-  H.-D. Alber. Asymptotics and numerical efficiency of the Allen-Cahn model for phase interfaces with low energy in solids. arXive, 1505.05442, 2015.
-  H.-D. Alber and P. Zhu. Comparison of a rapidely converging phase field model for interfaces in solids with the Allen-Cahn model. J. Elasticity, 111(2):153–221, 2013.
-  S. M. Allen and J. W. Cahn. A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening. Acta Metall., 27:1085–1095, 1979.
-  H. Amann. Linear and quasilinear parabolic problems. Vol. I, volume 89 of Monographs in Mathematics. Birkhäuser Boston, Inc., Boston, MA, 1995.
-  S. Bartels. A posteriori error analysis for time-dependent Ginzburg-Landau type equations. Numer. Math., 99:557–583, 2005.
-  J. F. Blowey and C. M. Elliott. Curvature dependent phase boundary motion and parabolic double obstacle problems. In Degenerate diffusions (Minneapolis, MN, 1991), volume 47 of IMA Vol. Math. Appl., pages 19–60. Springer, New York, 1993.
-  X. Chen, C. M. Elliott, A. Gardiner, and J. J. Zhao. Convergence of numerical solutions to the Allen-Cahn equation. Appl. Anal., 69(1-2):47–56, 1998.
-  P. G. Ciarlet. The finite element method for elliptic problems. North-Holland Publishing Co., Amsterdam-New York-Oxford, 1978.
-  U. Clarenz, F. Haußer, M. Rumpf, A. Voigt, and U. Weikard. On level set formulations for anisotropic mean curvature flow and surface diffusion. In Multiscale modeling in epitaxial growth, volume 149 of Internat. Ser. Numer. Math., pages 227–237. Birkhäuser, Basel, 2005.
-  E. De Giorgi. New problems in -convergence and -convergence. In Free boundary problems, Vol. II (Pavia, 1979), pages 183–194. Ist. Naz. Alta Mat. Francesco Severi, Rome, 1980.
-  K. Deckelnick and G. Dziuk. Convergence of a finite element method for non-parametric mean curvature flow. Numer. Math., 72:197–222, 1995.
-  K. Deckelnick, G. Dziuk, and C. M. Elliott. Computation of geometric partial differential equations and mean curvature flow. Acta Numer., 14:139–232, 2005.
-  M. Droske and M. Rumpf. A level set formulation for Willmore flow. Interfaces Free Bound., 6:361–378, 2004.
-  Q. Du and R. A. Nicolaides. Numerical analysis of a continuum model of phase transition. SIAM J. Numer. Anal., 28:1310–1322, 1991.
-  A. Ern and J.-L. Guermond. Theory and practice of finite elements, volume 159 of Applied Mathematical Sciences. Springer-Verlag, New York, 2004.
-  J. D. Eshelby. The elastic field outside an ellipsoidal inclusion. Proc. Roy. Soc. London. Ser. A, 252:561–569, 1959.
-  J. D. Eshelby. Elastic inclusions and inhomogeneities. In Progress in Solid Mechanics, Vol. II, pages 87–140. North-Holland, Amsterdam, 1961.
-  L. C. Evans, H. M. Soner, and P. E. Souganidis. Phase transitions and generalized motion by mean curvature. Comm. Pure Appl. Math., 45:1097–1123, 1992.
-  X. Feng and A. Prohl. Numerical analysis of the Allen-Cahn equation and approximation for mean curvature flows. Numer. Math., 94:33–65, 2003.
-  X. Feng, H. Song, T. Tang, and J. Yang. Nonlinear stability of the implicit-explicit methods for the Allen-Cahn equation. Inverse Probl. Imaging, 7:679–695, 2013.
-  X. Feng, T. Tang, and J. Yang. Stabilized Crank-Nicolson/Adams-Bashforth schemes for phase field models. East Asian J. Appl. Math., 3:59–80, 2013.
-  X. Feng and H.-j. Wu. A posteriori error estimates and an adaptive finite element method for the Allen-Cahn equation and the mean curvature flow. J. Sci. Comput., 24:121–146, 2005.
-  E. Fried and M. E. Gurtin. Dynamic solid-solid transitions with phase characterized by an order parameter. Phys. D, 72:287–308, 1994.
-  H. Garcke. On mathematical models for phase separation in elastically stressed solids. PhD thesis, 2000. Habilitation thesis.
-  H. Garcke. Curvature driven interface evolution. Jahresber. Dtsch. Math.-Ver., 115:63–100, 2013.
-  T. Ilmanen. Convergence of the Allen-Cahn equation to Brakke’s motion by mean curvature. J. Differential Geom., 38:417–461, 1993.
-  H. Y. Jian. A relation between -convergence of functionals and their associated gradient flows. Science in China, Ser. A, 42, 1999.
-  D. Kessler, R. H. Nochetto, and A. Schmidt. A posteriori error control for the Allen-Cahn problem: circumventing Gronwall’s inequality. M2AN Math. Model. Numer. Anal., 38(1):129–142, 2004.
-  O. A. Ladyženskaja, V. A. Solonnikov, and N. N. Ural’ceva. Linear and quasilinear equations of parabolic type. Translated from the Russian by S. Smith. Translations of Mathematical Monographs, Vol. 23. American Mathematical Society, Providence, R.I., 1968.
-  R. H. Nochetto and C. Verdi. Convergence past singularities for a fully discrete approximation of curvature-driven interfaces. SIAM J. Numer. Anal., 34:490–512, 1997.
-  S. Osher and N. Paragios. Geometric Levelset Methods in Imaging, Vision and Graphics. Springer, 2003.
-  A. Pazy. Semi-groups of linear operators and applications to partial differential equations. Department of Mathematics, University of Maryland, College Park, Md., 1974.
-  J. Rubinstein, P. Sternberg, and J. B. Keller. Fast reaction, slow diffusion, and curve shortening. SIAM J. Appl. Math., 49:116–133, 1989.
-  J. Shen, T. Tang, and J. Yang. On the maximum principle preserving schemes for the generalized Allen-Cahn equation. Commun. Math. Sci., 14:1517–1534, 2016.
-  J. Shen and X. Yang. Numerical approximations of Allen-Cahn and Cahn-Hilliard equations. Discrete Contin. Dyn. Syst., 28:1669–1691, 2010.
-  A. Visintin. Models of phase transitions. Progress in Nonlinear Differential Equations and their Applications, 28. Birkhäuser Boston, Inc., Boston, MA, 1996.
-  N. J. Walkington. Algorithms for computing motion by mean curvature. SIAM J. Numer. Anal., 33(6):2215–2238, 1996.
-  J. Zhang and Q. Du. Numerical studies of discrete approximations to the Allen-Cahn equation in the sharp interface limit. SIAM J. Sci. Comput., 31:3042–3063, 2009.