Numerical analysis of hierarchical methods for phase field problems

Significance of phase field models

The Cahn–Hilliard equation, for example, includes fourth order space derivatives and is used for modelling membranes in biological cells [5], phase-separation in alloys [3], but also in image processing to reconstruct binary images with missing pieces [2].

In physics, phase fields are described by virtue of the Ginzburg–Landau equations, a representative of the general class of non-linear Schrödinger equations. They describe the time evolution of supercurrent and also superfluids like Bose–Einstein condensates.

While the main practical interest lies in efficient computational methods for integrating these phase field equations in time, there is also a strong interest in numerical methods for studying the transient behavior and finding the steady states of a given system. These are domain configurations that remain stable over time; they often correspond to free energy optimal arrangements of the subdomains. The problem of tracking the steady states in the parameter spaces leads to questions in the field of bifurcation analysis. An automated way to locate and track the steady and transient configurations is highly desirable in many applications such as the design of superconducting nano-devices, systems biology, microstructure in material science.

The Ginzburg–Landau problem

To turn a superconductor into practical use, it is necessary to be able to control its structure and properties on a nano- and micrometer scale (i.e. the mesoscale), to maximize the critical current and to understand and exploit the concurrent electronic phases.

The state of a mesoscopic superconductor theoretically described as a wave functions that is the solution of a non-linear Schrödinger equation, known as the Ginzburg–Landau equation. This particular system treats the phase field as a complex valued function and has a magnetic field as an additional unknown (see ().

The equations

Let the domain of the superconductor Ω ⊆ ⁿ be a bounded and open.

Any state of a superconductor is completely described by the two observables supercurrent density ρ(x) and magnetic field B(x). The actual equations, however, merely include the so-called order parameter ψ C_ℂ²(Ω) and the magnetic vector potential A C^_n2(Ω). They are related to the observables by

()

where n is the normal unit vector, and κ a given parameter that influences the behavior of the superconductor. Besides the strong non-linearity, the Ginzburg–Landau equations feature a coupling between the magnetic potential field and the actual phase field.

Special properties of ()

As it is immediately seen, the pair (ψ,A) ≡ 0 is a solution to the system (). This corresponds to the physical notion of the system having no reason to diverge away from non-superconductive state without external excitation.

On the other hand, there are clearly more solutions to the problem which is typical for this kind of non-linear PDEs. It is hence of particular importance to get a concise idea of the – rather complex, as it turns out – solution landscape.

Simplification of ()

In the case of κ ≫ 1, it is justified to treat the magnetic field as constant B(x) = (0,0,H₀)^T , and to adapt the vector potential accordingly; for example

The equations () then break down to

(s)

Typical solutions for ψ on a square-shaped domain with different values for H₀ are shown in the figures 2 and 3.

Figure 2: Typical solution ψ of the Ginzburg–Landau equation (s) with H0 = 0.5 on a square-shaped domain with edge length 10. Left: The absolute value |ψ|. Right: The angle arg ψ. It is clearly seen that the application of the external magnetic field with H0 has the superconductor pierced at four vortices where there is (almost) no superconduction. Going around these vortices, the argument arg ψ changes by 2π.

Figure 3: Typical solution ψ of the Ginzburg–Landau equation (s) with H0 = 0.8 on a square-shaped domain with edge length 10. Left: The absolute value |ψ|. Right: The angle arg ψ. As in figure re 2, the sample is pierced at pronounced vortices. Because the applied magnetic field is somewhat stronger here, is comes as little surprise that there are more vortices and that there is less superconduction in general on the domain.

A chosen numerical obstacle

As mentioned above, the equations () (and also (_s)) have the effectively equal solutions (exp(iη)ψ₀,A₀) with arbitrary η once a solution (ψ₀,A₀) is given. This means that the Jacobian _ψ0 at ψ₀ has a non-trivial null space; it is more precisely

This poses a serious limitation to how one can treat the non-linear problem as standard methods (such as Newton’s) rely on equation systems with the Jacobian to be solved in each step. Approaching a solution, this will be harder and harder to do as the condition number increases with a dropping smallest eigenvalue. The respective effect is displayed in figure 4.

Figure 4: Newton sequence when solving (s). Left: 2-norm of F(ψ) at steps k. Right: Condition numbers at steps k. It is observed that the method works as expected when not too close to a solution. When approaching = 0, the Jacobian gets more and more ill-conditioned and finally, the system is solved with so large errors that update takes the iterates away from a solution. Then, the process is restarted without ever actually reaching a solution.

There are different ways to tackle the problem, a prominent one being described [4]. The key is basically to identify additional (redundant) variables and to come up with suitable extra conditions which are then incorporated into the system (see, for an application, figure 5).

Figure 5: Newton sequence when solving (s) with a bordering strategy applied to the Jacobian. Left: 2-norm of F(ψ) at steps k. Right: Condition numbers at steps k. As opposed to the situation in figure 4, the condition number remains bounded even when getting close to a solution. The system can be solved without much ado and Newton’s method assures quick convergence.

Notes on the implementation

To date, there implementations of the solver of (_s) have been developed for both Fortran (free form) and MATLAB.

For the Fortran code in particular, routines of the packages

• LAPACK,
• BLAS (in particular Kazushige Goto’s excellent version),
• ARPACK (for calculating the k smallest eigenvalues of a matrix),
• Higham’s norm estimator (plus adapted versions), and
• SuperLU (for solving the appearing equations systems)

are in use.

For numerical continuation, it is planned to make use of MATCONT or Friedman’s modified version of it on the MATLAB side; and (at least) one of the libraries

• AUTO and
• LOCA.

on the Fortran side.

As of now, the compiled Fortran code delivers about 20 times faster performance than the equivalent MATLAB code which is mostly due to MATLAB’s being an interpreted (and not compiled) language.

References