Solving the fixed gravimetric boundary value problem by the finite element method using mapped infinite elements.

The numerical approach for solving the fixed gravimetric boundary value problem (FGBVP) based on the finite element method (FEM) with mapped infinite elements is developed and implemented. In this approach, the 3D semi-infinite domain outside the Earth is bounded by the triangular discretization of the whole Earth’s surface and extends to infinity. Then the FGBVP consists of the Laplace equation for unknown disturbing potential which holds in the domain, the oblique derivative boundary condition (BC) given directly at computational nodes on the Earth’s surface, and regularity of the disturbing potential at infinity. In this way, it differs from previous FEM approaches, since the numerical solution is not fixed by the Dirichlet BC on some part of the boundary of the computational domain. As a numerical method, the FEM with finite and mapped infinite triangular prisms has been derived and implemented. In experiments, at first, a convergence of the proposed numerical scheme to the exact solution is tested. Afterwards, a numerical study is focused on a reconstruction of the harmonic function (EGM2008) above the Earth’s topography. Here, a special discretization of the Earth’s surface which is able to fulfil the conditions that arise from correct geometrical properties of finite elements, and it is suitable for parallel computing is implemented. The obtained solutions at nodes on the Earth’s surface as well as nodes that lie approximately at the altitude of the GOCE satellite mission have been tested.


Introduction
In mathematical modelling, there are some problems which are obviously unbounded like the object of our study -gravity of the Earth. However, also many other engineering problems, where, e.g., one concentrates only on a small subdomain of an extremely large domain may be also termed as infinite domain problems, since in these cases, the total area or volume can be for practical purposes taken as infinite. These mathematical models then occur in various engineering applications, e.g. in elasticity [21,29], potential problems [10], fluid flow [17,39,42], etc. In the finite element analysis these problems are solved in various ways, all of which have advantages and disadvantages. The simplest way is to truncate the domain by adding some artificial boundary with the prescribed artificial BC at some large but finite distance, see e.g. [12,14]. In our approach, we will use the so-called mapped infinite elements (MIE), which were originally pioneered by Bettess in [2]. We will show that such an approach provides an efficient and effective alternative to other wellknown numerical approaches.
Gravity field modelling is usually performed as solving the so-called geodetic boundary value problems (BVPs). In our study, we will focus on one of them, namely the fixed gravimetric boundary value problem (FGBVP) that repre-sents an exterior oblique derivative BVP for the Laplace equation, cf. Koch and Pope [20], Freeden and Kersten [15], Bjerhammar and Svensson [4], Holota [16]. As the oblique derivative boundary conditions (BC) given on the discretized Earth's surface will be taken surface gravity disturbances. The detailed overview of various procedures for solving the oblique derivative BVP can be found, e.g., in Minarechová et al. [34].
With an expansion of high-performance computing (HPC) facilities in the last decades, approaches based on numerical methods have become a powerful and efficient tool finding applications also in gravity field modelling. The pioneering studies on numerical methods applied to gravity field modelling were based on the finite element method (FEM), cf. Meissl [32] or Shaofeng and Dingbo [38]. Later, the finite difference method (FDM) was applied by Keller [18], and the indirect boundary element method (BEM) approach was developed by Klees [19] and Lehmann and Klees [22]. Afterwards, the direct BEM approach was introduced by Cunderlík et al. [7] andČunderlík and Mikula [8]. The oblique derivative BVP treated by BEM was discussed iň Cunderlík et al. [9]. At that time, also new studies on FEM were published, cf. Fašková et al. [12,14] and Mráz et al. [35]. The FEM with MIE for BVP with Neumann BC has been introduced by Šprlák et al. [41]. The first application of the finite volume method (FVM) to gravity field modelling was introduced by Fašková [13] and its parallel implementation by Minarechová et al. [33]. The FVM applied to the oblique derivative BVP has been studied in Macák et al. [24,26,27]. Lately, Medl'a et al. [31] presented the FVM for solving the oblique derivative BVP on 3D unstructured meshes above the real Earth's topography and Droniou et al. [11] analysis of FVM for elliptic equations with oblique derivatives. Recently, Yin and Sneeuw in [43] published a new approach to the gravitational field modelling by using CFD (Computational fluid dynamics) techniques based on FEM, and Macák et al. in [28] and Minarechová et al. in [34] applied the FEM for local solutions of the oblique derivative BVP.
In previous numerical approaches based on FEM or FVM, see for example [12-14, 24, 26-28, 31, 33, 34, 41], a condition of the regularity at infinity has been abandoned by introducing an artificial upper boundary approximately at altitudes of the GOCE satellite orbits. Here the FEM or FVM numerical solutions have been fixed by the prescribed Dirichlet BC in terms of the disturbing potential. In this paper we avoid this drawback and the presented FEM approach is developed to solve the original FGBVP including the condition of the regularity at infinity. This is a main contribution of this study in comparison with the previous ones.
The paper is organized as follows. In Section 2, we formulate the FGBVP with the oblique derivative BC. In Section 3, we derive a numerical scheme of the FEM with MIE for the FGBVP. Then numerical experiments are presented in Sec-tion 4 that aim to test accuracy of the developed approach.
Here we also describe a special triangular discretization of the Earth's surface that is suitable for FEM. The paper ends with Conclusion and summary.

Formulation of the FGBVP with the oblique derivative BC
Let us consider the FGBVP, cf. [4,16,20]: where the computational domain , see Fig. 1, is bounded by the boundary , i.e., the Earth's surface, from the bottom and extends to infinity. T (x) stands for the disturbing potential defined as a difference between the real and normal gravity potential at any point x = (x, y, z), δg(x) is the gravity disturbance, and the vector s(x) = −∇U (x)/|∇U (x)| is the unit vector normal to the equipotential surface of the normal potential U (x) at any point x.

Solution to the FGBVP with the oblique derivative BC by the FEM with MIE
In our approach, we follow the fundamental principles of FEM presented by Reddy in [37] with the theory of MIE Fig. 1 The semi-infinite computational domain for global gravity field modelling. It is bounded from the bottom by the boundary discretizing the Earth's surface, and extends to infinity in the vertical direction published by Bettes in [2,3], Marques and Owen in [29] or Zienkiewicz et al. in [44,45].

Discretization of the computational domain
The FEM is a numerical method that assumes discretization of the whole computational domain by a union of elements e , e = 1, ..., N , where N denotes the number of elements in the domain . In our FEM approach, we divide the semiinfinite computational domain into two parts, see Fig. 2, where the lower one, denoted by F E , is meshed with finite elements and the upper one, M I E , is meshed with one layer of infinite elements.
For our purpose -to fit irregular boundary , we have chosen triangular prisms, i.e., finite pentahedral elements with six nodes and five faces, and corresponding mapped infinite pentahedral elements with nine nodes and five faces, see Fig. 3. In this way, we divide the computational domain into n 1 , n 2 , n 3 elements in the latitudinal, longitudinal and vertical direction, respectively, and to specify the position of an element e we use indexes k, l, m, where k = 1, ..., n 1 , l = 1, ..., n 2 and m = 1, ..., n 3 .

Derivation of the weak formulation on the element
Let us consider an arbitrary element e from our finite element discretization with indexes k = 1, ..., n 1 , l = 1, ..., n 2 and m = 2, ..., n 3 . We multiply the differential equation (1) by a weight function w and using Green's identity (we omit (x) to simplify the notation in the following equations) we obtain the weak formulation (WF) of (1) over an arbitrary above defined element e e ∇T · ∇w dxdydz = ∂ e ∇T · n w dσ, where n denotes the unit normal to ∂ e . Since on the bottom boundary the oblique derivative BC (2) is prescribed, for the row of elements that lie on this boundary (see Fig. 4, elements depicted by yellow), i.e., k = 1, ..., n 1 , l = 1, ..., n 2 and m = 1, we modify (4) in such a way that has been also presented, for instance, in [28] or [34]. We split the oblique vector s into one normal and two tangential components where n is the normal vector and t 1 , t 2 are tangent vectors to e ⊂ ∂ e ⊂ R 3 , where e denotes the bottom boundary of an element e .
where we assume that c 1 = 0, since in practical experiments we always use a 'nonzero' horizontal resolution of the grid points that discretized the real Earth's surface. In such manner, s is never perpendicular to n, so the assumption c 1 = 0 is always fulfilled. Now, we insert (7) to (4) to get After some rearrangement, we have ∇T · n w dσ.
Thus we have obtained the weak formulation (4) or (9) of the BVP (1) -(3) on every element e of our finite ele-ment discretization. The study of weak solution of the oblique derivative BVP is included in the book by Lieberman [23].

Solution by the Finite Element Method
For a finite pentahedral element e with six nodes, see Fig. 3 a), we can write i. e., we take an approximation of the unknown value T as T e , a linear combination of basis functions ψ j with coefficients T e j , j = 1, ..., 6. Then we substitute it into the weak formulation (4), namely for elements e with indexes k = 1, ..., n 1 , l = 1, ..., n 2 and m = 2, ..., n 3 − 1 (see Fig. 5, elements depicted by green), and consider ψ i for weight function w. We obtain the i th equation in the form where q n = ∇T · n denotes the projection of the vector ∇T along the unit normal n.
For the row of elements e given by indexes k = 1, ..., n 1 , l = 1, ..., n 2 and m = 1 (see Fig. 5, elements depicted by yellow), we follow the same way and after inserting (10) into (9) and considering w = ψ i , we obtain the i th equation in the form where index j = 1, ..., 3 refers to nodes of the element e that lie on the bottom boundary of the computational domain . As we can see in Eq. (12), there are 2 tangent vectors corresponding to the nodes lying on the bottom boundary of the element belonging to .
Finally, for the mapped infinite pentahedral element e with nine nodes, see Fig. 3

b), we can write
We substitute (13) for elements e with indexes k = 1, ..., n 1 , l = 1, ..., n 2 and m = n 3 (see Fig. 5, elements depicted by grey) into (4), consider ψ i for weight function w and we obtain the i th equation in the form where q n = ∇T · n again denotes the projection of the vector ∇T along the unit normal n.

Shape and mapped functions
The basis function ψ i is the piecewise quadratic function and it is uniquely determined by choosing value 1 at N i and 0 at every N j , i = j, cf. [5]. In our approach, we will work with isoparametric coordinates ξ, η, ζ in local coordinate system, so we rewrite (10) and (13) in the form Table 1 The shape functions ψ i (ξ, η, ζ ) for the finite pentahedral element with six nodes N i defined by isoparametric coordinates ξ , η and ζ where the transformation between local coordinates ξ, η, ζ and global coordinates x, y, z is given by Then the shape functions ψ i (ξ, η, ζ ) for finite pentahedral element with six nodes (see Fig.3 a)) in local coordinate system that is used to mesh F E are defined in Table 1.
Since in FGBVP (1) - (3), the computational domain tends to infinity only in the vertical direction, we start by explaining MIE on one-dimensional problem and then we will apply this idea on our FGBVP.
Let us have the one-dimensional element e which extends from node N 1 with coordinate x 1 through N 2 with coordinate x 2 to the point N 3 infinity, see Fig 6. Then let us have the so-called Pole that is the point O with coordinate x 0 and lies on the line that passes through points N 1 , N 2 and N 3 . The relationship between coordinates x 0 , x 1 and x 2 is given by It is obvious that the position of the Pole has an impact on the results and has to be chosen very carefully, since an incorrect Pole position can lead to a distortion of the map- Fig. 6 Brief sketch of one-dimensional MIE e ping, cf. [2,3,44,45]. However, due to geometry of our computational domain , see Fig. 1, it is natural to have O = [0, 0, 0], so we will not discuss the problem of incorrect Pole position in detail. Then this element is mapped onto the parent element defined in the local coordinate system in the range −1 < ζ < 1 using formula where It is obvious from (19)- (21) that ζ = −1, 0, 1 correspond to the global positions of x 1 , x 2 , ∞, respectively. Then for infinite pentahedral element with nine nodes N i , we obtain the mapping functions M i (ξ, η, ζ ) by multiplying (21) with the shape functions ψ i (ξ, η, ζ ) where directions ξ and η are finite. These mapping functions M i (ξ, η, ζ ) can be seen in Table 2.
We remain that for more details about MIE, see e.g. [2,3] or [44,45]. Now we can write (11) and (12) in a compact matrix form where K e = [K i j ] stands for an element stiffness matrix, T e = (T 1 , ..., T 6 ) is a column vector of unknowns and Q e denotes the right-hand side vector.
To evaluate element matrices and vectors we proceed the following way. We choose one basis function ψ i per vertex N e i and we differentiate the basis functions with respect to a position of each node. To calculate two integrals over a boundary in Eq. (11) which include a tangential derivative, Table 2 The mapping functions M i (ξ, η, ζ ) for the infinite pentahedral element with nine nodes N i which extends to infinity in the vertical direction. Values ξ , η and ζ denote isoparametric coordinates of nodes we approximate derivatives in tangential direction like in the FDM, i.e., using values of basis functions at nodes N i of element e we have where d denotes the distance between two neighbouring nodes. The similar idea of approximating derivatives in tangential direction like in the FDM, however for hexahedral elements, has been presented in [34].

Assembly of element equations
Finally, we assemble all element equations using two principles, see Reddy [37]: (i) continuity of primary variables at the interelement nodes. It means that nodal values T e j and T e+1 j of two adjacent elements e and e+1 at the connecting nodes have to be the same.
(ii) balance of secondary variables in a weighted-integral sense, see Reddy [37].
In this manner, we obtain the global linear system of equations with a column vector of unknown global nodal values T where the matrix K is sparse, since most of its entries are zero, and positive definite, and Q is the column vector whose entries are also almost zero except that for nodes with the prescribed oblique derivative BC (2).

Numerical experiments
To test the proposed numerical scheme, we present several numerical experiments. In them, results of the proposed method are compared to either exact, in testing experiments, or EGM2008, in experiments with gravity data, solutions. In testing experiments, the numerical scheme is qualified according to the value of the experimental order of convergence (EOC), see e.g. [28], that can be determined by comparing numerical solutions and exact solutions on subsequently refined grids. In experiments with the reconstruction of EGM2008 data, the obtained numerical solution is compared with the disturbing potential generated from EGM2008 directly and is qualified according to the statistics of residuals. To solve the nonsymmetric linear system, BiCGSTAB(l) linear solver with l = 8, cf. [40], has been implemented.
To save the memory, only nonzero coefficients were stored. These large-scale parallel computations were performed on 4 nodes of our cluster with 1.0 TB of distributed memory (each node consists of four 8-core CPUs with 256 GB RAM). Thanks to the NUMA (Non-Uniform Memory Access) architecture of each node, we have implemented an MPI and OpenMP parallelization, see [25].

Discretization of the bottom boundary
As we have mentioned in the beginning, the FEM requires elements to meet some geometrical specifications, as well as detailed global gravity modelling requires parallel computing. In this section, we describe the generation of the bottom boundary mesh, which is suitable for the presented numerical approach. The final triangle mesh is geometrically a standard octahedron sphere, usually created by subdividing the faces of an octahedron. However, we use a numbering of vertices and faces, which is suitable for a parallelisation. The mesh generation is based on mapping of a rectangle to regular octahedron, see Fig. 7. The process of mesh generation is then performed in several steps described below.
where N j,1 , N j,2 denote x and y coordinate, respectively. Then this rectangle is divided into 8 auxiliary triangles T 1 , . . . , T 8 , see Fig. 7 a). 2. Mapping to octahedron. Then we map the rectangle onto the regular octahedron by a piecewise affine mapping φ : R 3 → R 3 . It means that each triangle T i with corresponding nodes N j is mapped onto the face T i of the octahedron, see Fig. 7, by the affine mapping φ i : R 3 → R 3 defined as where G i , G i denote the centroids of the triangles T i , T i , respectively, and M i is a 3 × 3 matrix. Therefore, the mapping φ is defined as φ(X ) := φ i (X ) for X ∈ T i . The nodes N j on rectangle are mapped to the nodes N j := φ(N j ) on the octahedron. The construction of the matrix M i is done as follows. Let X 1 X 2 X 3 denote a domain triangle with centroid G i and X 1 X 2 X 3 a range triangle with centroid G i . ThenX k := X k − G i andX k := X k − G i denote the coordinates of Final splitting and numbering of nodes and mesh triangles is depicted on rectangle. Duplicate nodes are highlighted by green vertices w.r.t. centroid, see Fig. 8. If we plug the vertices X k and X k = φ i (X k ) into the equation (26), we obtain a linear system M iXk =X k , k = 1, 2, 3 with 9 equations and 9 unknown elements of the matrix M i . Unfortunately, the equations are not linearly independent, becauseX 1 + X 2 +X 3 = 0 andX 1 +X 2 +X 3 = 0. To obtain an independent system of equations, we can include, e.g., normal vectorsX 1 ×X 2 ,X 1 ×X 2 instead ofX 3 ,X 3 , respectively. Therefore, we have to solve the system If we define 3 × 3 matrices we see that the matrix M i is given by M i =X X −1 . 3. Gluing. The mapping φ glues the triangle edges on the boundary of the rectangle. There are 5 gluings together of pairs of triangle edges, 4 of them are indicated by grey arrows in Fig. 7, left. The last one glues the left and right side of the rectangle. As a result, corresponding nodes N j need to be identified, see Fig. 9. 4. Normalisation to unit sphere. The nodes N j on the octahedron are mapped to the unit sphere simply by nor-malisationN j := N j / N j . 5. Scaling to the WGS84 ellipsoid. Finally, we scale the unit sphere to the ellipsoid to obtain final nodes V j := diag(a, a, b)N j , where diag stands for diagonal matrix, and a is the semi-major axis and b is the polar semi-minor axis of the WGS84 ellipsoid. 6. Splitting quads. In the last step, we split each quad diagonally to create two triangles. The direction of splitting as well as numbering of the triangles and nodes is shown in Fig. 9.
In numerical computations we use meshes with n = 2 l , where l is the level of refinement of the mesh. Examples of two coarse meshes (level l = 2, 4) are visualised in Fig. 10. Mesh statistics for levels l = 1, . . . , 10 is presented in Table 3.

Testing numerical experiments
As usual in numerical mathematics, at first we test the stability and behaviour of the numerical scheme derived in Section 3 by investigating its EOC. In these testing numerical experiments we have chosen the well-known example, gravitational potential generated by a homogeneous sphere, Table 3 Mesh statistics for meshes with the refinement level l. The number of divisions of rectangle in y direction is calculated as n = 2 l , the number of quads as 4n 2 , the number of triangles as 8n 2 , the number of vertices (nodes) on the rectangle as (4n + 1)(n + 1) and the number of duplicate vertices on ellipsoid as 5n − 1

Testing experiment 1
In the first testing experiment, the height of the finite domain F E has been 500 [km], so the finite/infinite element interface has been at height 6 871 [km]. Then the center of the infinite elements (which corresponds to the point x 2 in Fig. 6) has been at height 13 742 [km]. We have started with the mesh made up of 64 × 16 × 4 nodes, see Table 4, and then we have performed four successive refinements. Results have been compared with the exact solution in the form G M/R.
The statistics of residuals [m 2 s −2 ] on the bottom boundary computed between the obtained numerical solution and the exact solution can be found in Table 4. One can observe that by refining the mesh, the obtained numerical solution converges to the exact one and that the proposed method is second order accurate, see the EOC column.

Testing experiment 2
For the second testing experiment we have chosen the domain and mesh consisting of 256 × 64 × 16 nodes from Testing experiment 1, see Section 4.2.1. Then we have fixed the size of elements while redoubling the radius of the finite element domain. It means that with a doubling of the height of the finite element domain, the number of elements in the vertical direction, see in Table 5, redoubled as well to remain the size of elements. Again, we have compared obtained results with the exact solution in the form G M/R.
The statistics of residuals on the bottom boundary can be seen in Table 5. We can observe that in case of homogeneous sphere, the height of the domain doesn't have impact on the solution on the bottom boundary (Table 6).

Global gravity field modelling with the EGM2008 data as the reconstruction of the harmonic function
The bottom boundary has been the discretized Earth's surface created with Becker et al. [1] data. The horizontal resolution has been set in accordance with Table 3 presented in Section 4.1, namely 512×128, see Table 7, and 1024×256, see Table 8. The input surface gravity disturbances as BC (2) have been generated from Earth Gravitational Model 2008 (EGM2008) [36].  Tables 7 and 8). Then the obtained numerical solution on the bottom boundary and at altitude H=200 [km] has been compared to the disturbing potential generated from EGM2008 directly. Statistics of residuals for successive refinements can be found in Tables 7 and 8. In them, we can observe an initial rapid decrease of the stan-   . We can also observe an improvement of mean value of residuals when refining the mesh in the vertical direction. Also, it is not surprising that the standard deviation of residuals is better at H=200 [km] than at , since the disturbing potential solution smooths with increasing heights. When we look closely at the  Table 6) and then we have doubled it. The obtained residuals are visualised in Fig. 11 and their corresponding statistics is presented in Table 6. We can observe an improvement of all statistical characteristics when refining the computational grid, which proves the validity of the presented FEM approach to the global gravity field modelling.

Conclusion and summary
In the presented paper, a novel numerical approach for solving the original FGBVP including the condition of the regularity at infinity was developed. Here, the numerical scheme based on the FEM with finite and mapped infinite elements was derived and implemented. In numerical experiments, at first, a convergence of the proposed numerical scheme to the exact solution in different situations was studied. It was shown that in the case of potential generated by a homogeneous sphere, the proposed numerical approach is second order accurate. Then a numerical study focused on a reconstruction of the harmonic function (EGM2008) above the discretized Earth's topography was done. The obtained numerical solutions were tested at nodes on the Earth's surface as well as nodes lying approximately at the altitude of the GOCE satellite mission. Gained results showed a practical contribution of the presented approach to global gravity field modelling.