Summary
The finite element method is a powerful tool for 3 D DC resistivity modelling and inversion. The solution accuracy and computational efficiency are critical factors in using the method in 3 D resistivity imaging. This paper investigates the solution accuracy and the computational efficiency of two common element type schemes: trilinear interpolation within a regular 8 node solid parallelepiped, and linear interpolations within six tetrahedral bricks within the same 8 node solid block. Four iterative solvers based on the pre conditioned conjugate gradient method (SCG, TRIDCG, SORCG and ICCG), and one elimination solver called the banded Choleski factorization are employed for the solutions. The comparisons of the element schemes and solvers were made by means of numerical experiments using three synthetic models. The results show that the tetrahedron element scheme is far superior to the parallelepiped element scheme, both in accuracy and computational efficiency. The tetrahedron element scheme may save 43 per cent storage for an iterative solver, and achieve an accuracy of the maximum relative error of < 1 per cent with an appropriate element size. The two iterative solvers, SORCG and ICCG, are suitable options for 3 D resistivity computations on a PC, and both perform comparably in terms of convergence speed in the two element schemes. ICCG achieves the best convergence rate, but nearly doubles the total storage size of the computation. Simple programming codes for the two iterative solvers are presented. We also show that a fine grid, which doubles the density of a coarse grid, will require at least 27 = 128 times as much computing time when using the banded Choleski factorization. Such an increase, especially for 3 D resistivity inversion, should be compared with SORCG and ICCG solvers in order to find the computationally most efficient method when dealing with a large number of electrodes.
Introduction
The DC resistivity method is a well established and cost effective tool in geophysical exploration. It is widely used in mining exploration, civil and hydrologic engineering, and environmental investigations. Over the past decade, much effort has been put into understanding and interpreting 3 D resistivity data, with new modelling and inversion techniques being developed (Park & Van 1991; Shima 1992; Li & Oldenburg 1992; Ellis & Oldenburg 1994; Sasaki 1994; Dabas et al. 1994; Spitzer 1995; Zhang et al. 1995; Loke & Barker 1995, 1996; Zhao & Yedlin 1996).
Techniques for the numerical forward modelling of 3 D resistivity may be classified into four types in terms of the principles used: the integral equation method, transmission network method, finite difference method (FDM) and finite element method (FEM). The integral equation methods (Dieter et al. 1969; Lee 1975; Okabe 1981; Das & Parnasis 1987; Xu et al. 1988) use the Green's function of a uniform medium directly and are most suitable for simple model geometries. The transmission network method approximates the governing equations by a network of lumped impedances, where the impedances are proportional to the model resistivities (Zhang et al. 1995). The FDM approximates the differential equations with difference equations (Mufti 1976; Dey & Morrison 1979; Mundry 1984; James 1985; Lowry et al. 1989; Spitzer 1995; Zhao & Yedlin 1996), and the FEM converts the boundary value problem into the minimization of an integral functional (Coggon 1971; Fox et al. 1980; Pridmore et al. 1981; Holcombe & Jiracek 1984; Queralt & Marcuello 1991). The last three methods enable one to investigate the complicated potential field of a 3 D arbitrary medium. In particular, the FEM has a solid theoretical foundation that gives added reliability. Its main advantage as compared with other numerical methods is that complicated geometries, general boundary conditions, and spatially variable or non linear material properties can be handled relatively easily (Pridmore et al. 1981; Sasaki 1994). Furthermore, it does not suffer from a singularity problem at the source point in the numerical resistivity modelling, as the source singularity is naturally and effectively smoothed out by finding the weighted residual solution in an integral sense, or minimizing an integral formed functional (a refinement of the FDM is required in order to remove the source singularity, see Lowry et al. 1989 and Zhao & Yedlin 1996).
To implement the FEM, two simple choices for a 3 D discretization—parallelepiped or tetrahedron elements—are used. They have different shape functions for approximating the solution (Schwarz 1988). Each element may involve many nodes or a high order polynomial for improvement of the accuracy, which leads to a number of different finite element schemes. For 3 D resistivity modelling or inversion, the simplest scheme is to employ trilinear interpolation within a solid parallelepiped or linear interpolation within a tetrahedron based on the basic nodes, so as easily to minimize the dimension of 3 D computations. Pridmore et al. (1981) and Sasaki (1994), for example, using the basic nodes, presented a tetrahedron element scheme for 3 D DC resistivity modelling and inversion. In general, different element schemes will yield different accuracies and computing efficiencies of the FEM.
Another important factor is the type of equation solver used for the linear equation system of 3 D FEM. The different solvers have quite different computing efficiencies and consumption of computer resources. Two approaches, an iterative solver and an elimination solver, have been recognized in FEM theory (Schwarz 1988). In fact, the two kinds of solvers are complementary. The former requires much less computer storage than the latter, but the solution procedure has to be repeated for different source vectors. The latter often needs huge amounts of memory, but it is only implemented once for different source vectors. The choice of the solver depends on the available computer resources and the actual computing task. However, no matter what kind of solver is employed, the convergence rate should be rapid. The elimination solver uses the minimum computer memory and requires the least computer time. Most researchers have preferred an iterative solver so as to utilize the sparsity of the coefficient matrix in 3 D simulations. In particular, it has been recognized in structural analysis that the pre conditioned conjugate gradient method (PCG) has better convergence features than other classical iterative solvers (David 1978; Ajiz & Jennings 1984; Manolis & Michael 1991). In 3 D DC resistivity modelling, Spitzer (1995) applied one algorithm (SSORCG) of the PCG to the FDM, and Zhang et al. (1995) applied a different algorithm (ICCG) of the PCG to the transmission network method.
In order to achieve an accurate solution and the maximum computing efficiency, this paper investigates the performance of the two different element schemes—trilinear interpolation within the 8 node solid parallelepiped and linear interpolations within six tetrahedral bricks within the same 8 node solid block—for 3 D DC resistivity modelling. In addition, we have applied four iterative algorithms of the pre conditioned conjugate gradient method and one elimination solver to the FEM in order to compare their computing accuracies, efficiencies and computer storage requirements. By means of numerical experiments on three synthetic models, it was found that the tetrahedron element scheme is far superior to the parallelepiped element scheme, both in accuracy and computing efficiency, and two iterative PCG solvers are suitable options for 3 D FEM resistivity modelling on a PC.
Basic equations
3 D DC resistivity modelling reduces to solving the following mixed boundary value problem:
where G is the Green's function of the electric potential due to a unit current source I = 1 (magnitude) at position rc. The quantity σ is the conductivity of the medium whose spatial range and boundary are denoted by Ω and ∂Ω. In general, the conductivity is a function of the spatial coordinates r = (x, y, z). Here ν is a coefficient due to the artificial boundary (Dey & Morrison 1979), and it may be written in general form when considering a buried source:
where r = |
| and r′ = |
′| are the distances from the current source and its image position (due to the Earth's surface) to the boundary; and θ and θ′ are the angles between the normal vector of the boundary and the vectors 
and 
′, respectively. On the Earth's surface, we set ν = 0 because of ∂G/∂n = 0.
Accordingly, we have the following expressions for the potential, apparent resistivity and Fréchet derivative based on the Green's function (Zhou & Greenhalgh 1999):
where δG(r) = G(r, r) − G(r, r), K is the geometry factor, which can be calculated from the positions of the four electrodes rA, rB, rM and rN in a configuration (i.e. pole–pole, pole–bipole, bipole–pole and bipole–bipole), and σe is the conductivity of a model element Ωe. In eqs (4) and (5), if a remote electrode (infinite location r∞) is involved in the computation, we set G(r, r∞) = 0. Eqs (4) and (5) are thus available for any electrode configuration. Eqs (3), (4) and (5) show that in theory the potential, the apparent resistivity and the Fréchet derivative can be calculated by the Green's function, which can be obtained by solving eq. (1).
The FEM can be used in solving eq. (1). According to finite element theory, two principles, namely the variational principle (Pridmore et al. 1981) and Galerkin's criterion (Zienkiewicz 1971), are available for various applications. It has been shown that, if the same shape function is employed in the two schemes, they lead to the same linear equation system. Galerkin's criterion, however, gives a simpler formulation and is easier to understand. In the following section, we adopt Galerkin's criterion for the FEM formulation.
Finite element solution
According to Galerkin's criterion, the solution of eq. (1) reduces to solving the following integral equation:
where Wj is a weighting function that belongs to Hilbert space H1(Ω). After carrying out integration by parts and substituting for the boundary condition in eq. (1), eq. (6) becomes
By discretizing the 3 D region into finite elements, for example Ω = ΣeΩe and ∂Ω = ΣeΓe, in which the Green's function can be expressed as a combination of the shape functions Nie(r), i.e.
where n is the total number of shape functions, one obtains the discrete form of eq. (7):
The weighting function in Ωe can be chosen to be the same as the shape function, that is Wje(r) = Nje(r), j = 1, 2,…, n, so that eq. (9) can be written in the matrix form
Here, G˜ is the vector of Green's functions, bs = (0, 0,…, 1,…, 0) is the source vector with only one non zero component, at the current injection position, and M is an n × n square matrix assembled from the individual element matrices, namely
which depends on the element integral over Ωe and the boundary integral on Γe. Obviously, the element matrix Mije is symmetric. Once the shape functions [Ni(r), i = 1, 2,…, n] are chosen and the conductivity σ is given, the element matrix can be calculated using eq. (11), and one can obtain the Green's function vector G˜ by solving eq. (10). The electric potential or the apparent resistivity at any observational point can be calculated with the Green's function with expressions (3) or (4). Therefore, the FEM for resistivity modelling reduces to calculating the matrix M and solving the linear equation system (10).
Element matrix calculations
In order to calculate the element matrix (see eq. 11), one must choose the elements {Ωe}, which together compose the whole model range Ω. For 3 D problems, two common element types—parallelepiped and tetrahedron elements—are often employed for the computation of the element matrix.
Parallelepiped element
A parallelepiped element is shown in Fig. 1(a), simply represented by eight nodes {I1, I2,…, I8}. If trilinear interpolation is used within the element, each node has the following shape function (Schwarz 1988):
Two common element specifications in the 3 D finite element method: (a) parallelepiped element; (b) tetrahedron element.
where
and
The coordinates (x1, x2), (y1, y2) and (z1, z2) represent the range of the parallelepiped element.
Substituting eq. (12) into eq. (11), the element integral becomes
from which one can see that, due to the symmetry between i and j, only 36 independent non zero components need to be calculated for each parallelepiped element.
To calculate the boundary integral in eq. (11), one may suppose that the global boundary consists of the Earth's surface and five subsurface planes—the left, the right, the front, the back and the bottom boundaries. On the Earth's surface the integral vanishes due to ν = 0. On the others it gives rise to five boundary integrals to be calculated. As an example, the calculation for the left boundary is given. Setting η = 0 in eq. (12) and substituting into the boundary integral of eq. (11), we have
from which one can see that only 10 components need to be calculated, as a result of the symmetry. The contributions from the other boundaries can be obtained in a similar way.
Tetrahedron element
A tetrahedron element is shown in Fig. 1(b), simply represented by four nodes (i, j, l, m). Pridmore et al. (1981) and Sasaki (1994) presented a tetrahedron element scheme in which the 3 D region is divided into a large number of bricks and each brick is assembled from five tetrahedron elements. Zhou & Zhong (1984) presented another tetrahedron element scheme in which each of the bricks is assembled from six tetrahedron elements. Both use a linear shape function within the tetrahedron for 3 D resistivity modelling. Here, we use the six tetrahedron element scheme and briefly outline the procedure. The six tetrahedrons can be denoted by their nodes (see Fig. 1): (I1, I2, I3, I5), (I2, I3, I5, I6), (I3, I5, I6, I7), (I2, I3, I4, I6), (I3, I4, I6, I7) and (I4, I6, I7, I8). So, from eq. (11) the integral over the brick is calculated by the summation of the six tetrahedron integrals:
where ΓeΔ represents the side of the tetrahedron element that coincides with the artificial boundary. For simplicity, in the following text we use (i, j, l, m) to represent the four nodes of a tetrahedron element, and their shape functions can be written as follows:
which satisfies Np(xq, yq, zq) = δpq (p, q = i, j, l, m). Here Ve is the bulk of the tetrahedron element, and ap, bp, cp and dp are constants which can be calculated from the node coordinates {(xp, yp, zp), p = i, j, l, m}. Substituting (16) into the bulk integral of eq. (15), we have
which is much simpler than eq. (13).
The boundary integral in eq. (15) can be calculated in a similar way to in the parallelepiped element scheme, but using the shape function given by eq. (16) (see Zhou & Zhang 1984).
Linear equation system solvers
The final step of the FEM is to solve the linear equation system (10). From the previous section, the coefficient matrix M is a large, sparse, symmetric and positive definite matrix. One can employ a direct matrix method or an iterative solver, depending on the available computer resources and the task of the computation. This section is concerned with two kinds of solvers applied to our case. One is called diagonal band elimination, specifically the banded Choleski factorization method, and the other is the pre conditioned conjugate gradient method (PCG) with different pre conditioners. They each use different amounts of computer resources and offer different efficiencies for 3 D resistivity modelling.
Diagonal band elimination
The attractive feature of diagonal band elimination is that the solving procedure is conducted only once for all the source vectors (different vectors bs in eq. 10). This advantage suggests that one should apply a band elimination method to solve eq. (10) for a large number of source vectors. In our case, it may be an option for 3 D resistivity inversions, because the Green's functions of all the current and potential electrodes in a 3 D measurement must be calculated for the Fréchet derivative (see eq. 5). A 3 D resistivity measurement may employ hundreds of electrodes. This means that, if applying an iterative solver to eq. (10), the iterative procedure has to be repeated hundreds of times for all the Green's functions. It may be that the computer time required for implementing the band elimination once is less than that required when employing an iterative solver for a large number of electrodes.
Considering the various algorithms of band elimination methods, the banded Choleski factorization algorithm is a suitable choice for this case, because it can deal with only the lower half of the symmetric band matrix M, thus saving a significant amount of memory (Schwarz 1988). The banded Choleski factorization algorithm implements the factorization M = LLT, where L is a lower triangular matrix but has the same bandwidth m as M (Mij = 0, Lij = 0 for |i − j| > m). Therefore, only one array, with size (n, m + 1), is required for storage, and then by performing forward and back substitution one can obtain all the solutions for the different vectors bs. The amount of calculation necessary for the banded Choleski factorization is directly proportional to the order n and to the square of the bandwidth m (Schwarz 1988); that is,
It is therefore particularly important that the bandwidth m should be kept as small as possible for lower storage requirements and minimal computation.
2010 Pre-conditioned conjugate gradient method (PCG)
The basic idea of the PCG is to multiply a matrix M˜−1 (M˜ is called a pre conditioner) with the linear equation system (10) so that the resultant coefficient matrix is close to the identity matrix M˜−1M≈I or M˜≈M and the CG algorithm has a fast convergence rate for the linear equation system M˜−1Mx=M˜−1bs. Applying this idea to the CG algorithm, we write the PCG algorithm as follows:
- 1
initialization: solving M˜r0=b−Mx0 and setting p0 =r0;
- 2
loop for i = 0, 1, 2, 3,….
19
from which one can find that the pre conditioner M˜ should be close to M and has some specific properties for efficiently solving the linear equation M˜qi = Mpi. Simple choices for M˜ may be the diagonal matrix of M: M˜ = diag(M11, M22,…, Mnn) or the tridiagonal partition of M: M˜ = {M11, M12, M21, M22, M23,…, M(n−1)n, Mnn}. The former is actually the scaled conjugate gradient method (SCG), and the latter (called TRIDCG) has more components than the diagonal choice; with both of them it is easy to obtain the qi in the algorithm (19). A further two choices for M˜ are the symmetric successive over relaxation matrix method (SSORCG; see Axelsson 1984; Spitzer 1995) and the incomplete Cholesky decomposition of M (ICCG; see David 1978; Ajiz & Jennings 1984; Manolis & Michael 1991; Zhang et al. 1995).
Following the SSORCG, we rewrite the pre conditioner in the following factorization form:
where D and E are the diagonal and the lower triangular matrices that consist of the diagonal and the lower triangular components of M, respectively (M = D + E +ET), w is called the relaxation parameter, and F(w) represents the two w dependent terms on the right hand side of eq. (20). Eq. (20) shows that the conditioner M˜ is not necessarily symmetric (in the SSORCG M˜ is made symmetric by pre scaling), and the factors H1 = D + wE and H2 = I + wD−1ET are directly assembled from the lower and upper triangular matrices M. With such factors, the vector qi can be obtained just by forward and back substitution in the PCG. This method may be called SORCG. Moreover, eq. (20) shows that the optimal relaxation parameter w is the value that minimizes the norm of the matrix F(w). It has been shown that w depends on the eigenvalue of the multiplication of the matrices M˜−1M (Axelsson 1984). In general, it is difficult to find the optimal value, but it may be chosen by trial and error in different cases.
The ICCG solver, based on the factorization M = LLT, defines the pre conditioner with an incomplete Choleski factorization—performing the Choleski factorization with rejection by magnitude or by position (David 1978). Manolis & Michael (1991) pointed out that the incomplete factorization based on rejection by magnitude has the disadvantage that the storage size and pattern are not known at the start and the cost in computer space may become high since some auxiliary vectors will be needed for storage of the pre conditioning matrix, but the factorization with the rejection by position (M˜ keeps the same sparsity pattern as matrix M) prescribes the additional storage beforehand and the computing times are very competitive. Following Manolis and Michael's method (1991), we rewrite ICCG based on rejection by position as follows:
where P(Mij) represents the sparsity pattern of M and cii(k) is calculated by the following equations:
The diagonal modification for the case Mii*≤0 in eq. (21) was suggested by Ajiz & Jennings (1984) to retain the stability of ICCG.
Computer implementation
As presented above, the parallelepiped and tetrahedron element schemes lead to different coefficient matrices M, but they result in the same number of nodes in a 3 D discretization. The different solvers mentioned above can be applied to eq. (10) for the 3 D solution. It will be seen that with the two element schemes one can obtain different solution accuracies and with the different solvers one can obtain different computing efficiencies, because these schemes and solvers have quite different storage size requirements and computer time consumption.
Storage size of matrix M
With the parallelepiped approximation each node has a maximum of 26 adjacent nodes around it (see Fig. 2). This means that the coefficient matrix M generally has 27 non zero components in each row. Considering the symmetric band, the computation and the compact storage may only be that required to calculate and store the lower or upper band; that is, 14 non zero elements are calculated in each row of M. As an example, if NX, NY and NZ (n = NX × NY × NZ) are the total numbers of the nodes in the x , y and z directions, and NXY (= NX × NY) is the number of the nodes in the xy plane, and the condition NXY≤NXZ≤NYZ (NXZ = NX × NZ, NYZ = NY × NZ) holds, the 14 non zero elements and their node numbers are shown in Fig. 2. These node numbers are also indicators of their position in M. So, the parallelepiped element scheme requires a storage size of at least M(n, 14) for any iterative solver and M(n, NXY + NX + 2) for the banded Choleski factorization (Jennings & Mckeown 1992). With the tetrahedron element scheme, however, the storage size may reduce to M(n, 8) for any iterative solver and M(n, NXY + 1) for the banded Choleski factorization (see Fig. 2). A great deal of computer memory (about 43 per cent for an iterative method) is saved. In general, the tetrahedron element scheme needs a storage size M(n, min(NXY, NXZ, NYZ) + 1) for the band elimination algorithm. Obviously, for 3 D modelling the storage size is still significant: a 128 MB PC can only handle a 3 D grid size not larger than 31 × 31 × 31 with the banded Choleski factorization, while with an iterative method one may have a wide range of grid sizes for different applications.
The adjacent nodes of the ith node in a 3 D discretization. The nodes with a number are the non zero components of the ith row in the coefficient matrix M when the tetrahedron element is used. The nodes with a subscript yield the non zero entries in the ith row in the coefficient matrix M when the parallelepiped element is employed.
Programming of the PCG
From algorithm (19), one finds that the only difference between the PCG and the standard CG than solving the linear equation system M˜qi = Mpi (if M˜ =I, the PCG becomes the normal CG). The crucial point is the choice of the pre conditioner M˜ according to M. Matrix M˜ should (1) achieve rapid convergence, (2) have an additional storage size as small as possible, and (3) obtain qi efficiently. In fact, the four schemes of the PCG (SCG, TRIDCG, SORCG and ICCG) referred to above have no problem with the third criterion, because these pre conditioners are diagonal, tridiagonal, lower and upper triangular matrices by which qi can be efficiently obtained by just performing the forward and back substitution. Therefore, only the convergence and the storage size of the pre conditioners need to be examined against the parallelepiped and tetrahedron element schemes for 3 D resistivity modelling.
Based on the above compact storage of M, it is not difficult to write the segments of the program to obtain qi for SCG, TRIDCG and SORCG. In fact, these three pre conditioners are implicitly contained in M, so it is unnecessary to prescribe additional storage for M˜. As an example, with the tetrahedron element the segment of the FORTRAN code for the SORCG may be written as follows according to eq. (20):
where CMJ(I,MJ) is a subroutine to find the node numbers MJ(15) for the 15 adjacent nodes of the ith node (they indicate the positions of the non zero elements of the ith row in M), w is the relaxation parameter, M(n, 8) is the compact storage of M, and B(n)=Mpi. Part 1 is the implementation of the forward substitution with H1 = D + wE, and Part 2 obtains the vector qi through back substitution with H2 = I + wD−1ET. From this code, one finds that only one vector V(n) needs to be prescribed, and the other array MJ(15) is very small in size (15 elements).
For ICCG, Ajiz & Jennings (1984) gave the code for the incomplete Choleski factorization based on rejection by magnitude. Manolis & Michael (1991) improved the scheme and gave another program for general structural analysis. In our case, we found that with the compact storage M(n, 14) (parallelepiped element) or M(n, 8) (tetrahedron element), the algorithm based on rejection by position becomes very simple; that is, using the tetrahedron element scheme, the major computation of Σk = 1i − 1lkilkj in eq. (21) is fulfilled by the following segment:
where MI(15) and MJ(15) are two integer arrays for the positions of the non zero elements of the ith and jth rows, and L(n, 8) is the incomplete Choleski factorization. From this algorithm one can see that the main additional storage is L(n, 8), which has the same size as M(n, 8). After the factorization, the same code as for SORCG can be used for ICCG by replacing M(n, 8) with L(n, 8) and setting w = 1.
Numerical experiments
In order to examine the efficiency and accuracy of the FEM schemes (parallelepiped and tetrahedron elements) and the above linear equation solvers, we conducted numerical experiments with three synthetic models, namely homogeneous, two layered and vertical contact, whose analytic solutions are easily obtained for comparison. The computer used for the experiments was a Pentium II, 128 MB PC running at 450 MHz. We recorded the storage sizes and computing times of these solvers and the accuracies of the parallelepiped and tetrahedron element schemes. One can compare these records and determine an accurate and efficient element scheme and solver for 3 D DC resistivity solutions.
Figs 3, 4 and 5 show the solutions with the tetrahedron element scheme for the three models. They depict equipotential contours in the horizontal and vertical planes, as well as the voltage versus depth profile for model 1 (see Fig. 3), apparent resistivity sounding with a dipole–dipole array for model 2 (see Fig. 4), and voltage change along the x axis for model 3 (see Fig. 5). The grid size of the FEM was 119 × 119 × 80 (Δx = Δy = Δz = 0.5 m), and the SORCG was employed for the solutions. The time required was about 45 min for each model. Fig. 6 is a comparison of the accuracy of the parallelepiped and tetrahedron element schemes. These results show that the tetrahedron element scheme is much better than the parallelepiped element scheme, not only in the storage size (savings of 43 per cent—see the previous section) but also in the accuracy (the maximum relative errors are ≤4 per cent and ≤1 per cent for element sizes of 1 m and 0.5 m respectively).
FEM result of a buried current electrode in a homogenous half space. The left diagram is the equipotential contour volume and the right diagram is the voltage versus depth profile in the borehole sketched in the left diagram. The exact analytic solution is also shown for comparison.
FEM result of a dipole–dipole survey in a two layer model. A and B are current electrodes, M and N are potential electrodes. The left diagram shows equipotential contours in the horizontal and vertical planes. The right diagram is the apparent resistivity graph for the FEM. The exact (analytic) solution is shown for comparison.
FEM result of a pole–pole survey near a vertical contact. A and M represent the current injection and potential point. The equipotential diagram and the voltage versus AM separation curve are shown.
Comparison of solution accuracy between the parallelepiped and tetrahedron element schemes. (a) 3 D grid size is 57 × 57 × 49 (element size=1 m); (b) 3 D grid size is 119×119 × 80 (element size=0.5 m).
As mentioned above, with the SORCG solver one has to predetermine the relaxation parameter w in applications (see eq. 20). In order to detect the sensitivity of w to the convergence of SORCG, we repeated the above three simulations with four relaxation parameters (w = 0.5, 1.0, 1.5, 2.0) and recorded the computer times for the same degree of convergence (10−10). We found that with both the parallelepiped and tetrahedron element schemes the optimal value of the relaxation parameter was w = 1.5 for the three models (see Fig. 7). This value is very close to Spitzer's (1995) result in finite difference experiments. These experiments show that the optimal value is relatively stable with the two element schemes and the three models. Therefore, it can be confidently applied to 3 D resistivity computations.
Effect of the relaxation parameter of SORCG for three synthetic models. Computing time is least for a value of w = 1.5. (a) Parallelepiped; (b) tetrahedron.
A further comparison of the computing efficiency of the four PCG algorithms (SCG, TRIDCG, SORCG and ICCG) was made with the two element schemes and three models. Table 1 gives the details of the comparison using the tetrahedron element scheme, and shows the differences in additional work, additional storage size, number of iterations and computer times to obtain the Green's function of one source vector. The parallelepiped element scheme has similar results to those in Table 1. From this table one can see that (1) SCG and TRIDCG do not require any additional storage, SORCG has one more vector, but ICCG needs much more computer memory [nearly double the storage size of others, because L(n, 8) is the same size as M]; and (2) TRIDCG has little improvement in convergence rate over SCG, SORCG converges much faster than SCG and TRIDCG, and ICCG has the best convergence speed of the four solvers. These results indicate that SORCG is a suitable option for large scale computations, but if the computer memory is available for additional storage L(n, 8), ICCG is the best choice of the iterative solvers. Both of them may save on computer time by at least ∼50–70 per cent.
Comparison of PCG algorithms for solving MG̃=bs.
To investigate the computing efficiency of a band elimination solver, we applied the banded Choleski factorization method to the three models. Owing to the limitation of the PC memory, only the tetrahedron element scheme and two grid sizes, n = 17 × 17 × 14 and n = 33 × 33 × 27, could be employed for the experiments. The second grid has nearly double the density of the first one. Although the two grid sizes are not large enough for a real case, their storage sizes and the computer times may help us to estimate the computing efficiency for a larger size simulation and inversion. According to the storage requirement of the tetrahedron element scheme (M[n, min{NX*NY, NX*NZ, NY*NZ} + 1]), the second grid size requires at least 128 MB of memory for the banded Choleski factorization. This is the maximum size we can handle with our PC. The results of the experiments (omitted here) showed that the two grid sizes hardly attain reasonable levels of accuracy for the three models; the first grid size took about 20 s and the second grid size spent over one hour solving eq. (10). The experiments proved that a coarse grid cannot give satisfactory solutions, and a fine grid, for example just doubling the density in the three directions of the coarse grid, will cost over 150 times as much computer time as a coarse grid. In fact, the time consumption of the factorzation method can be estimated using the relation (18), from which we have the following computer time for the Choleski elimination method for two grid sizes:
where T1 and T2 are assumed to be the computing times for the banded Choleski factorization using grid 1 and grid 2, whose dimensions in the x , y and z directions are (NX1, NY1, NZ1) and (NX2, NY2, NZ2), respectively. From this expression one can see that if grid 2 has double the size of grid 1, the computing time is at least 27 = 128 times that of grid 1. Our experiments showed that the actual computing time is greater than the estimate, because eq. (23) is just for the banded Choleski factorization. Such an increase in the computing time should be compared with the iterative solver SORCG or ICCG, even for many source vectors, otherwise the procedure still costs more computing time than the iterative solvers. Consequently, the efficiency of the elimination method is only achieved when the following inequality is true:
where TCF, TSORG andTICCG stand for the computing times required to solve eq. (10) using the Choleski factorization, SORCG and ICCG solvers, respectively. NRD is the total number of electrodes.
Conclusions
For 3 D DC resistivity modelling and inversion, trilinear interpolation within an 8 node solid parallelepiped and linear interpolations within six tetrahedral bricks within the same 8 node solid block may be employed in the FEM, but the latter is far superior to the former in both accuracy and computational efficiency. With the compact storage, the tetrahedron element scheme saves at least 43 per cent in computer memory for an iterative solver and can produce much better simulating accuracy than the parallelepiped element scheme with the same grid size.
The pre conditioned conjugate gradient method is a suitable option as an iterative solver for 3 D FEM DC resistivity solutions, specifically the versions SORCG and ICCG. Our numerical experiments show that both have a much better convergence rate than SCG and TRICG. Applying the tetrahedron element scheme and ICCG, one can obtain the best convergence speed of all the iterative solvers. With the compact storage M(n, 8) in the tetrahedron element scheme, the main programming segments for SORCG and ICCG become very simple and it is shown that SORCG has only one additional vector storage V(n), but that ICCG needs nearly double the storage size of the others for the incomplete Choleski factorization L(n, 8).
The banded Choleski factorization algorithm may be an alternative to a standard elimination solver for a 3 D resistivity inversion that involves hundreds of electrodes for the measurement, but it must be noted that the minimum storage size needs a 2 D array M(n, min(NXY, NXZ, NYZ) + 1), and for a fine grid size which is half of a coarse grid size in each direction, the computing time becomes over 27 = 128 times that required for the coarse grid size. Such an increase in the computing time should be compared with SORCG or ICCG so as to achieve the maximum computing efficiency.
Acknowledgments
The authors are grateful to the Department of Primary Industries and Resources, South Australia and to the CSIRO Center for Groundwater Studies for funding the project. They also thank the anonymous reviewers for their comments.
References
