## Servicios Personalizados

## Revista

## Articulo

## Indicadores

- Citado por SciELO
- Accesos

## Links relacionados

- Similares en SciELO

## Compartir

## Geofísica internacional

##
*versión On-line* ISSN 2954-436X*versión impresa* ISSN 0016-7169

### Geofís. Intl vol.52 no.3 Ciudad de México jul./sep. 2013

Original paper

**Parallel Algorithms for Computational Models of Geophysical Systems**

**Antonio Carrillo-Ledesma, Ismael Herrera* and Luis M. de la Cruz**

*Instituto de Geofísica, Universidad Nacional Autónoma de México, Ciudad Universitaria, Delegación Coyoacán, 04510, México D.F., México.* *Corresponding author: iherrera@geofisica.unam.mx

Received: October 10, 2012

Accepted: February 05, 2013

Published on line: June 28, 2013

**Resumen**

Los modelos matemáticos de muchos sistemas geofísicos requieren el procesamiento de sistemas algebraicos de gran escala. Las herramientas computacionales más avanzadas están masivamente paralelizadas. El software más efectivo para resolver ecuaciones diferenciales parciales en paralelo intenta alcanzar el *paradigma* *de los métodos de descomposición de dominio*, que hasta ahora se había mantenido como un anhelo no alcanzado. Sin embargo, un grupo de cuatro algoritmos –los *algoritmos DVS*- que lo alcanzan y que tiene aplicabilidad muy general se ha desarrollado recientemente. Este artículo está dedicado a presentarlos y a ilustrar su aplicación a problemas que se presentan frecuentemente en la investigación y el estudio de la Geofísica.

**Palabras clave:** computational-geophysics, computational-PDEs, non-overlapping DDM, BDDC; FETI-DP.

**Abstract**

Mathematical models of many geophysical systems are based on the computational processing of large-scale algebraic systems. The most advanced computational tools are based on massively parallel processors. The most effective software for solving partial differential equations in parallel intends to achieve the *DDM-paradigm*. A set of four algorithms, the *DVS-algorithms*, which achieve it, and of very general applicability, has recently been developed and here they are explained. Also, their application to problems that frequently occur in Geophysics is illustrated.

**Key words:** computational-geophysics, computational-PDEs, non-overlapping DDM, BDDC, FETI-DP.

**1. Introduction**

Mathematical models of many systems of interest, including very important continuous systems of Earth Sciences and Engineering, lead to a great variety of partial differential equations (PDEs) whose solution methods are based on the computational processing of large-scale algebraic systems. Furthermore, the incredible expansion experienced by the existing computational hardware and software has made amenable to effective treatment problems of an ever increasing diversity and complexity, posed by scientific and engineering applications [PITAC, 2006].

Parallel computing is outstanding among the new computational tools and, in order to effectively use the most advanced computers available today, massively parallel software is required. Domain decomposition methods (DDMs) have been developed precisely for effectively treating PDEs in parallel [DDM Organization, 2012]. Ideally, the main objective of domain decomposition research is to produce algorithms capable of *'obtaining the* *global solution by* *exclusively solving* *local problems'*, but up-to-now this has only been an aspiration; that is, a strong desire for achieving such a property and so we call it *'the DDM-paradigm'*. In recent times, numerically competitive DDM-algorithms are *non-overlapping*, *preconditioned* and necessarily incorporate *constraints* [Dohrmann, 2003; Farhat *et al*., 1991; Farhat *et al.*, 2000; Farhat *et al.*, 2001; Mandel, 1993; Mandel *et al.*, 1996; Mandel and Tezaur, 1996; Mandel *et al.*, 2001; Mandel *et* *al.*, 2003; Mandel *et al*., 2005; J. Li *et al*., 2005; Toselli *et al*., 2005], which pose an additional challenge for achieving the *DDM-paradigm*.

Recently a group of four algorithms, referred to as the *'DVS-algorithms'*, which fulfill the *DDM-paradigm*, was developed [Herrera *et al*., 2012; L.M. de la Cruz *et al*., 2012; Herrera and L.M. de la Cruz *et al*., 2012; Herrera and Carrillo-Ledesma *et al*., 2012]. To derive them a new discretization method, which uses a non-overlapping system of nodes (the *derived-nodes*), was introduced. This discretization procedure can be applied to any boundary-value problem, or system of such equations. In turn, the resulting system of discrete equations can be treated using any available DDM-algorithm. In particular, two of the four *DVS-algorithms* mentioned above were obtained by application of the well-known and very effective algorithms BDDC and FETI-DP [Dohrmann, 2003; Farhat *et al*., 1991; Farhat *et al*., 2000; Farhat *et al.*, 2001; Mandel *et al.*, 1993; Mandel *et al.*, 1996; Mandel and Tezaur, 1996; Mandel *et al.*, 2001; Mandel *et al.*, 2003; Mandel *et al*., 2005; J. Li *et al.*, 2005; Toselli *et* *al.*, 2005]; these will be referred to as the *DVS-BDDC* and *DVS-FETI-DP* algorithms. The other two, which will be referred to as the *DVS-PRIMAL* and *DVS-DUAL* algorithms, were obtained by application of two new algorithms that had not been previously reported in the literature [Herrera *et al*., 2011; Herrera *et al*., 2010; Herrera *et* *al*., 2009; Herrera *et al.*, 2009; Herrera, 2008; Herrera, 2007]. As said before, the four *DVS-algorithms* constitute a group of preconditioned and constrained algorithms that, for the first time, fulfill the *DDM-paradigm* [Herrera *et al*., 2013; L.M. de la Cruz *et al.*, 2012].

Both, BDDC and FETI-DP, are very well-known [Dohrmann, 2003; Farhat *et al*., 1991; Farhat *et al*., 2000; Farhat *et al*., 2001; Mandel *et al.*, 1993; Mandel *et al.*, 1996; Mandel and Tezaur, 1996; Mandel *et al*., 2001]; and both are highly efficient. Recently, it was established that these two methods are closely related and its numerical performance is quite similar [Mandel *et al*., 2003; Mandel *et al*., 2005]. On the other hand, through numerical experiments, we have established that the numerical performances of each one of the members of *DVS-algorithms* group (*DVS-BDDC*, *DVS-FETI-DP*, *DVS-PRIMAL* and *DVS-DUAL*) are very similar too. Furthermore, we have carried out comparisons of the performances of the standard versions of BDDC and FETI-DP with *DVS-BDDC* and *DVS-FETI-DP*, and in all such numerical experiments the DVS algorithms have performed significantly better.

Each *DVS-algorithm* possesses the following conspicuous features:

• It fulfills the

DDM-paradigm;• It is applicable to symmetric, non-symmetric and indefinite matrices (i.e., neither positive, nor negative definite); and

• It is preconditioned and constrained, and has update numerical efficiency.

Furthermore, the uniformity of the algebraic structure of the matrix-formulas that define each one of them is remarkable.

This article is organized as follows. In Section 2 the basic definitions for the DVS framework are given; here we define the set of 'derived-nodes', internal, interface, primal and dual nodes, the 'derived-vector-space', among others. Section 3 is devoted to define the new set of vector spaces that conforms the DVS framework; the Euclidean inner product, is also defined here. In Section 4 the 'transformed-problem' on the derived-nodes is explained in detail, and this is our starting point to define the DVS algorithms. Section 5 presents a summary of the four DVS-algorithms: DVS-BDDC, DVS-FETI-DP, DVS-PRIMAL and DVS-DUAL. In Section 6 we give the numerical procedures we use to fulfilling the DDM-paradigm, and we explain in detail the implementation issues. Finally, in Section 7 we show some numerical results obtained after the application of the DVS-algorithms in the solution of several boundary values problems of interest in Geophysics. We studied examples for a single-equation, for the cases of symmetric, non-symmetric and indefinite problems. We also present results for an elasticity problem, where a system of PDE equations is solved.

**2. DVS Framework: A Summary**

The *'derived-vector-space framework (DVS-framework)'* is applied to the discrete system of equations that is obtained after the partial differential equation, or system of such equations, has been discretized. The procedure is independent of the method of discretization that is used. Thus, the DVS-framework's starting point is a system of linear algebraic equations that is referred to as the *'original problem'*:

However, in the *DVS* setting one does not work with the set of nodes originally used for discretizing the problem the *original-nodes'* (Figure 1). Instead, one uses an auxiliary set of nodes: the *'derived-nodes'*. Each one of such nodes has the property that it belongs to one and only one subdomain of the *coarse mesh*.

Indeed, generally after a *coarse-mesh* has been introduced, some *original-nodes* belong to more than one subdomain of the *coarse-mesh* (Figure 2), which is inconvenient for achieving the *DDM-paradigm*. Therefore, in the* DVS-framework*, each *original-node* that belongs to more than one subdomain is divided into as many new nodes –the *derived-nodes* (Figure 3) - as subdomains it belongs to. Then, the *derived-nodes* so obtained are distributed into the *coarse-mesh* subdomains so that each *derived-node* is assigned to one and only one subdomain of the *coarse-mesh* (Figure 4). Once this has been done, a convenient notation is to label each* derived-node* by a pair of natural numbers: the first one indicating the *original-node* from which it derives and the second one, the subdomain to which it is assigned.

The real-valued functions defined in the set of *derived-nodes* constitute a vector-space: the *'derived-vector-space'*, *W*. This space becomes a finite-dimensional Hilbert-space when it is supplied with the inner-product that is usually introduced when dealing with real-valued functions defined in a set of nodes; this is referred to as the *Euclidean* *inner-product*.

Afterwards, a new problem (referred to as the *'transformed problem'*) is defined in the *derived-vector-space*, which is equivalent to the original system of discrete equations. Thereafter, all the numerical and computational work is carried out in the *DVS-space*.

Before leaving this Section, we dwell a little further on the meaning of a *coarse-mesh*. By it, we mean a partition of Ω into a set of non-overlapping subdomains {Ω_{1},..., Ω_{E}}, such that for each *α*=1, ..., *E*, Ω_{α}, is open and:

Where stands for the closure of Ω_{α}. The set of *'subdomain-indices'* will be

, *α*=1,..., *E*, will be used for the subset of *original-nodes* that correspond to nodes pertaining to . As usual, nodes will be classified into *'internal'* and *'interface-nodes'*: a node is *internal* if it belongs to only one partition-subdomain closure and it is an *interface-node*, when it belongs to more than one. For the application of *dual-primal* methods, *interface-nodes* are classified into *'primal'* and *'dual'* nodes. We define:

as the set of

internal-nodes;as the set of

interface-nodes;as the set of

primal-nodes^{1}; andas the set of

dual-nodes.

The set of *primal-nodes* is required to be a subset of and, in principle, could be otherwise chosen arbitrarily. However, the algorithms considered by *domain decomposition methods* are iterative-algorithms and their rate of convergence depends crucially on the selection of the set . Thus, criteria for selecting have been studied extensively (see [Toselli *et al.*, 2005], for detailed discussions of this topic). Each one of the following two families of node-subsets is disjoint: and . Furthermore, these node subsets fulfill the relations:

Throughout our developments the *original matrix* is assumed to be non-singular (i.e., it defines a bijection of into itself). The following assumption (*'axiom'*) is also adopted in throughout the *DVS-framework*: "When the indices and are *internal original-nodes*, while *α* ≠ *β*, then and are unconnected". We recall that unconnected means:

**3. The Derived-Vector Space (DVS)**

In order to have at hand a sufficiently general framework, we consider functions defined on the set X of *derived-nodes* whose value at each *derived-node* is a *dD*−*Vector*. The numerical applications that will be discussed in this paper correspond to two possible choices of *d*: when the application refers to a single partial differential equation (PDE), *d*=1, and for the problems of elasticity that will be considered, which are governed by a three-equations system, *d*=3.

Independently of the chosen value for *d*, the set of such functions constitute a vector space, *W*, referred to as the *'derived-vector space'*. When , we write *u*(*p*, *α*) for the value of *u* at the *derived-node* (*p*, *α*). We observe that, in general, *u*(*p*, *α*) itself is a *d*−*Vector* and we adopt the notation *u*(*p*, *α*, *i*), *i*=1, ..., *d*. For the *i*−*th* component of *u*(*p*, *α*). When *d*=1 the index *i* is irrelevant and, in such a case, will deleted throughout.

For every pair of functions, and , the *'Euclidean inner product'* is defined to be

Here, *u*(*p*, *α*) *w*(*p*, *α*) stands for the inner-product of the *dD*−*Vectors* involved; thus,

A fundamental property of the *derived-vector* *space* *W*, is that it constitutes a finite dimensional *Hilbert-space* with respect to the *Euclidean inner-product*.

Let *W'* ⊂ *W* be a linear subspace and assume M ⊂ X is a subset of *derived-nodes*. Then, the notation *W'*(M) will be used to represent the vector subspace of *W'*, whose elements vanish at every *derived-node* that does not belong to M. Furthermore, corresponding to each *local subset* *of derived-nodes*, X^{α}, there is a *'local subspace of* *derived-vectors'*, *W*^{α}, which is defined by

Clearly, when whenever *β* ≠ *α*. We observe that

*W* =*W*^{1}⊕... ⊕*W*^{E} (3.4)

A derived-vector is said to be *continuous* when is independent of *α*. The set of *continuous vectors* constitute the linear subspace, *W*_{12}.

The orthogonal complement (with respect to the Euclidean inner-product) of *W*_{12} ⊂ *W* is *W*_{11} ⊂ *W*. Then *W*= *W*_{11}⊕*W*_{12}. Two projection-matrices and are here introduced; they are the projection-operators, with respect to the *Euclidean inner-product* on *W*_{12} and *W*_{11}, respectively. When , one has

the vectors and are said to be the *'jump'* and the *'average'* of , respectively. Therefore, *W*_{11} is the *'zero-average'* subspace, while *W*_{12} is the *'zero-jump'* subspace.

*Original-nodes* are classified into *'internal'* and *'interface-nodes'*: a node is *internal* if it belongs to only one subdomain-closure of the *coarse-mesh*, and it is an *interface-node* when it belongs to more than one of such closure-subdomains. Some subspaces, significant for our developments, are listed next:

At present, numerically competitive algorithms need to incorporate *restrictions* and to this end, in the *DVS-framework*, a *'restricted subspace'* *W*_{r} ⊂ *W* is selected. In the developments that follow, it is assumed that:

The matrix will be the projection-operator on *W*_{r}. We observe that when , one has . We also notice that

*W* =*W*_{Ι}⊕*W*_{Γ}=*W*_{Ι}⊕*W*_{ϖ}⊕*W*_{Δ} (3.7)

**4. The Transformed Problem**

The *transformed-problem* consists in finding such that

Where:

and

together with

The function *m* (*p, q*) is said to be the *'multiplicity'* of the pair (*p, q*). The *'derived-nodes'* are created after a *coarse-mesh* has been introduced, by dividing the *original-nodes* as explained in the Overview (Section 2), and then with each *'derived-node'* we associate a unique pair of numbers (*p,* *α*) such that and . In what follows, we identify *derived-nodes* with such pairs.

Then, in order to incorporate the constraints, we define

then, the matrix defined by

has the property that

Hence, Eq. (4.1) is replaced by

For matrices and vectors the following notation is adopted:

where the matrices

furthermore,

The matrix will be referred to as the *'transformed-matrix'*. We observe that when .

In turn, the *transformed problem* of (4.8) can be reduced, see [Herrera *et al*., 2010; Herrera *et* *al.*, 2009; Herrera, 2008; Herrera, 2007; Farhat *et* *al.*, 2000] for details, into the following problem, which is expressed in terms of the values of the solution at *dual-nodes*, exclusively: "Find (Δ) that satisfies

Here, and the *'Schur-complement* *matrix with constraints'* are defined by

and

respectively.

**5. The DVS-Algorithms**

Generally two kinds of approaches are distinguished: primal –these are direct approaches, which do not resort to Lagrange multipliers- and dual –indirect approaches that use Lagrange multipliers-. However, when DDMs are formulated using a setting as general as that supplied by the *DVS-framework*, such a distinction is irrelevant. The feature that is conspicuous for different options is the information that the algorithm seeks. Indeed, four algorithms will be obtained by seeking successively for the vectors: , , , and . However, in the presentation that follows we stick to the *'primal vs. dual-algorithms'* classification.

**5.1 Primal Formulations**

**The DVS Version of BDDC**

This is a primal algorithm which seeks directly for . Pre-multiplying Eq. (4.12) by , one gets:

In [Farhat *et al.*, 2000], it was shown that Eq. (5.1) is equivalent to Eq. (4.12). This equation is the DVS-version of BDDC.

**The DVS-Primal Algorithm**

For this algorithm, the *sought-information* is:

Applying to Eq. (5.2) it is seen that . Furthermore,

Therefore

Eq. (5.4) does not define an iterative algorithm. In order to obtain such an algorithm, we project on , to obtain:

This algorithm is referred to as the *'DVS-primal* *algorithm'*. The solution is given by

We observe that we could have written instead of Eq. (5.6). However, the application of the projection operator is important when and are not computed with exact arithmetic, as it is the case when using numerical methods, because when it is applied it replaces by the continuous-vector closest (with respect to the Euclidean distance) to it.

**5.2 Dual Formulations**

**The DVS Version of FETI-DP**

For this algorithm the *sought-information* is defined to be: . This algorithm can be easily derived from the *DVS-primal* formulation that has just been presented. We observe that , , in view of Eq. (5.2), and . This permits transforming Eq. (5.5) into

Applying to the first of these equations, it is obtained:

As for Eq. (5.6), it becomes:

**The DVS-Dual Algorithm**

In this algorithm, the *sought-information* is: . Then, . Replacing this in Eq. (5.1), one gets:

Finally, multiplying by the first of these equalities, it is obtained:

When is known, can be recovered by means of

A comment similar to that made immediately after Eq. (5.6), goes here: we have applied the projection matrix , in Eq. (5.12) because we are assuming that exact arithmetic generally will not be used.

**6. Numerical Procedures Fulfilling the DDM-Paradigm**

Summarizing, the preconditioned *DVS-algorithms* with constraints are:

**6.1 Comment on the DVS Numerical Procedures**

The outstanding uniformity of the formulas given in Eqs. (6.1) to (6.4) yields clear advantages for code development, especially when such codes are built using object-oriented programming techniques. Such advantages include:

I. The construction of very robust codes. This is an advantage of the

DVS-algorithms, which stems from the fact the definitions of such algorithms exclusively depend on the discretized system of equations, obtained after discretization of the partial differential equations considered (referred to as theoriginal problem), but which is otherwise independent of the problem that motivated it. In this manner, for example, essentially the same code was applied to treat 2-D and 3-D problems; indeed, only the part defining the geometry had to be changed, and that was a very small part of it;II. The codes may use different local solvers, which can be direct or iterative solvers;

III. Minimal modifications are required for transforming sequential codes into parallel ones; and

IV. Such formulas also permit developing codes which fulfill the

DDM-paradigm; i.e., in which "the solution of the global problem is obtained by exclusively solving local problems".

This last property makes the DVS-algorithms very suitable as a tool to be used in the construction of massively-parallelized software, so much needed for efficiently programming the most powerful parallel computers available at present. In the next Subsection, procedures for constructing codes possessing Property IV are explained with some detail.

All the DVS-algorithms of Eqs. (6.1) to (6.4) are iterative and can be implemented with recourse to Conjugate Gradient Method (CGM), when the matrix is definite and symmetric, or some other iterative procedure such as GMRES, when that is not the case. At each iteration step, depending on the *DVS-algorithm* that is applied, one has to compute the action on a *derived-vector* of one of the following matrices: , , , or . Such matrices in turn are different permutations of the matrices , , , and . Thus, to implement any of the preconditioned *DVS-algorithms*, one only needs to separately develop codes capable of computing the action of each one of the matrices , , or on an arbitrary *derived-vector*, of *W*.

Therefore, next we present numerical procedures for computing the application of each one of the matrices , , and , which fulfill the *DDM-paradigm*. It will be seen that only requires exchange of information between derived-nodes belonging to different subdomains; actually, between *derived-nodes* that are descendants of the same *original-node* (the exchange of information is minimal). As for , once the action of has been computed, no further exchange of information is required.

**6.2 Application of **

From Eq. (4.13), we recall the definition of the matrix . In order to evaluate the action of on any *derived-vector*, we need to successively evaluate the action of the following matrices , , and . Nothing special is required except for . A procedure for evaluating the action of this matrix, which fulfills the *DDM-paradigm* is explained next.

We have

Let , be an arbitrary *derived-vector*, and write

Then, is characterized by

and can obtained iteratively. Here,

and, with as the projection-matrix into ,

We observe that fulfilling the *DDM-paradigm* when computing the action of is straightforward because

is parallelizable. Once has been obtained, to derive one can apply:

this completes the evaluation of .

**6.3 Application of **

We define

and observe that

Therefore, the matrix can be written as:

Furthermore, fulfills

Another property that is relevant for the following discussion is:

*W _{r}*(Σ) =

*W*(Σ) (6.15)

for any , let us write

then, fulfills

Here, , where the matrix is the projection operator on *W*_{r}, while

Furthermore, we observe that

In order to use Eq. (6.19) as a means of parallelizing the DVS-algorithms, however, the detailed discussion of such procedures will be presented separately [Herrera *et al*., 2013; L.M. de la Cruz *et al.*, 2013]. It is necessary that the local matrices, , be invertible. This is granted when invertible in *W*_{r}, which generally is achieved by taking a sufficiently large number of *primal-nodes*.

Eq. (6.17) is solved iteratively. Once has been obtained, we apply:

This procedure permits obtaining in full; however, we only need . We observe that

The vector can be obtained by the general procedure presented above. Thus, take and

Therefore,

**6.4 Application of and .**

We use the notation

then [Herrera *et al.*, 2010]:

while therefore,

Therefore, only the evaluation of requires exchange of information between subdomains. In general, such numbers are very small; for example in application to single-equation problem, when an orthogonal grid is used, they are at most: 4, for problems in 2D, and 8 for problems in 3D.

As for the right hand-sides of Eqs. (4.14), all they can be obtained by successively applying to some of the operators that have already been discussed. Recalling Eq. (4.14), we have

The computation of does not present any difficulty and the evaluation of the actions of and were already analyzed.

**7. Numerical Results**

Taking into account the general description of the DVS-framework given of Section 2, it can be seen that each one of the DVS-algorithms is uniquely defined by:

1. The original-matrix;

2. The partition of the set of original-nodes, which is induced by the

coarse-meshthat is applied; and3. The set of constraints.

In turn, the original-matrix is determined by the partial differential equation, or system of such equations, the discretization method chosen and the *fine-mesh* adopted. As explained in Section 2, the partition of the set of original-nodes depends when the *fine-mesh* has already been defined, on the *coarse-mesh* (i.e., the domain decomposition) used. The *coarse-mesh* is constituted by a family of non-overlapping subdomains {Ω_{1},..., Ω_{E}} of Ω, the domain of definition of the boundary-value problem to be solved. In all the examples that are presented in this article, the *constraints* are fully determined by the *primal-nodes* and consist in requiring continuity of the *derived-vectors* at them.

Several codes were developed to treat the examples, which were written in C++ language, using the MPI library for the communications. In the computational implementations, the methods of solution used to treat the *original-problems* are: CGM, when such a linear system is symmetric and positive-definite and GMRES when the discrete system is non-symmetric or indefinite. Both are applied with a tolerance of 10^{-6}. Each *DVS-algorithm* was applied to each one of the examples considered, except for that referring to elasticity.

The results obtained for Examples 1 to 5 are summarized in Table 1, Table 2, Table 3, Table 4 and Table 5, respectively. In them, the acronym *dof* stands for to the number of degrees of freedom of the *original problem*, but it should be mentioned that the procedures used to treat such examples are such that the nodes that lie on the external boundary do not contribute to the *dof*. The notation to indicate the meshes that were adopted is as follows: In 2D cases, we use *(nxm)x(qxr),* where *(nxm)* refers to the *coarse-mesh*, while *(qxr)* to the *fine-mesh*; and similarly, in 3D cases, we use *(nxmxp)x(qxrxs),* where *(nxmxp)* define the *coarse-mesh* and *(qxrxs)* the *fine-mesh.* The constrains are imposed on the primal nodes, in all of our experiments the primal nodes were located at vertex in 2D and at edges in 3D of the subdomains, this coinciding with the algorithm "D" in [Toselli *et al.*, 2005].

Each Table contains at most ten columns. The first four indicate respectively: 1) the meshes used, 2) the number of subdomains of the *coarse-mesh*, 3) the *dof*, and 4) the number of *primal-nodes* used. The figures appearing in columns 5 to 9 correspond to the number of iterations that were required for convergence of each one of the algorithms applied. Columns 9 and 10 were only included in Table 3. For Example 3, in order to cover a wide range of values of the Peclet-number, the diffusion coefficient in Eq. (7.3), *v*, was varied and the tenth column in Table 3 indicates the different values of *v* for which the corresponding boundary-value problem was solved. Furthermore, the results obtained when the DVS-algorithms were applied were compared with those obtained in [Da Conceição *et al*., 2006] for the same problem, using the standard version of BDDC.

**7.1 Application of the DVS-algorithms to a Single-Equation**

The applicability of the *DVS-algorithms* is wide, as previously said it can be applied to general equation systems. In Section 3, it was announced that in this paper we present examples for which *d*, the number of equations of the system, is one and three. In this Subsection the examples for which *d*=1 will be discussed, leaving for the next Subsection the treatment of static-elasticity models, for which *d*=3.

Four boundary value problems corresponding to a single-equation will be presented. The first two are symmetric and positive definite boundary-value problems, whose definition involves the Laplace differential operator. The other two correspond to advection-diffusion transport, and the corresponding boundary-value problems are non-symmetric and indefinite. The discretization methods used in this Subsection are based on central finite differences (CFD), which are directly applicable to the symmetric problems. To apply CFD to the advection-diffusion problems it was necessary to stabilize the advection-diffusion differential-operator and to this end artificial diffusion was incorporated.

Despite the simplicity of the examples presented in this Subsection, they are very important because a wide range of geophysical systems give rise to similar problems [Herrera and Pinder, 2012]. The diversity of physical interpretations of the boundary-value problems here discussed is enormous. All the differential operators involved can be classified as advection-diffusion operators, since Laplace operator is obtained from the general advection-diffusion differential-operator when the transport-velocity vanishes. Transport processes of heat and solutes occur in a great diversity of geophysical systems. However, the physical processes governed by such differential-equations go far beyond transport phenomena.

Example 1. Poisson equation in two-dimensions.

We can see from Table 1, that the four algorithms perform very well as the number of subdomains and the degrees of freedom (dof) are increased. In this example, the DVS-DUAL algorithm presents the best performance, requiring only 11 iterations from 12x12 until 30x30 subdomains, and the same number of dof. All other algorithms show similar behavior. The numerical solution of this example can be seen in the Figure 5.

Example 2. Similar to Example 1, but it is formulated in a 3D domain.

In Table 2, we observe a similar performance of the algorithms as in the two-dimensional case. One more time the DVS-DUAL algorithm presents a little better behavior with respect all others.

Example 3. The boundary-value problem treated is:

This is an advection-diffusion transport problem in 2D, for which the differential operator is not self-adjoint.

This example is very interesting because it contains diffusion and advection terms, which are common in several complex geophysics phenomena. In this example, the Péclet number is defined as , where L is a characteristic length (in this case L = 1). We also define a local Péclet number as . Using these definitions, fixing the global partition to h=1/512, and the varying the viscosity from 0.01 to 0.0001, we have that the Péclet number varies from 316 to 316,227, and the local Péclet number varies from 0.617 to 617. In this case the linear system is non-symmetric, therefore we choose the GMRES method with a tolerance of 10^{-6}.

In Table 3 presents the results that the *DVS-algorithms* yielded and compares them with those obtained in [Da Conceição *et al.*, 2006]. We observe that, with fixed *coarse* and *fine* *meshes*, as the viscosity coefficient is reduced, so that the Péclet number increases, generally the iterations required for convergence reduce. Increasing the Péclet number implies that the effect of the advection term enlarges, and the numerical solution generally becomes unstable. However, the performance of the discretization strategy based on CFD combined with stabilization of the numerical-scheme by means of artificial viscosity is resilient to Péclet-number variations. For comparison purposes, the examples presented here were chosen to be the same as those presented in [Da Conceição *et al*., 2006], where the standard BDDC algorithm was applied with the same set of constraints; namely, the same set of subdomains and vertex nodes were chosen to be *primal*. As can be seen in Table 3, when the comparison criterion is based on the number of iterations required for convergence, the observed performance of the *DVS-algorithms* in these examples is slightly better than that of the standard BDDC algorithm. Finally, an illustration of the kind of numerical solution obtained is shown in Figure 6 and Figure 7.

The relative-residual decay for a coarse mesh (16X16) and several fine meshes is presented in Figure 8. We consider in these computations b=(1,3) and *v*=0.00001, in such a way that P*e*=3.16*e*+5. We observe that the best convergence is obtained when the fine mesh is increased, and the convergence slows when the dof occurring in the subdomains is reduced.

Example 4. The boundary-value problem treated is:

This is an advection-diffusion transport problem in 3D, for which the differential operator is not self-adjoint.

The diffusion and advection-diffusion differential-operator appears in the equations of the examples presented above. They are very important in natural and industrial phenomena. For example, the flow and transport of solutes in subsurface groundwater, the movement of aerosol and trace gases in the atmosphere, mixing of fluids in processes of crystal growth, among many other important applications [Tood, 1980; Pinder *et al*., 2006; Herrera *et al*., 1969; Herrera *et al*., 1973; Herrera *et al*., 1977; Herrera G.S. *et al*., 2005; L.M. de la Cruz *et al*., 2006]. In all our examples, we have shown that the *DVS-algorithms* obtain the numerical solution efficiently on parallel machines. In this respect, we remark that for advection-diffusion problems the matrices of the discrete linear systems are non-symmetric.

**7.2 Application to a System-Equations**

We use the *DVS-framework* to solve a Dirichlet boundary value problem, where displacements are zero over the boundary of the elastic body that occupies the domain Ω of the physical space. Over each one of such subdomains is solved a local problem by FEM, using linear functions as basis. On each node a of the mesh is defined a vector valued function with each component identified as *u*_{ai} for *i*=1, 2, 3.

Because our operators are symmetric and positive definite, we use CGM as an iterative procedure to solve those linear systems of equations that we have defined in the DVS framework.

The code used in the previous section, which was originally developed to solve a single equation using finite differences, was adapted for solving systems of equations with FEM. We added the corresponding functionality in order to be able to solve systems of equations, in this case the elasticity problem.

Example 5. A system of partial differential equations in three-dimensions has also been treated. This is the system of differential equations of static elasticity; namely:

which was subject to the following Dirichlet boundary conditions:

The domain of study for our numerical experiments is a homogeneous isotropic linearly elastic unitary cube. In all of our experiments the primal nodes were located at edges of the subdomains, which is enough for not being singular.

We consider constant coefficients l and m equal to one. With these conditions we have a problem that has analytical solution, and is written as follows:

= (sin*π*x sin*π*y sin*π*z, sin*π*x sin*π*y sin*π*z, sin*π*x sin*π*y sin*π*z) (7.7)

The Tables 5, summarizes the numerical results obtained using the DVS methods with a tolerance of 10^{-7}.

**8. Conclusions**

Mathematical models of many geophysical systems lead to a great variety of partial differential equations (PDEs) whose solution methods are based on the computational processing of large-scale algebraic systems [Herrera and Pinder, 2012]. Parallel computing is outstanding among the new computational tools and, in order to effectively use the most advanced computers available today, massively parallel software is required. Domain decomposition methods (DDMs) have been developed precisely for effectively treating PDEs in parallel [DDM Organization, 2012]. What domain decomposition methods ideally intend to do has been summarized in this paper in the *"DDM-paradigm"*: to develop algorithms that *'obtain the* __global__ solution by*exclusively solving* * local problems'*.

In conclusion, in this paper:

1. We have presented a *non-overlapping* *discretization* method (the *DVS-discretization*) -in the sense that it uses a system of nodes such that each one of them belongs to one and only one subdomain of the *coarse-mesh*- applicable to a wide class of well-posed boundary problems associated with elliptic systems of equations. In particular, the differential operators may be symmetric, non-symmetric or indefinite (non-positive-definite);

2. Four algorithms –the *DVS-algorithms* [Herrera *et al.*, 2011]*-*, which were derived using the *DVS-discretization* and achieve the *DDM-paradigm* have been explained. Two of them are the result of using the BDDC and FETI-DP algorithms after applying *DVS-discretization* to the boundary value problem considered. The other two are obtained when two new algorithms, which had not been reported previously in the literature, were used instead;

3. Numerical procedures that permit achieving the *DDM-paradigm* with each one of the *DVS-algorithms* have been also presented;

4. Codes were developed and applied to several boundary values problems that occur in the modeling of certain geophysical phenomena, such as transport of solutes by both, free-fluids and fluids in a porous medium. We also present results for a static elasticity problem, which thereby illustrates the application of the algorithms to systems of differential equations; and

5. Besides their attractive parallelization properties, in the numerical examples the *DVS-algorithms* exhibited significantly improved numerical performance with respect to standard versions of BDDC and FETI-DP.

**Acknowledgement**

The authors express their gratitude to Alberto Rosas-Medina e Iván Contreras-Trejo, both PhD students of the Earth-Sciences Graduate Program at UNAM, for having permitted us to reproduce some numerical results of their research work.

Luis M. de la Cruz wishes to acknowledge the support from the project PAPIIT-UNAM TB100112 to develop this research.

**Bibliography**

Da Conceição D.T. Jr., 2006, Balancing domain decomposition preconditioners for non-symmetric problems, Instituto Nacional de Matemática Pura e Aplicada, Agência Nacional do Petróleo PRH-32, Rio de Janeiro, May. 9. [ Links ]

DDM Organization, 2012, Proceedings of 21 International Conferences on Domain Decomposition Methods. http://www.ddm.org. [ Links ]

De la Cruz L.M., Herrera I., 2013, Generic and Parallel Software based on DVS algorithms for engineering, Submitted at Advances in Engineering Software. [ Links ]

De la Cruz L.M., Ramos E., 2006, Mixing with time dependent natural convection, Int. Comm. in Heat and Mass Transfer, 33, 2, pp 191-198. [ Links ]

Dohrmann C.R., 2003, A preconditioner for substructuring based on constrained energy minimization, SIAM *J. Sci. Comput*. 25, 1, 246-258. [ Links ]

Farhat Ch., Roux F., 1991, A method of finite element tearing and interconnecting and its parallel solution algorithm, *Internat. J.* *Numer. Methods Engrg.* 32:1205-1227. [ Links ]

Farhat C., Lessoinne M. LeTallec P., Pierson K., Rixen D., 2001, FETI-DP a dual-primal unified FETI method, Part I: A faster alternative to the two-level FETI method, *Int. J. Numer.* *Methods Engrg.* 50, pp 1523-1544. [ Links ]

Farhat C., Lessoinne M., Pierson K., 2000, A scalable dual-primal domain decomposition method, *Numer. Linear Algebra Appl.,* 7, pp 687-714. [ Links ]

Herrera I., Rosas-Medina A., 2013, The Derived-Vector Space Framework and Four General Purposes Massively Parallel DDM Algorithms, Engineering Analysis with Boundary Elements, in press. [ Links ]

Herrera I., De la Cruz L.M., Carrillo-Ledesma A., Rosas-Medina A., Contreras I., 2012, Foundations of the DVS-framework: Theory and algorithms, Memoria No. 8 del Grupo de Modelación Matemática y Computacional del Instituto de Geofísica, UNAM. [ Links ]

Herrera I., Carrillo-Ledesma A., Rosas-Medina A., 2012, Four general purposes massively parallel DDM algorithms, available as Memoria No.7 del Grupo de Modelación Matemática y Computacional del Instituto de Geofísica, UNAM. [ Links ]

Herrera I., Carrillo-Ledesma A., Rosas-Medina A., 2011, A brief overview of non-overlapping domain decomposition methods, *Geofísica* *Internacional*, 50, 4, pp 445-463. [ Links ]

Herrera I., Yates R.A., 2009, The multipliers-free dual primal domain decomposition methods for non-symmetric Matrices, NUMER. METH. PART D. E. 2009, DOI 10.1002/Num. 20581. [ Links ]

Herrera I., Yates R.A., 2010, The multipliers-free domain decomposition methods, NUMER. METH. PART D. E. 26: 874-905 July 2010, DOI 10.1002/num. 20462. [ Links ]

Herrera I., Yates R., 2009, Unified multipliers-free theory of dual-primal domain decomposition methods, NUMER. METH. PART D. E. Eq. 25:552-581, May 2009, (Published on line May 13, 08) DOI 10.1002/num. 20359. [ Links ]

Herrera I., 2008, New formulation of iterative substructuring methods without Lagrange Multipliers: Neumann-Neumann and FETI, NUMER METH PART D E 24, 3, pp 845-878, DOI 10.1002 NO. 20293. [ Links ]

Herrera I., 2007, Theory of differential equations in discontinuous piecewise-defined-functions, NUMER METH PART D E, 23, 3, pp 597-639, DOI 10.1002, 20182. [ Links ]

Herrera I., Pinder G.F., 2012, Mathematical modeling in science and engineering: An axiomatic approach, Wiley, 243 p. [ Links ]

Herrera I., Figueroa V. G.E., 1969, A correspondence principle for the theory of leaky aquifers, *Water Resources Research*, 5, 4, P. 900. [ Links ]

Herrera I., Rodarte L., 1973, Integrodifferential equations for systems of leaky aquifers and applications, Part 1: The nature of approximate Theories, *Water Resources Research*, 9, 4, pp. 995-1005. [ Links ]

Herrera I., Yates R., 1977, Integrodifferential equations for systems of leaky aquifers. Part 3. A numerical method of unlimited applicability, *Water Resources Research*, 13, 4, pp. 725-732. [ Links ]

Herrera G.S., Pinder G.F., 2005, Space-time optimization of groundwater quality sampling networks, *Water Resources Research*, 41, W12407, 15 pp. [ Links ]

Li J., Widlund O., 2005, FETI-DP, BDDC and block Cholesky methods, *Int. J. Numer. Methods* *Engrg.,* 66, 250-271. [ Links ]

Mandel J., 1993, Balancing domain decomposition, Commun. *Numer.* *Methods Engrg.* 233-241. [ Links ]

Mandel J., Brezina M., 1996, Balancing domain decomposition for problems with large jumps in coefficients, *Math. Comput*. 65, pp 1387-1401. [ Links ]

Mandel J., Tezaur R., 1996, Convergence of a substructuring method with Lagrange multipliers, Numer. Math 73, 4, 473-487. [ Links ]

Mandel J., Tezaur R., 2001, On the convergence of a dual-primal substructuring method, SIAM *J. Sci. Comput*., 25, pp 246-258, 2001. [ Links ]

Mandel J., Dohrmann C.R., 2003, Convergence of a balancing domain decomposition by constraints and energy minimization, Numer. *Linear Algebra Appl*., 10, 7, 639-659, 2003. [ Links ]

Mandel J., Dohrmann C.R., Tezaur R., 2005, An algebraic theory for primal and dual substructuring methods by constraints, *Appl.* *Numer. Math*., 54: 167-193. [ Links ]

Pinder G.F., Celia M.A., 2006, Subsurface hydrology, Wiley, 468 p. [ Links ]

PITAC, 2005, Computational Science: Ensuring america's competiveness, Report to the President of the United States, President Information, Technology Advisory Committee, Executive Office of the President of the United States, June. [ Links ]

Todd D.K., 1980, Groundwater hydrology, 2^{nd} ed. Wiley. [ Links ]

Toselli A., Widlund O., 2005, Domain decomposition methods-algorithms and Theory, Springer Series in Computational Mathematics, Springer-Verlag, Berlin, 450 p. [ Links ]

**Note**

^{1}In order to mimic standard notations, we should have used *Π* instead of the *low-case* *ϖ*. However, the modified definitions given here yield some convenient algebraic properties.