{'fraction_non_alphanumeric': 0.03790203432208206, 'fraction_numerical': 0.013126491646778043, 'mean_word_length': 4.478985056039851, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
First-Order System Least Squares and the Energetic Variational Approach for Two-Phase Flow

This paper develops a first-order system least-squares (FOSLS) formulation for equations of two-phase flow. The main goal is to show that this discretization, along with numerical techniques such as nested iteration, algebraic multigrid, and adaptive local refinement, can be used to solve these types of complex fluid flow problems. In addition, from an energetic variational approach, it can be shown that an important quantity to preserve in a given simulation is the energy law. We discuss the energy law and inherent structure for two-phase flow using the Allen-Cahn interface model and indicate how it is related to other complex fluid models, such as magnetohydrodynamics. Finally, we show that, using the FOSLS framework, one can still satisfy the appropriate energy law globally while using well-known numerical techniques.


Introduction
Complex fluids involve multi-physics and multi-scale phenomena that require advanced techniques in scientific computing in order to be solved efficiently. The goal of this paper is to present numerical techniques that have been well tested and shown to give accurate solutions with a reasonable amount of computational cost, while also preserving the global energy of the system. While the systems can be described with partial differential equations (PDEs), in this case a time-dependent and nonlinear system of equations, an underlying energy law also describes the physics. It is from these laws that the PDE is derived. Therefore, approximate solutions to complex fluid problems should adequately approximate this energy law.
Here, we look at one type of complex fluid problem, two-phase flow. The conventional sharp interface description of the mixture of two Newtonian fluids with density, ρ i , and viscosity, µ i , assumes the mixture occupies the overall domain, Ω = Ω 1 ∪ Ω 2 , as in figure 1. The interface between two materials is the free moving interface, Γ t . The overall system includes the Navier-Stokes equation in each Ω i , where i = 1 or 2: Here, u i , p i , ρ i , and µ i are the velocities, pressures, densities, and viscosities of each fluid, respectively. Two types of boundary conditions are possible. Kinematic boundary conditions yield u i ·n = V n , where n is the "directional" normal of Γ t and V n is the normal speed of Γ t . Traction free boundary conditions (also known as force balance conditions) satisfy [τ i ] · n = surface force on Γ t , where [τ i ] is the jump of the stress tensor, Remark Assume u 2 | ∂Ω = 0. With the traction free boundary condition, this implies Γ t ∩ ∂Ω = ∅.
In case the surface force is attributed to the surface tension, i.e., F surface = kHn, where H is the mean curvature of Γ t and k is the surface tension parameter, then the whole system possesses the underlying energy law, Ωi In this paper, we take a different approach, called the diffusive interface method (phase-field method). The label function, φ(x, t), is introduced such that φ(x, t) = +1 in fluid 1, −1 in fluid 2.
The interaction between the two materials is reflected through the "mixing" energy, or the Ginzburg-Landau type of energy: where Ω 1 2 |∇φ| 2 dx represents the tendency of the mixture to be homogeneous. In other words, it is the "phillic" interaction. In addition, Ω 1 4 2 φ 2 − 1 2 dx represents the "phobic" interaction, or the tendency of the mixture to be separated. The parameter reflects the competition between the two opposite tendencies. Due to the nature of this empirical energy, it is present in many theories in physics, such as superconductivity [1,2] and phase transitions [3,4]. It can be shown that φ → ±1 almost everywhere as → 0, thus representing the separation of the two immiscible materials. Moreover, the thermodynamics and patterns of the interface, Γ t , can be captured by the level sets of φ, for instance, {x|φ(x, t) = 0}. In addition, the energy Ω W (φ) → c 0 Area(Γ t ), which is equivalent to the surface tension case in the sharp interface formulation [5]. It is worth pointing out that, while one can view the diffusive interface formula as being the regularizing approximation of the sharp interface formulation, it is more physical to view the sharp interface formulation as being the ideal approximation of the diffusive interface method. Therefore, it can be understood that the phase field description really captures the microscopic interaction of the mixtures.
In the examples presented here, we consider problems where we have one or more bubbles of one phase interacting in a region of a fluid of the other phase. From this system, the continuous energy laws are derived as is the continuous system of PDEs. We then show that these laws are adequately preserved numerically.
To solve the system, we use the first-order system least-squares (FOSLS) finite element approach. This method has been well-tested on many different types of PDEs and, specifically, on several complex fluid systems [6,7,8,9,10], including those of magnetohydrodynamics (MHD) [11,12]. We show that use of the FOSLS method allows both an accurate solution of the system and preservation of the energy laws. Moreover, because we are using this method, we have the advantage of being able to use other numerical techniques, such as nested iteration and adaptive local refinement, which greatly improve the efficiency of solving these problems. This paper is structured as follows. In Section 2, we describe the two-phase model and present the equations to be solved. In section 3, we discuss the FOSLS formulation that is used and show how it presents a well-posed discrete problem. Section 4 introduces the energy laws and discusses what quantities are preserved with the numerical schemes being used. Numerical results for certain bubble test problems and an MHD simulation are presented in section 5, where it is shown that this approach is accomplishing the intended goal. Finally, a discussion of the results and summary is given in section 6.

The Two-Phase Equations
The dynamics of mixing between two fluids plays an important role in many physical applications. Here, we use a model that couples the Navier-Stokes equations with a phasefield transport equation. We use the Allen-Cahn-type equations, which were originally introduced in [13] to describe the motion of anti-phase boundaries in crystalline solids. 3 The Allen-Cahn equations have since been used in many interface problems in fluid dynamics through a phase-field approach [3,4,5,14,15]. They are Here, u is the fluid velocity, p is the fluid pressure, and φ is the phase-field function. The physical parameters are the mixing energy parameter, λ, the viscosity of the two fluids, µ, the elastic relaxation time of the system, γ, and the mixing layer width, . We assume that both fluids have the same viscosity and mass density (which is assumed to be 1). It should be noted that as → 0, the parameter λ is associated with the surface tension at the boundary. With the "mixing" energy, W (φ), described above, and following the general approach of Onsager [16,17,18,19], we can prescribe the whole system as a dissipative system: where Ω |φ| 2 dx represents the microscopic internal damping, with the kinematic assumption ofφ = ∂φ ∂t + u · ∇φ. We can then employ the least action principle (LAP) for the Hamiltonian part of the system and the maximum dissipation principle (MDP) for the dissipative part of the system [16,17,18,19,20]. In the case that µ 1 = µ 2 = µ and ρ i = 1, we obtain equations (1)-(3). In case γ = 0, i.e., we neglect the internal damping, we can show that (1)-(3) approaches the sharp interface formulation [15,21]. In particular, where λ ∇ 2 φ∇φ reflects the surface tension, kHn, modulating the pressure contributions.
Remark Later, we indicate the structural similarities between the diffusive interface formulation and that of resistive MHD systems. In this setting, the terms λ∇ · (∇φ ⊗ ∇φ) − ∇ |∇φ| 2 2 are equivalent to the Maxwell stress, and λ∇ 2 φ∇φ is equivalent to the Lorentz force.
Next, we discuss the numerical methods that are used to solve this problem.

FOSLS Formulation
To solve the nonlinear system of equations, which for discussion is denoted by L(u) = f , a first-order system least-squares discretization (FOSLS) is used [6,7,22]. First, consider a linear first-order system, denoted by Lu = f . Then, the linear problem is recast as the minimization of a functional constructed by taking the L 2 norm of the residual of each equation. This is written as where u * is the solution in an appropriate Hilbert space V.
Thus, u * satisfies [G (u * )]v = 0, which is the Fréchet derivative of G at u * applied to v ∈ V. This results in the weak form of the problem: find u * ∈ V such that where ·, · is the usual L 2 inner product. Desirable properties of the bilinear form, coercivity These properties imply the existence of a unique solution, u * ∈ V, for the weak problem (5).
Next, we approximate u * by restricting (4) to a finite-dimensional subspace, V h ⊆ V, which leads to (5) restricted to V h . Since V h is a subspace of V, the discrete problem is also well-posed. Choosing an appropriate basis, V h = span{Φ j }, yields an algebraic system of equations involving the matrix A with elements In general, if V is a subspace of H 1 (product space), and continuity and coercivity hold we say that the FOSLS functional is H 1 equivalent. In this case, the linear system is amenable to an iterative solution by multilevel techniques. In particular, AMG has been shown to work well on a wide range of such problems [6,7,8,9,23,24,25,26,27,28,29].
In addition, note that the functional yields an a posteriori error measure.
The last relation comes from the continuity and coercivity bounds in (6) and (7). Thus, the functional value is equivalent to the error measured in the Hilbert space norm. In general, the constant c in equation (8) depends on the continuity and coercivity constants, c 2 and c 1 . These constants, of course, depend on properties of the PDE as well as on boundary conditions and the computational domain. For this paper, we study simple domains and boundary conditions for which these constants are not large. In addition, our numerical results show that the problem parameters do not adversely affect our accuracy either. Addressing these issues analytically would be a topic of another paper 5 and has been studied in other problems such as in [8,9,10]. Locally, on any subset of Ω, such as an element, the functional yields an estimate of the error. This property of FOSLS helps make it possible to efficiently solve complex systems. At each step in the solution algorithm, a local measure of the functional is available. Since this is the norm we are minimizing, this allows judgements to be made based on estimates of the increase of accuracy that results from an increase in computational cost. As a result, methods such as nested iteration and adaptive local refinement can be used to improve the efficiency of solving such systems. This has been applied to various complex fluids problems, including magnetohydrodynamics [11,12,30,31]. A brief description of the algorithms used is described in section 5.

Linearization (Newton-FOSLS)
Since the two-phase system is nonlinear, we first linearize it before we put it into a FOSLS weak form. The system becomes where u 0 is the current approximation and [L (u 0 )] is the Fréchet derivative evaluated at u 0 . The functional then becomes Minimization of the linearized functional yieldsû * that satisfies the weak form: Onceû is found, it is added to the previous iterate to get the next approximation, Thus, in a nonlinear setting, the FOSLS approach can be applied and, if the linearized functional retains continuity and coercivity, (6) and (7), it also retains the desirable properties as described for linear systems. In addition, the nonlinear functional, that is, the L 2 norm of the residual of the nonlinear system, can be computed as well as the linear functional or the L 2 norm of the linearized system. This allows the two functional values to be compared after each linearization and helps determine if the Newton iterations are converging as expected.
A similar approach would be to create the FOSLS functional of the nonlinear problem and linearize this functional instead. This FOSLS-Newton approach generally involves more terms than the Newton-FOSLS method described here. As the approximation approaches the solution, these additional terms are higher-order and tend to zero faster than the overall functional and the two methods tend to be the same. Thus, when nested iteration is used, as it is here, there tends to be very little difference in the convergence behavior of these two approaches. FOSLS-Newton may be more robust in some especially challenging situations, but the Newton-FOSLS approach is simpler and has been successful in a number of applications [25,32], so we confine our presentation to it in this paper. 6

Velocity-Grad Formulation
To reformulate (1)-(3) as a first-order system, we introduce new variables that represent the gradient of the velocity vector, u, and the phase-field function, φ: This formulation has been shown to work well for the Navier-Stokes equations alone when velocity boundary conditions are given [8,9]. To make the system H 1 equivalent, a few auxiliary equations are added to the system. These are consistent with the original equations and make the first-order system well-posed. The resulting system is as follows: In practice, equation (18) is achieved by eliminating one of the diagonal components of V in terms of the other diagonal components. This comes from the divergence-free constraint on the velocity. Thus, in 2D, there are 13 equations and 9 unknowns and, in 3D, there are 29 equations and 16 unknowns.

Equivalence of the FOSLS Functional
As described above, if continuity, (6), and coercivity, (7), of the FOSLS functional holds, we say the functional is H 1 equivalent. This can be shown for the linearized system at each time step. We assume that an implicit backward difference time-stepping scheme is used and consider this semi-discrete system first. This says that, for each Newton step and time step, the FOSLS discretization gives a well-posed numerical system to solve. We state the main theorem here and give the proof. However, many of the details are omitted to save space and since they have been shown in previous works.
Theorem 3.1 (cf. [6,7,8,9,33]). Let L be the linearized first order two-phase flow operator as in equations (15)- (22) and let G be the FOSLS functional as defined in equation (10) of this linearized system. Let the domain, Ω, be a bounded convex polyhedron with connected boundary or a simply connected region with a C 1,1 boundary. Denote U = (u, V , p, φ, B) T and let Assume u n , V n , B n , and φ n are given members of our finite-element space and, therefore, ||u n || ∞ , ||V n || ∞ , ||B n || ∞ , ||∇ · B n || ∞ , and ||φ n || ∞ < ∞. Then there exists positive constants, c 1 and c 2 , that yield the following bounds: These bounds mean that the functional, G, is equivalent to the H 1 norm of the error.
Proof of 3.1. Proving continuity, (23), of the linearized first-order system involves several instances of the triangle inequality and is, therefore, easy to show. Coercivity, (24), is accomplished by treating the two coupled block systems, Navier-Stokes and the phase field equation, separately. A linearized system of (15)- (22) can be written in block form as follows: Here, A N S represents the linearized Navier-Stokes equations in first-order form and without the coupling term, (15)- (19). Similarly, A phase represents the linearized phase function equations without the velocity field coupling, (20)- (22). The coupling terms are as follows: The superscript n represents that the solution is from a previous Newton step. Since we are working on discrete finite element spaces, we can assume the L ∞ norms of all the previous iterations are bounded.
In [8,9], coercivity has been shown for the A N S block given here. For the phase field system, the coercivity bound can be shown in a similar fashion. In fact, ignoring the nonlinear terms, the system resembles the FOSLS formulation used on the heat equation [6,7]: Forming the functional of the linearized system for just the phase field block, we write Again, it is assumed that an implicit time-stepping scheme is used so α depends on the time step. With part of the nonlinear coupling included in this block, β depends on the parameters of the problem and the previous solutions. This in fact makes the problem resemble an advection-diffusion type system. Since this functional only adds zeroth-order terms (i.e. no derivatives) to this system, a standard compactness argument, [34], can be used to show that there exists a constant, C > 0, such that C ||B|| 2 1 + ||φ|| 2 1 ≤ G phase . It should be noted that this requires c||B|| 2 1 ≤ ||∇ · B|| 2 0 + ||∇ × B|| 2 0 , for some constant c > 0. However, this is true due to the assumptions on the boundary and the boundary condition, n × B = 0 [34]. Thus, the uncoupled system is H 1 equivalent: there exists positive constants,ĉ L andĉ U , such that Next, we consider adding the off-diagonal coupling blocks of the system. In [33], the following lemma was proved: Let A 1 and A 2 be invertible and let there exist a positive constant, C, such that Looking at the upper triangular portion of the two-phase flow system, we see that for some positive constants C 1 and C 2 . Adding a curl term at the end and using a Poincaré inequality, it is easily shown from the coercivity of the phase-field block that for some constant, C 3 > 0. Thus, from Lemma 3.2, we know that Finally, to show equivalence of the full linearized system, it should be noted that the lower left block, A P N , involves a term with no derivatives and another standard compactness argument can be used. Here, again, we must assume that the previous iterations of u and B are bounded. Thus, there exists positive constants c 1 and c 2 such that Therefore, for each time step and linearization, there exists a weak solution to our linear system, and our multilevel solvers will converge to that solution. The constants c 1 and c 2 depend on physical parameters such as γ, µ, and , as well as the previous Newton step iterations and time step approximations. However, they do not depend on the number of degrees of freedom or the size of the finite element grids being used. In the next section, we discuss the energetics of complex fluid flow.

Energetics
While most physical systems can be described by a set of PDEs, usually an underlying energy law also describes the system. In many cases, the PDE is actually derived from the energy laws via a calculus of variation (see e.g. [20]). The main idea is that the change in the total energy of a system over time must equal the total dissipation of the system. If the system is conservative, or there is no dissipation, then the total change is zero and it is a Hamiltonian system. Thus, the interaction or coupling between different scales or phases plays a crucial role in understanding complex fluids. Any set of equations that describe the system should then be a result of the energy laws. The energetic variational approach (EVA) takes the energy laws of a complex fluid model and, using the least action principle (LAP) and the maximum/minimum dissipation principle (MDP) [16,17,18,19], yields a weak form of the system, which we approximately preserve in our numerical simulations. The energy variation is based on the following energy dissipation law for the whole coupled system: where E total is the total energy of the system and D is the dissipation. The LAP, which is also referred to as the Hamiltonian principle, or principle of virtual work, gives us the Hamiltonian (reversible) part of the system related to the conservative force. At the same time, the MDP gives the dissipative (irreversible) part of the system related to the dissipative force.
In [4,5,13,14,15] and others, an energy law is developed for the Allen-Cahn phasefield model described above. The total energy is a combination of the kinetic energy (driven by the fluid) and the internal energy (driven by the interface): where W (φ) is the Ginzburg Landau mixing energy described above [4], thus representing a competition between two fluids with their hydrophilic and hydrophobic properties. The dissipation is a result of a diffusive term from the viscosity of the fluid and from the diffusion of the interface itself: Thus, any system describing this behavior, such as equations (1)-(3) or the FOSLS formulation (15)- (22), should globally produce results that approximately satisfy the following energy law: The overall energy law reflects the multi-scale, multi-physics nature of the system. As mentioned above, the variable φ represents the microscopic description of the mixture.
The kinematic assumption,φ = ∂φ ∂t + u · ∇φ, stands for the influence of macroscopic dynamics on the microscopic scale. The form of the internal dissipation is the consequence that we are looking at the long (or macroscopic) time dynamics. This brings us to the equation for φ, ∂φ ∂t This "gradient flow" dynamics is another formulation of the near equilibrium, linear response theory [18,19]. As for the macroscopic force balance, the LAP gives the following Hamiltonian system, with γ = 0: This describes the transient dynamics of the system. On the other hand, the MDP yields the dissipation of system, which stands for the long macroscopic time dynamics: shows the consistency between the LAP and the MDP [4].

System (1)-(3) is really the hybrid of these two systems.
Remark It is worth pointing out that due to the generality of the formulation, analogies of (1)-(3) can be found in many different physics. In particular, without the Ω 1 4 2 φ 2 − 1 2 dx term, the system is equivalent to the resistive MHD system. In fact, it is this understanding that allows us to arrive at a similar weak form for the FOSLS algorithms as in [11]. The choice of using the gradient of φ as an auxiliary variable in (14) was made for this reason.
Finally, we point out that the advantage of the above energetic variation approach enables us to derive the thermodynamic-consistent models involving different physics at different scales. Restricting to a finite-dimensional space, we show numerically that the FOSLS discretization method is capable of globally preserving this energy law.

Numerical results
In [11,12], an algorithm is devised to solve a system of nonlinear equations, L(u) = f . Starting on a coarse grid with an initial guess there, the system is linearized and the linearized FOSLS functional is then minimized on a finite element space. At this point, several algebraic multigrid (AMG) V-cycles are performed until there is little to gain in accuracy-per-computational-cost. The system is then relinearized and the minimum of the new linearized FOSLS functional is searched for in the same manner. After each set of linear solves, the relative difference between the computed linearized functional and the functional of the nonlinear system is checked. If they are close and near the minimum of the linearized functional, then it is concluded that we are close enough to the minimum of the functional of the nonlinear operator and, hence, we have a good approximation to the solution on the given grid. Next, the approximation is interpolated to a finer grid and the problem is solved on that grid. At this stage, the grid can be refined uniformly or locally. This process is repeated to yet finer grids until an acceptable error has been reached, or until we have run out of computational resources, such as memory. If, as in the case of the 2-phase flow equations, it is a time-dependent problem, the whole process is performed at each time step. This algorithm is summarized in the following pseudocode. The main goal behind this algorithm is to reduce the amount of work needed to solve the system of equations. Relatively cheap work is done on the coarse computational grid to get better starting guesses for the solution of the system on the finer finite element grids. Then, less iterations are needed at the resolution that would require the most work. The algorithm described above also takes into account the fact that resolving too much on coarse grids is unnecessary. Using the FOSLS functionals as estimates of the errors, stopping criteria is implemented for each stage of the algorithm, optimizing the amount of work per computational cost.
In addition to the nested iteration algorithm, adaptive local refinement is used. Again, using the FOSLS a posteriori error estimates, an efficiency-based adaptive local refinement scheme, known as ACE [30,35] can be used fairly easily. The idea behind this scheme is to predict the amount of work and the error reduction obtained from refining only certain portions of the computational domain. Then, the optimal grid is chosen. In simulations of magnetohydrodynamics, which have similar structure to the complex fluid problems described here, this scheme has been shown to reduce the amount of work by a factor of ten in the FOSLS framework [30,31].
For the following test problems, we use the FOSLS finite element method along with nested iteration, AMG, and adaptive local refinement. Here we show that, using these methods, we are able to solve the complex fluid problems well, while still adequately preserving the energetics of the system.

Coalescence
In this test problem, we start with two bubbles of a certain phase immersed in a fluid of a different phase on the unit square, Ω = [0, 1] × [0, 1]. The initial interface between 13 the two fluids is given such that the two bubbles are osculating or "kissing", and then the system is allowed to evolve. As a result, the two bubbles begin to coalesce into one and, eventually, due to the dissipation in the system, the newly formed bubble slowly shrinks and is absorbed by the outside fluid. The initial conditions for this system are similar to that in [4] and are as follows: The parameters used are η = 0.01, µ = 1, = 0.01, γ = 0.01, and λ = 0.0001. A second-order fully implicit backward differencing formula (BDF-2) is used and allows us to take a larger time step for the simulation than in previous work [4]. For the test problem shown, 100 time steps of size ∆t = 0.01 are taken. Plots of the phase field, φ, are given in figure 2. Here, one can see the coalescence of the two bubbles over time and the eventual dissipation. In addition, more refinement is done in the region of the interface. Therefore, most of the computation is being focussed on resolving this region, which is driving the physics in the system. To quantify how well the nested iteration algorithm with adaptive local refinement performs, the simulation was run using uniform refinement. An estimate of the work is computed by calculating the number of non-zeroes in each discrete system and determining how many equivalent fine-grid matrix-vector operations are performed. We call this 14 measure a work unit (WU), and it is a measure of the cost of performing one relaxation sweep, such as Gauss-Seidel or Jacobi, on the finest grid of the simulation. Table 1 shows that, while using adaptive refinement requires more work units (i.e. more V-cycles are needed to solve the problem), on average the number of non-zeroes on the finest grid is about one-tenth of that for uniform refinement. This corresponds to a work ratio of about 0.21. In other words, using adaptive local refinement saved about 80% of the computational cost. Therefore, on average, the above strategies on this test problem cost the equivalent of performing less than 49 Gauss-Seidel relaxation sweeps on a quadratic uniform grid with 65, 536 (256×256) elements. It should be noted here that the functional or estimate of the error for adaptive refinement is about half of the functional obtained from using uniform refinement on average. Both schemes got below the prescribed tolerance, but, on many time steps, the adaptive algorithm got to a much lower error (occasionally by an order of magnitude). Therefore, as far as accuracy-per-computational cost, the adaptive refinement scheme may be even better than the numbers indicate. What we can say is that it used less than 9% of the elements and took almost one-fifth of the computational time to solve the problem.  Table 1: Average estimates of work and accuracy using uniform refinement versus adaptive local refinement In addition, Table 2 shows that, for each refinement level, the number of Newton steps needed to solve the nonlinear problem decreases from almost 5 on the coarse grids to 1 on the finer grids. For all time steps, on the finest-grid, only one Newton step was needed. This greatly reduces the cost of the algorithm compared to if one were to only solve on the finest grid.
Next, we see in figure 3 that the energy law of this system is satisfied during the simulation. The total conservative energy over the whole domain is decreasing and the change in energy coincides with the dissipation in the system. Thus, the numerical simulation is resolving the expected physics of the system.

Square Bubble
In this test problem, we start with a single bubble immersed in a fluid of a different phase on the unit square, Ω = [0, 1] × [0, 1]. The initial interface between the two fluids is given such that the bubble is in the shape of a square. As time evolves, the interface   forces cause the bubble to transform to a circular shape before diffusing away. In the process, oscillations between a square and diamond shape are formed. This reflects the elastic short-time behavior of the system. The initial conditions for this system are similar to that in [14] and are as follows: otherwise .
The parameters used are µ = 0.1, = 0.02, γ = 0.01, and λ = 0.1. Again, a second-order fully implicit backward differencing formula (BDF-2) is used. For the test problem shown, 100 time steps of size ∆t = 0.01 are taken. The phase 16 field is shown in figure 4. As in the previous example, the refinement is being placed in the appropriate regions and requires about 20% of the work that would be required if uniform refinement were used. Due to the bigger jump at the interface, though, more Newton steps are needed in this problem than in the first. However, using nested iteration results in many of these linearizations being performed on the coarser grids. As an example, at time step 3, when the bubble begins to lose its square shape, the number of Newton steps goes from 10 to 1 as we move up through the grids. See table 3. Steps  9  16  10  8  64  4  7  112  5  6  208  5  5  463  5  4  1,423  4  3  3,766  3  2  12,793  3  1 40,030 1 In figure 5, one can see that the energy laws are being satisfied once again. An initial increase in energy occurs before the dissipation in the system takes over, causing the energy to decay. This is most likely due to the fact that the initial condition of the phase field is discontinuous at the interface. Once the simulation is started, this function is immediately smoothed out and the energy behaves as expected. By increasing γ, the elastic relaxation time in the system, the dissipation is increased and actually causes the bubble to dissipate away much faster. In fact, the entire bubble vanishes after time step 30 and the total energy goes to zero. This can be seen in the energy plots in figure 6. 17

Grid Elements Newton
The tiny jump in energy appears to occur when the region containing the bubble finally disappears and the system is changed to that of one fluid. With a smaller time step, this transition occurs more smoothly.

MHD
As described above, the characteristics of the two-phase system is similar to that of resistive MHD. The main results and description of this problem using the FOSLS and nested iteration framework can be found in [11,12,30]. However, here, we briefly mention the energetics of the system. For this system, the energy law must satisfy the following equation: Here, u represents the fluid velocity and B represents the magnetic field. R e and S L are the fluid Reynolds number and magnetic Lundquist number, respectively. A magnetic reconnection problem is simulated using the above numerical algorithm. Figure 7 shows that the discrete energy law is being satisfied. The total conservative energy decays as a direct result of the dissipation in the system. The oscillations in the change in energy plot are due to a lower order differentiation being used to compute dE dt after the simulation was complete. Similar to the first two examples, it took on average of 10 WU to solve the problem using nested iteration and adaptive local refinement.

Discussion
This paper shows that using the Newton-FOSLS finite-element formulation, along with nested iteration, algebraic multigrid, and adaptive local refinement, accurately resolves two-phase flow simulations while still preserving the energetics of the system. The numerical algorithms used were designed in a general setting and have the advantage of reducing the amount of computational cost. For the two-phase flow problems, it is possible to resolve the physics in the equivalent of less than 50 fine-grid relaxation sweeps for each time step. This is accomplished by reducing the amount of work done on highly resolved grids and by reducing the number of elements needed to solve the problem via adaptive local refinement. As a result, complicated systems of PDEs can be solved very efficiently without losing the ability to measure certain physical quantities. This paper shows that, when applied to two-phase flow simulations, all the advantages of using implicit finite element methods are obtained, and the FOSLS methodology is capable of capturing the important aspects of the model, such as the energy laws.
Despite the success of the above algorithms in capturing the correct energy laws efficiently, there has been one bottleneck. In all the simulations, classical algebraic multigrid (AMG) was used to solve the discrete linear systems. Classical AMG was designed for scalar elliptic problems and M-matrices. As can be seen in figure 8, the average convergence factors for AMG over each time step can be quite poor when it is used on systems for which it was not designed. While it is always bounded away from 1, one can see that the worse the AMG convergence factor gets, the more computational work is needed to solve the problem. Therefore, future work will involve investigating various multigrid methods that will be more well-suited to the test problems described above. Algebraic methods that are designed more specifically for systems of PDEs will be studied as well as geometric multigrid. All the domains have been well structured and, therefore, it is possible to use more of the geometric information to create a better solver. If it is possible to improve the linear solver using these methods, then the amount of work can be reduced even further.