Add Journal to My Library
Geophysical Journal International
, Volume 212 (3) – Mar 1, 2018

35 pages

/lp/ou_press/finite-element-solution-to-multidimensional-multisource-vNc37b80J0

- Publisher
- The Royal Astronomical Society
- Copyright
- © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.
- ISSN
- 0956-540X
- eISSN
- 1365-246X
- D.O.I.
- 10.1093/gji/ggx530
- Publisher site
- See Article on Publisher Site

Summary We propose an approach to solving multisource induction logging problems in multidimensional media. According to the type of induction logging tools, the measurements are performed in the frequency range of 10 kHz to 14 MHz, transmitter–receiver offsets vary in the range of 0.5–8 m or more, and the trajectory length is up to 1 km. For calculating the total field, the primary-secondary field approach is used. The secondary field is calculated with the use of the finite-element method (FEM), irregular non-conforming meshes with local refinements and a direct solver. The approach to constructing basis functions with the continuous tangential components (from Hcurl(Ω)) on the non-conforming meshes from the standard shape vector functions is developed. On the basis of this method, the algorithm of generating global matrices and a vector of the finite-element equation system is proposed. We also propose the method of grouping the logging tool positions, which makes it possible to significantly increase the computational effectiveness. This is achieved due to the compromise between the possibility of using the 1-D background medium, which is very similar to the investigated multidimensional medium for a small group, and the decrease in the number of the finite-element matrix factorizations with the increasing number of tool positions in one group. For calculating the primary field, we propose the method based on the use of FEM. This method is highly effective when the 1-D field is required to be calculated at a great number of points. The use of this method significantly increases the effectiveness of the primary-secondary field approach. The proposed approach makes it possible to perform modelling both in the 2.5-D case (i.e. without taking into account a borehole and/or invasion zone effect) and the 3-D case (i.e. for models with a borehole and invasion zone). The accuracy of numerical results obtained with the use of the proposed approach is compared with the one obtained by other codes for 1-D and 3-D anisotropic models. The results of this comparison lend support to the validity of our code. We also present the numerical results proving greater effectiveness of the finite-element approach proposed for calculating the 1-D field in comparison with the known codes implementing the semi-analytical methods for the case in which the field is calculated at a large number of points. Additionally, we present the numerical results which confirm the accuracy advantages of the automatic choice of a background medium for calculating the 1-D field as well as the results of 2.5-D modelling for a geoelectrical model with anisotropic layers, a fault and long tool-movement trajectory with the varying dip angle. Electrical properties, Downhole methods, Numerical approximations and analysis, Numerical modelling 1 INTRODUCTION For processing the electromagnetic soundings data, the multidimensional modelling is increasingly used (e.g. Oldenburg et al.2013; Kamm & Pedersen 2014; Kordy et al.2016b). The inversions based on the 2-D or 3-D modelling make it possible to recover the structure of the investigated medium with greater adequacy (e.g. Persova et al.2013; Schwarzbach & Haber 2013; Haber & Schwarzbach 2014; Usui 2015; Persova et al.2016, etc.). However, these inversions require high computational costs. Therefore in practice, until recently, quite simple and low-cost approaches, though not always adequate, have often been used. They are based on calculating the electromagnetic fields in the 1-D medium or other approaches, which make use of simplified mathematical models of the electromagnetic field (among others Macnae 2015; Hunkeler et al.2016). In some cases, these approaches allow us to obtain acceptable results (e.g. Hennessy & Macnae 2010; Hunkeler et al.2016). Nevertheless, the researchers often note that the simplified and 1-D approaches not rarely give the distorted results, especially near local inhomogeneities (e.g. Ley-Cooper et al.2015; Persova et al.2015b, Ullmann et al.2016). In induction logging, this problem is important, e.g. when investigating the horizontal wells, and if the geometry and electrophysical properties of the borehole environment are locally and drastically changed or it is necessary to correctly define the intersection of the fault by the borehole. In these cases, for adequately recovering the medium, it is necessary to use the multidimensional approaches based on accurate (without significant simplifications) mathematical models of the electromagnetic field. Therefore, the increasingly important problem is the development of such methods of electromagnetic fields calculation that, on the one hand, would make it possible to model the field in a complex medium with the required adequacy (with the use of mathematical models without simplifications) and, on the other hand, would be computationally low-cost and convenient to use. For multidimensional electromagnetic fields modelling in complex 3-D media, among the mesh-based numerical methods the following are more frequently used: the finite-difference (FD) method (Grayver et al.2013; Jaysaval et al.2015; Shantsev & Maaø 2015, etc.), the finite volume method (Haber & Ruthotto 2014; Haber & Schwarzbach 2014; Jahandari & Farquharson 2014; Jahandari & Farquharson 2015, etc.) and the finite-element (FE) method (Farquharson & Miensopust 2011; Ansari & Farquharson 2014; Börner et al.2015; Persova et al.2015a; Kordy et al.2016a, etc.). Besides, for modelling the multidimensional electromagnetic fields including the solution of 3-D induction logging problems, the integral equation method is applied quite successfully (Fang et al.2003; Nie et al.2013; Yang et al.2017, etc.). Its computational effectiveness is especially high if the volume of 3-D inhomogeneities is small. When solving the multidimensional electromagnetic problems with the use of mesh-based numerical methods, a significant lowering of computational costs can be achieved by means of the following approaches. First, this is the primary-secondary field approach which uses the separate calculation of the horizontally layered medium field (that is called the primary field). Among the first papers, which considered this approach for FD and FE methods, one should mention the papers (Newman & Alumbaugh 1995; Soloveichik et al.1997; Soloveichik et al.1998; Badea et al.2001). Nowadays, this approach continues to be actively applied (e.g. Mukherjee & Everett 2011; Zaslavsky et al.2011; Da Silva et al.2012; Grayver et al.2013). The more similar the investigated medium to the horizontally layered one, which can be taken different for different logging tool positions, the higher the effectiveness of the primary-secondary field approach. In this case, for achieving the required accuracy of the summary (total) field calculation, when the anomalous (secondary) field is calculated, the number of nodes in the 3-D mesh can be decreased by several times, and the computational time of one 3-D problem solution can be decreased by an order or more. Second, a significant lowering of the computational costs can be achieved if we use the unstructured meshes (e.g. Lelièvre & Farquharson 2013; Wang et al.2013; Um et al.2015) including the non-conforming ones with local refinements (e.g. Grayver & Bürg 2014; Persova et al.2014b; Grayver 2015). In this case, in comparison with the regular meshes, the number of elements in the mesh can be decreased by several times (without the loss of the accuracy of the 3-D problem solution), and this allows the additional reduction of the computational costs by an order. Third, if the number of the logging tool positions is big, a significant lowering of the computational costs can be achieved if we use the direct solvers for the finite-difference or FE equations system (Da Silva et al.2012; Oldenburg et al.2013; Chung et al.2014; Haber & Schwarzbach 2014; Puzyrev et al.2016) (e.g. PARDISO solver (Schenk & G̈artner 2004)). In this case, the matrix factorization, which is the most computational cost procedure in calculating a 3-D field, is performed only once for several logging tool positions simultaneously, and the computational time of 3-D problems for all logging tool positions joined in one relatively small group is close to the computational time of the 3-D problem for one tool position. However, it is necessary to take into account that the computational time for one tool position when the direct solver is used can be much higher in comparison with the time when the iteration solver is used. In this paper, we consider the approach to calculating the 3-D field in the induction logging problems, in which all of the three above methods are used for lowering the computational costs. Additionally, we propose: a new method of FE calculation of the primary field in a horizontally layered medium, which is effective when it is necessary to calculate the primary field at many points simultaneously (that is necessary for calculating the secondary field); the method of 3-D electromagnetic field FE calculation on irregular non-conforming meshes with the construction of conforming basis functions from the non-conforming ones and corresponding procedure of assembling matrices of FEs; the method of grouping the logging tool positions (in order to use the direct solver when calculating the 3-D field for all positions in one group) with the automatic choice of the 1-D background medium for calculating the primary fields and automatic generation of the irregular 3-D mesh with hanging nodes for one group of logging tool positions. 2 MATHEMATICAL MODELS 2.1 The model describing the electromagnetic field in a 3-D anisotropic medium We represent a total electric field in the form $${{{\bf \vec{E}}}^t} = {{{\bf \vec{E}}}^p} + {{{\bf \vec{E}}}^s}$$ where $${{{\bf \vec{E}}}^p}$$ is the primary electric field induced by a source in a horizontally layered (1-D background) medium, and $${{{\bf \vec{E}}}^s}$$ is the 3-D secondary electric field determined by the impact of 3-D inhomogeneities such as layers pinchout, faults, non-parallelism of layers, etc. We can seek $${{{\bf \vec{E}}}^s}$$ by solving the following vector differential equation: \begin{equation} \frac{1}{{{\mu _0}}}\nabla \times \nabla \times {{{\bf \vec{E}}}^s} + i\omega \sigma {{{\bf \vec{E}}}^s} - {\omega ^2}\varepsilon {{{\bf \vec{E}}}^s} = - i\omega \left( {\sigma - {\sigma ^p}} \right){{{\bf \vec{E}}}^p} + {\omega ^2}\varepsilon {{{\bf \vec{E}}}^p}, \end{equation} (1)where μ0 denotes the vacuum magnetic permeability, σ denotes the 3 × 3 conductivity tensor in a 3-D medium, σp denotes the 3 × 3 conductivity tensor in a horizontally layered (1-D background) medium (for which $${{{\bf \vec{E}}}^p}$$ is calculated), ε(x, y, z) denotes the dielectric permeability of the 3-D medium, ω = 2πν denotes the angular frequency (ν is the electromagnetic field frequency), and $$i = \sqrt { - 1} $$. The components σij(x, y, z) of the tensor σ and components $$\sigma _{ii}^p( z )$$ of the tensor σp are the piecewise constant functions of coordinates. The components σii are not equal to $$\sigma _{ii}^p$$ only in the subdomains where 3-D inhomogeneities are located. When i ≠ j, σij are not equal to zero only in the subdomains where 3-D anisotropic inhomogeneities are located, and the anisotropic axes of these inhomogeneities do not coincide with the coordinate axes. Dielectric permeability ε can also be defined as a tensor. If thin layers in the anisotropic medium are not parallel to Cartesian planes (i.e. anisotropic axes do not coincide with coordinate axes), the conductivity tensor can be obtained by the relation \begin{equation} \sigma = {{{\bf M}}^{{\bf T}}}\tilde{\sigma }{{\bf M}}, \end{equation} (2)where M is the rotation matrix, $$\tilde{\sigma }$$ is the diagonal tensor with the components corresponding to the lateral and vertical conductivities (in the local coordinate system related with anisotropic axes), the superscript ‘T’ denotes transposition. In a similar way, the tensor of dielectric permeability can be defined in the anisotropic medium. Eq. (1) defines the secondary field for the case when the primary field is calculated without the account of displacement currents. In general, in the induction logging problems, when electromagnetic field is calculated for the frequency ν < 200 kHz, the displacement currents effects are minor, and the terms with ε can be excluded from eq. (1) (i.e. ε is assumed to be zero). If the frequency ν > 200 kHz, the displacement currents effect on the calculated electromagnetic field can be significant. In this case, the displacement currents are to be taken into account in calculating the secondary field. In induction logging measurements, the transmitters (induction coils) with different orientation can be used. Most often the coils are coaxial with the borehole. But sometimes the tools with several (two or three) orthogonal coils (the axis of one of these coils coincides with the tool axis) are considered. For computing $${{{\bf \vec{E}}}^p}$$ induced by an arbitrary oriented coil in the horizontally layered background medium, it is required to compute two fields for two coil orientations: in the first case the coil axis is perpendicular to the layers, in the second one the coil axis is parallel to the layers. In general, for calculating $${{{\bf \vec{E}}}^p}$$, analytical or semi-analytical methods can be used. The first and the simplest is the approach where the field of homogeneous space is taken as $${{{\bf \vec{E}}}^p}$$. This field can be calculated by the analytical method. The disadvantage of this approach is the following. Even if the conductivity of the homogeneous space is taken equal to the conductivity of the layer in which the transmitter is located, it is impossible to reach a comparable effectiveness that can be achieved if we take the horizontally layered medium field as $${{{\bf \vec{E}}}^p}$$. This is because, in the first case, the layers (as 3-D inhomogeneities) are the sources of the secondary field. Therefore, for reaching the required accuracy, it is necessary to use sufficiently fine meshes (that will lead to much higher computational costs). Until recently, semi-analytical methods have been developed quite intensively (Key 2009). They are very effective if there is a need to compute the field at a relatively moderate number of points (corresponding to receiver locations). But if it is required to compute the field values at a great number of points, the computational efficiency of these methods is not so high. For implementing the method of calculating the 3-D electromagnetic field when the horizontally layered medium is taken as a background medium, it is required to compute the primary field $${{{\bf \vec{E}}}^p}$$ at a great number of points (cells) in subdomains (of the computational domain) where 3-D inhomogeneities are located. Even if the optimized irregular non-conforming meshes are used, the primary field needs to be calculated at tens of thousands of points. Consequently, in this case the approaches which allow us to compute the horizontally layered medium field in all the investigated volume simultaneously are more advantageous (Persova et al.2011; Mogilatov et al.2016). Therefore, for calculating the primary field $${{{\bf \vec{E}}}^p}$$, we propose another approach based on the use of the FE method. 2.2 Calculating $${{{\bf \vec{E}}}^p}$$ in a horizontally layered background medium We represent $${{{\bf \vec{E}}}^p}$$ as a sum of solutions of several axisymmetrical problems. When describing mathematical models in the cylindrical coordinate system, we use the diagonal 2 × 2 tensor σp, rz(z) with the components $$\sigma _r^{p,rz} = \sigma _{11}^p = \sigma _{22}^p$$ and $$\sigma _z^{p,rz} = \sigma _{33}^p$$ (where $$\sigma _{ii}^p$$ are the diagonal components of the tensor σp in eq. 1). At first, we consider the mathematical model for computing the electrical field $${{{\bf \vec{E}}}^{p,1}}$$ which is induced by a ‘horizontal’ coil. This coil lies in plane which is parallel to the layers of a horizontally layered background medium (and, obviously, the coil axis is perpendicular to the layers borders). $${{{\bf \vec{E}}}^{p,1}}$$ induced by this coil is axisymmetrical and has only one non-zero component $$E_\varphi ^{}( {r,z} )$$ in cylindrical coordinates {r, φ, z}. This component ($$E_\varphi ^{}( {r,z} )$$) can be found by solving the equation \begin{equation} \! - \frac{1}{{{\mu _0}}}\Delta E_\varphi ^{} + \frac{1}{{{\mu _0}{r^2}}}E_\varphi ^{} + i\omega \sigma _r^{p,rz}E_\varphi ^{} = - i\omega I{\delta ^{{J_\varphi }}}, \end{equation} (3)where I is the value of the coil current, and $${\delta ^{{J_\varphi }}}$$ is the delta-function concentrated on the circle r = rδ, z = zδ corresponding to the coil contour ($$I{\delta ^{{J_\varphi }}}$$ is the source current density Jφ). The components of the vectors $${{{\bf \vec{E}}}^{p,1}}$$ and $${{{\bf \vec{B}}}^{p,1}}$$ can be obtained using $$E_\varphi ^{}$$ by the following relations: \begin{equation} E_x^{p,1}\left( {x,y,z} \right) = - E_\varphi ^{}\left( {r,z} \right)\frac{{y - {y^c}}}{r},\quad E_y^{p,1}\left( {x,y,z} \right) = E_\varphi ^{}\left( {r,z} \right)\frac{{x - {x^c}}}{r},\quad E_z^p = 0, \end{equation} (4) \begin{equation*} B_x^{p,1}\left( {x,y,z} \right) = \frac{1}{{i\omega }}\frac{{\partial {E_\varphi }\left( {r,z} \right)}}{{\partial z}}\frac{{x - {x^c}}}{r},\quad B_y^{p,1}\left( {x,y,z} \right) = \frac{1}{{i\omega }}\frac{{\partial {E_\varphi }\left( {r,z} \right)}}{{\partial z}}\frac{{y - {y^c}}}{r}, \end{equation*} \begin{equation} B_z^{p,1}\left( {x,y,z} \right) = - \frac{1}{{i\omega }}\left( {\frac{{{E_\varphi }\left( {r,z} \right)}}{r} + \frac{{\partial {E_\varphi }\left( {r,z} \right)}}{{\partial r}}} \right), \end{equation} (5)where xc and yc are the coordinates of the coil centre, and $$r = \sqrt {{{( {x - {x^c}} )}^2} + {{( {y - {y^c}} )}^2}} $$. The electrical field $${{{\bf \vec{E}}}^{p,2}}$$ induced by a ‘vertical’ coil (its axis is parallel to the borders of the layers of the horizontally layered medium) is 3-D even for a 1-D (horizontally layered) medium, that is, the components of $${{{\bf \vec{E}}}^{p,2}}$$ depend on all three space coordinates. In order to compute this field, we apply the approach similar to the one considered by Persova et al. (2011) for modelling the field induced by the horizontal electrical dipole. It is based on the representation of the field, which is in fact a 3-D field in the 1-D medium, as a sum of the fields each of which can be obtained by the solution of the problem of a less dimension. However, when we use the approach considered by Persova et al. (2011) there appears a problem (this is particularly true for harmonic modes) which is related to the fact that the total field is formed by the components with the opposite signs the moduli of which are close to each other. This makes it necessary to increase significantly the accuracy of the components calculation that, in its turn, leads to the significant increase in computational costs. Therefore, in this paper we propose the modification of the approach described by Persova et al. (2011) that makes it possible to obtain $${{{\bf \vec{E}}}^{p,2}}$$ (the field induced by the ‘vertical’ coil) as a sum of several axisymmetrical fields without mutually subtracting components (Soloveichik et al.2014; Persova et al.2014a). Each of these fields (forming $${{{\bf \vec{E}}}^{p,2}}$$) is induced by a special source and can be obtained by solving the corresponding axisymmetrical problem. At first we note that, in the problem considered, the sizes of the transmitter coil are small relative to the distance to the receivers. Therefore, we can take any geometrical form of coil. We represent the ‘vertical’ coil as a rectangular (square) frame. In this case the total source (frame) can be represented as a sum of sources shown in Fig. 1. This representation allows us to obtain the electromagnetic field in a horizontally layered background medium by means of solutions of several axisymmetrical problems. Figure 1. View largeDownload slide The view of (a) non-axisymmetric source ‘I1’ and (b) two axisymmetrical sources ‘I2’; (c) the current of the upper part of the source ‘I1’ in plane XY; (d) the currents of one of NGR sources in plane XY; (e) currents configuration which is obtained when the negative source ‘I2’ is added to the source ‘I1’; (f) currents configuration which is obtained when the positive source ‘I2’ is added to the source shown in panel (e) (as a result, the ‘vertical’ rectangular current frame is obtained); (g) the approximation of the horizontal ‘current lines’ by the point sources located in the centres of subintervals of these ‘current lines’ (the length of subintervals is $$h_L^{}$$). Dashed light-green arrows show the direction of the surface radial currents flowing in plane XY. Solid black lines show the currents flowing in the ‘vertical’ current frame. Red points show the positions of the point sources approximating the horizontal ‘current lines’. Figure 1. View largeDownload slide The view of (a) non-axisymmetric source ‘I1’ and (b) two axisymmetrical sources ‘I2’; (c) the current of the upper part of the source ‘I1’ in plane XY; (d) the currents of one of NGR sources in plane XY; (e) currents configuration which is obtained when the negative source ‘I2’ is added to the source ‘I1’; (f) currents configuration which is obtained when the positive source ‘I2’ is added to the source shown in panel (e) (as a result, the ‘vertical’ rectangular current frame is obtained); (g) the approximation of the horizontal ‘current lines’ by the point sources located in the centres of subintervals of these ‘current lines’ (the length of subintervals is $$h_L^{}$$). Dashed light-green arrows show the direction of the surface radial currents flowing in plane XY. Solid black lines show the currents flowing in the ‘vertical’ current frame. Red points show the positions of the point sources approximating the horizontal ‘current lines’. The rectangular ‘vertical’ current frame is represented as a sum of source ‘I1’ shown in Fig. 1(a) (we consider the top and bottom currents as the parts of one source ‘I1’) and two sources ‘I2’ shown in Fig. 1(b). The source ‘I1’ consists of four surface radial currents, which flow in plane XY (in Fig. 1, their direction is shown by dashed light-green arrows), and two horizontal ‘current lines’ (in Fig. 1, the currents in them are shown by solid black lines). We consider the left and right currents as different sources of the type ‘I2’, the former is negative and the latter is positive with regard to the direction of the axes shown in Fig. 1 (in this case we direct the z-axis down because in all examples presented in this paper, the direction of the coordinate axes is the same). Each of the sources ‘I2’ consists of the vertical ‘current line’ (the current in which shown by solid black line) and two surface radial currents (their direction is shown by dashed light-green arrows). Joining the top and bottom currents in one source ‘I1’ and separating the left and right currents into two sources of the type ‘I2’ are convenient for further reduction of the fields induced by these sources to the corresponding axisymmetrical fields. The Z-coordinate of the upper part of source ‘I1’ coincides with the Z-coordinate of the upper frame border, and the Z-coordinate of the lower part of the source ‘I1’ (with the opposite sign of the current) coincides with the Z-coordinate of the lower frame border. The axis of one of the sources ‘I2’ coincides with the left frame border, and the axis of another source ‘I2’ (with an opposite sign of the current) coincides with the right frame border. Convergent and divergent surface radial currents are the parts of these sources. They provide zero current divergence in sources ‘I1’ and ‘I2’. Further, we will refer to as ‘non-grounded radial sources’ (NGR sources) for each of these surface radial currents converging into the point Pi or diverging from it. Fig. 1(c) shows the current of the upper part of the source ‘I1’ in plane XY, and Fig. 1(d) shows the currents of one of NGR sources in plane XY. For each of NGR sources, the same current passes through any circle C the centre of which coincides with the NGR source centre (see Fig. 1d, the integrals along circles С1 and С2 of surface radial currents of the NGR source are equal). Figs 1(e) and (f) demonstrate how the current in the ‘vertical’ frame is composed from the source ‘I1’ and two sources ‘I2’. In Fig. 1(e) the negative source ‘I2’ with the axis, which coincides with the left side of the ‘vertical’ frame, is added to the source ‘I1’. The radial currents of the source ‘I2’ compensate the radial currents of the source ‘I1’ in the left top and left bottom angles (these currents are equal and are in the opposite directions). When the positive source ‘I2’ (with the axis coinciding with the right side of the ‘vertical’ frame) is added to the source shown in Fig. 1(e), the radial currents in the right top and right bottom angles are compensated. As a result, we obtain the source with the current in the ‘vertical’ frame shown in Fig. 1(f). We would like to reiterate that each of the sources ‘I2’ generates the axisymmetrical field (in the coordinate system in which the z-axis coincides with the vertical current of the source). The mathematical model for calculating the axisymmetrical field generated by the source ‘I2’ is presented below (see eq. 21–23). We first consider how the field induced by the non-axisymmetric source ‘I1’ can be calculated by solving several axisymmetrical problems. We introduce the vector potential $${{\bf \vec{A}}}$$ that is related with the magnetic induction and electric field by the relationships $${{\bf \vec{B}}} = {\rm{rot}}{{\bf \vec{A}}}$$ and $${{\bf \vec{E}}} = - \frac{{\partial {{\bf \vec{A}}}}}{{\partial t}} = - i\omega {{\bf \vec{A}}}$$. Then, the vector potential $${{{\bf \vec{A}}}^{{\rm{I1}}}}$$ defining the electromagnetic field induced by source ‘I1’ in a horizontally layered medium (with conductivity σp) can be found by solving the equation \begin{equation} \frac{1}{{{\mu _0}}}\nabla \times \nabla \times {{{\bf \vec{A}}}^{{\rm{I1}}}} + i\omega {\sigma ^p}{{{\bf \vec{A}}}^{{\rm{I1}}}} = {{{\bf \vec{J}}}^{{\rm{I1}}}}, \end{equation} (6)where $${{{\bf \vec{J}}}^{{\rm{I1}}}}$$ is the density of the currents defining source ‘I1’. We define this source in such a way that, first, its currents have only horizontal non-zero components, and, second, for any bounded subdomain Ωf, the value of the current entering Ωf is equal to the value of the outward current. Therefore, for the current density $${{{\bf \vec{J}}}^{{\rm{I1}}}}$$, the ratio \begin{equation} \!\nabla \cdot {{{\bf \vec{J}}}^{{\rm{I1}}}} = 0 \end{equation} (7)is true and $$J_z^{{\rm{I1}}} \equiv 0$$. It would now be easy to show that eq. (6) is satisfied for the vector potential $${{{\bf \vec{A}}}^{{\rm{I1}}}}$$ whose z-component is equal to zero and the components $$A_x^{{\rm{I1}}}$$ and $$A_y^{{\rm{I1}}}$$ satisfy the equations \begin{equation} \! - \frac{1}{{{\mu _0}}}\Delta A_x^{{\rm{I1}}} + i\omega \sigma _{11}^pA_x^{{\rm{I1}}} = J_x^{{\rm{I1}}}, \end{equation} (8) \begin{equation} \! - \frac{1}{{{\mu _0}}}\Delta A_y^{{\rm{I1}}} + i\omega \sigma _{22}^pA_y^{{\rm{I1}}} = J_y^{{\rm{I1}}}, \end{equation} (9)where Δ denotes the Laplacian, and $$\sigma _{11}^p = \sigma _{22}^p = \sigma _r^{p,rz}(z)$$. We will prove it. For this purpose, we differentiate eq. (8) with respect to x, and eq. (9) with respect to y (taking into account that $$\sigma _{11}^p$$ and $$\sigma _{22}^p$$ do not depend on x and y and equal to $$\sigma _r^{p,rz}$$). Changing the order of differentiation in the first summands of the left-hand sides of these equations and summing them, we obtain \begin{equation} \! - \frac{1}{{{\mu _0}}}\Delta \left( {\frac{{\partial\! A_x^{{\rm{I1}}}}}{{\partial x}} + \frac{{\partial\! A_y^{{\rm{I1}}}}}{{\partial y}}} \right) + i\omega \sigma _r^{p,rz}\left( {\frac{{\partial\! A_x^{{\rm{I1}}}}}{{\partial x}} + \frac{{\partial\! A_y^{{\rm{I1}}}}}{{\partial y}}} \right) = \left( {\frac{{\partial J_x^{{\rm{I1}}}}}{{\partial x}} + \frac{{\partial J_y^{{\rm{I1}}}}}{{\partial y}}} \right). \end{equation} (10)Taking into account that $$A_z^{{\rm{I1}}} \equiv 0$$ and $$J_z^{{\rm{I1}}} \equiv 0$$ (and, consequently, $$\frac{{\partial\! A_z^{{\rm{I1}}}}}{{\partial z}} = 0$$ and $$\frac{{\partial J_z^{{\rm{I1}}}}}{{\partial z}} = 0$$), we obtain \begin{equation} \! - \frac{1}{{{\mu _0}}}\Delta ( {\nabla \cdot {{{{\bf \vec{A}}}}^{{\rm{I1}}}}}) + i\omega \sigma _r^{p,rz}\left( {\nabla \cdot {{{{\bf \vec{A}}}}^{{\rm{I1}}}}} \right) = \nabla \cdot {{{\bf \vec{J}}}^{{\rm{I1}}}}. \end{equation} (11)As far as eq. (7) is true for $${{{\bf \vec{J}}}^{{\rm{I1}}}}$$, the function $$u = \nabla \cdot {{{\bf \vec{A}}}^{{\rm{I1}}}}$$ satisfies the homogeneous equation \begin{equation} \! - \frac{1}{{{\mu _0}}}\Delta u + i\omega \sigma _r^{p,rz}u = 0 \end{equation} (12)for which the function u = 0 is the solution. Thus, the divergence of vector potential $${{{\bf \vec{A}}}^{{\rm{I1}}}}$$ (whose x- and y-components satisfy eqs (8) and (9), and $$A_z^{{\rm{I1}}} \equiv 0$$) is equal to zero. Taking into account $$A_z^{{\rm{I1}}} = 0$$, the equation system (8), (9) can be written in the vector form \begin{equation} \! - \frac{1}{{{\mu _0}}}\Delta {{{\bf \vec{A}}}^{{\rm{I1}}}} + i\omega {\sigma ^p}{{{\bf \vec{A}}}^{{\rm{I1}}}} = {{{\bf \vec{J}}}^{{\rm{I1}}}}, \end{equation} (13)where Δ denotes the vector Laplacian. As far as $$\nabla \cdot {{{\bf \vec{A}}}^{{\rm{I1}}}} = 0$$ (we have just shown it), and taking into account the identical relation $$\nabla \times \nabla \times {{{\bf \vec{A}}}^{{\rm{I1}}}} \equiv - \Delta {{{\bf \vec{A}}}^{{\rm{I1}}}} + \nabla \nabla \cdot {{{\bf \vec{A}}}^{{\rm{I1}}}}$$, the vector potential $${{{\bf \vec{A}}}^{{\rm{I1}}}}$$, which is obtained by solving eqs (8) and (9) (or by solving vector eq. 13), satisfies eq. (6). Now we show how the solution of eq. (13) could be obtained by solving several axisymmetrical problems. For this purpose, we divide the source ‘I1’ into the following parts: four NGR sources and two ‘current lines’. By ‘current line’ we mean the source which is concentrated on one of the horizontal borders P1P2 or P3P4 of the ‘vertical’ current frame (see Fig. 1). The density of such sources can be defined via delta-function δL, which is concentrated on the line and defined in the following way: the integral of δL taken over an arbitrary volume Ωf is equal to the length of the L line segment, which lies inside the volume Ωf. Thus, we represent right-hand side of vector eq. (13) as a sum of three summands: \begin{equation} {{{\bf \vec{J}}}^{{\rm{I1}}}} = {{{\bf \vec{J}}}^{{\rm{I1,NGR1}}}} + {{{\bf \vec{J}}}^{{\rm{I1,NGR2}}}} + {{{\bf \vec{J}}}^{{\rm{I1,hor}}}}, \end{equation} (14)where $${{{\bf \vec{J}}}^{{\rm{I1,NGR1}}}}$$ is the density of surface currents in two (axial) NGR sources, while the currents of one of these sources converge into the point P1, and currents of another source diverge from the point P3 (see Fig. 1a); $${{{\bf \vec{J}}}^{{\rm{I1,NGR2}}}}$$ is the density of surface currents in other two (axial) NGR sources (currents of one of these sources diverge from the point P2, and currents of another source converge to the point P4); $${{{\bf \vec{J}}}^{{\rm{I1,hor}}}} = {{{\bf \vec{J}}}^{{\rm{I1,}}{{\rm{P}}_1}{{\rm{P}}_2}}} - {{{\bf \vec{J}}}^{{\rm{I1,}}{{\rm{P}}_3}{{\rm{P}}_4}}}$$ is the density of currents in the horizontal ‘current lines’ P1P2 and P3P4. Then, the solution of eq. (13) can be found as a sum of solutions of three equations (they are in fact independent) \begin{equation} \! - \frac{1}{{{\mu _0}}}\Delta {{{\bf \vec{A}}}^{{\rm{I1,NGR1}}}} + i\omega {\sigma ^p}{{{\bf \vec{A}}}^{{\rm{I1,NGR1}}}} = {{{\bf \vec{J}}}^{{\rm{I1,NGR1}}}}, \end{equation} (15) \begin{equation} \! - \frac{1}{{{\mu _0}}}\Delta {{{\bf \vec{A}}}^{{\rm{I1,NGR2}}}} + i\omega {\sigma ^p}{{{\bf \vec{A}}}^{{\rm{I1,NGR2}}}} = {{{\bf \vec{J}}}^{{\rm{I1,NGR2}}}}, \end{equation} (16) \begin{equation} \! - \frac{1}{{{\mu _0}}}\Delta {{{\bf \vec{A}}}^{{\rm{I1,hor}}}} + i\omega {\sigma ^p}{{{\bf \vec{A}}}^{{\rm{I1,hor}}}} = {{{\bf \vec{J}}}^{{\rm{I1,}}{{\rm{P}}_1}{{\rm{P}}_2}}} - {{{\bf \vec{J}}}^{{\rm{I1,}}{{\rm{P}}_3}{{\rm{P}}_4}}}. \end{equation} (17)Right-hand sides of eqs (15) and (16) are the functions with axial symmetry (the z-axis passes through the points P1,P3 or P2,P4, respectively). In the corresponding cylindrical coordinate system, the functions $${{{\bf \vec{J}}}^{{\rm{I1,NGR1}}}}$$ and $${{{\bf \vec{J}}}^{{\rm{I1,NGR2}}}}$$ of the right-hand sides of these equations have the only non-zero radial component: $$J_r^{{\rm{I1,NGR1}}} = \frac{I}{{2\pi r}}( {{\delta ^{{{\hat{\Gamma }}_2}}} - {\delta ^{{{\hat{\Gamma }}_1}}}} )$$ and $$J_r^{{\rm{I1,NGR2}}} = \frac{I}{{2\pi r}}( {{\delta ^{{{\hat{\Gamma }}_1}}} - {\delta ^{{{\hat{\Gamma }}_2}}}} )$$, where I is the value of current in the ‘vertical’ frame, $${\delta ^{{{\hat{\Gamma }}_1}}}$$ and $${\delta ^{{{\hat{\Gamma }}_2}}}$$ are the delta-functions concentrated on the surfaces $${\hat{\Gamma }_1}$$ and $${\hat{\Gamma }_2}$$ (i.e. on the planes z = z1 and z = z2 corresponding to top and bottom frame borders, see Fig. 1) where the radial currents of the top and bottom NGR sources flow. Thus, the solution of both eqs (15) and (16) can be found as a result of solving one axisymmetrical problem: find the solution $$A_r^{}( {r,z} )$$ of equation \begin{equation} \! - \frac{1}{{{\mu _0}}}\Delta A_r^{} + \frac{1}{{{\mu _0}{r^2}}}A_r^{} + i\omega \sigma _r^{p,rz}A_r^{} = \frac{I}{{2\pi r}}( {{\delta ^{{{\hat{\Gamma }}_1}}} - {\delta ^{{{\hat{\Gamma }}_2}}}} ), \end{equation} (18)where $$\sigma _r^{p,rz}$$ is the lateral conductivity of a horizontally layered medium. The field $${{{\bf \vec{A}}}^{{\rm{I1,hor}}}}$$ defined by eq. (17) is not axisymmetrical. The vector function $${{{\bf \vec{J}}}^{{\rm{I1,}}{{\rm{P}}_1}{{\rm{P}}_2}}} - {{{\bf \vec{J}}}^{{\rm{I1,}}{{\rm{P}}_3}{{\rm{P}}_4}}}$$ defining the right-hand side of this equation has only one non-zero component depending on the orientation of the ‘vertical’ current frame: if the frame is oriented in the way shown in Fig. 1, then the y-component is the only non-zero component of the right-hand side vector of eq. (17) and, consequently, the vector function $${{{\bf \vec{A}}}^{{\rm{I1,hor}}}}$$ has the only non-zero component $$A_y^{{\rm{hor}}}( {x,y,z} )$$. If the ‘vertical’ current frame is oriented so that the lines P1P2 and P3P4 are parallel to the x-axis, then the x-component is the only non-zero component of the right-hand side vector of eq. (17) and, therefore, the vector function $${{{\bf \vec{A}}}^{{\rm{I1,hor}}}}$$ has the only non-zero component $$A_x^{{\rm{hor}}}( {x,y,z} )$$. These components can be found by solving the 3-D scalar equation \begin{equation} \! - \frac{1}{{{\mu _0}}}\Delta A_\xi ^{{\rm{hor}}} + i\omega \sigma _r^{p,rz}A_\xi ^{{\rm{hor}}} = I\left( {{\delta ^{{{\rm{P}}_1}{{\rm{P}}_2}}} - {\delta ^{{{\rm{P}}_3}{{\rm{P}}_4}}}} \right), \end{equation} (19)where $$A_\xi ^{{\rm{hor}}}( {x,y,z} )$$ is the vector potential component corresponding to the horizontal line direction (x-direction or y-direction), $${\delta ^{{{\rm{P}}_1}{{\rm{P}}_2}}}$$ and $${\delta ^{{{\rm{P}}_3}{{\rm{P}}_4}}}$$ are delta-functions concentrated on the lines P1P2 and P3P4. The right-hand side of eq. (19) can be approximated with any required accuracy by the set of point sources positioned along the lines P1P2 and P3P4 (see Fig. 1g) and, therefore, the 3-D component $$A_\xi ^{{\rm{hor}}}( {x,y,z} )$$ of the field induced by the sources concentrated on the lines P1P2 and P3P4 can be computed as a sum of axisymmetrical fields $$A_\xi ^{}( {r,z} )$$ induced by the point sources positioned along these lines. The field induced by two point sources (lying one below (above) another) can be calculated by the solution of the equation \begin{equation} \! - \frac{1}{{{\mu _0}}}\Delta A_\xi ^{} + i\omega \sigma _r^{p,rz}A_\xi ^{} = I{h_L}\left( {{\delta ^{{z_1}}} - {\delta ^{{z_2}}}} \right), \end{equation} (20)where $$A_\xi ^{}( {r,z} )$$ is the vector potential component corresponding to the ‘current line’ direction (i.e. x- or y-component), $${\delta ^{{z_1}}}$$ and $${\delta ^{{z_2}}}$$ are the delta-functions concentrated in the points (centres of the ‘current line’ subintervals, in which the point sources are located), hL is the length of the ‘current line’ subinterval (see Fig. 1g), and, I as before, is the value of the current in the frame. By summing axisymmetrical fields $$A_\xi ^{}( {r,z} )$$ taking into account the translation of the cylindrical coordinate system centres into the corresponding points of one of the ‘current lines’ of source ‘I1’, we obtain the 3-D distribution of one component of the vector potential, which describes the field of the ‘current line’. Thus, we can calculate the field induced by the non-axisymmetric source ‘I1’ by summing axisymmetrical fields taking into account the translation of the cylindrical coordinate system centres into corresponding points. The formulae for calculating the electric and magnetic field will be presented below (we take into account the fact that, when summing the components of the fields induced by separate parts of the source, the radial currents are compensated, and, therefore, the fields induced by these radial currents are not added to the resulting electric and magnetic field). In order to compute the axisymmetrical field induced by the source of type ‘I2’, we use the following mathematical model. Taking into account that the divergence of current density $${{{\bf \vec{J}}}^{{\rm{I2}}}} = {( {{J_r},{J_z}} )^{\rm{T}}}$$ is equal to zero, the electromagnetic field induced by each source ‘I2’ can be obtained by the solution of the differential equation system \begin{equation} \! - \frac{{\rm{1}}}{{{\mu _0}}}\Delta A_r^{} + \frac{{\rm{1}}}{{{\mu _0}{r^2}}}A_r^{} + \sigma _r^{p,rz}\left( {i\omega A_r^{} + \frac{{\partial V}}{{\partial r}}} \right) = {J_r}, \end{equation} (21) \begin{equation} \! - \frac{{\rm{1}}}{{{\mu _0}}}\Delta A_z^{} + \sigma _z^{p,rz}\left( {i\omega A_z^{} + \frac{{\partial V}}{{\partial z}}} \right) = {J_z}, \end{equation} (22) \begin{equation} \! - \nabla \cdot \left( {{\sigma ^{p,rz}}\nabla V} \right) - i\omega \nabla \cdot \left( {{\sigma ^{p,rz}}{{{{\bf \vec{A}}}}^{{\rm{I2}}}}} \right) = 0, \end{equation} (23)where $${{{\bf \vec{A}}}^{{\rm{I2}}}} = {( {{A_r}( {r,z} ),{A_z}( {r,z} )} )^{\rm{T}}}$$ is the vector potential in the cylindrical coordinate system, V = V(r, z) is the scalar potential, and Δ = ∇ · ∇ is the scalar Laplacian. We represent the r-component of the vector potential $${{{\bf \vec{A}}}^{{\rm{I2}}}}$$ as $$A_r^{} = A_r^* + \tilde{A}_r^{}$$ where $$A_r^*$$ is the solution of equation in the form of (18). Then system (21)–(23) takes the form \begin{equation} \! - \frac{{\rm{1}}}{{{\mu _0}}}\Delta \tilde{A}_r^{} + \frac{{\rm{1}}}{{{\mu _0}{r^2}}}\tilde{A}_r^{} + \sigma _r^{p,rz}\left( {i\omega \tilde{A}_r^{} + \frac{{\partial V}}{{\partial r}}} \right) = 0, \end{equation} (24) \begin{equation} \! - \frac{{\rm{1}}}{{{\mu _0}}}\Delta A_z^{} + \sigma _z^{p,rz}\left( {i\omega A_z^{} + \frac{{\partial V}}{{\partial z}}} \right) = {J_z}, \end{equation} (25) \begin{equation} \! - \nabla \cdot \left( {{\sigma ^{p,rz}}\nabla V} \right) - i\omega \nabla \cdot \left( {{\sigma ^{p,rz}}{{( {\tilde{A}_r^{},A_z^{}} )}^{\rm{T}}}} \right) = i\omega \nabla \cdot \left( {{\sigma ^{p,rz}}{{( {A_r^*,0} )}^{\rm{T}}}} \right). \end{equation} (26)The field $$A_r^*$$ in the right-hand side of eq. (26) is fully identical with $$A_r^{}$$ calculated for the source ‘I1’ (with the opposite sign) that allows us to solve eq. (18) only once and (what is more important!) to exclude these two mutually subtracting fields when we form the total field of the ‘vertical’ frame. Thus, in order to compute the field of the ‘vertical’ frame, it is necessary to solve three problems: eq. (18), (20) and equation system (24)–(26). Besides, the proposed scheme can be further simplified (that will increase its computational effectiveness) because there is no need to solve eq. (18). We define function u as $$u = \nabla \cdot {(A_r^*,0)^{\rm{T}}}$$. Differentiating eq. (18) properly we obtain the equation \begin{equation} \! - \frac{1}{{{\mu _0}}}\Delta u + i\omega \sigma _r^{p,rz}u = I\left( {{\delta ^{P,1}} - {\delta ^{P,2}}} \right), \end{equation} (27)which is identical to eq. (20). Equation system (24)–(26) takes the form \begin{equation} \! - \frac{{\rm{1}}}{{{\mu _0}}}\Delta \tilde{A}_r^{} + \frac{{\rm{1}}}{{{\mu _0}{r^2}}}\tilde{A}_r^{} + \sigma _r^{p,rz}\left( {i\omega \tilde{A}_r^{} + \frac{{\partial V}}{{\partial r}}} \right) = 0, \end{equation} (28) \begin{equation} \! - \frac{{\rm{1}}}{{{\mu _0}}}\Delta A_z^{} + \sigma _z^{p,rz}\left( {i\omega A_z^{} + \frac{{\partial V}}{{\partial z}}} \right) = {J_z}, \end{equation} (29) \begin{equation} \! - \nabla \cdot \left( {{\sigma ^{p,rz}}\nabla V} \right) - i\omega \nabla \cdot \left( {{\sigma ^{p,rz}}{{( {\tilde{A}_r^{},A_z^{}} )}^{\rm{T}}}} \right) = i\omega \sigma _r^{p,rz}u. \end{equation} (30)Finally, in order to obtain the field of the ‘vertical’ frame, it is necessary to solve two problems: eq. (20) and equation system (28)–(30). Further, we consider how to calculate the components of $${{{\bf \vec{E}}}^{p,2}}$$ using the solution of eq. (20) and equation system (28)–(30). In the case when the ‘vertical’ frame axis is parallel to the x-axis (i.e. the coordinates of its vertexes are (xc, y1, z1), (xc, y2, z1), (xc, y2, z2), and (xc, y1, z2)), the values of the components of $${{{\bf \vec{E}}}^{p,2}}$$ at the point (x, y, z) can be calculated with the help of the following formulae: \begin{equation} E_x^{p,2}\left( {x,y,z} \right) = \left( {i\omega \tilde{A}_r^{}\left( {{r_1},z} \right) + \frac{{\partial V\left( {{r_1},z} \right)}}{{\partial r}}} \right) \cdot \frac{{\left( {x - {x_c}} \right)}}{{{r_1}}} - \left( {i\omega \tilde{A}_r^{}\left( {{r_2},z} \right) + \frac{{\partial V\left( {{r_2},z} \right)}}{{\partial r}}} \right) \cdot \frac{{\left( {x - {x_c}} \right)}}{{{r_2}}}, \end{equation} (31) \begin{equation} E_y^{p,2}\left( {x,y,z} \right) = - i\omega A_y^{{\rm{hor}}} + \left( {i\omega \tilde{A}_r^{}\left( {{r_1},z} \right) + \frac{{\partial V\left( {{r_1},z} \right)}}{{\partial r}}} \right) \cdot \frac{{\left( {y - {y_1}} \right)}}{{{r_1}}} - \left( {i\omega \tilde{A}_r^{}\left( {{r_2},z} \right) + \frac{{\partial V\left( {{r_2},z} \right)}}{{\partial r}}} \right) \cdot \frac{{\left( {y - {y_2}} \right)}}{{{r_2}}}, \end{equation} (32) \begin{equation} E_z^{p,2}\left( {x,y,z} \right) = \left( {i\omega \tilde{A}_z^{}\left( {{r_1},z} \right) + \frac{{\partial V\left( {{r_1},z} \right)}}{{\partial z}}} \right) - \left( {i\omega \tilde{A}_z^{}\left( {{r_2},z} \right) + \frac{{\partial V\left( {{r_2},z} \right)}}{{\partial z}}} \right), \end{equation} (33)where $${r_1} = \sqrt {{{( {x - {x_c}} )}^2} + {{( {y - {y_1}} )}^2}} $$, $${r_2} = \sqrt {{{( {x - {x_c}} )}^2} + {{( {y - {y_2}} )}^2}} $$ and $$A_y^{{\rm{hor}}}( {x,y,z} )$$ is obtained by solving eq. (27) as a result of summing the functions u(r, z) (after recalculating u(r, z) from the cylindrical coordinate system into the Cartesian one), which are calculated for the sources positioned in the centres of subintervals of the horizontal ‘current lines’ of the frame. The components of $${{{\bf \vec{B}}}^{p,2}}$$ at the point (x, y, z) can be calculated with the help of the following formulae: \begin{equation*} B_x^{p,2}\left( {x,y,z} \right) = - \frac{{\partial\! A_y^{{\rm{hor}}}\left( {r,z} \right)}}{{\partial z}} + \left( {\frac{{\partial {A_r}\left( {{r_1},z} \right)}}{{\partial z}} - \frac{{\partial {A_z}\left( {{r_1},z} \right)}}{{\partial r}}} \right) \cdot \frac{{\left( {y - {y_1}} \right)}}{{{r_1}}} - \left( {\frac{{\partial {A_r}\left( {{r_2},z} \right)}}{{\partial z}} - \frac{{\partial {A_z}\left( {{r_2},z} \right)}}{{\partial r}}} \right) \cdot \frac{{\left( {y - {y_2}} \right)}}{{{r_2}}}, \end{equation*} \begin{equation*} B_y^{p,2}\left( {x,y,z} \right) = - \left( {\frac{{\partial {A_r}\left( {{r_1},z} \right)}}{{\partial z}} - \frac{{\partial {A_z}\left( {{r_1},z} \right)}}{{\partial r}}} \right) \cdot \frac{{\left( {x - {x^c}} \right)}}{{{r_1}}} + \left( {\frac{{\partial {A_r}\left( {{r_2},z} \right)}}{{\partial z}} - \frac{{\partial {A_z}\left( {{r_2},z} \right)}}{{\partial r}}} \right) \cdot \frac{{\left( {x - {x^c}} \right)}}{{{r_2}}}, \end{equation*} \begin{equation*} B_z^{p,2}\left( {x,y,z} \right) = \frac{{\partial\! A_y^{{\rm{hor}}}\left( {r,z} \right)}}{{\partial x}}. \end{equation*} If the transmitter coil is oriented so that the angle between its axis and the axis z is equal to α (see Fig. 2a), then $${{{\bf \vec{E}}}^p}$$ can be calculated using $${{{\bf \vec{E}}}^{p,1}}$$ and $${{{\bf \vec{E}}}^{p,2}}$$ by formula \begin{equation*} {{{\bf \vec{E}}}^p} = {{{\bf \vec{E}}}^{p,1}}\cos \alpha + {{{\bf \vec{E}}}^{p,2}}{\sin \alpha}. \end{equation*} Figure 2. View largeDownload slide (a) The orientation of transmitter whose axis is perpendicular to the y-axis and (b) the rotation of the coordinate system when the transmitter orientation is arbitrary. Figure 2. View largeDownload slide (a) The orientation of transmitter whose axis is perpendicular to the y-axis and (b) the rotation of the coordinate system when the transmitter orientation is arbitrary. For arbitrary transmitter coil orientation, at first, we can calculate $${{{\bf \vec{\tilde{E}}}}^p}$$ by formulae (31)–(33) in the coordinate system $$\{ {\tilde{x},\tilde{y},\tilde{z}} \}$$, in which the coil axis is perpendicular to the $$\tilde{y}$$-axis, and then recalculate $${{{\bf \vec{\tilde{E}}}}^p}$$ into the coordinate system {x, y, z} by formulae $${{{\bf \vec{E}}}^p} = {{\bf M}}{{{\bf \vec{\tilde{E}}}}^p}$$ where M is the rotation matrix defined as follows: \begin{equation} {{\bf M}} = \left( {\begin{array}{c@{\quad}c@{\quad}c} {\cos \beta }&{ - \sin \beta }&0\\ {\sin \beta }&{\cos \beta }&0\\ 0&0&1 \end{array}} \right). \end{equation} (34)In formula (34), β is the angle of rotation around the axis $$z = \tilde{z}$$ which transfers the coordinate system {x, y} to the coordinate system $$\{ {\tilde{x},\tilde{y}} \}$$(see Fig. 2b). The components of the magnetic field induction are calculated in a similar way. Thus, the proposed approach together with the possibility to take the field of the source in the horizontally layered medium as $${{{\bf \vec{E}}}^p}$$, which is much similar to the realistic 3-D medium in the vicinity of the logging tool location, makes it possible to exclude the appearance of the opposite sign components of the field, the moduli of which are close to each other. 3 FINITE-ELEMENT APPROXIMATION FE solutions of all the equations considered above and describing the secondary field $${{{\bf \vec{E}}}^s}$$ and the components of the primary field (which are axisymmetrical, if the coordinate system is properly translated) are found in the limited 3-D domain Ω (for $${{{\bf \vec{E}}}^s}$$) or 2-D domain Ωrz (for components of $${{{\bf \vec{E}}}^p}$$). On the outer boundary of Ω and Ωrz, we set homogeneous essential boundary conditions. For $${{{\bf \vec{E}}}^s}$$, these conditions represent the equality of $${{{\bf \vec{E}}}^s}$$ components to the zero, which are tangential to the boundary of Ω: $${{{\bf \vec{E}}}^s} \times {{\bf \vec{n}}} = {{\bf \vec{0}}}$$ where $${{\bf \vec{n}}}$$ is the normal to the boundary of Ω (these essential boundary conditions are set when the vector FEM is used). If the field $${{{\bf \vec{E}}}^s}$$ has a symmetry plane, then on this plane, we set either the homogeneous natural boundary condition $$( {\nabla \times {{{{\bf \vec{E}}}}^s}} ) \times {{\bf \vec{n}}} = {{\bf \vec{0}}}$$ (where $${{\bf \vec{n}}}$$ is the normal to the symmetry plane), if the current frame lies in this plane, or the homogeneous essential boundary condition $${{{\bf \vec{E}}}^s} \times {{\bf \vec{n}}} = {{\bf \vec{0}}}$$, if the current frame is sectioned by the symmetry plane. For calculating the axisymmetrical components of the field $${{{\bf \vec{E}}}^p}$$, on the outer boundaries of the computational domain Ωrz, we set homogeneous Dirichlet boundary conditions for all potentials describing the corresponding field. The size of computational domains Ω and Ωrz depends on the frequency and ‘average’ conductivity of the medium. Besides, the size of Ω is much smaller than the size of Ωrz. For example, for the frequency of 200 kHz, Ω = [−20,20] × [−15,15] × [−20,20] m3 (or Ω = [−20,20] × [0,15] × [−20,20] m3 if y = 0 is the symmetry plane of the field), and Ωrz = [0,70] × [−70,70] m2 (if the ‘average’ conductivity of the medium is about 0.05–1 Sm m−1). 3.1 Finite-element approximation for mathematical models used for $${{{\bf \vec{E}}}^p}$$ calculation The equivalent variational formulation for eq. (3), which defines the field of ‘horizontal’ coil (loop), takes the following form: \begin{equation} \frac{1}{{{\mu _0}}}\left( {\int_{{{\Omega ^{rz}}}}{{\nabla E_\varphi ^{} \cdot \nabla \nu r{\rm{d}}r{\rm{d}}z}} + \int_{{{\Omega ^{rz}}}}{{\frac{1}{r}E_\varphi ^{}\nu {\rm{d}}r{\rm{d}}z}}} \right) + i\omega \int_{{{\Omega ^{rz}}}}{{\sigma _r^{p,rz}E_\varphi ^{}\nu r{\rm{d}}r{\rm{d}}z}} = - i\omega I\nu \left( {{r_\delta },{z_\delta }} \right){r_\delta }, \end{equation} (35)where ν(r, z) is a real test function, Ωrz is the computational domain in the cylindrical coordinate system. The expression in the right-hand side of eq. (35) is obtained as a result of the delta-function $${\delta ^{{J_\varphi }}}$$ (which concentrates on the circle r = rδ, z = zδ) operation on the test function ν(r, z). More specifically, the function $$E_\varphi ^{}$$ satisfying eqs (3) and (35) goes to infinity when approaching the circle r = rδ, z = zδ. The FE solution at the point (rδ, zδ) of the cylindrical domain Ωrz is finite, but together with mesh refinement around this point, the FE solution will grow infinitely. However, at some fixed (non-zero) distance from the point (rδ, zδ), the value of $$E_\varphi ^{}$$ can be calculated with any required accuracy by means of FEM. Therefore, when calculating $${{{\bf \vec{E}}}^{p,1}}$$, this fact must be taken into account, and when forming $${{{\bf \vec{E}}}^p}$$ for the solution of eq. (1), it is not allowed to come to the circle r = rδ, z = zδ closer than some pre-set distance. The solution $$E_\varphi ^{}$$ of eq. (35) is expressed as \begin{equation} E_\varphi ^{} = \sum\limits_{j = 1}^{{n^{rz}}} {q_j^{{E_\varphi }}{\psi _j}} , \end{equation} (36)where ψj(r, z) are (real) nodal basis functions, nrz is their number (which coincides with the number of nodes in the 2-D mesh). As a result of substituting (36) into (35) and replacing the test function ν with all basis functions ψl in turn, the following system of linear algebraic equations (SLAE) is obtained: \begin{equation} \!\left( {{{{\bf G}}^{{E_\varphi }}} + i\omega {{{\bf M}}^{{E_\varphi }}}} \right){{{\bf q}}^{{E_\varphi }}} = - i\omega {{{\bf f}}^{{E_\varphi }}}, \end{equation} (37)where the components of the matrices $${{{\bf G}}^{{E_\varphi }}}$$, $${{{\bf M}}^{{E_\varphi }}}$$ and vector $${{{\bf f}}^{{E_\varphi }}}$$ are calculated by formulae \begin{equation*} G_{lj}^{{E_\varphi }} = \frac{1}{{{\mu _0}}}\left( {\int_{{{\Omega ^{rz}}}}{}\nabla {\psi _l} \cdot \nabla {\psi _j}r{\rm{d}}r{\rm{d}}z + \int_{{{\Omega ^{rz}}}}{}\frac{1}{r}{\psi _l}{\psi _j}{\rm{d}}r{\rm{d}}z} \right),\quad M_{lj}^{{E_\varphi }} = \int_{{{\Omega ^{rz}}}}{}\sigma _r^{p,rz}{\psi _l}{\psi _j}r{\rm{d}}r{\rm{d}}z, \end{equation*} \begin{equation} f_l^{{E_\varphi }} = I{\psi _l}\left( {{r_\delta },{z_\delta }} \right){r_\delta }. \end{equation} (38)When FE approximation is used for this problem, we exclude a small subdomain in the vicinity of r = 0 from the computational domain Ωrz, on its boundary Γ0 the homogeneous Dirichlet boundary condition $${ {E_\varphi ^{}} |_{{\Gamma _0}}} = 0$$ is set (as on the distant boundaries of Ωrz). As mentioned above, for obtaining the field of the ‘vertical’ frame, it is necessary to solve two problems: eq. (20) and equation system (28)–(30). The equivalent variational formulation for eq. (20) takes the form \begin{equation} \frac{1}{{{\mu _0}}}\int_{{{\Omega ^{rz}}}}{{\nabla A_\xi ^{} \cdot \nabla \nu r{\rm{d}}r{\rm{d}}z}} + i\omega \int_{{{\Omega ^{rz}}}}{{\sigma _r^{p,rz}A_\xi ^{}\nu r{\rm{d}}r{\rm{d}}z}} = Ih_L^{}\left( {\nu \left( {0,{z_1}} \right) - \nu \left( {0,{z_2}} \right)} \right), \end{equation} (39)where ν(0, z1) and ν(0, z2) are the results of delta-functions $${\delta ^{{z_1}}}$$ and $${\delta ^{{z_2}}}$$ operations on the test function ν(r, z). When this problem is solved in the vicinities of the points (0, z1) and (0, z2), FE solution can grow infinitely, and, therefore, when calculating $${{{\bf \vec{E}}}^{p,2}}$$, it is not allowed to come to these points closer than some pre-set distance. The solution $$A_\xi ^{}$$ of eq. (39) is also expressed as a linear combination of basis functions ψj(r, z) \begin{equation} A_\xi ^{} = \sum\limits_{j = 1}^{{n^{rz}}} {q_j^{A_\xi ^{}}{\psi _j}} . \end{equation} (40)As a result of substituting (40) into (39) and replacing the test function with all basis functions ψl in turn, the following SLAE is obtained: \begin{equation} \left( {{{{\bf G}}^{A_\xi ^{}}} + i\omega {{{\bf M}}^{A_\xi ^{}}}} \right){{{\bf q}}^{A_\xi ^{}}} = {{{\bf f}}^{A_\xi ^{}}}, \end{equation} (41)where the components of matrices $${{{\bf G}}^{A_\xi ^{}}}$$, $${{{\bf M}}^{A_\xi ^{}}}$$, and vector $${{{\bf f}}^{A_\xi ^{}}}$$ are calculated by formulae \begin{equation*} G_{lj}^{A_\xi ^{}} = \frac{1}{{{\mu _0}}}\int_{{{\Omega ^{rz}}}}{}\nabla {\psi _l} \cdot \nabla {\psi _j}r{\rm{d}}r{\rm{d}}z,\quad M_{lj}^{A_\xi ^{}} = \int_{{{\Omega ^{rz}}}}{}\sigma _r^{p,rz}{\psi _l}{\psi _j}r{\rm{d}}r{\rm{d}}z, \end{equation*} \begin{equation} f_l^{A_\xi ^{}} = Ih_L^{}\left( {{\psi _l}\left( {0,{z_1}} \right) - {\psi _l}\left( {0,{z_2}} \right)} \right). \end{equation} (42)The equivalent variational formulation for the equation system (28)–(30) takes the form \begin{equation} \frac{1}{{{\mu _0}}}\left( {\int_{{{\Omega ^{rz}}}}{{\nabla \tilde{A}_r^{} \cdot \nabla \nu r{\rm{d}}r{\rm{d}}z}} + \int_{{{\Omega ^{rz}}}}{{\frac{1}{r}\tilde{A}_r^{}\nu {\rm{d}}r{\rm{d}}z}}} \right) + i\omega \int_{{{\Omega ^{rz}}}}{{\sigma _r^{p,rz}\tilde{A}_r^{}\nu r{\rm{d}}r{\rm{d}}z}} + \int_{{{\Omega ^{rz}}}}{{\sigma _r^{p,rz}\frac{{\partial V}}{{\partial r}}\nu r{\rm{d}}r{\rm{d}}z}} = 0, \end{equation} (43) \begin{equation} \frac{1}{{{\mu _0}}}\int_{{{\Omega ^{rz}}}}{{\nabla A_z^{} \cdot \nabla \nu r{\rm{d}}r{\rm{d}}z}} + i\omega \int_{{{\Omega ^{rz}}}}{{\sigma _z^{p,rz}A_z^{}\nu r{\rm{d}}r{\rm{d}}z}} + \int_{{{\Omega ^{rz}}}}{{\sigma _z^{p,rz}\frac{{\partial V}}{{\partial z}}\nu r{\rm{d}}r{\rm{d}}z}} = \int_{{{\Omega ^{rz}}}}{{{J_z}\nu r{\rm{d}}r{\rm{d}}z}}, \end{equation} (44) \begin{equation} \int_{{{\Omega ^{rz}}}}{{{\sigma ^{p,rz}}\nabla V \cdot \nabla \nu \,\,r{\rm{d}}r{\rm{d}}z}} + i\omega \int_{{{\Omega ^{rz}}}}{{{\sigma ^{p,rz}}{{( {\tilde{A}_r^{},A_z^{}} )}^{\rm{T}}} \cdot \nabla \nu \,\,r{\rm{d}}r{\rm{d}}z}} = i\omega \int_{{{\Omega ^{rz}}}}{{\sigma _r^{p,rz}u\nu r{\rm{d}}r{\rm{d}}z}}, \end{equation} (45)where ν(r, z) is the test function as before. Functions $$\tilde{A}_r^{}$$, $$A_z^{}$$ and V are represented in the form \begin{equation} \tilde{A}_r^{} = \sum\limits_{i = 1}^{{n^{rz}}} {q_i^{\tilde{A}_r^{}}{\psi _i}} ,\quad A_z^{} = \sum\limits_{i = 1}^{{n^{rz}}} {q_i^{{A_z}}{\psi _i}} ,\quad V = \sum\limits_{i = 1}^{{n^{rz}}} {q_i^V{\psi _i}} . \end{equation} (46)Values $$q_i^{\tilde{A}_r^{}}$$, $$q_i^{{A_z}}$$ and $$q_i^V$$ are found from the SLAE which is obtained by substituting (46) into equation system (43)–(45) and replacing the test function ν with the basis function ψi in each equation. The obtained SLAE takes the following form: \begin{equation} {{{\bf C}}^{AV}}{{{\bf q}}^{AV}} = {{{\bf f}}^{AV}}. \end{equation} (47)The size of system (47) is equal to 3nrz × 3nrz. The matrix CAV has a block structure with 3 × 3 block elements, vectors qAV and fAV also have a block structure with 3 × 1 block elements. The vector qAV has the following form: $${{{\bf q}}^{AV}} = {\{ {q_1^{\tilde{A}_r^{}},q_1^{A_z^{}},q_1^V,q_2^{\tilde{A}_r^{}},q_2^{A_z^{}},q_2^V,...} \}^{\rm{T}}}$$. The block elements of the matrix CAV take the following form: \begin{equation} {{\bf C}}_{lj}^{AV} = \left[ {\begin{array}{c@{\quad}c@{\quad}c} {G_{ij}^{{{\tilde{A}}_r}} + i\omega M_{lj}^r}&0&{D_{lj}^r}\\ 0&{G_{ij}^{{A_z}} + i\omega M_{lj}^z}&{D_{lj}^z}\\ {i\omega D_{jl}^r}&{i\omega D_{jl}^z}&{G_{ij}^V} \end{array}} \right], \end{equation} (48)where \begin{equation} G_{lj}^{{{\tilde{A}}_r}} = \frac{1}{{{\mu _0}}}\left( {\int_{{{\Omega ^{rz}}}}{}\nabla {\psi _l} \cdot \nabla {\psi _j}r{\rm{d}}r{\rm{d}}z + \int_{{{\Omega ^{rz}}}}{}\frac{1}{r}{\psi _l}{\psi _j}{\rm{d}}r{\rm{d}}z} \right),\quad G_{lj}^{{A_z}} = \frac{1}{{{\mu _0}}}\int_{{{\Omega ^{rz}}}}{}\nabla {\psi _l} \cdot \nabla {\psi _j}r{\rm{d}}r{\rm{d}}z, \end{equation} (49) \begin{equation} G_{lj}^V = \int_{{{\Omega ^{rz}}}}{{{\sigma ^{p,rz}}\nabla {\psi _l} \cdot \nabla {\psi _j}\,\,r{\rm{d}}r{\rm{d}}z}},\quad D_{lj}^r = \int_{{{\Omega ^{rz}}}}{}\sigma _r^{p,rz}{\psi _l}\frac{{\partial {\psi _j}}}{{\partial r}}r{\rm{d}}r{\rm{d}}z,\quad D_{lj}^z = \int_{{{\Omega ^{rz}}}}{}\sigma _z^{p,rz}{\psi _l}\frac{{\partial {\psi _j}}}{{\partial z}}r{\rm{d}}r{\rm{d}}z, \end{equation} (50) \begin{equation} M_{lj}^r = \int_{{{\Omega ^{rz}}}}{}\sigma _r^{p,rz}{\psi _l}{\psi _j}r{\rm{d}}r{\rm{d}}z,\quad M_{lj}^z = \int_{{{\Omega ^{rz}}}}{}\sigma _z^{p,rz}{\psi _l}{\psi _j}r{\rm{d}}r{\rm{d}}z. \end{equation} (51)The block elements of the vector fAV take the following form: \begin{equation} {{\bf f}}_l^{AV} = \left[ {\begin{array}{@{}*{1}{c}@{}} 0\\ {f_l^z}\\ {i\omega f_l^u} \end{array}} \right],\quad f_l^z = \int_{{{\Omega ^{rz}}}}{{{J_z}{\psi _l}r{\rm{d}}r{\rm{d}}z}},\quad f_l^u = \int_{{{\Omega ^{rz}}}}{{\sigma _r^{p,rz}u{\psi _l}r{\rm{d}}r{\rm{d}}z}}. \end{equation} (52)Note that when calculating $$f_l^u$$ in (52), the function u is found as $$u = \sum_{j = 1}^{{n^{rz}}} {q_j^u{\psi _j}} $$ by solving the equation, which is analogous to eq. (20) (or 39), that is, with the use of eqs (41) and (42) when $$h_L^{} = 1$$. Therefore, the component $$f_l^u$$ in (52) can be calculated as $$f_l^u = \sum_j {M_{lj}^rq_j^u} $$ where $$M_{lj}^r$$ is defined by (51). There is a further point to be made here. In order to obtain the complex symmetrical SLAE, every third row of the system must be divided by iω. 3.2 Finite-element approximation for mathematical models used for calculating the 3-D field in anisotropic medium The equivalent variational formulation for eq. (1) takes the following form: \begin{eqnarray} \int_{\Omega }{{\frac{1}{{{\mu _0}}}}}( {\nabla \times {{{{\bf \vec{E}}}}^s}} ) \cdot ( {\nabla \times {{\bf \vec{\boldsymbol\psi }}}} ){\rm{d}}\Omega + i\omega \int_{\Omega }{{\sigma {{{{\bf \vec{E}}}}^s}}} \cdot {{\bf \vec{\boldsymbol\psi }}}{\rm{d}}\Omega - {\omega ^2}\int_{\Omega }{{\varepsilon {{{{\bf \vec{E}}}}^s}}} \cdot {{\bf \vec{\boldsymbol\psi }}}{\rm{d}}\Omega &=& - i\omega \int_{\Omega }{{\left( {\sigma - {\sigma ^p}} \right){{{{\bf \vec{E}}}}^p}}} \cdot {{\bf \vec{\boldsymbol\psi }}}{\rm{d}}\Omega + {\omega ^2}\int_{\Omega }{{\varepsilon {{{{\bf \vec{E}}}}^p}}} \cdot {{\bf \vec{\boldsymbol\psi }}}{\rm{d}}\Omega \,\,, \end{eqnarray} (53)where Ω is the computational domain (and dΩ = dxdydz), the test vector function $${{\bf \vec{\boldsymbol\psi }}}$$ and the unknown vector function $${{{\bf \vec{E}}}^s}$$ are the members of the space Hcurl(Ω) of the vector functions $${{\bf \vec{\Phi }}}$$ with continuous tangential components, which satisfy conditions $$\int_{\Omega }{{{{\bf \vec{\Phi }}} \cdot {{\bf \vec{\Phi }}}{\rm{d}}\Omega }} < \infty $$ and $$\int_{\Omega }{{( {\nabla \times {{\bf \vec{\Phi }}}} ) \cdot ( {\nabla \times {{\bf \vec{\Phi }}}} ){\rm{d}}\Omega }} < \infty $$. The function $${{{\bf \vec{E}}}^s}$$ is expressed as \begin{equation} {{{\bf \vec{E}}}^s} = \sum\limits_{j = 1}^n {q_j^s{{\vec{\boldsymbol\psi }}_j}} , \end{equation} (54)where $${\vec{\boldsymbol\psi }_j}$$ are the basis vector functions, n is their number (which coincides with the number of the edges in the 3-D mesh). As a result of substituting (54) into (53) and replacing the test function $${{\bf \vec{\boldsymbol\psi }}}$$ with the basis functions $${\vec{\boldsymbol\psi }_l}$$ in turn, we obtain the following SLAE: \begin{equation} \left( {{{\bf G}} + i\omega {{{\bf M}}^\sigma } - {\omega ^2}{{{\bf M}}^\varepsilon }} \right){{{\bf q}}^s} = {{\bf f}}, \end{equation} (55)where the components of matrices G, Mσ and Mε, and the vector f are defined by relations \begin{equation} {G_{lj}} = \frac{1}{{{\mu _0}}}\int_{\Omega }{}\left( {\nabla \times {{\vec{\boldsymbol\psi }}_l}} \right) \cdot \left( {\nabla \times {{\vec{\boldsymbol\psi }}_j}} \right){\rm{d}}\Omega ,\quad M_{lj}^\sigma = \int_{\Omega }{}\sigma {\vec{\boldsymbol\psi }_l} \cdot {\vec{\boldsymbol\psi }_j}{\rm{d}}\Omega ,\quad M_{lj}^\varepsilon = \int_{\Omega }{}\varepsilon {\vec{\boldsymbol\psi }_l} \cdot {\vec{\boldsymbol\psi }_j}{\rm{d}}\Omega , \end{equation} (56) \begin{equation} {f_l} = - i\omega \int_{\Omega }{}\left( {\sigma - {\sigma ^p}} \right){{{\bf \vec{E}}}^p} \cdot {\vec{\boldsymbol\psi }_l}{\rm{d}}\Omega + {\omega ^2}\int_{\Omega }{{\varepsilon {{{{\bf \vec{E}}}}^p}}} \cdot {\vec{\boldsymbol\psi }_l}{\rm{d}}\Omega . \end{equation} (57)The function $${{{\bf \vec{E}}}^p}$$ is expressed as a linear combination of the basis functions $${\vec{\boldsymbol\psi }_j}$$ \begin{equation} {{{\bf \vec{E}}}^p} = \sum\limits_{j = 1}^n {q_j^p{{\vec{\boldsymbol\psi }}_j}} , \end{equation} (58)where expansion coefficients $$q_j^p$$ in basis $$\{ {{{\vec{\boldsymbol\psi }}_j}} \}$$ are the values of tangential components of $${{{\bf \vec{E}}}^p}$$ on the FEs’ edges. In this case, the vector f in the right-hand side of FE equation system (55) can be calculated by the formula \begin{equation} {{\bf f}} = \left( { - i\omega {{{{\bf \bar{M}}}}^\sigma } + {\omega ^2}{{{\bf M}}^\varepsilon }} \right){{\bf q}}_{}^p, \end{equation} (59)where components of the matrix $${{{\bf \bar{M}}}^\sigma }$$ are defined by the formula \begin{equation} \bar{M}_{lj}^\sigma = \int_{\Omega }{}\left( {\sigma - {\sigma ^p}} \right){\vec{\boldsymbol\psi }_l} \cdot {\vec{\boldsymbol\psi }_j}{\rm{d}}\Omega . \end{equation} (60) 4 SPECIAL ASPECTS OF IMPLEMENTATION For the solution of axisymmetrical problems (when we calculate the primary field $${{{\bf \vec{E}}}^p}$$), we use the regular rectangular meshes with smoothly increasing of the step when moving away from the source, which is located on the z-axis. The step in the vicinity of the source is set about 0.01 m, the factor χ of increase of the mesh step when moving away from the source in the coordinates R and Z is set in the range of 1.05–1.08 (i.e. the value of the next mesh step hnext in each of the coordinates is defined by the value of the previous step hprev in the corresponding coordinate as hnext = χ ⋅ hprev). The approach proposed in the previous sections for modelling the primary field requires much lower computational costs in comparison with computational costs required for calculating the 3-D secondary field. Therefore, we use the special algorithms of optimization only for the 3-D meshes. These algorithms make it possible to reduce the computational time required for calculating the secondary field by approximately an order in comparison with the calculations on regular meshes (which are generated according to the principle considered above for 2-D meshes). Therefore, in this Section, we consider the algorithms of generating the optimized irregular non-conforming 3-D meshes and principles of constructing the FE approximation on them, as well as another aspects (grouping of logging tool positions and automatic choice of the background medium), which allow additional significant reduction of the computational costs. Note that we do not use the approaches with adaptive mesh generation, which allow the convergence towards the solution using the sequence of meshes. These approaches are based on various error estimators (Schwarzbach et al.2011; Grayver & Kolev 2015; Yin et al.2016, etc.) and require the additional costs caused by the calculations on several meshes in order to obtain the adaptive mesh and for each transmitter position. At the same time, these approaches do not always guarantee the required accuracy when the signals are calculated in the receivers. This fact is implicitly confirmed by Grayver & Kolev (2015) where, after generating the adaptive mesh (and satisfying the corresponding criteria), the additional refinements are performed in the vicinities of the receivers. 4.1 Generation of the optimized finite-element meshes For the FE solution, we use the optimized meshes with local refinements. Note that the trajectory length in the logging problems is usually much bigger than the size of the computational domain, in which the field for one position of the source can be computed with the required accuracy. For increasing the numerical solution accuracy it is necessary to refine the mesh in the vicinity of the transmitters and receivers. Therefore, during the calculation of the signals received by the logging tool when it is moving along the trajectory, it is necessary to change not only the positions of the local refinements of the mesh but also the computational domain itself (the analogous approach has been considered in (Yang et al.2014), but it uses regular meshes). The accuracy of approximation of the computational domain geometry is the important aspect for obtaining the high accuracy solution. On the one hand, this approximation should be quite good in the vicinity of the logging tool (i.e. the error caused by the inaccurate boundaries approximation should not exceed the error of the numerical solution approximation). On the other hand, the excessive mesh refinement caused by good approximation of the boundaries, which are distant from the logging tool, can lead to the significant increase in the number of FEs and, consequently, computational costs. The numerical results have shown that high accuracy geometry approximation is required only in some subdomain around the logging tool, and rough geometry approximation can be used at some distance from it. Note that though we use the irregular non-conforming meshes, we do not use FEs with curvilinear or inclined boundaries. This is related to the fact that, if we use the FEs with curvilinear or inclined boundaries, which smoothly approximate the corresponding boundaries of the domain, the application of the primary-secondary field approach can lead to the intersection of FEs of the 3-D mesh by the layers interfaces of the background medium. In this case, the secondary field $${{{\bf \vec{E}}}^s}$$ should be discontinuous inside those FEs which are intersected by the interfaces of background medium layers with different conductivity (because $${{{\bf \vec{E}}}^t} = {{{\bf \vec{E}}}^p} + {{{\bf \vec{E}}}^s}$$ inside the FE should be continuous, but $${{{\bf \vec{E}}}^p}$$ will be discontinuous if the FE is intersected by the layers with a different conductivity). Therefore, we use another approach (stair-like approximation) in which the FEs are never intersected by the layers interfaces of the background medium. Mesh generation algorithm consists of four main stages. At the first stage, the size and structure of the computational domain are defined depending on the logging tool position (or several positions; it will be described in Section 4.3). The second stage is the stage of generating the so-called base mesh, which does not include inner borders of the computational domain and is refined in the subdomain containing the logging tool. At the third stage, the geoelectrical inhomogeneities are embedded into the mesh taking into account the requirements to the accuracy of their boundaries approximation depending on the distance from the logging tool. And, finally, at the fourth stage, the local refinements in the vicinity of receivers and transmitters are performed. Thus, when generating the base mesh, we select the parallelepiped subdomain, which contains the transmitters and receivers locations for one tool position (or several tool positions included in one group). Inside this subdomain, the base mesh is generated with a constant step in all directions, and outside this subdomain, the cells size increases smoothly with moving away from the selected subdomain boundaries. Thus, the base mesh for a group of tool positions on the trajectory is defined by three parameters. They are the following: the step inside the selected subdomain, the distance from the tool to boundaries of the computational domain, and the factor of increase of the cell size (when moving away from the selected subdomain boundaries). Based on multiple numerical experiments performed for different geoelectrical models, we have obtained the following values of these parameters. For the frequency of 100 kHz, which we take as a ‘basic’ one, the size of the minimal step in the base mesh is set equal to 0.5 m, and the distance from the tool to the boundaries of the computational domain is set equal to 20 m in all directions. For other frequencies, these sizes are scaled by multiplying by the square root of the frequencies ratio. The factor of increase of the cells size (when moving away from the boundaries of the subdomain containing all locations of transmitters and receivers included in a group) is taken equal to 1.2. Before the third stage of the algorithm, all geoelectrical inhomogeneities included in the computational domain (or the whole computational domain) are approximated by parallelepipeds the size of which is defined by the step of geometry approximation. The step defining the approximation of geometry of embedded geoelectrical inhomogeneities in the vicinity of receivers’ locations is about 0.2 m for the problems without a borehole and invasion zone, and it is about 0.05 m for the problems in which the borehole and invasion zone need to be taken into account. Outside the vicinity of receivers’ locations, the step of boundaries approximation increases proportionally to the increase in the step size of the base mesh. The local refinements in the vicinity of receivers (at the fourth stage of the algorithm) are defined by the following parameters: the number of refinements inside the FEs into which the receivers get, the number of neighbouring levels of FEs, and the number of refinements for each neighbouring levels. For the examples considered in this paper, the values of these parameters are taken equal to 2, 3 and 1, respectively. As mentioned above, for solving the problems considered here, we use the mathematical models based on the representation of the total field as a sum of the primary (corresponding to the field induced by the source in the horizontally layered background medium) and secondary (which is the impact field of 3-D and 2-D geoelectrical inhomogeneities) components. Note that the local refinements in the vicinity of the transmitter should be performed only in the case if the transmitter is located in the subdomain (of the computational domain) in which σp ≠ σ. As an example, consider the geoelectrical model the vertical section of which is shown in Fig. 3. The logging tool trajectory is shown by black solid line. Figure 3. View largeDownload slide The geoelectrical model specific for induction logging problems. The model contains five layers (interfaces of which are slightly sloped) and the fault. The layers of equivalent conductivity are shown by the same colour. The conductivity of the middle layer is anisotropic: lateral conductivity is 0.256 Sm m−1 and vertical conductivity is 0.025 Sm m−1. Figure 3. View largeDownload slide The geoelectrical model specific for induction logging problems. The model contains five layers (interfaces of which are slightly sloped) and the fault. The layers of equivalent conductivity are shown by the same colour. The conductivity of the middle layer is anisotropic: lateral conductivity is 0.256 Sm m−1 and vertical conductivity is 0.025 Sm m−1. We consider the logging tool that includes one transmitter and two receivers the offsets of which are 4.5 and 9.5 m. The distance between two adjacent logging tool positions is 1 m. The value of the frequency is equal to 100 kHz. Fig. 4(a) shows the base mesh, which is obtained after the first and second stages of the algorithm, for one group of three neighbouring logging tool positions on the trajectory. Fig. 4(b) shows the mesh obtained after the third stage of the algorithm, that is, after embedding 2-D geoelectrical inhomogeneities into the mesh. Fig. 5(a) shows the mesh obtained after the fourth stage of the algorithm, that is, after local refinement in the vicinity of the receivers and transmitters. Fig. 5(b) shows the 3-D view of the sectioned FE mesh. In Fig. 6, the FE mesh is shown for other logging tool positions on the trajectory. Figure 4. View largeDownload slide (a) The vertical section of the base mesh obtained after the first and second algorithm stages and (b) the vertical section of the mesh obtained after embedding 2-D geoelectrical inhomogeneities into the mesh. Ti denotes the location of the ith transmitter; Ri,j denotes the location of the jth receiver of the ith transmitter. Figure 4. View largeDownload slide (a) The vertical section of the base mesh obtained after the first and second algorithm stages and (b) the vertical section of the mesh obtained after embedding 2-D geoelectrical inhomogeneities into the mesh. Ti denotes the location of the ith transmitter; Ri,j denotes the location of the jth receiver of the ith transmitter. Figure 5. View largeDownload slide (a) The vertical section of the mesh obtained after local refinements in the vicinities of receivers and transmitters and (b) 3-D view of the sectioned finite-element mesh. Ti denotes the location of the ith transmitter; Ri,j denotes the location of the jth receiver of the ith transmitter. Figure 5. View largeDownload slide (a) The vertical section of the mesh obtained after local refinements in the vicinities of receivers and transmitters and (b) 3-D view of the sectioned finite-element mesh. Ti denotes the location of the ith transmitter; Ri,j denotes the location of the jth receiver of the ith transmitter. Figure 6. View largeDownload slide The vertical sections of the meshes generated for different groups of the adjacent tool positions. Ti denotes the location of the ith transmitter; Ri,j denotes the location of the jth receiver of the jth transmitter. Figure 6. View largeDownload slide The vertical sections of the meshes generated for different groups of the adjacent tool positions. Ti denotes the location of the ith transmitter; Ri,j denotes the location of the jth receiver of the jth transmitter. 4.2 Finite-element approximation on the meshes with hanging nodes 4.2.1 ‘Non-conforming’ and ‘conforming’ bases on the meshes with hanging nodes The meshes presented in the previous section contain the non-conforming FEs. On these meshes, we introduce two sets of the basis functions: $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$ and $$\{ {\vec{\boldsymbol\psi }_k^c} \}$$, where the superscript ‘nc’ denotes non-conforming, and superscript ‘c’ denotes conforming. The basis functions $$\vec{\boldsymbol\psi }_k^{nc}$$ are associated with all edges of the mesh with hanging nodes. On each FE Ωe containing the corresponding (kth) edge, $$\vec{\boldsymbol\psi }_k^{nc}$$ is defined in the standard way: at this edge, the basis function $$\vec{\boldsymbol\psi }_k^{nc}$$ has the only non-zero (equal to 1) component, which corresponds to the coordinate (x, y or z) coinciding with the edge direction, and, on Ωe, this component decreases to zero as a bilinear function of two other coordinates when coming to three other edges with the same direction. Obviously, the functions $$\vec{\boldsymbol\psi }_k^{nc} \notin {{{\bf H}}^{{\rm{curl}}}}( \Omega )$$ occur in the basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$. As an example, consider the fragment of the irregular non-conforming mesh (with hanging nodes), which is shown in Fig. 7. Denote the basis function, which is associated with the edge connecting nodes with numbers 5 and 26, as $$\vec{\boldsymbol\psi }_5^{nc}$$. This function is non-zero on the FEs Ω1, Ω2, and Ω3, and zero on all other FEs. Further, we will denote the edges, FEs, and their faces by the sets of nodes, which are listed in braces. Thus, the basis function $$\vec{\boldsymbol\psi }_5^{nc}$$ is associated with the edge {5, 26}, FE Ω1 = {1, 2, 4, 5, 22, 23, 25, 26}, FE Ω1 has the faces {1, 2, 4, 5}, {1, 2, 22, 23}, etc. As an example, consider the basis functions $$\vec{\boldsymbol\psi }_{36}^{nc}$$, $$\vec{\boldsymbol\psi }_{37}^{nc}$$ and $$\vec{\boldsymbol\psi }_{38}^{nc}$$, which are associated with the edges {5, 10}, {10, 16} and {16, 26}, respectively. These edges are the sub-edges of the edge {5, 26}. The basis function $$\vec{\boldsymbol\psi }_{36}^{nc}$$ is non-zero only on FE Ω4 = {5, 6, 8, 9, 10, 12, 13, 15}, basis function $$\vec{\boldsymbol\psi }_{37}^{nc}$$ is non-zero only on FE Ω5 = {10, 11, 13, 14, 16, 17, 19, 20}, basis function $$\vec{\boldsymbol\psi }_{38}^{nc}$$ is non-zero only Ω7 = {16, 17, 19, 20, 26, 27, 30, 31}, and on all other FEs, these basis functions are equal to zero. Figure 7. View largeDownload slide The fragment of the irregular non-conforming finite-element mesh with hanging nodes. Figure 7. View largeDownload slide The fragment of the irregular non-conforming finite-element mesh with hanging nodes. Besides the basis functions $$\vec{\boldsymbol\psi }_k^{nc}$$ associated with the sub-edges of ‘big’ edges (which are the edges of other FEs), in the ‘non-conforming’ basis, there are functions, which are associated with the FEs’ edges lying within the faces of other FEs (naturally, such edges do not coincide with the boundaries of these faces). In Fig. 7, these are the functions $$\vec{\boldsymbol\psi }_{39}^{nc}$$ and $$\vec{\boldsymbol\psi }_{40}^{nc}$$, which are associated with the edges {11, 17} and {17, 27}, respectively. The function $$\vec{\boldsymbol\psi }_{39}^{nc}$$ is non-zero only on FE Ω5 (see above its nodes) and Ω6 = {11, 12, 14, 15, 17, 18, 20, 21}, and the function $$\vec{\boldsymbol\psi }_{40}^{nc}$$ is non-zero only on FE Ω7 and Ω8 = {17, 18, 20, 21, 27, 28, 31, 32}. In general, by the basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$ we mean the basis functions associated with all edges (including sub-edges). These basis functions are defined in the standard form on the FEs including these edges. On the non-conforming FE mesh with hanging nodes, we introduce the basis functions $$\vec{\boldsymbol\psi }_j^c$$ with continuous tangential components in the computational domain Ω, that is, $$\vec{\boldsymbol\psi }_j^c \in {{{\bf H}}^{{\rm{curl}}}}( \Omega )$$. These basis functions are associated only with the ‘main’ edges of the mesh. By the ‘main’ edges we mean such edges that are not the sub-edges and do not lie within the faces of any FEs. Define the basis functions $$\vec{\boldsymbol\psi }_j^c$$ in the form of linear combination of the basis functions $$\vec{\boldsymbol\psi }_k^{nc}$$ \begin{equation} \vec{\boldsymbol\psi }_j^c = \sum\limits_k {{T_{jk}}\vec{\boldsymbol\psi }_k^{nc}} , \end{equation} (61)that is, the basis $$\{ {\vec{\boldsymbol\psi }_j^c} \}$$ is constructed from the basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$ with the help of the transition matrix T = {Tjk}. The size of this matrix is nc × n, where n is the total number of the edges in the mesh (including sub-edges and edges, which lie within the faces of any FEs), and nc is the number of ‘main’ edges only. Thus, specifically with the help of the transition matrix T, we construct the basis $$\{ {\vec{\boldsymbol\psi }_j^c} \}$$ in which all vector functions have continuous tangential components in Ω (i.e. ∀j$$\vec{\boldsymbol\psi }_j^c \in {{{\bf H}}^{{\rm{curl}}}}( \Omega )$$). We note a very important aspect related to the ‘non-conforming’ basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$. Besides the basis functions with the discontinuous tangential components (which we have called ‘non-conforming’), the ‘conforming’ basis functions (with continuous tangential components) are also included into the basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$. They are associated with the edges, which we will name ‘regular’. These edges do not contain subedges, are not subedges themselves, and do not lie within any face. For example, in Fig. 7, the ‘regular’ edges are the edges {1, 22}, {4, 25}, {2, 23}, etc. (if we assume that either the faces of FE Ω1 lie on the computational domain boundary, or FE Ω1 and FEs adjacent to Ω1, which are not shown in this mesh fragment, are conforming). The superscript ‘nc’ added to these basis functions (associated with ‘regular’ edges) denotes only that they are included in the basis in which there are vector functions with discontinuous tangential components. Thus, the superscript ‘nc’ added to the function from the ‘non-conforming’ basis does not always denote that this function necessarily have discontinuous tangential components, that is, among the functions from $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$, there are functions with continuous tangential components, which are associated with the ‘regular’ edges. Therefore, further, the superscript ‘nc’ added to the basis function $$\vec{\boldsymbol\psi }_k^{nc}$$ will not rigorously denote that this function necessarily has discontinuous tangential components (i.e. the term ‘conforming’ function $$\vec{\boldsymbol\psi }_k^{nc}$$ is applicable and denotes that such a function from the basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$ is associated with the ‘regular’ edge). However, all ‘regular’ edges are the ‘main’, but unlike ‘regular’, the ‘main’ edges can contain the sub-edges (which are the edges of some FEs). For example, in Fig. 7, the edge {5, 26} is ‘main’, but not ‘regular’, because it contains subedges {5, 10} (the edge of FE Ω4), {10, 16} (the edge of FE Ω5), and {16, 26} (the edge of FE Ω7). 4.2.2 The construction of the ‘conforming’ basis from the ‘non-conforming’ Consider how the ‘conforming’ basis function $$\vec{\boldsymbol\psi }_j^c$$ can be constructed from the ‘non-conforming’ functions $$\vec{\boldsymbol\psi }_k^{nc}$$. In so doing, we really develop the algorithm of generating the matrix T. Let the face of one FE (this face will be called ‘big’) be adjacent to the several faces of other FEs (these faces will be called ‘small’). Then, in order to ensure on this ‘big’ face the continuity of the tangential component of the basis function $$\vec{\boldsymbol\psi }_j^c$$ associated with one of its edges, it is necessary to define the weights Tjk in the expansion (61) in basis functions $$\vec{\boldsymbol\psi }_k^{nc}$$. The weights Tjk are equal to the values of the tangential component $$\vec{\boldsymbol\psi }_j^c$$ on the corresponding edges of the ‘small’ faces. We illustrate this by the example of the mesh shown in Fig. 7. Consider the face {5, 6, 26, 28} of FE Ω2 (this is the ‘big’ face) adjacent to the faces {5, 6, 10, 12} of FE Ω4, {10, 11, 16, 17} of FE Ω5, {11, 12, 17, 18} of FE Ω6, {16, 17, 26, 27} of FE Ω7 and {17, 18, 27, 28} of FE Ω8 (these are the ‘small’ faces). Let, as before, the vector function associated with the edge {5, 26} have number 5 in the ‘non-conforming’ basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$, and vector functions associated with the edges {5, 10}, {10, 16}, {16, 26}, {11, 17} and {17, 27} have numbers 36, 37, 38, 39 and 40, respectively. In the ‘conforming’ basis $$\{ {\vec{\boldsymbol\psi }_j^c} \}$$, the vector function associated with the edge {5, 26} has number 5 too. Then, the ‘conforming’ vector function associated with the edge {5, 26} takes the form \begin{equation} \vec{\boldsymbol\psi }_5^c = \vec{\boldsymbol\psi }_5^{nc} + \vec{\boldsymbol\psi }_{36}^{nc} + \vec{\boldsymbol\psi }_{37}^{nc} + \vec{\boldsymbol\psi }_{38}^{nc} + \frac{{{h_{11,12}}}}{{{h_{10,12}}}}\left( {\vec{\boldsymbol\psi }_{39}^{nc} + \vec{\boldsymbol\psi }_{40}^{nc}} \right), \end{equation} (62)where hν, η is the length of the edge {ν,η}. Coefficients T5, 5 = 1, T5, 36 = 1, T5, 37 = 1, T5, 38 = 1, $${T_{5,39}} = {{{h_{11,12}}} / {{h_{10,12}}}}$$, and $${T_{5,40}} = {{{h_{11,12}}} / {{h_{10,12}}}}$$ are the values of the function $$\vec{\boldsymbol\psi }_5^c$$ at the edges {5, 10}, {10, 16}, {16, 26}, {11, 17} and {17, 27}, which lie on the ‘big’ face {5, 6, 26, 28} and with which the functions $$\vec{\boldsymbol\psi }_{36}^{nc}$$, $$\vec{\boldsymbol\psi }_{37}^{nc}$$, $$\vec{\boldsymbol\psi }_{38}^{nc}$$, $$\vec{\boldsymbol\psi }_{39}^{nc}$$ and $$\vec{\boldsymbol\psi }_{40}^{nc}$$, are associated. Note that the function $$\vec{\boldsymbol\psi }_5^c$$ identically coincides with the function $$\vec{\boldsymbol\psi }_5^{nc}$$ on the ‘big’ face, because on FE Ω2, both $$\vec{\boldsymbol\psi }_5^c$$ and $$\vec{\boldsymbol\psi }_5^{nc}$$ are defined by the same shape function. We also note that formula (62) defines the function $$\vec{\boldsymbol\psi }_5^c$$ only inside the FEs, and its values on the faces are defined as continuous extensions of its values from the FEs. On the face {5, 6, 26, 28} itself, formula (62) should be interpreted only in the sense of continuous extensions of the basis functions $$\vec{\boldsymbol\psi }_k^{nc}$$ onto this face from the FEs, which contain the edges corresponding to $$\vec{\boldsymbol\psi }_k^{nc}$$, because all functions $$\vec{\boldsymbol\psi }_k^{nc}$$ in the right-hand side of (62) are discontinuous on this face. As a result, we really obtain the function $$\vec{\boldsymbol\psi }_5^c$$ with continuous tangential components on the face {5, 6, 26, 28}. Thus, we propose the following algorithm of calculating the matrix T, which makes it possible to obtain the ‘conforming’ basis functions $$\vec{\boldsymbol\psi }_j^c$$ from the ‘non-conforming’ (but standard for the vector FEM) basis functions $$\vec{\boldsymbol\psi }_k^{nc}$$. To begin with, we assume that in the ‘non-conforming’ basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$, the functions associated with the ‘main’ edges are numbered first. The same numbering of the ‘main’ edges is used for the functions from the ‘conforming’ basis $$\{ {\vec{\boldsymbol\psi }_j^c} \}$$. Then, the left block of the matrix T with the size nc × nc is the identity matrix (i.e. Tjj = 1, and Tjk = 0 if k ≠ j). This is convenient both for the storage of the matrix T (i.e. in fact, there is no need to store this block) and for assembling the matrices of FEs. Note that the basis functions $$\vec{\boldsymbol\psi }_j^c$$ associated with the ‘regular’ edges are identically equal to $$\vec{\boldsymbol\psi }_j^{nc}$$, and for these edges, all components of the corresponding (jth) rows in the right block of the matrix T are equal to zero (i.e. Tjk = 0 if k > nc). Non-zero components Tjk if k > nc are calculated in the following way. If the edge numbered k is the sub-edge of the edge numbered α ≤ nc, then Tαk = 1, and all other components of the kth column of the matrix T are zero (i.e. ∀β ≠ α Tβk = 0). If the edge numbered k lies within the ‘big’ face, on which the basis functions $$\vec{\boldsymbol\psi }_\alpha ^{nc}$$ and $$\vec{\boldsymbol\psi }_\beta ^{nc}$$ for α, β ≤ nc have non-zero component tangential to this edge, then only two components Tαk and Tβk are non-zero in the kth column of the matrix T, and their values are equal to the values of functions $$\vec{\boldsymbol\psi }_\alpha ^{nc}$$ and $$\vec{\boldsymbol\psi }_\beta ^{nc}$$ at the kth edge, respectively. Thus, for the mesh shown in Fig. 7, if the edge {11, 17}, which lies within the face {5, 6, 26, 28}, has number 39, and edges {5, 26} and {6, 28} have numbers 5 and 6, respectively, then only the components $${T_{5,39}} = {{{h_{11,12}}} / {{h_{10,12}}}}$$ and $${T_{6,39}} = {{{h_{10,11}}} / {{h_{10,12}}}}$$ (where, as before, hν, η is the length of the edge {ν,η}) are non-zero in the 39th column of the matrix T. Finally, we consider the situation when the edge numbered k > nc lies within the ‘big’ face and even if one of the edges of this ‘big’ face, which are parallel to the k-th edge, is not the ‘main’, that is, its number is α > nc. This situation is shown in Fig. 8. The edge {10, 16} lies within the face {4, 5, 24, 26} and, therefore, this face is ‘big’ over the edge {10, 16} and FEs {9, 10, 12, 13, 15, 16, 18, 19} and {10, 11, 13, 14, 16, 17, 19, 20}, which contain it. But the edge {4, 24} of this ‘big’ face lies within the face {2, 7, 22, 28}, which is ‘big’ over this edge and FE {2, 3, 4, 5, 22, 23, 24, 26}, which contains this edge. Therefore, the edge {4, 24} is not ‘main’ and its number is α > nc. In this case, the non-zero components of the kth column of the matrix T are calculated with the use of a ‘chain’. First, the values $${\tilde{T}_{\alpha k}}$$ are calculated as values of tangential components of the basis function $$\vec{\boldsymbol\psi }_\alpha ^{nc}$$ at the kth edge. For the example shown in Fig. 8 and edges {10, 16} (its number is equal to k) and {4, 24} (its number is equal to α), $${\tilde{T}_{\alpha k}}$$ is equal to $${{{h_{10,11}}} / {{h_{9,11}}}}$$. The tilde ‘∼’ over Tαk have been added because if α > nc, there are no such components in the matrix T (because its size is nc × n), and the values $${\tilde{T}_{\alpha k}}$$ if α > nc are really fictitious and used only in the ‘chain’ for calculating the non-fictitious components of the matrix T (i.e. components Tjk if j ≤ nc). Figure 8. View largeDownload slide The fragment of the irregular non-conforming finite-element mesh with hanging nodes such that the ‘chains’ appear when the transition matrix is generated. Figure 8. View largeDownload slide The fragment of the irregular non-conforming finite-element mesh with hanging nodes such that the ‘chains’ appear when the transition matrix is generated. Further, the ‘big’ face within which the edge numbered α lies is processed. For the edges of this face, which are parallel to the edge numbered α, the values of basis functions $$\vec{\boldsymbol\psi }_\gamma ^{nc}$$ and $$\vec{\boldsymbol\psi }_\nu ^{nc}$$ associated with these edges are calculated at the edge numbered α. If γ ≤ nc and ν ≤ nc, that is, the edges numbered γ and ν are the ‘main’ ones, then the components $${T_{\gamma k}} = {T_{\gamma \alpha }} \cdot {\tilde{T}_{\alpha k}}$$ and $${T_{\nu k}} = {T_{\nu \alpha }} \cdot {\tilde{T}_{\alpha k}}$$ are non-zero components in the kth column of the matrix T. Thus, for the mesh shown in Fig. 8, if the edges {10, 16}, {2, 22} and {7, 28} are numbered k, γ and ν, respectively, then $${T_{\gamma ,k}} = \frac{{{h_{10,11}}}}{{{h_{9,11}}}} \cdot \frac{{{h_{4,7}}}}{{{h_{2,7}}}}$$ and $${T_{\nu ,k}} = \frac{{{h_{10,11}}}}{{{h_{9,11}}}} \cdot \frac{{{h_{2,4}}}}{{{h_{2,7}}}}$$. In so doing, the generation of the kth column of the matrix T is completed. If at least one of the edges is not the ‘main’ one (i.e. γ > nc and/or ν > nc), then for this edge (or these two edges), instead of Tγα and/or Tνα the ‘fictitious’ values $${\tilde{T}_{\gamma \alpha }}$$ and/or $${\tilde{T}_{\nu \alpha }}$$ are calculated, and the corresponding ‘chains’ (or one of them) must be continued in a similar manner. Thus, the generation of the kth column of the matrix T will be completed when all ‘chains’, which have already been begun, are finished (i.e. the ‘main’ edge should be at the end of each ‘chain’). Note that a number of non-zero components of the matrix T (even without taking into account the identity submatrix with the size of nc × nc) is small. Therefore, for its storage it is reasonable to use the sparse format as for the global matrix of the FE equations system. The use of the matrix T allows sufficiently convenient generation of matrices G, Mσ, and Mε of FE equations system (55) for the ‘conforming’ basis $$\{ {\vec{\boldsymbol\psi }_j^c} \}$$ by assembling the standard FE matrices corresponding to the ‘non-conforming’ basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$. 4.2.3 Generation of global matrices of finite-element equations system We find the solution $${{{\bf \vec{E}}}^s}$$ of eq. (53) as a linear combination of basis functions $$\vec{\boldsymbol\psi }_j^c$$: \begin{equation*} {{{\bf \vec{E}}}^s} = \sum\limits_{j = 1}^{{n_c}} {q_j^s\vec{\boldsymbol\psi }_j^c} . \end{equation*} Note once again that basis functions $$\vec{\boldsymbol\psi }_j^c$$ have continuous tangential components in the whole computational domain Ω ($$\vec{\boldsymbol\psi }_j^c \in {{{\bf H}}^{{\rm{curl}}}}( \Omega )$$) and, therefore, are suitable for variational equation (53). We show how the matrices G, Mσ, and Mε, which define the matrix of FE equations system (55), can be calculated using the FE matrices, the components of which are obtained with the use of the basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$. The basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$ is defined by the standard shape vector functions on FEs, including the ‘non-conforming’ ones, that is, on each FE the basis functions $$\vec{\boldsymbol\psi }_k^{nc}$$, which are non-zero on it, are the shape functions. Thus, shape functions are the local basis functions (on FE) of the basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$. We obtain the expressions for calculating the components of the matrix G (see eq. 56) with the use of the basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$. Considering that tangential components of basis functions $$\vec{\boldsymbol\psi }_k^{nc}$$ can be discontinuous on the faces of the FEs Ωe, before turning to the basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$, we represent the integral over $$\Omega = \bigcup\limits_e {{\Omega _e}} $$ in (56) as a sum of integrals over FEs Ωe, and after that we use expression (61): \begin{eqnarray*} {G_{lj}} &=& \frac{1}{{{\mu _0}}}\int_{\Omega }{}\left( {\nabla \times \vec{\boldsymbol\psi }_l^c} \right) \cdot \left( {\nabla \times \vec{\boldsymbol\psi }_j^c} \right){\rm{d}}\Omega = \frac{1}{{{\mu _0}}}\sum\limits_e {\int_{{{\Omega _e}}}{}\left( {\nabla \times \vec{\boldsymbol\psi }_l^c} \right) \cdot \left( {\nabla \times \vec{\boldsymbol\psi }_j^c} \right){\rm{d}}\Omega }\nonumber \\ &=& \frac{1}{{{\mu _0}}}\sum\limits_e {\int_{{{\Omega _e}}}{}\left( {\nabla \times \sum\limits_k {{T_{lk}}\vec{\boldsymbol\psi }_k^{nc}} } \right) \cdot \left( {\nabla \times \sum\limits_m {{T_{jm}}\vec{\boldsymbol\psi }_m^{nc}} } \right){\rm{d}}\Omega}. \end{eqnarray*} Now the summation over k and m can be taken off the sign of integration: \begin{equation} {G_{lj}} = \frac{1}{{{\mu _0}}}\sum\limits_e {\sum\limits_k {\sum\limits_m {{T_{lk}}{T_{jm}}} } \int_{{{\Omega _e}}}{}\left( {\nabla \times \vec{\boldsymbol\psi }_k^{nc}} \right) \cdot \left( {\nabla \times \vec{\boldsymbol\psi }_m^{nc}} \right){\rm{d}}\Omega } . \end{equation} (63)Obviously, the integral in (63) defines the matrix $${{{\bf \hat{G}}}^e}$$ of the FE Ωe for basis functions $$\vec{\boldsymbol\psi }_k^{nc}$$ which are the standard shape vector functions on Ωe. Thus, formula (63) defines the rule of assembling matrices $${{{\bf \hat{G}}}^e}$$ when the ‘conforming’ basis functions $$\vec{\boldsymbol\psi }_j^c$$ are used on the non-conforming meshes, and shape functions on FEs (and, consequently, their matrices) are defined in a standard way for the vector FEM. This rule is as follows. The numbers k and m of global basis functions from the basis $$\{ {\vec{\boldsymbol\psi }_i^{nc}} \}$$ are defined by the numbers of local basis functions, and the considered component of the local matrix is added with the factor TlkTjm only to those components Glj of the global matrix G, for which both components Tlk and Tjm are non-zero. The generation of matrices Mσ and Mε of FE equation system (56) is performed similarly, that is, with the use of standard shape vector functions, which correspond to the ‘non-conforming’ basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$. Because the basis functions $$\vec{\boldsymbol\psi }_j^c$$ are the elements of the space Hcurl(Ω), all theoretical estimations of the convergence of FE solutions obtained with the use of these functions (Monk 2003) are true. If we perform the local mesh refinements in the needed subdomains (in the vicinity of the receivers and anomalies with significant influence), FE solutions obtained on such non-conforming meshes have essentially the same accuracy as FE solutions obtained on regular meshes, which correspond to these irregular non-conforming meshes. However, when such irregular meshes are used, the computational costs are lowered by an order or more (frequently, even more than 100 times) in comparison with the use of corresponding (i.e. with the cells refinement in the same subdomains) regular meshes. This will be shown below, in Section 5. 4.3 The grouping of logging tool positions and using the direct solver The application of the primary-secondary field approach with a separate calculation of the horizontally layered medium field and application of the unstructured optimized meshes with local refinements of the cells for calculating the secondary field allow us to significantly reduce the size of matrix of the FE system. In this case, the direct solvers become sufficiently effective. Because the solution of logging problems implies the computation of the large number of fields in the same medium for the different positions of transmitters and receivers, the application efficiency of the direct solvers can be highly increased. If the medium does not change and only the transmitter position changes, then the matrix of the FE system will remain the same for all these problems when the same mesh is used. In this case, only the right-hand side vectors of the FE system will be different (in accordance with the transmitter position). Thus, when the direct solvers for SLAE are used, the matrix factorization into triangular matrices can be performed only once for all transmitter positions, and when the field for each transmitter position is computed, one should solve only the systems with these triangular matrices (Da Silva et al.2012; Oldenburg et al.2013; Haber & Schwarzbach 2014). Note that in our code, we use Intel MKL PARDISO solver (Schenk & G̈artner 2004). However, simultaneous field calculation for all logging tool positions on the trajectory is inappropriate. First, the trajectories in the logging problems usually significantly exceed the size of the computational domain, in which the field for one transmitter position can be computed with a required accuracy. Second, in order to obtain the numerical solution with the required accuracy, it will suffice to refine the mesh only in a reasonably small vicinity of the tool position. Therefore, in order to achieve the required accuracy of the solution on the same mesh for several logging tool positions, it is required to refine the mesh in quite a large subdomain of the computational domain, and the number of nodes in the mesh (and, consequently, the SLAE size) can increase very much with the increase in the number of tool positions for which the same matrix of the FE system will be generated. As a result, the size of the FE system can increase so much that even the time of the solution of two systems with two triangular matrices (i.e. factors obtained due to the factorization) can exceed the total time of calculating the field for one tool position on the mesh with a small number of FEs (not to mention other computational costs: factorization of the big size matrix, recalculation of the primary field onto a fine mesh, and etc.). Thus, both the approach in which the field calculation is performed for each logging tool position on the individual small mesh and the approach in which the fields for all logging tool positions on the trajectory are calculated simultaneously (on the corresponding mesh with a great number of elements) are not optimal. Therefore, the most appropriate is the grouping of problems, which correspond to different tool positions on the trajectory, such that, on the one hand, there are more positions in one group, and, on the other hand, the size of the matrix of the FE system does not increase too much, because with the increase in the size of SLAE, the computational costs for matrix factorization (in order to construct the triangular factors) and other overhead costs grow very fast. In the case of increasing the number of positions in the group, the difference between the 3-D medium and the horizontally layered background medium taken for this group can also increase. In this case, in order to achieve the required accuracy, it is necessary to use a more fine mesh for calculating the secondary field that can lead to higher computational costs. We also note that the field calculation for groups can be performed in parallel that makes it possible to achieve the additional acceleration. In our code, the grouping of tool positions is performed automatically. In this case, the number of positions in a group is essentially defined by the number of positions getting into the sequent segment of the trajectory. The length of the trajectory segment (as well as the step sizes in the mesh) depends on the frequency. Thus, for the frequency of 100 kHz, the optimal length of the trajectory segment varies from 2 m (in the case, if the distance between adjacent positions is about 0.1 m) to 6 m (in the case, if the distance between adjacent positions is about 1–2 m). For other frequencies, the length of the trajectory segment is scaled by multiplying by the square root of the frequencies ratio. Thus, the number of positions in one group varies on the average from 3 to 20. 4.4 Automatic choice of the horizontally layered background medium for calculating the primary field The highest accuracy of the FE solution can be achieved without increasing the computational costs in the situation when the logging tool (i.e. its transmitter and receivers) is entirely located in such a subdomain of the 3-D computational domain the conductivity of which coincides with the conductivity of the background medium (i.e. the horizontally layered medium taken as a background for calculating the primary field). Moreover, the structures of geoelectrical models and logging tool motion trajectories, which are specific for the problems considered here (see e.g. Fig. 3), as a rule, make practically impossible the choice (pre-setting) of the single (for the entire trajectory) background medium, which would satisfy the requirements mentioned above. Therefore, we propose to perform the automatic choice of the background medium for each logging tool position (or group of positions) on the trajectory taking into account the location of the tool and the properties of the medium in close vicinity of it. For calculating the field for one group of logging tool positions, we propose the following algorithm of the automatic choice of the background medium. In the 3-D FE mesh, we select the subdomain ΩT-R, which in plane contains all transmitters and receivers of the considered group and some vicinity surrounding them (i.e. this subdomain is limited only in x- and y-directions, but in the z-direction, it coincides with the original 3-D computational domain). Using all FEs, which are inside this subdomain, the array ZT-R of z-coordinates of these FEs boundaries is generated. This array is ordered, and the numbers of FEs, the z-coordinates of which are located inside the interval ($$Z_i^{{\rm{T - R}}}$$, $$Z_{i + 1}^{{\rm{T - R}}}$$), are assigned to this interval. This information is formed as a list by one looping over all FEs from ΩT-R. Then, by means of looping over intervals the boundaries of which are defined by the array ZT-R, the horizontally layered (background) medium is formed by the following algorithm. For the current interval, we define the maximal conductivity in this interval by analysing all FEs associated with this interval and taken from the corresponding list. If this conductivity coincides with the maximal conductivity obtained for the previous interval, then we pass to the next interval. If we obtain a different value of maximal conductivity, then the interface between the current and previous intervals becomes the interface of the horizontally layered (background) medium to be formed. As conductivities of this horizontally layered (background) medium we take maximal conductivity values, which are obtained during the determination of the z-coordinates of layers interfaces. Because the size of the 3-D computational domain is smaller than the size of the computational domain for calculating the primary field, outside ΩT-R, the background horizontally layered medium is formed using the data describing the original problem. At the end of the algorithm, we check the thickness of the layers. If the layers with a very small thickness are presented in the formed medium, then these layers are joined with the adjacent ones (from two adjacent layers we choose the layer which is thicker), that is, in fact, a small thickness layer is excluded from the medium, and the layer adjacent to it is properly extended. We note that the 3-D mesh, which is used for the automatic choice of the background medium, can partially describe the geometry of the 3-D problem to be solved. Before automatically choosing, it is better not to include small size objects into this mesh (e.g. a borehole if its effect on the measured signal is rather small). We include these objects into the 3-D mesh after the layers of the background horizontally layered medium have been formed, and their interfaces have been included into the mesh. Fig. 9 (right) shows the sections of the background mediums, which are formed with the use of the algorithm considered above for three groups of tool positions, including the situation when the trajectory crosses the fault. Figure 9. View largeDownload slide The vertical sections of meshes (left) generated for different groups of adjacent tool positions for the case when the background medium (right) is chosen automatically. Ti denotes the location of the ith transmitter; Ri,j denotes the location of the jth receiver of the ith transmitter. Figure 9. View largeDownload slide The vertical sections of meshes (left) generated for different groups of adjacent tool positions for the case when the background medium (right) is chosen automatically. Ti denotes the location of the ith transmitter; Ri,j denotes the location of the jth receiver of the ith transmitter. Fig. 9 shows the FE meshes generated for the same tool positions on the trajectory as shown in Figs 5 and 6. But in Fig. 9, these meshes have been generated for the case when the background medium is chosen automatically for each group of the tool positions. The comparison of the numerical results for the cases when the background medium is pre-set as single and when the background medium is chosen automatically (for each group of the adjacent tool positions) will be presented below, in Section 6. 5 VALIDATION In this section, we present the results of the verification of the mathematical models and numerical schemes proposed in this paper and described in the previous sections. We will call the field calculation performed with the use of mathematical models (3), (20), (28)–(30) the 2-D approximation because, although the medium is horizontally layered (i.e. really 1-D), in order to calculate the electromagnetic field induced by the current coil, we use the solution of one or several axisymmetrical (i.e. 2-D) problems in the 2-D computational domain in the (r, z)-coordinates. We will call the field calculation performed with the use of mathematical model (1) the 3-D approximation even in the case when the 2.5-D problem is solved, because in order to obtain its solution, we find the FE solution in the 3-D computational domain. In order to verify the method proposed in this paper for calculating the primary field in the horizontally layered medium, we use the Dipole1D code (Key 2009). Because the Dipole1D code allows the calculation of the field induced only by the electrical dipole, then in order to obtain the field induced by the current frame (loop), we sum the fields of four dipoles making up the current frame. As a testing model, we use the geoelectrical model shown in Fig. 10. Figure 10. View largeDownload slide Testing geoelectrical model 1 of the horizontally layered medium and locations of the lines (Line1…Line5) along which the non-zero components of $${{\bf \vec{B}}}$$ and $${{\bf \vec{E}}}$$ are calculated with the use of the proposed approach and Dipole1D code (Key 2009). Figure 10. View largeDownload slide Testing geoelectrical model 1 of the horizontally layered medium and locations of the lines (Line1…Line5) along which the non-zero components of $${{\bf \vec{B}}}$$ and $${{\bf \vec{E}}}$$ are calculated with the use of the proposed approach and Dipole1D code (Key 2009). The source in the form of a square frame with the side of 0.1 m is fixed, and its centre is located at the point (0, 0, 13). The current is 100 A. The frequency is 100 kHz. In order to verify the field values in the receivers, we will output the non-zero components of the magnetic induction $${{{\bf \vec{B}}}^p}$$ along Line 1 and Line 2 (see Fig. 10). Besides, as we have described in the previous sections, when calculating the 3-D field, the right-hand side of the equations (see eq. 1) is defined by the electrical field induced by the source in the horizontally layered medium, which has been taken as a background. Therefore, we will also verify the accuracy of calculating the non-zero components of the electrical field $${{{\bf \vec{E}}}^p}$$. We will analyse the accuracy of calculating $${{{\bf \vec{E}}}^p}$$ along Line 3, Line 4, and Line 5 (see Fig. 10), and also along Line 6 and Line 7, which are parallel to Line 3 and have the same z-coordinate. The distance along the y-axis between these lines (Line 6 and Line 7) and Line 3 is 2 m and 4 m, respectively. The plots placed at the top of Fig. 11 show the real and imaginary x- and z-components of $${{{\bf \vec{B}}}^p}$$ along Line 1 and Line 2 for the case when the field is induced by the ‘vertical’ frame and ‘horizontal’ loop. The plots placed at the bottom of Fig. 11 show the relative (in per cent) deviations of $${{{\bf \vec{B}}}^p}$$, which are calculated using the approach proposed in this paper, from the values of $${{{\bf \vec{B}}}^p}$$, which are calculated with the use of the Dipole1D code (Key 2009). The results presented in Fig. 11 show that within the distance of 10 m, the relative deviations for all components do not exceed 0.5 per cent, and within the distance of 20 m, the deviations do not exceed 1–1.5 per cent. The error growth takes place in the subintervals of low values of the field. We do not consider the field values farther than 20 m, because this value (20 m) defines the maximal size of the 3-D computational domain. Figure 11. View largeDownload slide $$B_x^{{\rm{Re}}}$$ (dark blue curve), $$B_x^{{\mathop{\rm Im}\nolimits} }$$ (green curve), $$B_z^{{\rm{Re}}}$$ (red curve) and $$B_z^{{\mathop{\rm Im}\nolimits} }$$ (light bl