Research paperGEMM3D: An Edge Finite Element program for 3D modeling of electromagnetic fields and sensitivities for geophysical applications
Introduction
The Finite Element (FE) method has been used to simulate geophysical electromagnetic data for a long time. It generates a large sparse system of linear equations in which the coefficient matrix does not depend on the space distribution of the source. In formulations that use the separation of the field into primary and secondary terms, the primary field enters only at the right end of the equations. This means that to simulate multiple source locations, for example, it is necessary to factor the matrix only once per frequency, and to solve multiple problems by forward and back substitutions.
When the variables in an FE implementation are associated with the mesh nodes, the method is termed Nodal FE. An alternative, more recent, approach defines the variables as the components of the electric or the magnetic field in the directions of the edges of the 3D elements. In this case, the method is called Vector or Edge Finite Element (EFE) (Jin, 2002). This formulation has the advantages of allowing the presence of discontinuous field components, as when the electric field crosses interfaces between regions of discontinuous conductivity, and of ensuring the null divergence condition for the electric field, which is not guaranteed in the nodal formulation (Gerrit, 1994, Cai et al., 2014a). Edge Finite Elements have been applied to model widely different problems, like induction well logging (Heagy and Oldenburg, 2019), time-domain problems (Cai et al., 2017b, Qi et al., 2019), marine controlled source electromagnetics (Cai et al., 2014b, Um et al., 2017, Cai et al., 2017a, Castillo-Reyes et al., 2018). The method has been used in open source (Heagy et al., 2017) as well as in commercial (Butler and Zhang, 2016) software.
To implement a Finite Element (FE) formulation for an electromagnetic problem, an efficient approach (Wannamaker et al., 1987, Hohmann, 1987) is to build the models as a group of 3D finite bodies placed inside a (1D) layered earth background, and to write the EM field as the sum of two parts: the field generated by the geophysical source inside the background model, called primary field, and the so called secondary field, which is the difference between the total field in the 3D medium and the primary field. The differential equations are written in terms of the secondary field, and the primary field values go only in the function that builds the source term on the right hand side of the equations. As opposed to formulations that solve for the total field, this approach has the advantages of avoiding the representation of the singular sources in the mesh, and of generating a numerically stable solution to a problem where the secondary fields are orders of magnitudes less intense than the primary ones. The edge FE formulation using this field separation into primary and secondary components has been applied by Cai et al. (2014a) to simulate marine Controlled Source ElectroMagnetic (mCSEM) data in a grid of parallelepiped elements, and by Christoph et al. (2011), using a non-structured mesh with tetrahedral elements.
This paper presents a software package named Geophysical ElectroMagnetic Modeling 3D — GEMM3D that implements an EFE algorithm to simulate 3D frequency domain electromagnetic data from different geophysical sources: plane waves, electric and magnetic dipoles, and large loops.
An important feature of this program is the efficient implementation of the Adjoint State (AS) method (McGillivray et al., 1994) to generate sensitivities for a grid of homogeneous cells, for use in inverse problems. The AS method calculates the electric field from auxiliary dipole sources located at the measurement positions, together with the field from the simulated geophysical source, in a formulation that generates sensitivities by calculating the volume integrals of the inner product of these two electric fields over each cell. The discretization with hexahedral elements yields the components of the electric field in the directions of the three coordinate axes, which makes the calculation of inner products an easy and quick task.
By implementing the electric field formulation of the edge FE method, we are able to generate sensitivities efficiently, with only the cost of additional forward and back-substitutions for each receiver position.
In the following sections we detail the methods implemented in the program and present a few validation examples.
Section snippets
Methods
The GEMM3D program performs the forward field calculations in the frequency domain by using primary fields in the source term of a 3D edge finite element implementation to calculate the electric field components in the directions of the edges of the finite elements. The discretization mesh is a non-uniform grid of parallelepipeds, whose edges are aligned with the directions of the coordinate system. This choice of discretization allows an easy and efficient implementation of the Adjoint State
Some implementation details
For the plane wave simulations, the primary fields need only be calculated at the depths of the central points of the element edges in a single column. The values of the 1D fields at an edge in this column are simply repeated for the remaining edges of the grid, at the same depth.
The plane wave has analytical solutions in 1D layered models, with no frequency limitations to simulate MT data. Any frequency used in a real survey can be used to generate synthetic data. All the other sources require
Validation of the data simulation
To validate the code, we have simulated data from models found in the geophysics literature. Here we present results from two different field sources: plane waves to generate magnetotelluric data, and the horizontal electric dipole for the marine CSEM method. All results were generated in a workstation with a 6 core Intel Xeon CPU running at 3.5 GHz, and 64 GB of RAM.
Validation of sensitivities
The discretization of the medium with a structured hexahedral element grid has the advantage of producing the electric field components in the directions of the coordinate axes. To calculate sensitivities with the Adjoint State method, as described by McGillivray et al. (1994), the program makes use of these components inside each inversion cell to evaluate the volume integrals expressed by Eq. (36).
This section shows validation results for the sensitivity calculations, by comparing the results
Final remarks
This paper presented the GEMM3D program for modeling 3D electromagnetic fields from geophysical sources, using the Edge Finite Element method. The program is shown to be effective for 3D simulations of data of a wide variety of frequency domain geophysical methods. It uses a parallel implementation of a direct solver, and takes advantage of the possibility of factoring the coefficient matrix once and solving multiple systems in parallel. The implementation of the electric field formulation of
CRediT authorship contribution statement
Carlos Mateus Barriga Nunes: Methodology, Software. Cícero Régis: Supervision, Writing - original draft, Writing - review & editing.
Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Acknowledgments
The authors thank the Federal University of Pará, Brazil and the National Institute of Science and Technology of Petroleum Geophysics, Brazil, for the support. Nunes thanks CAPES, Brazil for the doctorate scholarship that funded the research for which the code was developed. We thank the editor Prof. Dario Grana and the anonimous reviewers that contributed to improve the quality of the paper.
References (39)
- et al.
Forward modeling of geophysical electromagnetic methods using comsol
Comput. Geosci.
(2016) - et al.
Parallelized 3d csem modeling using edge-based finite element with total field formulation and unstructured mesh
Comput. Geosci.
(2017) - et al.
Finite-element time-domain modeling of electromagnetic data in general dispersive medium using adaptive pad series
Comput. Geosci.
(2017) - et al.
3d controlled-source electromagnetic modeling in anisotropic medium using edge-based finite element method
Comput. Geosci.
(2014) - et al.
3d controlled-source electromagnetic modeling in anisotropic medium using edge-based finite element method
Comput. Geosci.
(2014) - et al.
Petgem, A parallel code for 3d csem forward modeling using edge finite elements
Comput. Geosci.
(2018) - et al.
A framework for simulation and inversion in electromagnetics
Comput. Geosci.
(2017) - et al.
Modeling electromagnetics on cylindrical meshes with applications to steel-cased wells
Comput. Geosci.
(2019) - et al.
Evaluation of parallel direct sparse linear solvers in electromagnetic geophysical problems
Comput. Geosci.
(2016) - et al.
Three-dimensional modeling of frequency- and time-domain electromagnetic methods with induced polarization effects
Comput. Geosci.
(2019)
Pardiso: a high-performance serial and parallel sparse linear solver in semiconductor device simulation
Future Gener. Comput. Syst.
A tetrahedral mesh generation approach for 3d marine controlled-source electromagnetic modeling
Comput. Geosci.
Methods for modelling electromagnetic fields
J. Appl. Geophys.
Numerical integration of related Hankel transforms of orders 0 and 1 by adaptive digital filtering
Geophysics
Optimized fast hankel transform filters
Geophys. Prospect.
Three-dimensional adaptive higher order finite element simulation for geo-electromagnetics—a marine csem example
Geophys. J. Int.
An introduction to marine controlled-source electromagnetic methods for hydrocarbon exploration
Geophysics
Fast marine csem inversion in the cmp domain using analytical derivatives
Geophysics
Cited by (9)
Iterative solver with folded preconditioner for finite element simulation of magnetotelluric fields
2022, Computers and GeosciencesHigh order edge-based finite elements for 3D magnetotelluric modeling with unstructured meshes
2022, Computers and GeosciencesCitation Excerpt :In contrast to nodal elements, edge elements only impose continuity of tangent fields, hence allowing a discontinuous normal field component. EFEM has been applied on applications related to hexahedral (e.g. Nam et al., 2007; Rivera-Rios et al., 2019; Nunes and Régis, 2020) and tetrahedral (e.g. Liu et al., 2008) elements. Tetrahedral elements have the advantage of providing greater flexibility in fitting complicated domains (cf. Sack and Urrutia, 2000 p. 319).
3D edge-based and nodal finite element modeling of magnetotelluric in general anisotropic media
2022, Computers and GeosciencesCitation Excerpt :However, to simulate a subsurface electrical structure that is sufficiently realistic, three-dimensional (3D) numerical modeling of the arbitrary anisotropy MT problem is needed. The forward modeling of 3D MT anisotropy is consistent with the isotropic numerical method, including the integral equation (IE) method (Avdeev and Knizhnik, 2009; Kruglyakov and Kuvshinov, 2018), FD (Newman and Alumbaugh, 1995; Wang et al., 2021), finite element (FE) (Jin, 2002; Nunes and Regis, 2020; da Piedade et al., 2021; Zhang et al., 2021), and finite volume (FV) methods (Weiss and Constable, 2006; Kruglyakov and Kuvshinov, 2018). At present, researchers have developed several 3D forward modeling algorithms for anisotropic media.
Computational cost comparison between nodal and vector finite elements in the modeling of controlled source electromagnetic data using a direct solver
2021, Computers and GeosciencesCitation Excerpt :Numerical derivatives are necessary only for calculating the remaining field, when it is needed. Since its inception, the VFE method has been implemented to simulate MCSEM data in 3D models using structured hexahedral grids (da Silva et al., 2012; Nunes and Régis, 2020), isoparametric grids (Cai et al., 2014), and non-structured tetrahedral meshes (Cai et al., 2017; Rochlitz et al., 2019), which is the chosen kind of mesh in this study. The successes in the application of the VFE method prompted us to try to measure quantitatively the differences in performance of the two variants of the Finite Element method, in terms of computer resources, specifically memory and processing time, to achieve the same level of accuracy in the modeling of geophysical EM data.