3D induction log modelling with integral equation method and domain decomposition preconditioning

The deployment of electromagnetic (EM) induction tools while drilling is one of the standard routines for assisting the geosteering decision-making process. The conductivity distribution obtained through the inversion of the EM induction log can provide important information about the geological structure around the borehole. To image the 3D geological structure in the subsurface, 3D inversion of the EM induction log is required. Because the inversion process is mainly dependent on forward modelling, the use of fast and accurate forward modelling is essential. In this paper, we present an improved version of the integral equation (IE) based modelling technique for general anisotropic media with domain decomposition preconditioning. The discretised IE after domain decomposition equals a fixed-point equation that is solved iteratively with either the block Gauss-Seidel or Jacobi preconditioning. Within each iteration, the inverse of the block matrix is computed using a Krylov subspace method instead of a direct solver. An additional reduction in computational time is obtained by using an adaptive relative residual stopping criterion in the iterative solver. Numerical experiments show a maximum reduction in computational time of 35 per cent compared to solving the full-domain IE with a conventional GMRES solver. Additionally, the reduction of memory requirement for covering a large area of the induction tool sensitivity enables acceleration with limited GPU memory. Hence, we conclude that the domain decomposition method is improving the efficiency of the IE method by reducing the computation time and memory requirement.


INTRODUCTION
State-of-the-art tools for electromagnetic (EM) induction logging-while-drilling (LWD) enable real-time mapping of formation boundaries tens of meters away from the borehole (Sinha et al. 2022).These tools typically consist of multiple antenna configurations that have different sensitivities to the electrical resistivity distribution in the medium around the borehole.The distribution of the electrical properties is quantified through an inversion process and provides structural information and characteristics of the surrounding medium.The studies in real-time geosteering inversion usually employ 1D or 2D approximations (Bakr et al. 2017;Noh et al. 2022;Pardo & Torres-Verdín 2015;Puzyrev 2019).However, for imaging complex geological structures, it is important to capture the 3D variability of the resistivity change around the borehole through 3D inversion methods (Puzyrev et al. 2019;Sinha et al. 2022).The work of Wilson et al. (2019) shows that it is possible to perform 3D inversion in real-time, however, it is challenging due to the large computational cost required for the 3D forward modelling, especially when quantification of the uncertainties in the inversion is required.Therefore, the study of a fast 3D forward solver that accurately models induction logs remains essential for the development and testing of new imaging methods.
The integral equation (IE) method is one of the most widely applied numerical methods for the 3D modelling of EM data (Avdeev 2005;Wang et al. 2021) alongside the finite difference (Newman & Alumbaugh 2002;Hou et al. 2006) and finite element methods (Puzyrev et al. 2013;Ren et al. 2014).One of the main advantages of using the IE method is that it has the accuracy of a semi-analytical solution (Wang et al. 2021).Without introducing many specific approximations, the EM fields around the borehole are obtained by solving the linear system arising from the discretization of the integral equations.As the linear system is dense, the computational memory and time required can be large compared to other numerical methods (Yoon et al. 2016; Integral equation method with domain decomposition 3 Zaslavsky et al. 2011).To overcome this challenge, the linear system can be efficiently solved using an iterative solver based on the Krylov subspace method in combination with the utilization of the FFTs to accelerate the convolution integral operations in the linear system (Fang et al. 2006).A faster convergence rate can be achieved by implementing the contraction IE formulation (Hursan & Zhdanov 2002) which works especially well in the presence of a high contrast or a high degree of anisotropy.Additionally, the application of GPUs further decreases computation times because GPUs enable the acceleration of mathematical operations that can be straightforwardly parallelized (Dyatlov et al. 2015;Saputera et al. 2022).
In the work of Zhdanov et al. (2006), the formulation of the IE method is extended by decomposing the region of interest into several sub-domains.The field in the entire domain is obtained by sequentially solving the linear system in each sub-domain and updating the interaction between the sub-domains iteratively until convergence.With this formulation, it becomes feasible to conduct large-scale modelling of surface EM data in heterogeneous media as the computational operation can be reduced to one sub-domain at a time.It is possible to obtain an additional reduction in computational costs by only considering sub-domains that contain an anomaly with respect to the background medium.This leads to a smaller number of discretization blocks required for the 3D modelling while still enabling FFT implementation (Endo et al. 2009) and an improved iterative solver convergence rate (Van Dongen et al. 2007).Typically, a horizontally layered model is chosen as the background medium as the theory of Green's functions for layered one-dimensional (1D) models is very well developed (Zhdanov et al. 2006).Hence, the IE method can be very efficient when the resistivity model only deviates in some areas from the 1D model.However, in our application, the subsurface structure can vary in all directions.The sub-domains containing an anomaly can be everywhere around the EM tools and it may not be possible to achieve a reduction in the number of discretizations by the domain decomposition.Additionally, the sub-domains from the decomposition can be adjacent to each other such that the interactions between neighbouring sub-domains are not negligible.
The domain decomposition method can lead to an efficient way of solving the linear system of the IE method (Jakobsen & Tveit 2018;Wang et al. 2017).In the work of Jakobsen & Tveit (2018), the domain decomposition method is used to efficiently compute the T-matrix for the inversion of controlled source EM data.It is also shown that the domain decomposition method opens up the possibility to compute the T-matrix in parallel.
In this paper, we demonstrate that the formulation of IE with domain decomposition (IE-DD) can be interpreted as a preconditioned linear system, offering a computational advantage.We illustrate that the IE-DD method can be represented as a fixed-point equation, which is iteratively solved using block Gauss-Seidel or Jacobi preconditioners (Saad 2003).In particular, we will use a Krylov subspace method to invert the block matrices that are present in the formulation.Instead of expressing the decomposition formulation in terms of the contrast source in each sub-domain as described in Zhdanov et al. (2006) and Endo et al. (2009), we present the domain decomposition formulation in terms of the electric field in each sub-domain and a different perspective on the derivation of the IE-DD formulation.Additionally, we propose the use of an inexact iterative solver when solving the IE linear system for each sub-domain where the target tolerance is adapted based on the full domain residual.
The outline of this paper is described as follows.In section 2 called theory, we give an overview of the theory and implementation of the conventional IE method and the IE-DD.In section 3 called numerical results and discussion, we present three numerical examples to show the performance of the IE-DD method and discuss the computational aspect of our implementation.First, we show an example with isolated sub-domains to verify if the domain decomposition formulation will produce the same numerical results as the conventional full-domain formulation.In the second example, we show a numerical experiment with three different IE-DD schemes and compare the performance of these schemes with each other and the full-domain IE as a reference.In the last example, we simulate a logging scenario across a large complex 3D model.Furthermore, we showcase the ability of the domain decomposition method to reduce the memory requirement for dealing with a large number of grid blocks in the last example.This feature lets us cover more portion of the subsurface receivers while keeping a fine grid size, which may not be straightforward to implement in our currently available computer without the domain decomposition method.In section 4, we provide a compact evaluation of the IE-DD implementation in this study and also Integral equation method with domain decomposition 5 some possible improvements for future research.This paper contains appendices with more indepth details of the IE-DD derivation and implementation.We also include the comparison of our conventional IE code and existing code as a benchmark of our work in the appendix.

The Integral Equation Method for 3D Induction Logs Modelling
Maxwell's equations for heterogeneous media (Wannamaker & Zhdanov 2002) are the basic theory for modelling the induction tools' response within the frequency domain: where E (r) and H (r) are the total electric and magnetic fields, respectively, at location r, J H (r) denotes the magnetic source term, ω is the angular frequency, µ is the magnetic permeability, σ (r) = σ (r) − iωε (r) is the complex electric conductivity, ε is the dielectric permittivity, and i = √ −1.We assume that the magnetic permeability is constant and it is set equal to the magnetic permeability of the vacuum µ 0 .Additionally, the imaginary part of the complex conductivity can be ignored in the diffusion regime, which is a typical assumption for the operating conditions of induction tools.
The total electric and magnetic fields can be formulated using the following integral equations (Fang et al. 2006) where the Ω indicates the domain of integration where anomalies in the conductivity relative to the background conductivity are present.The integral terms in equations ( 3) and (4) represent the scattered electric and magnetic fields, respectively, due to the presence of these anomalies.The (0) superscripts indicate the fields defined for a homogeneous isotropic background medium with conductivity σ 0 which are referred to as the background fields.We choose a homogeneous isotropic background medium for simplicity and efficiency when calculating Green's tensor (Fang et al. 2006), and we assume that the tool is not always surrounded by a horizontally layered medium.
The tensor ∆σ (r) = σ (r)−σ 0 I, denotes the conductivity contrast between the actual anisotropic medium and the background medium, and with I the identity tensor.The electric Green's tensor G E (r, r ′ ) and its relation to the magnetic Green's tensor G H (r, r ′ ) for a homogenous isotropic medium are (Fang et al. 2006) where g (r, r ′ ) is the scalar Green's function and k 0 = √ iωµ 0 σ 0 .To calculate the total magnetic fields, the total electric fields need to be obtained first by solving equation (3).Afterward, the calculation of the total magnetic fields is a straightforward addition of the background magnetic fields and the integral term as shown in equation (4).Therefore, the main computational challenge of the IE method is to solve the integral equation (3), which is classified as the Fredholm integral equation of the second kind (Fang et al. 2006).

Numerical Implementation of the Integral Equation Method
A numerical solution of the volume integral in equation ( 3) can be obtained using the method of moments (Gibson 2021).The subsurface model around the induction tool is discretized into a set of grid blocks with centroids r j and volume of ∆v, where j indicates the j-th grid block.
The discretization of equation (3) leads to a linear system of equations that can be expressed in operator form as where G is the operator that represents the discrete convolution integral of the electric Green's tensor G E (r, r ′ ) with the contrast source ∆σE in equation ( 3).The Green's function in equation ( 5) can be discretized by separating the non-singular part of the Green's function and dealing Integral equation method with domain decomposition 7 with the singularity by integrating the Green's function of a single grid block over a spherical domain with an equivalent volume (Gao et al. 2005;Jakobsen & Tveit 2018).The linear system in equation ( 8) can be efficiently solved using a Krylov subspace method because it does not require the matrix of the linear system to be formed explicitly.The desired accuracy of the iterative method is quantified by the relative residual e which is calculated as where • is the L 2 -norm.In this study, we use the generalized minimum residual or GMRES (Saad & Schultz 1986) as the linear system solver.
Green's tensor operator exhibits a convolution structure in each of the tensor components.This property enables the use of FFT to convolve a Green's tensor component G E pq and a component of the contrast source (∆σE) q efficiently (Fang et al. 2006).The p and q indices indicate the component of Green's tensor and the contrast source vector with p and q = x,y,z.At each step of the iterative solver, the convolution integral can be efficiently calculated by where F is the FFT operator and ⊙ denotes elementwise multiplication.This operation reduces the convolution computation complexity from O(N 2 ) to O(Nlog 2 N) with N the number of grid blocks.
It should be noted that the size of the discretized contrast source ∆σE needs to be padded by zeros such that the padded ∆σE has twice the original number of points in all directions to avoid the periodicity in the FFT convolution result.The FFT of Green's tensor can be pre-calculated before calling the iterative solvers to save computational time during the iterative process.

Domain Decomposition
The domain decomposition method attempts to solve the problem for the entire from solutions of the different sub-domains (Saad 2003).In our case, the spatial domain Ω is decomposed into M non-overlapping rectangular sub-domains Ω j , hence see Fig. 1.Adapting the domain decomposition formulation described in (Endo et al. 2009), the convolution integral term or the scattered electric field term in equation ( 3) can be expressed as a sum of scattered electric fields from each of the sub-domains.Subsequently, equation ( 3) can be written as where Ω j indicates the sub-domains with the conductivity anomaly.From equation ( 12), we obtain the following set of integral equations evaluated in each sub-domain: The terms E (i,0) , E (i) , and ∆σ (i) are the background electric field, total electric field, and the conductivity contrast defined at the sub-domain Ω i , respectively.The terms G (ij) ∆σ (j) E (j) in equation ( 13) are the discrete representations of the convolution integral in equation ( 12) which denote the scattered electric fields in the sub-domain Ω i due to the contrast source in the sub-domain Ω j .It can be seen in equation ( 13) that the region without a conductivity anomaly does not contribute to the sum and hence can be omitted from the discretization when calculating the electric field.
By collecting the scattered field terms into the left-hand side of the equations, the linear system of equations in (13) can be expressed with a block-matrix representation, viz.
where A is the block matrix of the rearranged linear system according to the domain decomposition Integral equation method with domain decomposition 9 with E and E (0) are the block vectors containing the total and background electric fields in different sub-domains, respectively.These terms are defined as . . . ,0) . . .
Each block in the matrix A indicates interaction terms between the sub-domains.The diagonal terms I−G (ii) ∆σ (i) in equation ( 15) can be interpreted as the intra-domain interaction within a sub-domain while the off-diagonal terms −G (ij) ∆σ (j) represent the inter-domain interaction terms.Since the sub-domains are rectangular, the convolution integrals with Green's tensor in the intra-and inter-domain interaction terms can still be calculated using the FFT.
To solve the rearranged linear system of equation with domain decomposition in equation ( 14), the matrix A is preconditioned by splitting the matrix into a strictly-lower-triangular (L), strictly-upper-triangular (U ), and diagonal (D) part (Barrett et al. 1994;Saad 2003): where the matrices L, U , and D are defined by , and respectively.By substituting the matrix splitting in (17) into the equation ( 14) and some simple algebra, we obtain which can be solved by choosing an initial guess of E and iteratively calculating the following with k the iteration number.The iteration described in equation ( 20) corresponds to the block Gauss-Seidel iterative method (Barrett et al. 1994;Saad 2003).The matrix (D + L) has a lower triangular form where the inverse can be obtained using forward substitution (Venkateshan & Swaminathan 2014).The forward substitution process to compute equation ( 20) is outlined in Appendix A. With the forward substitution, the total electric field update in each sub-domain according to equation ( 20) can be expressed in the simple form as: where i = 1,2,. . ., M denotes the number of inner iterations where the IE is solved for one subdomain and the number k denotes the number of the total domain sweeps where the electric field is updated for the entire domain.The inverse operation of the block intra-domain term in equation ( 21) is not calculated using the direct solver, but instead by using a Krylov subspace method to solve the following linear system of equations within each sub-domain: The domain sweep is carried out until the relative residual on the whole domain reaches a desired threshold.The resulting operation of solving equation ( 22) iteratively is equivalent to the formulation described in Zhdanov et al. (2006) and Endo et al. (2009).However, in our derivation, we can see the link between the original formulation to a block-preconditioned iterative method, which is the block Gauss-Seidel iterative method in this case.The convergence of the Gauss-Seidel iterative method depends on the diagonal dominance of the linear system matrix (Saad 2003).In this case, if the sum of the inter-domain terms' norm is small compared to the norm of the intradomain terms in equation ( 22), then this scheme is guaranteed to converge.Since the magnitude of Green's tensor elements depends on the distance between sub-domains, the interaction terms are small when the sub-domains are isolated from each other.When a sub-domain has small contrasts, the interaction is one-sided from the sub-domain with high contrast.These properties should be considered when designing the domain decomposition settings.
Instead of the Gauss-Seidel iterative method, one can also choose the Jacobi iterative method by taking only the diagonal part of the matrix A as the preconditioner of the fixed-point equation instead of its lower triangular part.The fixed-point equation that corresponds to the Jacobi iterative method can be written as which leads to the following linear system of equations to be solved in each sub-domain: Since the right-hand side of equation ( 24) only depends on the solutions at the k-th iteration, the Jacobi iterative method is more straightforward to be implemented in parallel computing environments (Barrett et al. 1994).In this case, the linear system of equations at each sub-domain can be solved with the Krylov solver in parallel and the interaction terms are updated after the Krylov solver computations are done in all sub-domains.The main drawback is that the Gauss-Seidel method generally has better convergence properties than the Jacobi method (Barrett et al. 1994).
To further improve the computation speed, we propose to use a Krylov solver with adaptive target residual when solving the IE linear system of a sub-domain.The main idea is the relative residual of the Krylov solver in a sub-domain only needs to be an order of magnitude less than the full-domain relative residual to achieve the convergence of the Gauss-Seidel or Jacobi iteration.
Inaccurate approximate solutions from the Krylov solver are acceptable at the beginning of the iteration and the relative residual target of the Krylov solver is lowered as the full-domain relative solver is decreasing during the Gauss-Seidel or Jacobi iteration.Additionally, the initial guess for the Krylov solver in the current outer iteration is updated from the result of the previous outer iteration.Detailed implementation of this strategy is shown in Algorithm 1.
To further accelerate the computation, it is possible to use the contraction integral form which accelerates the Krylov solver convergence rate (Endo et al. 2009;Zhdanov et al. 2006).However, in this study, we use the conventional integral equation formulation in the Krylov solver to reduce the complexity of evaluating the performance of the IE with domain decomposition.

NUMERICAL RESULTS & DISCUSSION
In this section, we present three numerical cases to demonstrate the effectiveness of the domain decomposition preconditioning of the IE method.The first case is a model with two anomalous sub-domains separated by an isotropic medium with conductivity equal to the background conductivity.In the second case, we present a model where the anomalous isotropic conductivity is surrounded by an anisotropic medium.Lastly, we simulate a logging scenario across a faulted sand formation surrounded by anisotropic shale layers.We use the conventional IE formulation as described in section 2 in the GMRES solver for both full-domain IE and IE with domain decomposition (IE-DD) method.All numerical experiments presented in this paper are performed on a laptop with an AMD Ryzen 7 4800H processor and NVIDIA GeForce RTX 3060 Laptop GPU using MATLAB with GPU support enabled.We have compared our full-domain IE code with existing 1D semi-analytical solution (Shahriari et al. 2018) and 3D finite-volume method (Hou et al. 2006).This comparison is shown in Appendix B and our results show a good agreement with less than one per cent average absolute difference.

Isolated Sub-domains Example
We consider two isolated anomalous sub-domains embedded in an isotropic medium background as shown in Fig. 2. The background conductivity σ 0 is equal to 0.1 S m −1 and the conductivity in the anomalous sub-domain is equal to 0.01 S m −1 .A transmitter with 24 KHz frequency is located at the origin (x = 0, y = 0, z = 0 m) and it is oriented in the x-direction.The whole domain is discretized into 128 × 128 × 128 grid blocks with a uniform grid size of 0.25 × 0.25 × 0.25 m 3 .
The two anomalous sub-domains are set to have an equal size of 30 × 30 × 7.5 m 3 .The distance between the closest edges of the two sub-domains is 10 m which is approximately equal to the skin depth of the background medium given the transmitter frequency.
To solve the conventional full-domain IE and the inverse operation in IE-DD formulation with Gauss-Seidel iterative method, we use restarted GMRES method with 10 restart iterations.We set the relative residual e = 10 -6 for both the conventional full-domain IE and IE-DD Gauss-Seidel iteration stopping criterion.For this case, we did not implement the adaptive relative residual scheme and set e = 10 -6 as the relative residual target for the GMRES solver stopping criterion in the IE-DD to analyze the convergence behavior of the method.Because the medium without contrast does not contribute to the scattering field, the relative residuals are only computed within the anomalous sub-domains in both cases.
The sub-domains without the conductivity anomalies are excluded from the discretization in the IE-DD iterations.This results in only half of the total number of grid blocks of the full-domain IE being discretized in the IE-DD iterations.Fig. 3 shows the convergence plot and the total number of GMRES iterations taken to reach the target residual.The number of GMRES iterations required is decreasing in each of the Gauss-Seidel iterations as the relative residual is decreasing.
This indicates that the changes in the electric fields due to the interaction terms become smaller as the initial guess for the GMRES solver is updated in each of the Gauss-Seidel iterations.Overall, it took only four Gauss-Seidel iterations for the method to converge below the threshold level with a total of 278 GMRES iterations and the total computational time is 23 s.The full-domain IE took 86 GMRES iterations to converge on the same error level, but the computation time is 53 s.Even though there is more GMRES iteration in IE-DD, the total computational time is approximately twice faster compared to the full-domain IE because the operation at each of the GMRES iterations in IE-DD works on a smaller domain with the number of blocks equal to a quarter of the full-domain grid blocks in each domain.
The magnetic field comparison between both methods is shown in Fig. 4. Qualitatively, there are no differences observed because both methods show similar numerical results within less than 0.01 per cent average normalized magnitude difference.Therefore, the IE-DD method will give the same result within the same relative residual level as the full-domain IE. iterative method and with the same adaptive GMRES relative residual as the second scheme.In the adaptive relative residual scheme, the relative residual stopping criterion is set to be one order of magnitude lower than the relative residual calculated on the whole domain or full-domain relative residual divided by ten.

Simple Anisotropic Medium Example
Fig. 6 displays the comparison between the total GMRES iteration in each Gauss-Seidel iteration for both IE-DD-GS schemes.The total number of GMRES iterations generally increases along with the Gauss-Seidel iteration in the adaptive relative residual scheme, while it is decreasing in the fixed relative residual scheme.In both cases, sub-domain 1, which contains the faulted resistive layer, took the greatest number of GMRES iterations.Since the number of GMRES iterations is proportional to the conductivity contrast, this indicates that sub-domain 1 has the largest conductivity contrast compared to the other two sub-domains.
Table 1 summarizes  shown in Fig. 7 indicates that the IE-DD-GS-A has a better convergence rate compared to the IE-DD-Jacobi-A.Besides having a smaller number of total outer iterations, the number of GMRES iterations taken in the IE-DD-GS-A is also fewer compared to the IE-DD-Jacobi-A.However, the computation time of the IE-DD-Jacobi-A can potentially be improved by further utilizing parallel computation to independently solve the linear system in each sub-domains.

Logging Simulation Across a Complex Formation
We simulated induction logs across the faulted anisotropic formation with an 85 • drilling angle as illustrated in Fig. 8a.This formation consists of anisotropic shale layers surrounding isotropic sand layers.The shale layers are indicated by the blue colours and the sand layers are indicated by the yellow colours in Fig. 8.The model has 2.5D main structural features with the addition of a simple 3D Gaussian perturbation only in the sand layers to imitate a fluid distribution in a reservoir.This perturbation is defined by where the subscripts sand denote the values located in the sand layers and the superscripts u indicate the defined unperturbed value; r c is the location of the maximum perturbation; α and γ are the factors that control the magnitude and range of the perturbation, respectively.In this example, we set the peak perturbation location r c at x = 500 m, y = 0 m, and z = 40 m; and define α = 4 and γ = 50 m.
We use a moving 3D forward modelling window to simulate a moving transmitter scenario.
The z-direction in the forward modelling window is directed to the drilling direction so it is consistent with the component direction of the induction tools (Pardo & Torres-Verdín 2015).Hence the coordinate system in the window is rotated from the cartesian coordinate according to the drilling direction as illustrated in Fig. 8a.Consequently, the conductivity tensor elements are transformed following the domain rotation (Gao 2006), see Appendix C for detail.In each of the forward modelling windows, we set a constant background conductivity σ 0 = 0.1 S m −1 .
Following the typical tool configurations described in Antonsen et al. ( 2022), we set a zoriented transmitter with a frequency of 24 KHz and three receivers with spacings of 7, 15, and of 32 × 32 × 32 m 3 may not be enough to capture the full sensitivity of all the receivers, especially the one with the largest spacing.Hence, we tested two different windows with different sizes of 32 × 32 × 64 m 3 and 64 × 64 × 64 m 3 to see different sensitivities of the receivers with the forward modelling domain size.We refer to the smaller window as window 1 and the larger one as window 2. In both windows, we keep a grid size of 0.25 × 0.25 × 0.25 m 3 resulting in a total of 128 × 128 × 256 and 256 × 256 × 256 grid blocks for window 1 and 2, respectively.Since our current computer GPU memory size can only handle a maximum of 128 3 grid blocks for the GMRES solver computation, we decompose the window 1 and 2 into two and eight sub-domains in the drilling direction, respectively, as illustrated in Fig. 10.In this way, the memory requirements to calculate the electric field with both windows are reduced to the memory requirement for the calculation using 128 3 grid blocks plus the memory for storing the electric fields.This allows us to fully take advantage of the acceleration with the GPU implementation.
The logging position starts at x = 0 m, y = 0 m, z = 0 m and ends at x = 900 m, y = 0 m, z = 78.74m.In each logging position, we use the IE-DD-GS-A scheme and set a full-domain tolerance of 10 -3 as the stopping criterion.Fig. 11 shows the magnitude of the z-component of the magnetic field |H zz | measured in the receivers at each transmitter position.Qualitatively, the differences in the results between the two window settings are increasing with the receiver spacings.This result shows different sensitivities of the transmitter-receiver configuration, and we can observe that the sensitivity range is proportional to the receiver spacing.
The computation time required to calculate the magnetic field for one logging position using the window 1 and 2 settings takes an average of approximately 1.5 and 15 minutes, respectively.
Updating the interaction terms is the most expensive part of the computation time, taking up around 80 per cent of the time at every iteration due to the operation acting on the entire domain that consists of a massive number of grid blocks.In every position, it took less than seven Gauss-Seidel iterations to reach the desired tolerance.

CONCLUSION
The linear system of equations arising in the IE method for 3D EM method modelling can be naturally decomposed into a set of linear systems of equations that correspond to the IE in different parts of the modelling domain.The IE-DD formulation is reducing the memory requirement to compute a large-scale problem as it provides the connection between each sub-domains while still maintaining the viability of using FFTs to calculate the convolution integral operation.By expressing these linear systems of equations in a block matrix representation where each block represents the interactions between the domains, we have made a link between the derivation in Zhdanov et al. ( 2006) and Endo et al. (2009) with a preconditioned fixed-point iteration using domain decomposition method.Depending on the choice of the preconditioner, the fixed-point iteration corresponds to the block Gauss-Seidel and Jacobi iterative method.In every Gauss-Seidel or Jacobi iteration, the inverse of the block intra-domain interaction term is calculated using the Krylov subspace method instead of a direct solver.
Our numerical experiment results show that a reduction in computation time can be achieved although the total number of GMRES solver iterations in IE-DD schemes is more than in the fulldomain IE.This speed-up is due to the GMRES solver in the decomposed domains being cheaper to compute and it is shown that it only takes five to eight IE-DD outer iterations to reach the desired tolerance.Additionally, specifying adaptive relative residual stopping improves the computation time of the IE-DD by reducing the total number of GMRES iterations required for reaching desired error tolerance.The Gauss-Seidel preconditioning with adaptive relative residual has the fastest computation time among the schemes that we tested in this study.This scheme reduces the computation time of the conventional IE by approximately 35 per cent.The scheme with Jacobi preconditioning takes longer computation time compared to the one with Gauss-Seidel.However, the form of the Jacobi iterative method is more suitable for parallel computation as the operation in each sub-domain can be computed independently, which is a subject for future implementation.
In this study, we have only implemented IE-DD with a simple iterative update corresponding to Gauss-Seidel and Jacobi iterative method.The Gauss-Seidel and Jacobi iterative methods are in general not very competitive in terms of convergence compared to the Krylov subspace method Integral equation method with domain decomposition 19 (Barrett et al. 1994).Therefore, further potential improvement of the IE-DD presented in this study is obtained by implementing the Krylov subspace as the outer iteration update instead of the Gauss-Seidel and Jacobi iteration update.Another interesting application of the domain decomposition in the IE method would be to incorporate a direct method that can be computed in parallel into the domain decomposition preconditioner, for example using the T-matrix method (Jakobsen & Tveit 2018;Sommer & Jakobsen 2018).

METHOD
The inverse of a lower-triangular matrix can be obtained through forward substitution (Venkateshan & Swaminathan, 2014).For simplicity, we demonstrate the forward-substitution process in the case of domain decomposition with three sub-domains.We can write the block-matrix representation and its splitting for three sub-domains: with the matrices D, L, and U are the block diagonal, strictly lower-triangular, and strictly upper-triangular of A, respectively.These terms are defined as

Integral equation method with domain decomposition 23
Here each of the blocks denotes the inter-domain and intra-domain operators described in section 2 as The fixed-point equations using matrix A that corresponds to the Gauss-Seidel iterative method By substituting the matrices D, L, and U for 3 sub-domains and calculating the inverse of (D + L) , we obtain: where and the terms R i denote the i-th row of the second term in the right-hand side of equation (A.5) By multiplying the matrix on the right-hand side of equation (A.6), we obtain the following equations: Notice that the term A −1 11 R 1 in the second right-hand side term of equations (A.10) and the third right-hand side term of (A.11) can be substituted by E (1),k+1 from equation (A.9).Also, components calculated using our IE code is less than one per cent compared to the 1D semianalytical result.

APPENDIX C: CONDUCTIVITY TENSOR TRANSFORMATION FOR TRANSVERSELY ANISOTROPIC FORMATION
The tensor structure of electrical conductivity σ for anisotropic media is generally expressed as with an off-diagonal symmetry σ ij = σ ji .For a vertical transverse isotropic medium with the z-axis as the vertical axis, the conductivity tensor is written as (Gao 2006;Jakobsen & Tveit 2018) where σ h and σ v are the conductivity in the horizontal and vertical directions, respectively.In the case where the angle between the formation layering and the drilling trajectory is not 90, it is necessary to rotate the conductivity tensor from the formation coordinate system to the induction tool coordinate system.For a transversely isotropic formation, the explicit expressions of the rotated conductivity tensor are as follows (Gao 2006) where θ and φ are the z-axis rotation and y-axis rotation angles from the formation coordinate system to the induction tool coordinate system, respectively.The resulting conductivity tensor has non-zero off-diagonal components and our method can deal with this complication without making any other changes in the implementation.b = E (i,0) + i−1 j=1 G (ij) ∆σ (j) E (j),k+1 + M j=i+1 G (ij) ∆σ (j) E (j),k else if Jacobi preconditioning b = E (i,0) + M j=1 G (ij) ∆σ (j) E (j),k end if set : initial guess = E (i),k , threshold = e/10 E (i),k+1 = GMRES A = I − G (ii) ∆σ (i) , b, initial guess, threshold end for k = k + 1 end while

Fig. 5
Fig.5shows an xz-plane view of a faulted resistive isotropic medium with a conductivity of 0.01 S m −1 surrounded by an anisotropic medium with vertical transverse isotropy.The conductivity tensor of the anisotropic medium consists of the conductivity in the horizontal and vertical direction with the value of σ h = 0.2 S m −1 and σ v = 0.1 S m −1 , respectively.The conductivity of the media does not vary in the y-direction.For the background medium, we choose a homogenous isotropic medium with the conductivity of σ 0 = 0.01 S m −1 .A transmitter with 24 KHz frequency is located at the origin and is oriented in the x-direction.We set 10 -6 as the relative residual target for both methods.The whole domain is discretized into 120 × 120 × 120 grid blocks with a grid size of 0.25 × 0.25 × 0.25 m 3 .
scheme compared to the other two methods.The IE-DD-GS-A scheme shows a faster computation time compared to the IE-DD-GS-F because there are fewer GMRES iterations in the IE-DD-GS-A scheme.Therefore, specifying the adaptive relative residual for the Krylov solver in the IE-DD improves the computation time of the original IE-DD formulation with the cost of going through more Gauss-Seidel iterations.The full-domain relative residual plot in each of the outer iterations

Figure 1 .
Figure 1.Schematic of the domain decomposition.The black arrows indicate the field coming from the transmitter Tx to the sub-domains Ω j .The double-headed green arrows indicate the scatterers' interaction between the sub-domains.The transmitter can also be located in the anomalous domain.

Figure 2 .Figure 3 .
Figure 2. (a) xz-plane view of the model at y = 0 m.(b) xy-and xz-slice of the model in 3D view.

Figure 4 .
Figure 4. xz-plane slice of magnetic fields at y = 0 m.(a) Real and (b) imaginary parts of the magnetic field obtained from the full-domain IE.(c) Real and (d) imaginary parts of the magnetic field obtained from the IE-DD.A transmitter oriented to the x-direction is located at x = 0 m, y = 0 m, and z = 0 m.

Figure 5 .
Figure 5. (a) xz-plane view of a faulted resistive layer surrounded by an anisotropic medium at y = 0 m and the domain decomposition setting.(b) 3D view with the anisotropic layers removed.A transmitter oriented to the x-direction is located at x = 0 m, y = 0 m, and z = 0 m.

Figure 6 .Figure 7 .
Figure 6.Total GMRES iterations within each Gauss-Seidel iteration for IE-DD-GS with (a) fixed and (b) adaptive relative residual stopping criterion.
Figure 8.(a) xz-plane view at y = 0 m.(b) xy-plane view at z = 40 m.(c) 3D view of the model with the shale layers removed.The magenta dashed lines indicate the drilling trajectory.The black box is a forward modelling window example with a size of 64 × 64 × 64 m 3 at one logging position.

Figure 9 .
Figure 9. Illustration of an induction tool with a single transmitter and three receivers.Tx and Rx stand for transmitter and receiver, respectively.

Figure 10 .Figure 11 .
Figure 10.Domain decomposition illustration of (a) window 1 and (b) window 2. z'-axis is the drilling direction.The range in the y' direction is equal to the range in the x' direction.

Figure B2 .
Figure B2.Comparison of the calculated magnetic field couplings with different numerical methods.
set : ǫ = threshold value, M = number of sub-domains, initialize : E 0 := E (0) , k : This work is part of the Center for Research-based Innovation DigiWells: Digital Well Center for Value Creation, Competitiveness and Minimum Environmental Footprint (NFR SFI project no.309589, https://DigiWells.no).The center is a cooperation of NORCE Norwegian Research Centre AS, the University of Stavanger, the Norwegian University of Science and Technology (NTNU), and the University of Bergen.It is funded by Aker BP, ConocoPhillips, Equinor, Lundin Energy, TotalEnergies, Vår Energi, Wintershall Dea, and the Research Council of Norway.

Table 1 .
Performance of IE with different schemes.suband e f ull are the relative residual of the sub-and full-domain, respectively.Algorithm 1. IE method with Domain Decomposition Preconditioning.