Phase field models, i.e., diffuse interface models, have been a successful tool in studying interfacial dynamics . Due to the important applications in physics, biology, material science and image processing [14, 28, 15, 25, 50, 7, 8], there are substantial interests in developing efficient numerical methods for phase-field models [39, 30, 29, 21, 34, 53, 55, 65, 52, 51].
From a modeling perspective, phase-field models can be classified into two categories, known as Allen-Cahn type and Cahn-Hilliard type . The Allen-Cahn type models are typical examples of gradient flows , while the Cahn-Hilliard type models, which are concerned with a conserved quantity, are examples of -diffusion . Although the numerical methods for both types of phase-field models are well developed [39, 53, 52, 24], most of them are Eulerian methods, which solve the equation of “phase” function in a fixed grid . In order to resolve the thin diffuse interface, one must have mesh sizes much smaller than the width of the thin diffuse interface [59, 45, 24], which requires huge computational efforts. This difficulty is often handled by using adaptive mesh techniques [49, 1] or moving mesh approaches [29, 21, 53, 55].
The main purpose of this paper is to propose a variational Lagrangian methods for Allen-Cahn type phase-field models. The approach used here can be extended to other -gradient flow systems. Compared with Eulerian methods, Lagrangian methods, which are often self-adaptive, have potential advantages for problems involving singularity, sharp interface and free boundary. Recently, there has been an increasing interest in applying Lagrangian schemes for generalized diffusions, such as the porous medium equation and nonlinear Fokker-Plank equations [12, 64, 13, 37, 11, 43, 10, 40]. However, it is more difficult to construct Lagrangian schemes for the -gradient flows, like Allen-Cahn equation. Unlike the generalized diffusions, which have natural variational structures on the Lagrangian maps [27, 37, 10, 40], the variational structures of
-gradient flows are on the physical variables defined on Eulerian coordinates. Moreover, as a drawback of all Lagrangian methods, the meshes of Lagrangian solutions may become too skew, which not only influence the accuracy of the solution, but also may result in premature termination of Lagrangian calculations.
In order to overcome these difficulties, we first propose an energy-dissipation law for a phase-field model, given by
where is a “phase” function satisfying a transport equation
is the effective velocity associated with the Lagrangian map, and is the mixture energy density. This model is inspired by phase-field model of mixture of two incompressible fluids [68, 36, 1]. By employing an energetic variational approach, we can obtain the corresponding PDE of this system, given by
subject to suitable initial and boundary conditions. Formally, it is to reformulate (1) and (3) in terms of a Lagrangian map and its time derivative, and we can construct a variational structure preserved Lagrangian methods by a discrete energetic variational approach  based on (1).
For , this model employs the same energy-dissipation law of Allen-Cahn type models. So one can view (1) as a modified Allen-Cahn type model. The additional term in the dissipation part of (1) can be viewed as a regularization term on Lagrangian maps, which plays an essential role in calculations. It is obvious that such modification would only change the dynamics of the system when approaching the equilibrium states, which are also equilibria of the original Allen-Cahn type model. Within the energy-dissipation law (1), we can employ a discrete energetic variational approach  to construct a variational Lagrangian scheme which inherits the variational structure.
The rest of this paper is organized as follows. We first give a detailed description to our phase-field model in the next section. Then we construct our variational Lagrangian scheme by the discrete energetic variational approach in Sect. 3. Several numerical results to validate our methods are shown in Sect. 4.
2 Model development
2.1 Energetic Variational Approach
An energetic variational approach, originated from pioneering work of Onsager [46, 47] and Rayleigh  provides a general framework to determine the dynamics of system from a prescribed energy-dissipation law through two distinct variational processes: Least Action Principle (LAP) and Maximum Dissipation Principle (MDP) [38, 32]. During the last decade, this approach has been successfully applied to build up many mathematical models in physics, chemistry and biochemistry [28, 38, 56, 26, 32, 61].
For an isothermal closed system, an energy-dissipation is given by
which is the consequence of the first and second laws of thermodynamics . Here is the total energy, which is the sum of the Helmholtz free energy and the kinetic energy , and is the rate of energy dissipation. The Least Action Principle states that the equation of motion for a Hamiltonian system can be derived from the variation of the action functional with respect to the flow maps (the trajectory in Lagrangian coordinates) if applicable [4, 32], i.e.,
It gives a unique procedure to derive the conservative forces for the system. On the other hand, for a dissipative system (), the dissipative force can be obtained by minimization of the dissipation functional with respect to the “rate” in the regime of linear response , known as Onsager’s Maximum Dissipation Principle (MDP), i.e.,
Then, the force balance condition () results in
which is the dynamics of the system. We refer the reader to  for more detailed descriptions of energetic variational approaches. We only consider systems without kinetic energy, i.e. , throughout this paper.
2.2 Energetic variational approaches to phase-field models
From an energetic variational viewpoint, Allen-Cahn and Cahn-Hilliard type of models, provide a dynamics to minimize the mixture energy
for in some admissible set subject to some boundary conditions on . Here,
is a “phase” function that introduced to identify the two phases, is the mixture energy density given by
where is the interfacial (potential) energy that often taken as a double-well potential
Different phase-field models can be derived by different choices of admissible sets and dissipation functional .
In Allen-Cahn type models, is often chosen to be with a suitable boundary conditions, and the energy-dissipation law is given by
where is the constant of the dissipation rate . We take in the following. The energy-dissipation law (12) can be viewed as a gradient flow dynamics on the phase function . According to the general framework of an energetic variational approach, the corresponding gradient flow equation can be derived by first performing LAP with respect to and MDP with respect to
where we assume all boundary terms vanish due to the given boundary condition. Then the force balance equation (7) leads to a Allen-Cahn type equation
The stationary solutions of the Allen-Cahn type equations, satisfy the Euler-Lagrangian equation of the functional (8)
The above derivation performs an energetic variational approach in terms of and . We call this is an Eulerian approach, in which we can view as generalized coordinate in the system . There is an alternative way to derive a dynamic of the system, known as the Lagrangian approach . Instead of studying the evolution of phase function directly, the Lagrangian approach study the evolution of a Lagrangian map, or flow map, for given initial condition . For fixed , is a diffeomorphism between the initial domain and the current domain , known as a deformation map [58, 33]. For fixed , is the trajectory of the particle labeled by . We can view are Lagrangian coordinates and are Eulerian coordinates.
We can define the “velocity” in Eulerian coordinate, as
for the given flow map . Another important quantity associated with is the deformation tenors , defined by
which carries all the information about how the physical quantity transport with the flow. Since is one-to-one map between and for fixed , we can enforce , which means the map is orientation-preserving for .
In order to get the equation of , we shall impose the kinematic relation to the physical quantity . Then the dynamics of will be totally determined by the dynamics of the flow map . For Allen-Cahn type model, it is often assumed that satisfies
where is the initial condition. One can view (18) as a composition between and inverse flow map at time , that is
From the kinematic equation (18), we have
Hence, satisfies scalar transport equation
in Eulerian coordinates.
The above transport relation (21) is the macroscopic transport on the microscopic variable , which might only be valid locally. The complicated phase evolution, such as interface merging or pinching off, which is the consequences of microscopic evolution of , cannot be described by this kinematic.
Within the kinematic (18), is determined by for given . Hence, we can propose a energy-dissipation law in terms of and to characterize the local dynamics of the system, that is
and is the rate of energy dissipation on , i.e., velocity . The energy-dissipation law (22) can be viewed as a generalized gradient flow of the flow map . Since we are only concerned with equilibria of the system, the choice of dissipation only effects the dynamics approaching to equilibria. We’ll discuss this later.
The evolution equation of the flow map can be derived by employing an energetic variational approach, that is
where [See Appendix for the detailed computation]
A stationary solution in the Lagrangian approach satisfies
The Lagrangian approach minimizes the mixture energy in the admissible set
which is smaller than the Eulerian approach. So it is subtle to choose a suitable to get a desired equilibrium. For classical Allen-Cahn equation, due to the maximum principle, we have , it is not difficult to choose a proper . In general, can be obtained by some Eulerian approach. We can also update during the evolution of the flow map.
If is a conserved quantity that satisfies
then the kinematic equation is given by
This is the kinematic for the Cahn-Hilliard type equation, which can be viewed as a generalized diffusion with the energy-dissipation law given by 
Both Allen-Cahn and Cahn-Hilliard equations types are driven by the same mixture energy 8, but the kinematic and dissipation mechanisms are different.
Although the equations for the stationary solution obtained by Eulerian approach (variation on the phase variable ) and the Lagrangian approach (variation with respect to the domain) looks different [(15) and (26) ], formally one can easily show that :
For a given energy functional (8), all smooth (regular enough) solutions of the Euler-Lagrangian equation:
also satisfy the equation
This result indicates connection between variation with respect to and the variation with respect to flow map through Legendre transform . Theoretically, all equilibria in the Eulerian approach can be obtained from the Lagrangian approach with a proper choice . However, for a given , the Lagrangian approaches and Eulerian approaches may not end up with the same equilibrium.
2.3 Dissipation Functional
In this subsection, we discuss the choice of dissipation for Lagrangian approach to phase-field models. The different choices of dissipation provide different dynamics approaching equilibria of the system. Since we may have multiple equilibria for the free energy like (8) , different dynamics may end up with different equilibria for given .
for given initial condition . The equation of the flow map can be obtained via a standard energetic variational approach, which is
In a recent work , the authors study numerical methods for equation (34) in one-dimension by discretizing directly. Their results show that the dynamics (34) can capture the thin interface of Allen-Cahn type equation with a small number of mesh points in 1D. However, the energy-dissipation law (33) may not be suitable for Lagrangian calculations, especially for high dimensions . Indeed, since is a rank one matrix,
is not a invertible matrix for, so is not well-defined everywhere. Moreover, even for one-dimensional cases, is almost zero in non-interfacial regions, which restricts the choice of .
The degeneracy of the equation (34) motivates us to consider a different dissipation functional by adding a new term, that is
where is a constant. By a direct computation, for such an energy-dissipation law, the dynamics of the system is given by
which gives us the equation of the flow map in Lagrangian coordinates. The energy-dissipation law (35) fix the degeneracy of . Moreover, from a computational perspective, can be viewed as a regularization term to the flow map , which control the quality of mesh generated by the flow map. In Lagrangian coordinates, (35) can be written as
in Lagrangian coordinates. In order to simplify the numerical numerical implementation, in the following, we replace by . Then the equation for flow map is given by
subject to the initial condition , where for and .
It is worth mentioning that the additional terms in both (35) and (37) are not physically unacceptable viscosity for compressible fluids. Both of them conflict with the frame-indifference. More specifically. let
According to the frame-indifference, we should have
hence, it is easy to show that the additional terms in both (35) and (37) do not satisfy the frame-indifference. For compressible flow, a physically acceptable viscosity in the dissipation is often taken as
At the end of this section, We should emphasize the above derivation is rather formal, in which we assume the flow map exists at least locally. The goal of this paper is designing some Lagrangian schemes that preserve the above variational structures in a discrete level. More analysis are certainly need to show the existence of flow map. We refer the interesting reader to [27, 18] for some theoretical results on some related but different systems.
3 Numerical Scheme
In this section, we construct our variational Lagrangian scheme for the phase-field model with the energy-dissipation law (1) by a discrete energetic variational approach . Instead of considering a particular weak form of the flow map equation (36), a discrete energetic variational approach, which performs energetic variational approaches in a semi-discrete level, derives a “semi-discrete equation” that preserves the variational structure from a discrete energy-dissipation law directly. By introducing a proper temporal discretization to the “semi-discrete equation”, we can construct an energy stable Lagrangian scheme to our phase-field model.
3.1 A discrete energetic variational approach
In general, for a system without kinetic energy, a discrete energy-dissipation law can be written as
where is the “discrete” state variable, is the discrete free energy and is the discrete dissipation. One can obtain a discrete energy-dissipation law (42) from the continuous energy-dissipation law by either discretizing physical quantity (Eulerian approach) or the flow map (Lagrangian approaches) in space.
Similar to an energetic variational approach in a continuous level, the governing equation of , a system of nonlinear ODEs, can be obtained from the force balance equation
where the right-hand side comes by performing LAP, taking variation of the discrete action functional with respect to , while the left-hand side comes by performing MDP, taking variation of the discrete dissipation functional with respect to .
The discrete energetic variational approach follows the strategy of “discrete-then-variation”, which has been a powerful tool to construct numerical schemes for complicated systems with variational structures [31, 17, 12, 13, 66, 40]. Compared with the traditional “variation-then-discrete” approach, the “semi-discrete” equation obtained by the discrete energetic variational approach can automatically inherit the variational structure from the continuous level. One may obtain the same “semi-discrete” equation in “variation-then-discrete” approach by choosing a particular weak form for the PDE.
In order to get a discrete energy-dissipation law for our phase-field model, we first introduce a piecewise linear approximation to the flow map . A piecewise linear approximation to the flow map can be constructed by a standard finite element method. In the following, we only discuss the two-dimensional case, the procedure can be easily extended to other dimensions. Let be a triangulation of domain , consists of a set of simplexes and a set of nodal points . Then the approximated flow map is given by
and is the hat function satisfying . Since , can be viewed as the coordinate of -th mesh points at , and defines the velocity of -the mesh point. Within the above spatial discretization, the discrete state variable of the system is defined by
where . For simplicity’s sake, we consider the natural boundary condition for the flow map through this section. The specific boundary conditions can be incorporated to numerical schemes easily by modifying velocity of the mesh points on the boundary.
The framework of finite element discretization enable us to compute the deformation matrix explicitly on each element [see the Appendix in  for the explicit form]. We denote the deformation matrix on each element , which is a constant matrix for fixed by . The admissible set of is defined by
It can be noticed that is not a convex set, which imposes difficulties in both simulation and numerical analysis.
Inserting (44) into the original energy-dissipation law, we obtain the discrete free energy
and the discrete dissipation functional
By a discrete energetic variation approach, we can derive a system of ordinary differential equations of, that is
We refer readers to the Appendix for the detailed computation of and the . Although the explicit forms of both and may not be available in a general mesh, both of them are easy to obtain during the numerical implementation by summing the results on each element over the mesh. It can be noticed that given by
Here () is the modified mass matrix defined by
where is the centroid of , and is the modified stiff matrix defined by
It is easy to show that and is a positive-define matrix if . Hence, the presence of term ensures that is positive-definite.
3.2 Temporal discretization
Now we discuss the temporal discretization. A numerical scheme can be obtained by introducing a suitable temporal discretization to the “semi-discrete equation” (50). An advantage of existing a variational structure in the semi-discrete level is that various of classical numerical schemes can be reformulated as optimization problems [24, 65, 44]. In the current study, we use implicit Euler for temporal discretization. It is not difficult to apply high-order temporal discretization, such as BDF2 or Crank-Nicolson  to our system, which will be studied the future work.
For given , the implicit Euler scheme for (50) is given by
Although (54) is a system of highly nonlinear equations that is often difficult to solve, by virtue of the variational structures in the semi-discrete level, we can reformulate (54) into an optimization problem, given by
There are various of advantages in solving optimization problem (55) instead of solving the original nonlinear system (54) directly. Since might not be a convex function, solving (54) by standard nonlinear solvers, such as fixed-point iterations or Newton-type methods, may only obtain a saddle point or a local minimizer of , which may not decrease the discrete energy. Moreover, for our problem, it is also difficult to guarantee the solution of (54) obtained by some nonlinear solver is in the admissible set . For the optimization problem (55), we can use some line-search based optimization method and set
. By doing this, we can guarantee the decay of the discrete energy and will be satisfied automatically, even though the exact global minimizer of may not be found.
For given , there exists a solution to numerical scheme (54) such that the following discrete energy dissipation law holds, i.e.,
Since for , , following the proof in the Lemma 3.1 in , the existence of a minimizer can be obtained by showing the set
is a non-empty compact subset of . Obviously, , so is non-empty. On the other hand, since is positive-definite, there exists such that
which indicates is bounded. So we only need to show is closed in . For any converged sequence , our goal is to show that the limit is in . Note for and all
where is the area of element . Since if , we can conclude that is uniformly bounded away from zero. So for all , which means .
If is a global minimizer of in , we have
which completes the proof.
Theorem 2 indicates that our scheme is energy stable if one can find a minimal solution of (55). However, since defined in (56) is a non-convex functional, we still need to choose a small value of and large value of such that the optimization problem can be handled by a standard optimization method. Indeed, the first term in (56) can be viewed as a regularization term, which restricts us to find a minimizer around . So it is important to have the positive-definite condition on , otherwise, may have infinite minimizer even around . In the numerical implementation. In all numerical experiments shown in the next section, we adopt L-BFGS with line search to find a minimizer in the admissible set that decreases the discrete energy. Within the energy stable property, we can prove the convergence of series for the discrete scheme for the given triangulation and fixed .
For the given triangulation and fixed , the series converges to a stationary solution of the discrete energy .
We first prove that there exist such that
Since , we only need to show, for , there exist such that
Indeed, note , we have
following the same argument in the proof of the previous theorem, we can show that (), which is uniformly bounded away from zero. Hence,
is the smallest eigenvalue of the stiff matrix.
In the numerical implementation, we can compute by
which is equivalent to set after each iteration as in . An advantage of this treatment that in each iteration, we only need to compute a close to identity map . So the optimization problem (55) is often easy to solve.
One can view this as a reinitialization procedure. More complicated reinitialization procedure can be incorporated in our numerical framework. Indeed, for given , we also obtain the numerical solution defined at mesh points, that is
When the mesh become too skew, we can interpolate the numerical solutioninto a more regular mesh, obtained by coarsening or refining the current mesh . More importantly, we can also apply Eulerian solver by using as the initial condition, to update the value at each mesh point. This is close to the idea in velocity-based moving mesh method , which update both positions and values of mesh points. Unlike the traditional velocity-based moving mesh methods, our solution is spontaneously updated when the mesh moves. We’ll explore reinitialization procedures in details in the future work.
4 Numerical validation and discussion
In this section, we apply our Lagrangian scheme to several Allen-Cahn type phase-field models. Most of numerical examples used here are widely studied by Eulerian methods previously [15, 29, 53, 21, 69]. Numerical results show that, with a suitable initial conditions, our methods can capture the thin diffuse interface with a small number of mesh points, and reach a desired equlibrium.
4.1 Quasi-1D example
First, we consider a quasi-1D problem, in which . We impose Dirichlet boundary condition on and , that is
and Neumann boundary condition on and , that is
The corresponding boundary condition for the flow map is
where . We take the initial condition as
Typical meshes and computed solution for and at and are shown Fig. 1 (a)-(c). We compare the obtain the equilibrium solution with the 1D exact solution in Fig. 1 (d), in which the circles represent the projection of mesh points in the x-z plane, and the red line is the exact solution. It can be noticed that the equilibrium numerical solution can capture the thin interface with small number of mesh points. Due to the present of term in the dissipation, the vertical velocity of all mesh points are almost zero, which is essential for a successful Lagrangian computation. Comparing with the 1D simulation results in , we can achieve a better approximation to the exact solution since we take rather than . The in the dissipation part enables us to use such an initial condition.
Fig. 1 (e) shows the the discrete energy as a function of time for different values of . One can noticed that our scheme is energy stable in all cases and all calculations go to the same equilibrium when is large enough. The convergence of numerical solution becomes slower when become larger. On the hand, numerical tests show that the optimization problem (55) in each iteration will be easier to solve when is larger. In general, the optimal value of is problem dependent, and the value of also effect the equality of the obtained mesh.
4.2 Shrinkage of a circular domain
As a numerical test, we consider shrinkage of a circular domain in two-dimension. It is a classical benchmark problem for Allen-Cahn equation [15, 29, 53], in which the circular interface governed by the Allen-Cahn equation will shrink and eventually disappear.
We take and impose the Dirichlet condition . The initial condition is taken as
so the Dirichlet condition satisfies numerically. It is worth pointing out that in our Lagrangian methods, it is crucial to choose a proper initial condition. For the phase model, it is often choose in a hyperbolic tangent form such that , and the width of initial interface shoud be larger than the mesh size, since we need enough mesh points in the region of interface.
Fig. 2 (a) shows the numerical results for with at various time in a uniform mesh (), while Fig. 2 (b) and (c) show the numerical results for with in uniform () and non-uniform meshes () respectively. The non-uniform mesh is generated by DistMesh . We choose larger for smaller to control the quality of the mesh. It can be noticed that in all three cases, the mesh points will be concentrated at thin interface after one time iteration and maintain concentrated at the moving interface all the time. The results in Fig. 2 (c) suggests that we can incorporated our Lagrangian method with adaptive mesh technique. Within the Lagrangian solver, we only need to adaptive the initial mesh. As a limitation, the Lagrangian calculation cannot reach the equilibrium, in which the circular domain is disappear. Such a problem can be handled easily by applying some Eulerian solver to the numerical solutions obtained by Lagrangian calculations at the late stage.
4.3 Phase-field model with the volume constraint
In this subsection, we consider an Allen-Cahn type phase-field model with the volume constraint. We impose the volume constraint by introducing a penalty term in the free energy. So the total free energy of the system is given by
We take , , and , and impose the Dirichlet boundary condition , throughout this section.
Fig. 3 (a) shows numerical results for initial condition
in which we use a non-uniform mesh generated by the DistMesh  As expected, due to the effect of surface tension and the volume constraints, the bubble deforms into a circular bubble, and the mesh points keep concentrated at the thin interface when the shape of interface changes. As a benefit of pure Lagrangian calculation, we can guarantee the numerical solution .
We also consider the initial condition