COMPUTATION OF EQUILIBRIUM VORTEX STRUCTURES FOR TYPE-II SUPERCONDUCTORS^{∗}

MARK T. JONES AND PAUL E. PLASSMANN^{†}

Abstract. A model of the vortex structure and configuration for homogeneous type-II superconductors can be determined by minimization of the Ginzburg-Landau free energy functional. A generalization of this formulation, the Lawrence-Doniach model, can be used for layered systems. A finite-dimensional approximation of this functional leads to a very large, sparse optimization problem. Recently, a number of algorithms have been used to attempt to solve this problem. The most successful of these methods is a damped, inexact Newton algorithm. The computational kernel of this algorithm is the iterative solution of a large, sparse, linear system that is nearly singular.

In this paper we present computational results for solving these three-dimensional problems on the Intel DELTA. Using a general-purpose, scalable, iterative solver based on the preconditioned conjugate gradient method, we obtained sustained computational rates of up to 4.26 gigaflops on 512 processors. We find improvements of over a factor of 100 in the total execution time when compared with the same problem run on the CRAY-2. This computational tool has given us the first computational view of three-dimensional vortex phenomena such as vortex locking.

1. Introduction. The recent discovery of superconductivity at liquid nitrogen temperatures (Wu et al., 1987) in a class of materials known as cuprate superconductors (Bednorz and Müller, 1986) has generated a renewed interest, both scientific and practical, in the magnetic properties of type-II superconductors. Potential applications for these new materials include high-performance superconducting magnets, efficient electric transmission lines, and superconducting magnetic energy storage (SMES) devices. The prospect of operating these devices at liquid nitrogen temperatures, instead of the more costly liquid helium temperatures, is particularly attractive from an economic perspective. However, despite large values for both the critical temperature and the upper critical magnetic field, the relatively small values for the critical current density in cuprate superconductors have greatly limited the scope of their applications thus far.

The current carrying capability of the cuprate superconductors is intricately related to the properties of the mixed (or vortex) state. Characteristic of the mixed state is the ability of the material to remain superconducting while allowing some magnetic flux into the interior of the material in the form of magnetic vortices. These phenomena are possible because of the formation of compensating vortex currents around the magnetic flux lines which shield the superconducting regions from the effect of the magnetic field. In Figure 1 we give a simplified phase diagram for a model type-II material.

These vortex phenomena are ultimately derived from microscopic effects; a microscopic theory is well developed for low-temperature superconductors. However, for high-temperature materials the microscopic mechanism responsible for the supercon-

∗ This work was supported by the Office of Scientific Computing, U.S. Department of Energy, under Contract W-31-109-Eng-38.

† Mathematics and Computer Science Division, Argonne National Laboratory, 9700 South Cass Avenue, Argonne, IL 60439

Fig. 1. Schematic phase diagram for a conventional type-II superconductor, showing the position of the mixed (vortex) state relative to the ideal-superconducting (Meissner) and nonsuperconducting (normal) states. The coordinates are the temperature T and the external magnetic field strength H.

ductivity is not understood. For this reason, the approach used in this paper is based on the phenomenological Ginzburg-Landau model. We use the Lawrence-Doniach model, an extension of the Ginzburg-Landau model, to model layered systems, the systems of most interest. Both the Ginzburg-Landau and the Lawrence-Doniach models have been remarkably successful in predicting the equilibrium vortex structures and configurations in these high-T_{c} superconductors. They are, however, analytically tractable in only simple configurations that ignore many of the important intrinsic features of high-temperature superconductors (such as arbitrary magnetic field orientations and extreme anisotropy). For practical applications, inclusion of these features is essential to adequately describe flux vortex behavior. Under these circumstances, progress can be made only through numerical simulations of realistic three-dimensional superconducting systems. Because of their three-dimensionality and complicated behavior, these models necessarily involve millions of degrees of freedom. The power of massively parallel computers, along with constant refinement of the models and the development of novel algorithms to solve these large models, has enabled us to consider a wide range of realistic physical scenarios.

In this paper we describe a parallel optimization code developed to solve for the three-dimensional equilibrium vortex configuration in layered systems. This code has been used by materials scientists in a production mode for a number of calculations: magnetization and detailed vortex structure computations, vortex locking measurements, and phase diagram investigations.

The remainder of this paper is organized as follows. In §2 we describe the Ginzburg-Landau approach to modeling these layered superconductors. The optimization algorithms used to solve these problems are described in §3. In §4 we describe the generalpurpose sparse iterative solver that we have developed and used in this particular application. In §5 we present computational results obtained on the Intel DELTA. In §6 we outline two paths of future research.

2. A massively parallel optimization problem. The Ginzburg-Landau (GL) model is a phenomenological model that expresses the local free energy, ΔF, of the physical system in terms of a complex-valued, pseudo-wave function ψ. The (phenomenological) superconducting electron density is given by |ψ|^{2}, and ψ is referred to as the order parameter. Magnetic field effects are included through the introduction of a vector potential, A. The local free energy can be written as the sum of three terms,

ΔF (ψ, A) = F_{cond}(ψ) + F_{kin}(ψ, A) + F_{field}(A),

where F_{cond} is the condensation energy, F_{kin} is the kinetic energy, and F_{field} is the field energy. Note that the order parameter and vector potential are coupled only through the kinetic energy terms. For homogeneous (nonlayered) systems these three terms are given (in cgs units) by the equations

F_{cond} = α|ψ|^{2} + \frac{β}{2} |ψ|^{4},
F_{kin} = \frac{1}{2m} |(\frac{h̄}{i}∇ − \frac{2e}{c}A) ψ|^{2} ,
F_{field} = \frac{1}{8π} |∇ × A|^{2} .

The total free energy over some volume Ω, the free energy functional F, is given by the integral

F = ∫_{Ω} ΔF (ψ, A) dΩ .

The physical configuration of the order parameter and vector potential is given by the global minimizer of this integral. For a homogeneous (nonlayered) material, the minimum energy solution was first predicted by Abrikosov (1988). In this case the vortices are packed in a hexagonal lattice. This lattice can be viewed in Figure 7. Since the system is isotropic, a two-dimensional computational model is sufficient to obtain these results. We note these vortices have the property that the magnitude of the order parameter goes to zero in the center of each of the vortices, where a unit quantum of magnetic flux penetrates the bulk of the material. In addition, the loop line integral of the phase of the complex order parameter is equal to 2πn, where n is the number of vortices enclosed by the path of the loop.

However, the physical systems of most interest are layered, not homogeneous. For layered systems, an extension of the homogeneous free energy developed by Lawrence and Doniach (1971) is appropriate. The material in Lawrence-Doniach (LD) model is composed of alternating layers of relatively thin superconducting sheets separated by insulating layers. An additional term is added to the free energy to account for the energy associated with a Josephson coupling between superconducting currents in adjacent layers. Let ψ_{j} and ψ_{j+1} be the order parameter in layers j and j + 1, respectively. Assume that the x-axis is perpendicular to the layers, and let A_{x} be

Fig. 2. The three-dimensional layered model partitioned in two-dimensions

the x-component of the vector potential. The energy, F_{Josephson}, associated with the Josephson coupling between adjacent layers j and j + 1 is given by the expression

(2.1) F_{Josephson} = ζ |ψ_{j+1} exp [\frac{2ie}{h̄c} ∫_{(j+1)_{s}}^{j_{s}} A_{x}dx] − ψ_{j}|^{2} .

The coupling strength is given by ζ, and (j)_{s} and (j + 1)_{s} are x-coordinates of the j and j + 1 superconducting sheets.

An important symmetry inherent in the above equations is that they are invariant under gauge transformations. That is, given any scalar function, χ, the free energy is unchanged with the replacement of variables

ψ ← ψ′e^{iχ},
A ← A′ + \frac{h̄c}{2ie}∇χ .

Thus, a discretization of the model in the variables ψ and A leads to an underdetermined system, (i.e., a minimizer of F is at best unique only up to a gauge transformation). This degeneracy poses a significant difficulty for any optimization algorithm; in §3 we discuss an effective approach for dealing with this problem.

Modeling this system with an external magnetic field applied at an arbitrary angle requires the full three-dimensional model. In our computational results, the free energy functional is defined on a three-dimensional rectangular mesh. The complex order parameter ψ is defined within the superconducting layers, and the vector potential A is defined within the insulating layers. The number of magnetic flux quanta and the angle of the applied magnetic field are enforced by the imposition of pseudo-periodic boundary conditions on the model.

We obtain numerical solutions for the order parameter and vector potential by first applying a finite-difference discretization method to the free energy expression for the Lawrence-Doniach model. Differencing is accomplished in a gauge invariant manner as described by Garner et al. (1992). Because of their three-dimensionality and complicated behavior, these models necessarily involve millions of degrees of freedom. The resulting discrete problem can then be solved by the large-scale optimization techniques discussed in §3. The solution of this problem yields an approximation to the equilibrium superconducting electron distribution and vector potential that minimizes the free energy functional.

To attack these computationally demanding problems, we are using the Intel DELTA both for its large memory (8 gigabytes^{1}) and its computing power (512 Intel i860 processors). Because the DELTA is a distributed-memory computer in which processors communicate via message passing, the algorithms and their implementations are very different from those used on a traditional sequential computer. We illustrate this difference by considering the parallel computation of the free energy function, its gradient, and Hessian (matrix of second partial derivatives).

Unlike on sequential computers where the entire domain resides on a single processor, we must partition the domain across the memories of all the processors. The performance of the algorithms in this paper depends on having a good partition of the domain. Fortunately, because our domain is regular and the level of discretization is constant across the domain, we know a good partitioning a priori. This partitioning (shown in Figure 2) divides the domain in only two dimensions because the number of points needed to resolve the layers is much smaller than that needed for the other two dimensions.^{2}

Because of the forward differencing approximations in the discrete approximation to the free energy, both local (to a processor) and nonlocal data are required to compute this approximation. The order parameter and vector potential data required to evaluate the function, gradient, and Hessian on the local grid assigned to a processor include both the local grid points and nonlocal, adjacent grid points. Therefore, prior to the evaluation of the function, gradient, and Hessian, each processor must gather the nonlocal data from other processors that it needs to compute these differences. Therefore, unlike a sequential program, we require a communication step prior to these evaluations.

Fortunately, we can simplify this communication step with the following observation. Because of the way we have partitioned the grid and because we are using forward differences, we need grid points only from the edge of each our immediate neighbors in the partitioning in order to do the evaluations. For example, in Figure 3 we show the “ghost” points required by processor 2 to evaluate its portion of the function, gradient, and Hessian. This regular communication pattern is required after each step computed by the optimization algorithm. In our implementation, the interprocessor

1 The problems described in the experimental section require all the available memory on the DELTA.

2 In this case, the physical partition can be used to generate a good vertex assignment to processors. When the determination of a partition is not straightforward, a partitioning heuristic would have to be used. Some possibilities exist; for example, recent advances in the automatic partitioning of threedimensional domains (Vavasis, 1991) or in spectral dissection methods (Pothen, Simon, and Liou, 1990) could be employed.

Fig. 3. The adjacent values, or “ghost” points, required for the function, gradient, and Hessian evaluation

communication required to obtain these nonlocal values is managed by the BlockComm communication package, developed by William Gropp (1992).

3. Large-scale optimization techniques. The original attempts at determining a minimizer of the Ginzburg-Landau free energy functional was done for a homogeneous (nonlayered), two-dimensional model (Doria, Gubernatis, and Rainer, 1990). A Monte-Carlo simulated annealing algorithm was used and was able to obtain solutions only for small grid sizes with at most two vortices per unit cell. These calculations required millions of iterations and many hours of Cray time; eventually, to try to make the problem tractable, the simplier London model was used. The simulated annealing method is most suitable for problems in which many local minima (i.e., metastable configurations) are present and for which derivatives of the function are difficult to obtain. When these conditions do not hold, local optimization methods are more appropriate.

A second group, Wang and Hu (1991), made this observation and used a first-order technique: a steepest descent algorithm with fixed step size. With this method they were able to obtain solutions, for the homogeneous model, for which the gradient norm was reduced to no less than 10^{−5} on two-dimensional grids of size up to 200 × 200.

More sophisticated techniques were introduced by Garner et al. (1992) with greatly improved results. Among these methods were three classes of alternative optimization methods: nonlinear conjugate gradient, limited memory quasi-Newton, and modified Newton methods. Even though they require the same kind of derivative information as the steepest descent method (that is, first derivatives), these algorithms have much better convergence properties. Typically, the first-order methods, limited-memory quasi-Newton and nonlinear conjugate gradients, require several thousand iterations to reduce the gradient norm to 10^{−7}. Unfortunately, a good local minimizer is not always obtained by the first-order methods. The inexact Newton method, a second-order method initially described in Garner et al. (1992), is much more robust. This method usually requires fewer than 20 iterations and can reduce the gradient norm to less than 10^{−12}. However, this method requires the approximate solution of a large sparse linear system at each iteration.

The optimization algorithms and their convergence properties. We briefly describe the three local minimization algorithms mentioned above, with reference to the canonical unconstrained optimization problem

min_{x} F (x) , (3.2)

where x is a real-valued vector of length N. In our case, x includes the values of the order parameter and vector potential the grid points on which they are defined. The iterates produced by each of the methods are denoted by x_{k}, k = 1, 2, ⋯, and we denote the local solution to which these algorithms converge by x^{∗}.

In steepest descent algorithms, x_{k+1} is obtained from x_{k} by choosing the search direction p_{k} = −∇F (x_{k}) and setting

x_{k+1} = x_{k} + α_{k}p_{k}, (3.3)

where α_{k} is a positive scalar chosen in such a way that F (x_{k+1}) is significantly smaller than F (x_{k}). (Choosing the “optimal” value for α_{k} is a one-dimensional optimization problem. In practice, it need not be found accurately.)

The full BFGS algorithm proceeds by maintaining an approximation to the inverse of the Hessian matrix ∇^{2}F (x_{k}), which is denoted by H_{k}. After x_{k+1} has been computed, this matrix is updated with a low-rank correction involving the local gradient and step information. The search direction p_{k} is computed from the formula

p_{k} = −H_{k}∇F (x_{k}), (3.4)

and x_{k+1} is obtained as in (3.3). Like all quasi-Newton algorithms, BFGS is superlinearly convergent under the “standard” assumptions on x^{∗}, that is, ∥x_{k+1} − x^{∗}∥/∥x_{k} − x^{∗}∥ → 0 as k → ∞, and it has proved to be very effective on small to medium sized practical problems. For large problems, it has the disadvantage that H_{k} is dense and therefore, in a naive implementation, requires N^{2}/2 storage locations even when the true Hessian ∇^{2}F (x_{k}) is sparse. When N is greater than about 10^{4}, this is beyond the range of most modern computers. Liu and Nocedal (1990) describe a “limited-memory BFGS” (L-BFGS) algorithm, in which only the last m differences of x_{k} and ∇F (x_{k}) are stored, where m is a small integer, rather than the entire time-history of the iterations. At each step m BFGS updates are applied using these difference vectors. When m = 5, a total of about 14N storage locations is required, while computation of p_{k} (which can be performed efficiently by using a recursive formula) requires about 33N floating-point operations.

In Newton’s method, the search direction p_{k} is obtained by solving the N × N linear system

∇^{2}F (x_{k})p_{k} = −∇F (x_{k}), (3.5) and the line search (to determine the α_{k} in equation (3.3)) is performed to obtain x_{k+1}. This algorithm has the advantage that, under the “standard” assumptions on x^{∗}, it converges quadratically:

lim_{k→∞} \frac{∥x_{k+1} − x^{∗}∥}{∥x_{k} − x^{∗}∥^{2}} < ∞. (3.6)

A number of additional complications over a first-order method are required in order to make Newton’s method practical for these problems. First, for the LD model, a significant effort is involved in development and verification of the code to calculate the Hessian matrix ∇^{2}F (x_{k}) analytically. Obtaining the Hessian by sparse finitedifferencing methods was found to be too computationally expensive.

A second issue is the efficient iterative solution of the large sparse linear system given by equation (3.5). This problem is discussed in more detail in the next section. This system is solved approximately; the iterative solver is terminated once the residual is small enough or too many iterations are required because of ill-conditioning of the damped Newton system.

A third issue arises from the fact that the “standard” assumptions usually made for the convergence analysis of optimization algorithms do not hold for these problems. As a results of the invariance of the free energy functional under gauge transformations, the solution is nonunique. Thus the Hessian, ∇^{2}F (x^{∗}), is positive semi-definite at the solution (i.e., it has a large number of zero eigenvalues), and as a consequence, the basic Newton iteration, equation (3.5), performs very poorly and may not even converge. We have found a “damped” variant to be very effective; instead of computing p_{k} from equation (3.5), we actually solve the damped system

(∇^{2}F (x_{k}) + γ_{k}I)p_{k} = −∇F (x_{k}), (3.7)

where γ_{k} is a positive scalar that is adjusted at each iteration. In practice, we have found an effective choice to be

γ_{k} = α_{k}∥∇F (x_{k})∥ , (3.8)

where α_{k} is a small constant whose value is adjusted according to the number of conjugate gradient iterations required on the previous iteration.

For the present problem, the storage requirement for the Newton method is comparable to that of L-BFGS and depends on problem configuration. Each iteration of the damped Newton method requires much more work than an iteration of L-BFGS, but of course many fewer iterations are required (typically, the iteration count is less than 25, as opposed to many thousands for L-BFGS). Because of the modifications described in the preceding paragraph, and the nonuniqueness of the solution, quadratic convergence can no longer be expected, but we still observe fast linear convergence.

4. Scalable iterative methods for the solution of sparse linear systems. In this section we discuss the algorithms used in the scalable libraries for solving sparse linear systems given in equation (3.7). To solve these systems, we have used software developed for the parallel iterative solution of such systems (Jones and Plassmann, 1992a). The iterative solver we have employed uses an incomplete matrix factorization as a preconditioner for the conjugate gradient algorithm (Manteuffel, 1980). This general-purpose preconditioner performs well for many problems, including the sparse problems arising in this application. An important advantage of this approach is that it is general purpose; it can be used for structured and unstructured problems.

The scalable implementation of this preconditioned conjugate gradient is straightforward with two exceptions which we will note later. The linear algebra software does not require the repartitioning of the matrix or unknowns; each processor is responsible for columns of the matrix and the unknowns that correspond to the grid points on that processor.

This implementation comprises three components: the conjugate gradient algorithm, the matrix by vector multiplication, and the incomplete matrix factorization and triangular matrix solution required for the preconditioner. The parallel implementation of the first two operations is relatively simple. The only communication required by the conjugate gradient algorithm is a global sum for the inner product computation.^{3} In the matrix by vector multiplication, a single communication step is needed to communicate the nonlocal elements of the vector to those processors that need them. Aside from these minor communication steps, both these operations are perfectly parallel.

However, two main obstacles impede the efficient implementation of a iterative solver based on this preconditioning on a high-performance, scalable architecture. First, the triangular system solves do not perform well for standard matrix orderings (Jones and Plassmann, 1991b). Second, it is not sufficient to achieve scalable performance; one must also achieve good computation rates on each processor. It is the nature of high-performance RISC chips (such as the i860) that the best performance is obtained by algorithms that exhibit good data locality and minimize indirect addressing. In the following two subsections, §4.1 and §4.2, we describe the methods used to overcome each of these obstacles and obtain scalable, efficient performance on high-performance, scalable architectures such as the Intel DELTA.

4.1. The scalable inversion of triangular systems. The triangular system solution is the central problem in the parallelization of the standard iterative methods. For example, it is involved in the application of an SOR or SSOR iteration, in addition to preconditioners derived from an incomplete factorization.^{4} The traditional serial approach to generate incomplete factors employs a “natural” ordering of the variables. Unfortunately, a scalable parallel implementation of this approach is impossible because the dependencies in the solution of triangular systems make this computation inherently sequential.

3 For a discussion of the conjugate gradient algorithm see Golub and Van Loan (1983).

4 The scalable computation of the incomplete factorization can be accomplished in the same fashion as the triangular system solution.

However, a reordering of the preconditioning matrix based on a coloring of the graph associated with the matrix does allow for its scalable solution. The reordered triangular system solution is scalable because the number of sequential communication steps is proportional to the chromatic number of the graph (Schreiber and Tang, 1982) which is essentially a function of the local graph structure, and independent of the size of the graph. In Figure 4 we give an example of a multicoloring ordering for a regular grid that requires four colors.

Fig. 4. The adjacency graph corresponding to a nine-point stencil requires four colors. An ordering of the variables of the corresponding linear system allows for the solution of a triangular system of the same structure to be solved in four major parallel steps: one step for the unknowns corresponding to each color, followed by interprocessor communication to update the right-hand side.

For less regular problems for which one does not have a priori knowledge of an optimal graph coloring, a graph coloring heuristic must be used. There are good, lineartime, sequential graph coloring heuristics (Coleman and Moré, 1983). However, because the graph coloring heuristic will be used on a massively parallel computer, it must also be scalable. We have also developed and implemented a scalable graph coloring heuristic based on finding a sequence of independent sets that generates colorings similar those found by sequential heuristics (Jones and Plassmann, 1993).

We have demonstrated that this combination of graph coloring heuristics and incomplete matrix factorization is effective for a range of structured and unstructured finite element and finite difference problems (Jones and Plassmann, 1991b). In addition, we note that recent theoretical results have shown that one does not see the dramatic increase of the number of iterations required for convergence with “many color” orderings that one sees with the red/black ordering for model problems (Jones and Plassmann, 1992b).

4.2. Graph reductions. As discussed previously, it is not sufficient to achieve scalable performance; one must also efficiently use each processing node. For example, a standard implementation of a sparse matrix-vector multiplication does not exhibit good data locality and uses a large amount of indirect addressing. To improve locality and minimize indirect addressing, one can take advantage of the special local structure inherent to many multicomponent problems. For example, large, dense cliques exist in these graphs and can be easily recognized. In all of the sparse operations that we perform, operations involving these cliques can utilize dense Level 2 and 3 BLAS. In addition, many rows of the sparse matrix have identical structure, but differing nonzero values. We can take advantage of this structure and significantly reduce the amount of indirect addressing. We note that these ideas have been used with dramatic effect in direct sparse factorization for several years.

It is often observed that the sparse systems arising in many applications have a great deal of special local structure, even if the systems are described as “unstructured.” We illustrate some of this local structure, and how it can be identified, in the following sequence of figures.

In Figure 5a we depict a subsection of a graph that would arise from a twodimensional, linear, multicomponent finite-element model with three degrees of freedom per node point. We use this example to illustrate the local graph characteristics used by the linear solver. These same basic features are present in the three-dimensional difference graph arising from the superconductor model, however it is clearer to illustrate these features in this two-dimensional example. We illustrate the three degrees of freedom by the three dots at each node point; the linear elements imply that the twelve degrees of freedom sharing the four node points of each face are completely connected. In the figure we show edges only between the nodes; these edges represent the complete interconnection of all the vertices on each element, or face.

The dashed lines in the figure represent a geometric partitioning of the grid; we assume that the vertices in the central region are all assigned to one processor. We make several observations about the local structure of this subgraph. First, we note that the adjacency structure of the vertices at the same geometric node (i.e., the nonzero structure of the associated variables) is identical, and we call such vertices identical vertices. It was noted by Schreiber and Tang (1982) that a coloring of the graph corresponding to the geometric nodes results in a system with small dense blocks, of order the number of degrees of freedom per node, along the diagonal. We note that this observation can also be used to decrease the storage required for indirect indexing of the matrix rows, because the structures are identical.

We also consider a further graph reduction based on the local clique structure of the graph. In Figure 5b the dotted lines show one possible way the vertices assigned

Fig. 5. (a) The top figure is a subgraph generated by a two-dimensional, linear finite element model with three degrees of freedom per node point. The geometric partition shown by the dotted lines yields an assignment of the vertices in the enclosed subregion to one processor. (b) The left figure is a partition of the vertices into cliques. (c) The right figure is a quotient graph given the clique partition in the left figure.

to the shown partition and its neighbors can be partitioned into cliques.^{5} Denote such a partition by Q. If we associate a super-vertex with each clique, the quotient graph G/Q can be constructed based on the rule that there exists an edge between two super-vertices v and w if and only if there exists an edge between two vertices of their respective partitions in G. The quotient graph constructed by the clique partition shown in Figure 5b is shown in Figure 5c.

When the matrix is reordered according to such a clique decomposition, the matrix has large, dense blocks along the diagonal that allow for the use of the higher-level dense BLAS. In addition, by coloring the quotient graph rather than the original graph, the number of colors needed is greatly reduced. Thus, rather than simultaneously solving the diagonal submatrices associated with each color, the processors now simultaneously solve block diagonal submatrices associated with each color.

Finally, we note that the efficient determination of identical nodes, and a local maximal clique decomposition, is straightforward. Since the adjacency structure of the vertices assigned to a processor is known locally, no interprocessor communication is required, and a greedy heuristic can be used to determine a clique partition.

5 Of course the quotient graph reduction is not limited to the choice of a maximal clique partition; any local partition of the subgraph assigned to a processor can be used to generate the reduced graph. Several alternatives are discussed in Jones and Plassmann (1991a).

5. Computational results. In this section we present experimental results obtained on the Intel DELTA for the superconductivity modeling problem. First, we describe the performance of the implementation for three representative problem test sets. Second, we present some recent computational measurements of the vortexlocking phenomenon in layered type II materials.

5.1. Performance on the Intel DELTA. In Table 1 we give the relevant parameters for three representative test problem sets. For these problems NX, NY, and NZ are the number of grid points used in the X, Y, and Z directions, respectively. NK is the discretization used within an insulating layer for the vector potential; thus, the number of layers is given by NX/NK. The number of vortices per unit cell induced by the pseudo-periodic boundary is given by VORNUM. The number of independent variables is given by N, and the number of nonzeros in the Hessian is given by NNZ.

Timings for the L-BFGS method and the damped, inexact Newton method described in §3 are given in Table 1. Shown are the number of iterations required, the norm of the gradient at termination, and the sustained computational rate for the inexact Newton method. For comparison, we have included the L-BFGS running time for PROBLEM-3 on a CRAY-2 (Garner, 1992).

PROBLEM-1 is representative of a problem where the resolution of the vector potential within the insulating layers is important. This resolution is required when the vortex coherence length is small relative to the interlayer spacing. For problems where coarser resolution of the layers can be allowed, as in PROBLEM-2 and PROBLEM-3, larger numbers of layers (NX/NK) can be investigated. PROBLEM-2 and PROBLEM-3 differ in the size and aspect ratio of the unit cell, the angle of the bulk magnetic field (determined by the aspect ratio of the unit cell), and the number of vortices used (set by the magnetic flux boundary conditions).

The difficulty of the problem solution greatly increases as the number of vortices and the problem size are increased. For the problems in Table 1, L-BFGS was at best able to obtain approximate solutions; for all these problems L-BFGS was not able to determine a solution that satisfied the convergence criteria, and thus was not always physically meaningful. For the results next to the symbol (**) the line search was unable to satisfy the Wolfe conditions, and for the results next to the symbol (*) the maximum number of L-BFGS iterations was exceeded. In general, L-BFGS is unable to adequately solve problems for which VORNUM is greater than two, whereas the inexact Newton code has been able to obtain solutions for hundreds of vortices per unit cell. Thus, to obtain accurate solutions for large problems or for problems involving larger numbers of vortices, the inexact Newton method must be used. Finally, in comparison with the performance of L-BFGS on the CRAY-2, we note that an improvement of approximately 95 times was obtained on the Intel DELTA for PROBLEM-3 (Garner, 1992), or an improvement of more than a 100 times when compared to the inexact Newton method. The L-BFGS code does not perform particularly well on the CRAY-2 because it executes primarily level-1 BLAS and scalar operations.

We emphasize that all the symbolic computation required to reorder the linear

Table 1 Problem set description.

PROBLEM-1 PROBLEM-2 PROBLEM-3
NX  24  64  20
NK  8  4  2
NY  80  64  150
NZ  96  96  150
VORNUM  4  4  1
N  6.0 × 10^{5}  1.6 × 10^{6}  1.8 × 10^{6}
NNZ  2.0 × 10^{8}  1.7 × 10^{8}  1.9 × 10^{8}
L-BFGS (CRAY-2)  ≈ 40 hr
L-BFGS (DELTA) 45.5 min (**) 61.3 min (*) 25.3 min. (*)
L-BFGS (∥g∥)  1.2 × 10^{−8}(∗∗) 5.0 × 10^{−7} (*) 1.0 × 10^{−6} (*)
L-BFGS (iters)  14,664 (**)  10,000 (*)  3,100 (*)
Newton (DELTA)  30.1 min  43.2 min  22.1 min
Newton (∥g∥)  3.0 × 10^{−12}  1.0 × 10^{−11}  1.0 × 10^{−12}
Newton (iters)  22  27  17
Newton (Gflops)  4.26  2.55  1.38

systems is done in parallel.^{6} Note that the coloring is not done for the original graph but, rather, for the quotient graph obtained for the original graph modulo the identical node structure and the local cliques. For example, in PROBLEM-3 the original graph size is 1.8×10^{6}, and the size of the clique graph is 7.3×10^{4}. Thus, the size of the graph colored (the clique graph) is actually quite modest. The time required to determine the clique graph on the DELTA for PROBLEM-3 is 3.359 seconds, and the time required to color the clique graph is 0.059 seconds. These operations are done only once, since the matrix structure does not change between iterations; these times are negligible with respect to the total solution times. Other characteristics of the clique graph for PROBLEM-3 are that the average degree of the graph is 105 and the number of colors required by the parallel heuristic is 11.

In Table 2 we note the effect of varying the layer discretization on the i860 processor performance during the solution of the linear systems. For these numbers we have used 128 processors and fixed the local problem size (the subproblem assigned to each processor) to be roughly equivalent. The second column shows the number of identical nodes found in the graph by the solver; the third column shows the average clique size found. The final column shows the average computational rate per processor during the solution of the linear systems.

We demonstrate the scalability of the implementation in Table 3 by fixing the problem size assigned to each processor and varying the number of processors. Note

6 The parallel coloring heuristic used is described by Jones and Plassmann (1993). The clique reduction algorithm is described by Jones and Plassmann (1991a).

Table 2

Effect of varying the layer discretization on the processor performance in solving the linear systems

I-node Avg. clique Avg. Mflops/
NK  size  size processor
2  8  32.0  2.97
4  14  44.8  5.42
6  20  60.0  6.71
8  26  78.0  8.96

that the average performance per processor appears to depend only slightly on the number of processors used. These numbers are close to the peak performance on the i860 that could be obtained with these block sizes using the Level 2 BLAS.

Table 3 Effect of varying the number of processors for fixed problem size per processor (NK = 8)

Number of Avg. Mflops/  Total
processors processor performance (Gflops)
64  8.75  0.56
128  8.96  1.15
256  8.75  2.24
512  8.33  4.26

5.2. Computational measurements of vortex phenomena. The model has been used to extract detailed information about the three-dimensional vortex structure. The advent of massively parallel computers and the development of efficient computational algorithms have made this first “computational” glimpse into a layered, type-II superconductor possible.

For example, in Figure 8 we show isosurfaces of the magnetic field magnitude for a system near H_{c2}. This solution was obtained by using the inexact, damped Newton method. The external magnetic field is applied at an angle fixed by the pseudo-periodic boundary conditions imposed on the model. This view is taken parallel to the layers of the material; we can see the vortices tilted with respect to the layers. One can observe the “pinching” of the vortices as the magnetic field isosurfaces pass through the superconducting sheets and then spread out again in the insulating layers. In addition, one can see the characteristic “staircase” structure of the vortices that has been observed experimentally when the applied magnetic field is tilted.

In Figure 6 we demonstrate computational measurements made of the “vortexlocking” phenomenon. These calculations were done using the Lawrence-Doniach model (with κ = 5) and required approximately 30,000 node-hours on the Intel DELTA. Corresponding experimental measurements have been made by Liu et al. (1991; 1992). Intuitively, this effect can explained by noting that it is energetically favorable for the

Fig. 6. Vortex locking measurements. The angle ΔΘ between the direction of the bulk magnetic field (Θ_{B}) and the direction of the applied magnetic field (Θ_{H}) is plotted as a function of the angle Θ_{B} between the direction of the bulk magnetic field and the direction normal to the superconducting layers. Three curves are shown for different values of the average magnetic field. These calculations were done using the Lawrence-Doniach model (κ = 5) and required approximately 30,000 node-hours on the Intel DELTA.

vortices to stay in the superconducting sheets. Thus, the angle of the magnetic field in the bulk of the material will lag behind the external field, with the magnetization of the material making up the difference.

6. Conclusions. We have described a methodology for the phenomenological modeling of high-temperature superconductors on massively parallel computers. This methodology relies on good numerical optimization algorithms as well as scalable sparse linear algebra algorithms. An implementation of this methodology on the Intel DELTA computer exhibited scalable performance of 4.26 gigaflops. This tool provides us with the capability to solve three-dimensional models of layered superconductors that were previously considered intractable on CRAY vector supercomputers. We described the first computational results for these three-dimensional models.

We will pursue two paths of further research. First, in conjunction with our collaborators we intend to continue to use the current program to investigate other aspects of high-temperature superconductors. Second, we will build a new tool that will allow us to model superconductors with large κ values. The large value of κ implies very different length scales for the vector potential and the order parameter; this implies that the current regular grid discretization would require billions of unknowns to model high κ superconductors. We plan to construct a finite element representation of the Lawrence-Doniach model that will allow an irregular mesh to be used to allow the modeling of high κ superconductors in a few million unknowns, rather than a few billion.

Acknowledgments. We acknowledge our scientific collaborators: James Garner (Physics Department, Bradley University); Stephen Wright and David Levine (Mathematics and Computer Science Division, Argonne); and Roy Benedek (Materials Science Division, Argonne).

REFERENCES

[1] A. Abrikosov, Fundamentals of the Theory of Metals, North-Holland Publ. Co., Amsterdam, 1988.

[2] G. Bednorz and A. Müller, Possible high T_{c} superconductivity in Ba-La-Cu-O system, Z. Phys. B, 64 (1986), pp. 189–193.

[3] T. F. Coleman and J. J. Moré, Estimation of sparse Jacobian matrices and graph coloring problems, SIAM Journal on Numerical Analysis, 20 (1983), pp. 187–209.

[4] M. M. Doria, J. E. Gubernatis, and D. Rainer, On solving the Ginzburg-Landau equations by simulated annealing, Physical Review B, 41 (1990), pp. 6335–6340.

[5] J. Garner, Unpublished information. Materials Science Division, Argonne National Laboratory, April 1992.

[6] J. Garner, M. Spanbauer, R. Benedek, K. Strandburg, S. Wright, and P. Plassmann, Critical fields of Josephson-coupled superconducting multilayers, Physical Review B, 45 (1992), pp. 7973–7983.

[7] G. H. Golub and C. F. V. Loan, Matrix Computations, The Johns Hopkins University Press, Baltimore, 1983.

[8] W. Gropp, (Unpublished information). Mathematics and Computer Sciences Division, Argonne National Laboratory, 1992.

[9] M. T. Jones and P. E. Plassmann, Scalable iterative solution of sparse linear systems, Preprint MCS-P277-1191, Mathematics and Computer Science Division, Argonne National Laboratory, Argonne, Ill., 1991.

[10] , BlockSolve v1.0: Scalable library software for the parallel solution of sparse linear systems, ANL Report ANL-92/46, Mathematics and Computer Science Division, Argonne National Laboratory, Argonne, Ill., 1992.

[11] , The effect of many-color orderings on the convergence of iterative methods, in Proceedings of the Copper Mountain Conference on Iterative Methods, SIAM LA-SIG, 1992.

[12] , An efficient parallel iterative solver for large sparse linear systems, in Graph Theory and Sparse Matrix Computation, A. George, J. Gilbert, and J. W. Liu, eds., vol. 56 of IMA Volumes in Mathematics and its Applications, Springer-Verlag, 1993.

[13] , A parallel graph coloring heuristic, SIAM Journal on Scientific Computing, 14 (1993), pp. 654–669.

[14] W. E. Lawrence and S. Doniach, Theory of layer structure superconductors, in Proc. 12th Int’l. Conf. on Low Temp. Physics, E. Kanda, ed., Kyoto, 1971, Academic Press of Japan, pp. 361–362.

[15] D. C. Liu and J. Nocedal, On the limited-memory BFGS method for large scale optimization, Mathematical Programming, 45 (1990), pp. 503–528.

[16] L. Liu, J. Kouvel, and T. Brun, Directional vortex-flus bifurcation in polycrystalline YBa_{2}Cu_{3}O_{7} rotated in a magnetic field, Physical Review B, 43 (1991), pp. 7859–7865.

[17] , Anisotropic-vortex-flux production and retention in grain-oriented YBa_{2}Cu_{3}O_{7}, Physical Review B, 45 (1992), pp. 3054–3058.

[18] T. A. Manteuffel, An incomplete factorization technique for positive definite linear systems, Mathematics of Computation, 34 (1980), pp. 473–497.

[19] A. Pothen, H. Simon, and K.-P. Liou, Partitioning sparse matrices with eigenvectors of graphs, SIAM Journal on Matrix Analysis, 11 (1990), pp. 430–452.

[20] R. Schreiber and W.-P. Tang, Vectorizing the conjugate gradient method. Unpublished

information, Department of Computer Science, Stanford University, 1982.

[21] S. Vavasis, Automatic domain partitioning in three dimensions, SIAM Journal on Scientific and Statistical Computing, 12 (1991), pp. 950–970.

[22] Z. D. Wang and C. R. Hu, Numerical relaxation approach for solving the general Ginzburg-Landau equations for type-II superconductors, Physical Review B, 44 (1991), pp. 11918– 11923.

[23] M. Wu, J. Ashburn, C. Torng, P. Hor, R. Meng, L. Gao, Z. Huang, Y. Wang, and C. Chu, Superconductivity at 93 K in a new mixed-phase Y-Ba-Cu-O compound system at ambient pressure, Physical Review Letters, 58 (1987), pp. 908–910.

Fig. 7. The magnitude of the order parameter shown for the two-dimensional homogeneous model. The plane of the image is perpendicular to the applied magnetic field; thus we see a two-dimensional slice through the vortex “tubes.” The global minimum for the free energy functional is obtained by a hexagonal packing of the vortices.

Fig. 8. View parallel to the layers of the vortex structure in the layered three-dimensional model. Pictured are isosurfaces of the magnetic field magnitude as the vortices cross superconducting layers.

Fig. 9. View looking down on the layers of the vortex structure in the layered three-dimensional model. Pictured are isosurfaces of the order parameter magnitude. Recall that the order parameter is defined only within the superconducting sheets.