- Email: [email protected]
A tetrahedral mesh generation approach for 3D marine controlled-source electromagnetic modeling
Computers & Geosciences 100 (2017) 1–9
Contents lists available at ScienceDirect
Computers & Geosciences journal homepage: www.elsevier.com/locate/cageo
crossmark
A tetrahedral mesh generation approach for 3D marine controlled-source electromagnetic modeling Evan Schankee Uma, Seung-Sep Kimb, Haohuan Fuc,d,
⁎
a
Affiliate of Earth and Environmental Sciences, Lawrence Berkeley National Laboratory, CA, USA Geology and Earth Environmental Sciences, Chungnam National University, Daejeon, South Korea Tsinghua University, Ministry of Education Key Laboratory for Earth System Modeling and Center for Earth System Science, Beijing, China d National Supercomputing Center in Wuxi, Wuxi, China b c
A R T I C L E I N F O
A BS T RAC T
Keywords: Marine electromagnetic geophysics Finite element modeling Parallel computing
3D finite-element (FE) mesh generation is a major hurdle for marine controlled-source electromagnetic (CSEM) modeling. In this paper, we present a FE discretization operator (FEDO) that automatically converts a 3D finitedifference (FD) model into reliable and efficient tetrahedral FE meshes for CSEM modeling. FEDO sets up wireframes of a background seabed model that precisely honors the seafloor topography. The wireframes are then partitioned into multiple regions. Outer regions of the wireframes are discretized with coarse tetrahedral elements whose maximum size is as large as a skin depth of the regions. We demonstrate that such coarse meshes can produce accurate FE solutions because numerical dispersion errors of tetrahedral meshes do not accumulate but oscillates. In contrast, central regions of the wireframes are discretized with fine tetrahedral elements to describe complex geology in detail. The conductivity distribution is mapped from FD to FE meshes in a volume-averaged sense. To avoid excessive mesh refinement around receivers, we introduce an effective receiver size. Major advantages of FEDO are summarized as follow. First, FEDO automatically generates reliable and economic tetrahedral FE meshes without adaptive meshing or interactive CAD workflows. Second, FEDO produces FE meshes that precisely honor the boundaries of the seafloor topography. Third, FEDO derives multiple sets of FE meshes from a given FD model. Each FE mesh is optimized for a different set of sources and receivers and is fed to a subgroup of processors on a parallel computer. This divide and conquer approach improves the parallel scalability of the FE solution. Both accuracy and effectiveness of FEDO are demonstrated with various CSEM examples.
1. Introduction In the past decade, 3D finite-element (FE) solutions have been widely used in marine controlled-source electromagnetic (CSEM) modeling (e.g. Key and Ovall, 2011; Schwarzbach et al., 2011). In contrast to finite-difference (FD) solutions that approximate noncoordinate-conforming structures with small rectangular stair steps, FE uses geometry-conforming tetrahedral meshes and precisely represents complex seafloor topography. However, the advantage of FE over FD comes with extra complication. It is considered difficult to iteratively solve a system of FE equations. A system matrix with tetrahedra is unstructured and not diagonally dominant. A simple Jacobi preconditioner used in FD solutions does not ensure convergence of FE Krylov solutions. FE solutions require numerically expensive preconditioners such as incomplete factorization (Um et al., 2013) and others. It is also difficult to
⁎
find a robust preconditioner suitable to a wide range of EM problems. Accordingly, direct solvers are often the method of choice for FE solutions despite their large memory requirements and low parallel scalabilities (Fu et al., 2015). It is also considered difficult and time-consuming to create 3D FE meshes. 3D FE mesh generation requires good literacy in computeraided design (CAD) software that may have a steep learning curve. To mitigate the difficulty, adaptive FE methods have been introduced (Li and Key, 2007; Key and Ovall, 2011; Schwarzbach et al., 2011). They start with coarse meshes and successively refine meshes until a required tolerance is met. However, in either use of interactive CAD or automated adaptive refinement methods, it is nontrivial to generate reliable and efficient 3D meshes for a complex multi-scale model. The goal of this paper is to present a 3D FE discretization operator (FEDO) that automatically generates reliable and efficient tetrahedral meshes for marine CSEM. We develop strategies for directly creating
Corresponding author. E-mail addresses: [email protected] (E.S. Um), [email protected] (S.-S. Kim), [email protected] (H. Fu).
http://dx.doi.org/10.1016/j.cageo.2016.11.007 Received 4 August 2016; Received in revised form 17 November 2016; Accepted 20 November 2016 Available online 21 November 2016 0098-3004/ © 2016 Elsevier Ltd. All rights reserved.
Computers & Geosciences 100 (2017) 1–9
E.S. Um et al.
complete FE meshes without adaptive refinement. The development requires understanding basic properties of tetrahedral meshes for diffusive EM modeling, which are often not obvious when CAD or adaptive refinement methods are employed. Ultimately, we aim to demonstrate that FEDO can quickly derive accurate and efficient FE meshes from a given seabed model without computational overhead associated with adaptive mesh refinement or interactive CAD works. The remainder of this paper is organized as follows. We first review the FE formulation for CSEM. Its numerical dispersion characteristics are examined to determine proper tetrahedral sizes. To avoid excessive mesh refinement around receivers, we introduce an effective receiver size. We follow this by presenting a matrix that maps conductivity from FD model to FE simulation meshes. The three components above are casted into FEDO. Finally, we apply FEDO to complex offshore models, compute their CSEM solutions and demonstrate its accuracy and efficiency.
Table 1 δ and λ in the seawater.
(1)
where e(r) is the electric field at position r, Js(r) is an electric source at an angular frequency ω , μo is the magnetic permeability of free space (4π×10-7 H/m), and σ is the conductivity. The development of the equivalent weak statement of Eq. (1) requires the multiplication of Eq. (1) by the edge basis function and the integration over the model domain of V, resulting in e
1 nie(r)⋅( ∇ × ∇ × ee(r) + l j ωσeee(r) + l j ω Js(r))dV = 0. μo th
(2)
nie(r)
The superscript e denotes the e tetrahedral element, with i varying from 1 to 6 is a set of edge basis functions. V e is the volume of the eth element. nie(r) is also chosen as the basis. Thus, the electric field at a point inside or on a given element is expanded as 6
ek(r) =
6
∑ ekj (r) = ∑ ujk nkj (r), j =1
j =1
(3)
is the unknown amplitude of the electric field on the jth edge where of the kth element. Accordingly, the FE method is the second-order accurate (Jin and Riley, 2008). Substituting Eq. (3) into Eq. (2) and using the homogeneous Dirichlet boundary conditions yield
ujk
(Ae + μo ωjl Be)ue = − μo ωiˆ se
(4)
where
(i , j ) element of Ae =
∭V
(i , j ) element of Be =
∭V
i element of se =
∭V
e
ue = [ u1e u 2e ... u6e].
e
e
∇ × nie(r)⋅∇ × nej(r)dV ,
(5)
nie(r)⋅σ e nej(r)dV ,
(6)
nie(r)⋅Js(r)dV ,
λ (m)
0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50
873.0 712.8 617.3 552.1 504.0 466.6 436.5 411.5 390.4
5484.9 4478.4 3878.4 3469.0 3166.7 2931.8 2742.5 2585.6 2452.9
To determine tetrahedral mesh density required for accurate solutions, we examine numerical dispersion characteristics for a homogeneous whole-space seawater (3.33 S/m) model. We are interested in efficient discretization of the seawater because the seawater can be a large portion of a CSEM model and requires the smallest elements due to its highest conductivity. The seawater model is 20×10×10 km in the x-, y- and z-direction, respectively. The lower left corner of the model is at (−5, −5, 5 km). A 250 m long x-oriented electric source is placed at (0, 0, 0 km). 20 m long x-oriented electric receivers are placed with 1 km spacing from x=1–10 km. We consider nine source frequencies from 0.1 to 0.5 Hz with 0.05 Hz interval. Skin depth (δ) and wavelength (λ) of the seawater are presented in Table 1. The model (Fig. 1) is discretized using tetrahedral elements whose vary from 10 m (near sources and receivers) to 400 m (most areas). This is slightly smaller than δ at 0.3 Hz. Note that the FE meshes are refined around not only the source but also the receivers. This is one of the major differences between FD and FE. The mesh refinement near the receivers will be discussed in the next section. The transition from small to large edges is controlled by the growth factor defined as the maximum rate at which the edge size can grow. The growth factor is typically 1.5–2.0 from one edge to the next. Amplitudes, relative amplitude errors and phases of FE solutions are plotted against analytic solutions in Fig. 2a-c, respectively. The relative errors are defined by ‖(numerical solution–analytic solution)/ analytic solution)‖. It is assumed that over 10 km source-receiver offset, 5% amplitude errors are accurate enough. Based on the criterion, the FE meshes support the lowest four frequencies: 0.10, 0.15, 0.20 and 0.25 Hz. At 0.10 Hz, the boundary effects appear at x=8 km. To highlight the numerical dispersion characteristics of the FE meshes, we also repeat the same experiments with a 2nd-order accurate FD method (Newman and Alumbaugh, 1995). The FD grid size is set to 400 m as same as that in the FE meshes. However, FD does not produce accurate solutions at all nine frequencies (not included here). We reduce the FD grid to 100 m and analyze amplitudes, their relative errors and phases (Fig. 2d–f). Due to their geometric simplicity, the FD grids are not presented here. The FD grids support the lowest four frequencies (Fig. 2e). At 0.10 Hz, the boundary
Since the discretization density required for accurate EM solutions is directly related with an FE formulation, we briefly describe the total field FE formulation (Um et al., 2013) used here. The electric-field diffusion equation is given by
∭V
δ (m)
3. Numerical dispersion analysis
2. Finite element formulation
1 ∇ × ∇ × e(r) + l j ωσe(r) + l j ω Js(r) = 0, μo
f (Hz)
(7) (8)
Eq. (4) is considered local as it results from integration over each element. The local systems from all elements are assembled into a single global system of equations. The global system is solved using a parallel direct solver.
Fig. 1. The FE meshes for the seawater model. Blue and red lines represent source and receivers, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
2
Computers & Geosciences 100 (2017) 1–9
E.S. Um et al.
0.10 Hz 0.15 Hz 0.20 Hz 0.25 Hz 0.30 Hz 0.35 Hz 0.40 Hz 0.45 Hz 0.50 Hz
0.10 Hz 0.15 Hz 0.20 Hz 0.25 Hz 0.30 Hz 0.35 Hz 0.40 Hz 0.45 Hz 0.50 Hz
) Fig. 2. The FE (the first row) and FD (second row) solutions to the seawater model and their comparison with analytic solutions. The solid lines with a square are phases of FE/FD solutions; the broken lines with a circle are those of analytic solutions.
capture rapidly-decaying EM fields. It is our experiences that some refinements are often required in both total and scattered field solutions. In FD, the fine grids propagate throughout the entire computational domain, resulting in a large number of cells. In contrast, FE allows local mesh refinements. However, a major drawback of FE is the fact that it also requires fine elements near receivers. By using the first-order edge basis functions, FE determines an average electric-field vector along each edge of tetrahedral elements. The averaged electricfield vectors are used to interpolate the electric field at a receiver through Eq. (3). If the edges are much longer than a finite-length receiver, the interpolation would be erroneous. Consequently, it is necessary to refine meshes near receivers. A larger number of receivers results in a larger FE problem size. This is a major disadvantage of the FE method presented here. In this section, we examine local mesh refinements around receivers for accurate FE solutions at 0.25 Hz and present a simple method for determining an economic mesh size around receivers. To do this, the seawater model is discretized with eight levels of receiver mesh refinement as summarized in Table 2. For the first seven levels, 20, 40, 80, 120, 160, 200 and 250 m long x-oriented edges are placed at the midpoint of 20 m long receivers. For the eighth level, the 250 m long edges are used again near receivers, but their directions are not aligned with the receiver direction, but are randomly oriented. Fig. 3 shows FE
effect appears at x=8 km. This distance is similar to that of the FE meshes. As demonstrated, FE meshes can be coarser than FD meshes. FD produces accurate solutions when δ is discretized with 4–7 nodes. As shown in Fig. 2e, the mesh density required for accurate FD solutions strongly depends on a source-receiver offset because FD numerical dispersion errors accumulate. In contrast, FE solutions are accurate enough if local element sizes are slightly smaller than local δ. Note that the FE numerical dispersion errors oscillate as the source-receiver offset increases. This is because the FE meshes are not structured but randomly oriented; the errors can be either positive or negative, canceling each other out rather than being accumulated (Jin and Riley, 2008). Thus, the required FE element size does not strongly depend on the maximum source-receiver offset. δ serves as a robust criterion for determining the FE element size. In this paper, the maximum local FE element size is chosen as 90% of local δ for safety. The element size determined by δ is often much larger than CSEM resolution. For example, δ for the seabed of 0.3 S/m is about 1.74 km at 0.25 Hz although CSEM aims to resolve several ten meters thick structures. This observation suggests a basic strategy of determining element sizes for CSEM. A CSEM model is partitioned into two conceptual regions. One region represents a part of seabed where detailed geological structures are described. In this region, the element sizes should be small enough such that fine seabed structures can be described. In the other regions (e.g. the other part of the seabed, the seawater and the air), the local δ criteria determine maximum local element sizes, leading to the economic discretization. Later in this paper, we will further discuss partitioning a CSEM model and its region-by-region discretization. In the next section, we examine the FE mesh refinement required around sources and receivers and present the relaxation of the requirement by introducing an effective receiver size.
Table 2 Mesh refinement levels according to a receiver edge size and resultant total number of elements required to discretize the computational domain.
4. Effective receiver size Around sources, both FD and FE require fine discretization to 3
Mesh refinement level
Edge size for receiver (m)
# of elements
1 2 3 4 5 6 7 8
20 40 80 120 160 200 250 250
322,076 310,831 299,571 294.913 293,244 290,568 287.914 286.749
Computers & Geosciences 100 (2017) 1–9
E.S. Um et al.
meshes at a few selected refinement levels. The FE solutions from each FE mesh are compared against the analytic solution where the actual receiver length is set to 20 m (Fig. 4a–c). As the receiver edge size increases, the relative errors oscillate with larger amplitude and eventually exceed an acceptable level (i.e. 5%). However, Fig. 4a–c suggests that FE solutions can be accurate enough although FE meshes do not precisely resolve the actual receiver length. For example, when the receiver edge size is set to 160 m, the relative errors are bounded within 5%. This relatively low error is due to the fact that λ (i.e. 3469 m) is large enough compared with the actual receiver size (20 m). This result is also verified by comparing the analytic solutions of different actual receiver sizes as shown in Fig. 4d–f where the 20 m long actual receiver can be replaced even with a 250 m long actual receiver. It is also important to place edges of elements at the centers of receivers and align them along the receiver direction. This alignment helps the receiver edge remove the unwanted influence from the other five edges via Eq. (3). For example, when the receiver edge size is set to 250 m, relative errors with and without the alignment are noticeably different (Fig. 4b). It is our experience that in CSEM modeling, it is not necessary to resolve the actual receiver size. We can find an ‘effective’ receiver size that is typically several factors larger than the actual size. As done in Fig. 4e, one can determine an effective receiver size using analytical solutions. Fig. 5 demonstrates the benefit of using an effective receiver size in terms of a problem size. When the FE meshes resolve the actual receiver size (20 m), the total number of elements quickly increases with the number of receivers. In contrast, when an effective receiver size (160 m) is used, the increase is modest. The use of an effective receiver size is particularly recommended when a number of receivers are placed in a grid pattern.
Fig. 3. FE meshes with different refinement levels around receivers. (a) Level 3, (b) Level 5. (c) Level 7.
Fig. 4. The first row: (a) The FE solutions to the seawater model with a receiver length varying from 20 to 250 m and the analytic solution (20 m long receiver, green broken line). (b) The relative amplitude errors of the FE solutions with respect to the analytic solution. (c) Phase plots of the FE and analytic solutions. The second row: (d) The analytic solutions to the seawater model with various receiver lengths. (e) The relative amplitude differences of (d) when the analytic solutions with various receiver lengths are compared with that of 20 m. (f) The phase plots of (d). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
4
Computers & Geosciences 100 (2017) 1–9
E.S. Um et al.
Fig. 5. The number of elements required for the seawater model is plotted as a function of number of receivers when a receiver is resolved by its true (20 m) and effective (160 m) length.
Fig. 8. The FE seafloor topography meshes. The red lines represent sources. The seawater-seabed interface is in blue. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
region 3 and 5. If sources or receivers are on the seafloor, their coordinates are augmented to the coordinates of the seafloor topography. Region 2 and 4 surround region 3 and 5 and provide ample space (e.g. λ) through which EM fields are attenuated. The δ criterion determines the maximum element sizes for region 1, 2 and 4. Element sizes are nearly uniform in these regions. In contrast, region 3 can have a range of element sizes. The bottom of region 3 consists of small elements and defines the irregular seafloor topography. Region 3 also includes a number of small line segments for sources and receivers. For the efficient discretization of region 3, element sizes gradually grow outward from the bottom and source/ receiver positions and are bounded by δ. Because geological structures are mapped from FD to FE meshes as will be described below, element sizes of region 5 should not only satisfy δ criterion but also be roughly equivalent to those of FD meshes in a volumetric sense. For example, if the size of a regular tetrahedral element is 200 m, its volume (9.43×105 m3) is similar to a 100 m FD cube. FEDO utilizes TetGen (Si, 2006), a tetrahedral mesh generation algorithm to discretize the wireframes. Note that the wireframes and their discretization honor only the seafloor interface. The seafloor is the most important interface to CSEM modeling because the receivers can be directly placed on the interface. By precisely discretizing the seafloor geometry, we ensure the correctness of local EM interface conditions and local EM interaction. In contrast to seismic modeling, it is our experiences that diffusive CSEM is much less sensitive to the detailed geometry of other interfaces below the seafloor. Thus, we do not consider other interfaces as critical constraints that need to be explicitly honored in tetrahedral meshing. Instead, they are realized by mapping conductivity values from FD to FE meshes in a volume average sense as described below. After the discretization, FEDO maps conductivity attributes from FD to FE meshes. Region 1, 2 and 3 are assumed to be homogeneous. Their conductivity attributes directly come from the FD model. In contrast, each tetrahedral element of region 4 and 5 is divided into N small sub-elements. Their centroidal conductivities are read from the FD model and are averaged in a volumetric sense. The average is assigned to the tetrahedral element. When some sub-elements of a seabed tetrahedral element reside in the seawater of the FD model, their centroidal conductivities are excluded from the averaging processes. When all sub-elements reside in the seawater of the FD model due to large discrepancies between FD and FE topography, conductivities of adjacent tetrahedral elements are used to interpolate the conductivity of the tetrahedral element. The mapping processes are casted into
Region 1: the air
Region 2: sea
Region 3: sea Source Receiver
Region 4: seabed
Region 5: seabed
Fig. 6. The 3D wireframes.
Fig. 7. The comparison of the FE and FD seafloor topography models.
5. FE discretization operator Based on the analysis above, we present an FEDO for 3D marine CSEM modeling. FEDO requires (1) information about CSEM configurations, (2) the description of seabed conductivities in an FD format and (3) coordinates of seafloor topography. The first two items are those used for FD modeling. The input conductivity model is in the FD format because porosity and fluid saturation models used for constructing a conductivity model are often in the FD format. In contrast, the third item is the new information. FEDO performs three tasks described below. First, FEDO sets up the wireframes of a background seabed model (Fig. 6). The wireframes consist of five internal regions. The coordinates of the seafloor topography is given in a grid pattern and defines the interface of 5
Computers & Geosciences 100 (2017) 1–9
E.S. Um et al.
Fig. 9. Comparisons of FE and FD solutions for the four seabed conductivities at selected sources. A number inside the parenthesis is a conductivity value for the seabed.
σ FE = Mσ FD ; (i , j ) element of Mij =
FD element, respectively. The intersection operator ∩ computes the overlapping volume of the FE and FD cell if they intersect. In summary, FEDO sets up wireframes that honor the given seafloor topography. Its central regions are discretized with fine elements and are surrounded with large elements. The use of effective receiver sizes lead to further efficient FE meshes. FEDO finally maps conductivity from FD to FE meshes in a volume-average sense. The
(9a)
1 × (viFE ∩ vjFD ). viFE
(9b)
where vectors σ FE and σ FD contain conductivity attributes of the FE and FD models, respectively. Matrix M is a mapping operator from FD to FE meshes. viFE and vjFD are the volume of the ith FE element and the jth 6
Computers & Geosciences 100 (2017) 1–9
E.S. Um et al.
sequence of the tasks automatically generates an FE model without user intervention. 6. Examples We demonstrate FEDO for modeling CSEM over three seabed models. The FE solutions to the three models are compared with reference FD solutions. For the three CSEM models, a source frequency is 0. 25 Hz. Receivers in the FD models are always horizontal, whereas receivers in the FE models are aligned with the true seafloor slope. 6.1. Complex seafloor topography model We consider a homogeneous seabed model that has the complex topography. For comparison, Fig. 7 shows the seafloor topography. By using fine rectangular cells (i.e. 100-by-100-by-20 m), the FD-based seafloor topography closely resembles the FE topography. 10 m long receivers (actual receiver length) are placed with 500 m spacing on the seafloor. The conductivity of the seawater is 3.33 S/m. We consider four conductivity values of the seabed: 3.33, 2.00, 1.43 and 1 S/m. The four FE models are discretized using tetrahedral elements whose maximum size is set to 90% of δ of the seabed. A central portion of the FE meshes for 1.43 S/m seabed is presented in Fig. 8. The maximum element sizes for the seawater (region 2 and 3) and the seabed (region 4 and 5) are set to 400 m and 750 m, respectively. The seawater-seabed interface along the survey line is discretized using elements whose maximum size is set to 50 m. This length is a conservatively determined effective receiver size. Thus, the mesh cost for the topography and receivers is shared. To attenuate the unwanted boundary effects, the size of the computational domain is set to 28-by-20–28 km in the x-, y- and zdirection, respectively. The FE model consists of 676,398 elements and 776,921 unknowns. Fig. 9 shows the comparison of the FE and FD solutions for the four seabed models at selected source positions (x=45, 50 and 55 km). The two solutions show reasonable agreements, demonstrating the accuracy of FEDO for the complex seafloor topography. However, Fig. 9 also reveals effects of the different seafloor geometry and the tilted seafloor receivers between the two models. For example, consider the FE and FD solutions with the source at x=45 km (Fig. 9a). As the seabed conductivity decreases, the differences of the two solutions gradually increases from x=48–55 km. In contrast, the two solutions for the four seabed models have nearly identical airwave responses at long sourcereceiver offsets (x=55–58 km) because the airwave responses are mainly controlled by the seawater. Thus, the effects of the different seafloor geometry are localized at the intermediate source-receiver
Fig. 10. Scalability of the FE method with MUMPS.
Fig. 11. The complex seabed model with the flat seafloor. The color bar represents the conductivity on the log scale. The log scale is used in all models in this paper. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Fig. 12. The FE complex seabed model with the flat seafloor.
Fig. 13. FD and FE solutions to the complex seabed model. (a) Amplitudes. (b) Phases. The arrows at x=−4 and 4 km indicates source positions.
7
Computers & Geosciences 100 (2017) 1–9
E.S. Um et al.
Fig. 14. (a) SEG salt model at y=5 km. (b) Its seafloor. The solid lines represent survey lines. Asterisks and circles represent receivers and sources, respectively. The broken line represents the boundary of the salt.
solved due to the limited scalability of the factorization algorithm. Instead, it is desirable to generate multiple small FE meshes that are optimized for a subset of sources and receivers and feed each FE model to a subset of processors on massively-parallel computer architecture. The effectiveness of this approach has been successfully demonstrated in FD modeling (Yang et al., 2014; Yang and Oldenburg, 2016). This can also improve the overall efficiency of the FE solution processes when many sources are employed. Such an example will be demonstrated later when survey lines are designed in a 3D grid pattern.
6.2. Complex seabed model with flat seafloor Fig. 15. The FE SEG salt model.
The second model has the complex distribution of conductivity attributes with the flat seafloor (Fig. 11). FEDO is applied to the FD model, generating the FE model (Fig. 12). The FD model consists of uniform 100 m cells in the x-, y- and z-direction. For the FE model, the maximum element sizes of region 3 (seawater) and region 5 (seabed) are set to 500 and 200 m, respectively. The maximum element size for region 4 (seabed) is set to 1700 m ( < δ of 4 Ω-m medium at 0.25 Hz). As tetrahedral elements of region 5 have nearly uniform 200 m edges, its volume is roughly equal to that of the 100 m FD cells. Thus, FD and FE get matched in the seabed in a volumetric sense. The FE model consists of 2,735,144 elements and 3,170,075 unknowns. The FE model illustrates the efficiency of partitioning the computa-
offsets. Fig. 10 shows the scalability of the FE solution with MUMPS 5.0 (Amestoy et al., 2001) that performs parallel factorization and backand-forward substitution. The computation is carried out using Intel 2.66 GHz processors in a distributed memory system. The factorization is the most expensive part of the FE solutions and shows only moderate scalability up to 96 cores. Accordingly, it is beneficial to share LU with multiple right-hand-side vectors (e.g. 15 sources). However, it is not recommended to form a single very large system of FE equations for too many sources because such a large system cannot be effectively
Fig. 16. CSEM responses to SEG salt model along y=5 km. (a) Amplitudes. (b) Phases.
8
Computers & Geosciences 100 (2017) 1–9
E.S. Um et al.
resulting FE meshes precisely honor coordinates of the seafloor topography. Second, the 3D distribution of seabed conductivity is mapped from rectangular FD cells to tetrahedral FE meshes in a volumetric average sense, avoiding interactive CAD works or adaptive meshing overhead. Third, the mesh refinement required around receivers is relaxed by introducing an effective receiver size. Fourth, by partitioning a seabed model into multiple regions and applying different element sizes, FEDO economically discretizes a large computational domain. Fifth, FEDO can generate multiple small FE problems that are optimized for a subset of sources and receivers and make it possible to compute them simultaneously using multiple subsets of processors. This divide and conquer approach helps overcome the limited scalability of the parallel direct solver.
tional domain and applying different element sizes to different regions. About 80% of the total elements are used to discretize the relatively small volume of region 5 where the complex distribution of conductivity attributes is mapped from FD to FE. A 200 m long source is placed 100 m above the seafloor at x= ± 4000 m. 10 m long receivers (actual receiver length) are placed with 500 m spacing on the seafloor. For the FE model, the effective receiver size is set to 100 m. The resulting FE and FD solutions are compared in Fig. 13. The two solutions show overall good agreements, demonstrating the accuracy of FEDO. 6.3. SEG salt model The third model is based on the FD-based SEG salt velocity model (Aminzadeh et al., 1997) that has both complex seabed and seafloor structures. The velocity model is translated into a conductivity model (Fig. 14) as done in Um et al. (2014). FEDO generates a corresponding FE model (Fig. 15). For region 3 and 5, the maximum element sizes are set to 500 and 200 m, respectively. For the other regions, the maximum element sizes are determined by δ. The FE seafloor topography is defined by sampling the original coordinates of the FD seafloor topography. We consider hundred sources and two hundred receivers in a grid pattern (Fig. 14b). Ten 10 km long survey lines are apart each other by 1 km from y=1–10 km. Along each line, ten x-oriented electric dipole sources are placed at 1 km interval from x=1–10 km. Twenty xoriented receivers are placed at 500 m interval from x=1–10 km. The lengths of the actual sources and receivers are 100 and 10 m, respectively. In the FE model, the effective receiver size is set to 50 m. FEDO constructs ten independent FE meshes such that each FE mesh can accommodate one survey line (i.e. ten sources and twenty receivers). The ten FE meshes are assigned to ten subsets of processors on a parallel computer. By assigning small problems to subsets of processors and solving them simultaneously, we overcome the limited scalability of the parallel direct solver. On average, one FE mesh consists of 1,206,850 elements and 1,381,135 unknowns. Using 1200 cores (ten groups of ten twelve-core AMD 2.1 GHz processors), the hundred FE solutions are completed in 319 s. The FE and FD solutions along a selected survey line are shown in Fig. 16. Around the sources, the two solutions show large differences because they are computed in the total field mode. Except them, the two solutions are in the right ballpark, but the comparison is not as good as those of the previous two examples. We believe that subtle differences between the two solutions come from the different seafloor topography and different receiver orientations, demonstrating the sensitivity of CSEM to the two factors.
Acknowledgements This work was supported in part by National Natural Science Foundation of China (Grant No. 61303003 and 41374113). S.S.K acknowledges support from Korea National Research Foundation of the Ministry of Science, ICT and Future Planning (Grant No. NRF2013R1A1A1076071 and NRF-2012M2A8A5007440). This work benefited from the constructive comments of Colin Farquharson and three anonymous reviewers. References Amestoy, P.R., Duff, I.S., Koster, J., L’Excellent, J.-Y., 2001. A fully asynchronous multifrontal solver using distributed dynamic scheduling. SIAM J. Matrix Anal. Appl. 23, 15–41. Aminzadeh, F., Brac, J., Kunz, T., 1997. 3-D Salt and Overthrust Models. SEG/EAGE, Tulsa, Oklahoma. Fu, H., Wang, Y., Um, E.S., Fang, J., Wei, T., Huang, X., Yang, G., 2015. A parallel finiteelement time-domain method for transient electromagnetic simulation. Geophysics 80, E213–E224. Jin, J., Riley, D., 2008. Finite Element Analysis of Antennas and Arrays. Wiley-IEEE Press, Hoboken, New Jersey. Key, K., Ovall, J., 2011. A parallel goal-oriented adaptive finite element method for 2.5-D electromagnetic modelling. Geophys. J. Int. 186, 137–154. Li, Y., Key, K., 2007. 2D marine controlled-source electromagnetic modeling, part 1: an adaptive finite element algorithm. Geophysics 72, WA51–WA62. Newman, G.A., Alumbaugh, D.L., 1995. Frequency domain modeling of airborne electromagnetic responses using staggered finite differences. Geophys. Prospect. 43, 1021–1042. Schwarzbach, C., Börner, R.-U., Spitzer, K., 2011. Three-dimensional adaptive higher order finite element simulation for geo-electromagnetics-a marine CSEM example. Geophys. J. Int. 187, 63–74. Si, H., 2006. TetGen: A Quality Tetrahedral Mesh Generator and 3D Delaunay Trianulator. Weierstrass Institute for Applied Analysis and Stochastics, Berlin, Germany. Um, E., Commer, M., Newman, G., 2013. Efficient pre-conditioned iterative solution strategies for the electromagnetic diffusion in the Earth: finite-element frequencydomain approach. Geophys. J. Int. 193, 1460–1473. Um, E., Commer, M., Newman, G., 2014. A strategy for coupled 3D imaging of largescale seismic and electromagnetic data sets: application to subsalt imaging. Geophysics 79, ID1–ID13. Yang, D., Oldenburg, D.W., Haber, E., 2014. 3-D inversion of airborne electromagnetic data parallelized and accelerated by local mesh and adaptive soundings. Geophys. J. Int. 196, (1942-1507). Yang, D., Oldenburg, D.W., 2016. Survey decomposition: a scalable framework for 3D controlled-source electromagnetic inversion. Geophysics 81, E69–E87.
7. Conclusions We have analyzed numerical characteristics of tetrahedral meshes for CSEM modeling. A major characteristic of tetrahedral meshes is that the tetrahedral mesh density required for accurate solutions does not strongly depend on the maximum source-receiver offset because its numerical dispersion errors oscillate. Based on the analysis, a sequence of mesh generation tasks is casted into FEDO and is automated. Benefits of FEDO are summarized as follows. First, unlike FD grids, the
9
Contents lists available at ScienceDirect
Computers & Geosciences journal homepage: www.elsevier.com/locate/cageo
crossmark
A tetrahedral mesh generation approach for 3D marine controlled-source electromagnetic modeling Evan Schankee Uma, Seung-Sep Kimb, Haohuan Fuc,d,
⁎
a
Affiliate of Earth and Environmental Sciences, Lawrence Berkeley National Laboratory, CA, USA Geology and Earth Environmental Sciences, Chungnam National University, Daejeon, South Korea Tsinghua University, Ministry of Education Key Laboratory for Earth System Modeling and Center for Earth System Science, Beijing, China d National Supercomputing Center in Wuxi, Wuxi, China b c
A R T I C L E I N F O
A BS T RAC T
Keywords: Marine electromagnetic geophysics Finite element modeling Parallel computing
3D finite-element (FE) mesh generation is a major hurdle for marine controlled-source electromagnetic (CSEM) modeling. In this paper, we present a FE discretization operator (FEDO) that automatically converts a 3D finitedifference (FD) model into reliable and efficient tetrahedral FE meshes for CSEM modeling. FEDO sets up wireframes of a background seabed model that precisely honors the seafloor topography. The wireframes are then partitioned into multiple regions. Outer regions of the wireframes are discretized with coarse tetrahedral elements whose maximum size is as large as a skin depth of the regions. We demonstrate that such coarse meshes can produce accurate FE solutions because numerical dispersion errors of tetrahedral meshes do not accumulate but oscillates. In contrast, central regions of the wireframes are discretized with fine tetrahedral elements to describe complex geology in detail. The conductivity distribution is mapped from FD to FE meshes in a volume-averaged sense. To avoid excessive mesh refinement around receivers, we introduce an effective receiver size. Major advantages of FEDO are summarized as follow. First, FEDO automatically generates reliable and economic tetrahedral FE meshes without adaptive meshing or interactive CAD workflows. Second, FEDO produces FE meshes that precisely honor the boundaries of the seafloor topography. Third, FEDO derives multiple sets of FE meshes from a given FD model. Each FE mesh is optimized for a different set of sources and receivers and is fed to a subgroup of processors on a parallel computer. This divide and conquer approach improves the parallel scalability of the FE solution. Both accuracy and effectiveness of FEDO are demonstrated with various CSEM examples.
1. Introduction In the past decade, 3D finite-element (FE) solutions have been widely used in marine controlled-source electromagnetic (CSEM) modeling (e.g. Key and Ovall, 2011; Schwarzbach et al., 2011). In contrast to finite-difference (FD) solutions that approximate noncoordinate-conforming structures with small rectangular stair steps, FE uses geometry-conforming tetrahedral meshes and precisely represents complex seafloor topography. However, the advantage of FE over FD comes with extra complication. It is considered difficult to iteratively solve a system of FE equations. A system matrix with tetrahedra is unstructured and not diagonally dominant. A simple Jacobi preconditioner used in FD solutions does not ensure convergence of FE Krylov solutions. FE solutions require numerically expensive preconditioners such as incomplete factorization (Um et al., 2013) and others. It is also difficult to
⁎
find a robust preconditioner suitable to a wide range of EM problems. Accordingly, direct solvers are often the method of choice for FE solutions despite their large memory requirements and low parallel scalabilities (Fu et al., 2015). It is also considered difficult and time-consuming to create 3D FE meshes. 3D FE mesh generation requires good literacy in computeraided design (CAD) software that may have a steep learning curve. To mitigate the difficulty, adaptive FE methods have been introduced (Li and Key, 2007; Key and Ovall, 2011; Schwarzbach et al., 2011). They start with coarse meshes and successively refine meshes until a required tolerance is met. However, in either use of interactive CAD or automated adaptive refinement methods, it is nontrivial to generate reliable and efficient 3D meshes for a complex multi-scale model. The goal of this paper is to present a 3D FE discretization operator (FEDO) that automatically generates reliable and efficient tetrahedral meshes for marine CSEM. We develop strategies for directly creating
Corresponding author. E-mail addresses: [email protected] (E.S. Um), [email protected] (S.-S. Kim), [email protected] (H. Fu).
http://dx.doi.org/10.1016/j.cageo.2016.11.007 Received 4 August 2016; Received in revised form 17 November 2016; Accepted 20 November 2016 Available online 21 November 2016 0098-3004/ © 2016 Elsevier Ltd. All rights reserved.
Computers & Geosciences 100 (2017) 1–9
E.S. Um et al.
complete FE meshes without adaptive refinement. The development requires understanding basic properties of tetrahedral meshes for diffusive EM modeling, which are often not obvious when CAD or adaptive refinement methods are employed. Ultimately, we aim to demonstrate that FEDO can quickly derive accurate and efficient FE meshes from a given seabed model without computational overhead associated with adaptive mesh refinement or interactive CAD works. The remainder of this paper is organized as follows. We first review the FE formulation for CSEM. Its numerical dispersion characteristics are examined to determine proper tetrahedral sizes. To avoid excessive mesh refinement around receivers, we introduce an effective receiver size. We follow this by presenting a matrix that maps conductivity from FD model to FE simulation meshes. The three components above are casted into FEDO. Finally, we apply FEDO to complex offshore models, compute their CSEM solutions and demonstrate its accuracy and efficiency.
Table 1 δ and λ in the seawater.
(1)
where e(r) is the electric field at position r, Js(r) is an electric source at an angular frequency ω , μo is the magnetic permeability of free space (4π×10-7 H/m), and σ is the conductivity. The development of the equivalent weak statement of Eq. (1) requires the multiplication of Eq. (1) by the edge basis function and the integration over the model domain of V, resulting in e
1 nie(r)⋅( ∇ × ∇ × ee(r) + l j ωσeee(r) + l j ω Js(r))dV = 0. μo th
(2)
nie(r)
The superscript e denotes the e tetrahedral element, with i varying from 1 to 6 is a set of edge basis functions. V e is the volume of the eth element. nie(r) is also chosen as the basis. Thus, the electric field at a point inside or on a given element is expanded as 6
ek(r) =
6
∑ ekj (r) = ∑ ujk nkj (r), j =1
j =1
(3)
is the unknown amplitude of the electric field on the jth edge where of the kth element. Accordingly, the FE method is the second-order accurate (Jin and Riley, 2008). Substituting Eq. (3) into Eq. (2) and using the homogeneous Dirichlet boundary conditions yield
ujk
(Ae + μo ωjl Be)ue = − μo ωiˆ se
(4)
where
(i , j ) element of Ae =
∭V
(i , j ) element of Be =
∭V
i element of se =
∭V
e
ue = [ u1e u 2e ... u6e].
e
e
∇ × nie(r)⋅∇ × nej(r)dV ,
(5)
nie(r)⋅σ e nej(r)dV ,
(6)
nie(r)⋅Js(r)dV ,
λ (m)
0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50
873.0 712.8 617.3 552.1 504.0 466.6 436.5 411.5 390.4
5484.9 4478.4 3878.4 3469.0 3166.7 2931.8 2742.5 2585.6 2452.9
To determine tetrahedral mesh density required for accurate solutions, we examine numerical dispersion characteristics for a homogeneous whole-space seawater (3.33 S/m) model. We are interested in efficient discretization of the seawater because the seawater can be a large portion of a CSEM model and requires the smallest elements due to its highest conductivity. The seawater model is 20×10×10 km in the x-, y- and z-direction, respectively. The lower left corner of the model is at (−5, −5, 5 km). A 250 m long x-oriented electric source is placed at (0, 0, 0 km). 20 m long x-oriented electric receivers are placed with 1 km spacing from x=1–10 km. We consider nine source frequencies from 0.1 to 0.5 Hz with 0.05 Hz interval. Skin depth (δ) and wavelength (λ) of the seawater are presented in Table 1. The model (Fig. 1) is discretized using tetrahedral elements whose vary from 10 m (near sources and receivers) to 400 m (most areas). This is slightly smaller than δ at 0.3 Hz. Note that the FE meshes are refined around not only the source but also the receivers. This is one of the major differences between FD and FE. The mesh refinement near the receivers will be discussed in the next section. The transition from small to large edges is controlled by the growth factor defined as the maximum rate at which the edge size can grow. The growth factor is typically 1.5–2.0 from one edge to the next. Amplitudes, relative amplitude errors and phases of FE solutions are plotted against analytic solutions in Fig. 2a-c, respectively. The relative errors are defined by ‖(numerical solution–analytic solution)/ analytic solution)‖. It is assumed that over 10 km source-receiver offset, 5% amplitude errors are accurate enough. Based on the criterion, the FE meshes support the lowest four frequencies: 0.10, 0.15, 0.20 and 0.25 Hz. At 0.10 Hz, the boundary effects appear at x=8 km. To highlight the numerical dispersion characteristics of the FE meshes, we also repeat the same experiments with a 2nd-order accurate FD method (Newman and Alumbaugh, 1995). The FD grid size is set to 400 m as same as that in the FE meshes. However, FD does not produce accurate solutions at all nine frequencies (not included here). We reduce the FD grid to 100 m and analyze amplitudes, their relative errors and phases (Fig. 2d–f). Due to their geometric simplicity, the FD grids are not presented here. The FD grids support the lowest four frequencies (Fig. 2e). At 0.10 Hz, the boundary
Since the discretization density required for accurate EM solutions is directly related with an FE formulation, we briefly describe the total field FE formulation (Um et al., 2013) used here. The electric-field diffusion equation is given by
∭V
δ (m)
3. Numerical dispersion analysis
2. Finite element formulation
1 ∇ × ∇ × e(r) + l j ωσe(r) + l j ω Js(r) = 0, μo
f (Hz)
(7) (8)
Eq. (4) is considered local as it results from integration over each element. The local systems from all elements are assembled into a single global system of equations. The global system is solved using a parallel direct solver.
Fig. 1. The FE meshes for the seawater model. Blue and red lines represent source and receivers, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
2
Computers & Geosciences 100 (2017) 1–9
E.S. Um et al.
0.10 Hz 0.15 Hz 0.20 Hz 0.25 Hz 0.30 Hz 0.35 Hz 0.40 Hz 0.45 Hz 0.50 Hz
0.10 Hz 0.15 Hz 0.20 Hz 0.25 Hz 0.30 Hz 0.35 Hz 0.40 Hz 0.45 Hz 0.50 Hz
) Fig. 2. The FE (the first row) and FD (second row) solutions to the seawater model and their comparison with analytic solutions. The solid lines with a square are phases of FE/FD solutions; the broken lines with a circle are those of analytic solutions.
capture rapidly-decaying EM fields. It is our experiences that some refinements are often required in both total and scattered field solutions. In FD, the fine grids propagate throughout the entire computational domain, resulting in a large number of cells. In contrast, FE allows local mesh refinements. However, a major drawback of FE is the fact that it also requires fine elements near receivers. By using the first-order edge basis functions, FE determines an average electric-field vector along each edge of tetrahedral elements. The averaged electricfield vectors are used to interpolate the electric field at a receiver through Eq. (3). If the edges are much longer than a finite-length receiver, the interpolation would be erroneous. Consequently, it is necessary to refine meshes near receivers. A larger number of receivers results in a larger FE problem size. This is a major disadvantage of the FE method presented here. In this section, we examine local mesh refinements around receivers for accurate FE solutions at 0.25 Hz and present a simple method for determining an economic mesh size around receivers. To do this, the seawater model is discretized with eight levels of receiver mesh refinement as summarized in Table 2. For the first seven levels, 20, 40, 80, 120, 160, 200 and 250 m long x-oriented edges are placed at the midpoint of 20 m long receivers. For the eighth level, the 250 m long edges are used again near receivers, but their directions are not aligned with the receiver direction, but are randomly oriented. Fig. 3 shows FE
effect appears at x=8 km. This distance is similar to that of the FE meshes. As demonstrated, FE meshes can be coarser than FD meshes. FD produces accurate solutions when δ is discretized with 4–7 nodes. As shown in Fig. 2e, the mesh density required for accurate FD solutions strongly depends on a source-receiver offset because FD numerical dispersion errors accumulate. In contrast, FE solutions are accurate enough if local element sizes are slightly smaller than local δ. Note that the FE numerical dispersion errors oscillate as the source-receiver offset increases. This is because the FE meshes are not structured but randomly oriented; the errors can be either positive or negative, canceling each other out rather than being accumulated (Jin and Riley, 2008). Thus, the required FE element size does not strongly depend on the maximum source-receiver offset. δ serves as a robust criterion for determining the FE element size. In this paper, the maximum local FE element size is chosen as 90% of local δ for safety. The element size determined by δ is often much larger than CSEM resolution. For example, δ for the seabed of 0.3 S/m is about 1.74 km at 0.25 Hz although CSEM aims to resolve several ten meters thick structures. This observation suggests a basic strategy of determining element sizes for CSEM. A CSEM model is partitioned into two conceptual regions. One region represents a part of seabed where detailed geological structures are described. In this region, the element sizes should be small enough such that fine seabed structures can be described. In the other regions (e.g. the other part of the seabed, the seawater and the air), the local δ criteria determine maximum local element sizes, leading to the economic discretization. Later in this paper, we will further discuss partitioning a CSEM model and its region-by-region discretization. In the next section, we examine the FE mesh refinement required around sources and receivers and present the relaxation of the requirement by introducing an effective receiver size.
Table 2 Mesh refinement levels according to a receiver edge size and resultant total number of elements required to discretize the computational domain.
4. Effective receiver size Around sources, both FD and FE require fine discretization to 3
Mesh refinement level
Edge size for receiver (m)
# of elements
1 2 3 4 5 6 7 8
20 40 80 120 160 200 250 250
322,076 310,831 299,571 294.913 293,244 290,568 287.914 286.749
Computers & Geosciences 100 (2017) 1–9
E.S. Um et al.
meshes at a few selected refinement levels. The FE solutions from each FE mesh are compared against the analytic solution where the actual receiver length is set to 20 m (Fig. 4a–c). As the receiver edge size increases, the relative errors oscillate with larger amplitude and eventually exceed an acceptable level (i.e. 5%). However, Fig. 4a–c suggests that FE solutions can be accurate enough although FE meshes do not precisely resolve the actual receiver length. For example, when the receiver edge size is set to 160 m, the relative errors are bounded within 5%. This relatively low error is due to the fact that λ (i.e. 3469 m) is large enough compared with the actual receiver size (20 m). This result is also verified by comparing the analytic solutions of different actual receiver sizes as shown in Fig. 4d–f where the 20 m long actual receiver can be replaced even with a 250 m long actual receiver. It is also important to place edges of elements at the centers of receivers and align them along the receiver direction. This alignment helps the receiver edge remove the unwanted influence from the other five edges via Eq. (3). For example, when the receiver edge size is set to 250 m, relative errors with and without the alignment are noticeably different (Fig. 4b). It is our experience that in CSEM modeling, it is not necessary to resolve the actual receiver size. We can find an ‘effective’ receiver size that is typically several factors larger than the actual size. As done in Fig. 4e, one can determine an effective receiver size using analytical solutions. Fig. 5 demonstrates the benefit of using an effective receiver size in terms of a problem size. When the FE meshes resolve the actual receiver size (20 m), the total number of elements quickly increases with the number of receivers. In contrast, when an effective receiver size (160 m) is used, the increase is modest. The use of an effective receiver size is particularly recommended when a number of receivers are placed in a grid pattern.
Fig. 3. FE meshes with different refinement levels around receivers. (a) Level 3, (b) Level 5. (c) Level 7.
Fig. 4. The first row: (a) The FE solutions to the seawater model with a receiver length varying from 20 to 250 m and the analytic solution (20 m long receiver, green broken line). (b) The relative amplitude errors of the FE solutions with respect to the analytic solution. (c) Phase plots of the FE and analytic solutions. The second row: (d) The analytic solutions to the seawater model with various receiver lengths. (e) The relative amplitude differences of (d) when the analytic solutions with various receiver lengths are compared with that of 20 m. (f) The phase plots of (d). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
4
Computers & Geosciences 100 (2017) 1–9
E.S. Um et al.
Fig. 5. The number of elements required for the seawater model is plotted as a function of number of receivers when a receiver is resolved by its true (20 m) and effective (160 m) length.
Fig. 8. The FE seafloor topography meshes. The red lines represent sources. The seawater-seabed interface is in blue. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
region 3 and 5. If sources or receivers are on the seafloor, their coordinates are augmented to the coordinates of the seafloor topography. Region 2 and 4 surround region 3 and 5 and provide ample space (e.g. λ) through which EM fields are attenuated. The δ criterion determines the maximum element sizes for region 1, 2 and 4. Element sizes are nearly uniform in these regions. In contrast, region 3 can have a range of element sizes. The bottom of region 3 consists of small elements and defines the irregular seafloor topography. Region 3 also includes a number of small line segments for sources and receivers. For the efficient discretization of region 3, element sizes gradually grow outward from the bottom and source/ receiver positions and are bounded by δ. Because geological structures are mapped from FD to FE meshes as will be described below, element sizes of region 5 should not only satisfy δ criterion but also be roughly equivalent to those of FD meshes in a volumetric sense. For example, if the size of a regular tetrahedral element is 200 m, its volume (9.43×105 m3) is similar to a 100 m FD cube. FEDO utilizes TetGen (Si, 2006), a tetrahedral mesh generation algorithm to discretize the wireframes. Note that the wireframes and their discretization honor only the seafloor interface. The seafloor is the most important interface to CSEM modeling because the receivers can be directly placed on the interface. By precisely discretizing the seafloor geometry, we ensure the correctness of local EM interface conditions and local EM interaction. In contrast to seismic modeling, it is our experiences that diffusive CSEM is much less sensitive to the detailed geometry of other interfaces below the seafloor. Thus, we do not consider other interfaces as critical constraints that need to be explicitly honored in tetrahedral meshing. Instead, they are realized by mapping conductivity values from FD to FE meshes in a volume average sense as described below. After the discretization, FEDO maps conductivity attributes from FD to FE meshes. Region 1, 2 and 3 are assumed to be homogeneous. Their conductivity attributes directly come from the FD model. In contrast, each tetrahedral element of region 4 and 5 is divided into N small sub-elements. Their centroidal conductivities are read from the FD model and are averaged in a volumetric sense. The average is assigned to the tetrahedral element. When some sub-elements of a seabed tetrahedral element reside in the seawater of the FD model, their centroidal conductivities are excluded from the averaging processes. When all sub-elements reside in the seawater of the FD model due to large discrepancies between FD and FE topography, conductivities of adjacent tetrahedral elements are used to interpolate the conductivity of the tetrahedral element. The mapping processes are casted into
Region 1: the air
Region 2: sea
Region 3: sea Source Receiver
Region 4: seabed
Region 5: seabed
Fig. 6. The 3D wireframes.
Fig. 7. The comparison of the FE and FD seafloor topography models.
5. FE discretization operator Based on the analysis above, we present an FEDO for 3D marine CSEM modeling. FEDO requires (1) information about CSEM configurations, (2) the description of seabed conductivities in an FD format and (3) coordinates of seafloor topography. The first two items are those used for FD modeling. The input conductivity model is in the FD format because porosity and fluid saturation models used for constructing a conductivity model are often in the FD format. In contrast, the third item is the new information. FEDO performs three tasks described below. First, FEDO sets up the wireframes of a background seabed model (Fig. 6). The wireframes consist of five internal regions. The coordinates of the seafloor topography is given in a grid pattern and defines the interface of 5
Computers & Geosciences 100 (2017) 1–9
E.S. Um et al.
Fig. 9. Comparisons of FE and FD solutions for the four seabed conductivities at selected sources. A number inside the parenthesis is a conductivity value for the seabed.
σ FE = Mσ FD ; (i , j ) element of Mij =
FD element, respectively. The intersection operator ∩ computes the overlapping volume of the FE and FD cell if they intersect. In summary, FEDO sets up wireframes that honor the given seafloor topography. Its central regions are discretized with fine elements and are surrounded with large elements. The use of effective receiver sizes lead to further efficient FE meshes. FEDO finally maps conductivity from FD to FE meshes in a volume-average sense. The
(9a)
1 × (viFE ∩ vjFD ). viFE
(9b)
where vectors σ FE and σ FD contain conductivity attributes of the FE and FD models, respectively. Matrix M is a mapping operator from FD to FE meshes. viFE and vjFD are the volume of the ith FE element and the jth 6
Computers & Geosciences 100 (2017) 1–9
E.S. Um et al.
sequence of the tasks automatically generates an FE model without user intervention. 6. Examples We demonstrate FEDO for modeling CSEM over three seabed models. The FE solutions to the three models are compared with reference FD solutions. For the three CSEM models, a source frequency is 0. 25 Hz. Receivers in the FD models are always horizontal, whereas receivers in the FE models are aligned with the true seafloor slope. 6.1. Complex seafloor topography model We consider a homogeneous seabed model that has the complex topography. For comparison, Fig. 7 shows the seafloor topography. By using fine rectangular cells (i.e. 100-by-100-by-20 m), the FD-based seafloor topography closely resembles the FE topography. 10 m long receivers (actual receiver length) are placed with 500 m spacing on the seafloor. The conductivity of the seawater is 3.33 S/m. We consider four conductivity values of the seabed: 3.33, 2.00, 1.43 and 1 S/m. The four FE models are discretized using tetrahedral elements whose maximum size is set to 90% of δ of the seabed. A central portion of the FE meshes for 1.43 S/m seabed is presented in Fig. 8. The maximum element sizes for the seawater (region 2 and 3) and the seabed (region 4 and 5) are set to 400 m and 750 m, respectively. The seawater-seabed interface along the survey line is discretized using elements whose maximum size is set to 50 m. This length is a conservatively determined effective receiver size. Thus, the mesh cost for the topography and receivers is shared. To attenuate the unwanted boundary effects, the size of the computational domain is set to 28-by-20–28 km in the x-, y- and zdirection, respectively. The FE model consists of 676,398 elements and 776,921 unknowns. Fig. 9 shows the comparison of the FE and FD solutions for the four seabed models at selected source positions (x=45, 50 and 55 km). The two solutions show reasonable agreements, demonstrating the accuracy of FEDO for the complex seafloor topography. However, Fig. 9 also reveals effects of the different seafloor geometry and the tilted seafloor receivers between the two models. For example, consider the FE and FD solutions with the source at x=45 km (Fig. 9a). As the seabed conductivity decreases, the differences of the two solutions gradually increases from x=48–55 km. In contrast, the two solutions for the four seabed models have nearly identical airwave responses at long sourcereceiver offsets (x=55–58 km) because the airwave responses are mainly controlled by the seawater. Thus, the effects of the different seafloor geometry are localized at the intermediate source-receiver
Fig. 10. Scalability of the FE method with MUMPS.
Fig. 11. The complex seabed model with the flat seafloor. The color bar represents the conductivity on the log scale. The log scale is used in all models in this paper. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Fig. 12. The FE complex seabed model with the flat seafloor.
Fig. 13. FD and FE solutions to the complex seabed model. (a) Amplitudes. (b) Phases. The arrows at x=−4 and 4 km indicates source positions.
7
Computers & Geosciences 100 (2017) 1–9
E.S. Um et al.
Fig. 14. (a) SEG salt model at y=5 km. (b) Its seafloor. The solid lines represent survey lines. Asterisks and circles represent receivers and sources, respectively. The broken line represents the boundary of the salt.
solved due to the limited scalability of the factorization algorithm. Instead, it is desirable to generate multiple small FE meshes that are optimized for a subset of sources and receivers and feed each FE model to a subset of processors on massively-parallel computer architecture. The effectiveness of this approach has been successfully demonstrated in FD modeling (Yang et al., 2014; Yang and Oldenburg, 2016). This can also improve the overall efficiency of the FE solution processes when many sources are employed. Such an example will be demonstrated later when survey lines are designed in a 3D grid pattern.
6.2. Complex seabed model with flat seafloor Fig. 15. The FE SEG salt model.
The second model has the complex distribution of conductivity attributes with the flat seafloor (Fig. 11). FEDO is applied to the FD model, generating the FE model (Fig. 12). The FD model consists of uniform 100 m cells in the x-, y- and z-direction. For the FE model, the maximum element sizes of region 3 (seawater) and region 5 (seabed) are set to 500 and 200 m, respectively. The maximum element size for region 4 (seabed) is set to 1700 m ( < δ of 4 Ω-m medium at 0.25 Hz). As tetrahedral elements of region 5 have nearly uniform 200 m edges, its volume is roughly equal to that of the 100 m FD cells. Thus, FD and FE get matched in the seabed in a volumetric sense. The FE model consists of 2,735,144 elements and 3,170,075 unknowns. The FE model illustrates the efficiency of partitioning the computa-
offsets. Fig. 10 shows the scalability of the FE solution with MUMPS 5.0 (Amestoy et al., 2001) that performs parallel factorization and backand-forward substitution. The computation is carried out using Intel 2.66 GHz processors in a distributed memory system. The factorization is the most expensive part of the FE solutions and shows only moderate scalability up to 96 cores. Accordingly, it is beneficial to share LU with multiple right-hand-side vectors (e.g. 15 sources). However, it is not recommended to form a single very large system of FE equations for too many sources because such a large system cannot be effectively
Fig. 16. CSEM responses to SEG salt model along y=5 km. (a) Amplitudes. (b) Phases.
8
Computers & Geosciences 100 (2017) 1–9
E.S. Um et al.
resulting FE meshes precisely honor coordinates of the seafloor topography. Second, the 3D distribution of seabed conductivity is mapped from rectangular FD cells to tetrahedral FE meshes in a volumetric average sense, avoiding interactive CAD works or adaptive meshing overhead. Third, the mesh refinement required around receivers is relaxed by introducing an effective receiver size. Fourth, by partitioning a seabed model into multiple regions and applying different element sizes, FEDO economically discretizes a large computational domain. Fifth, FEDO can generate multiple small FE problems that are optimized for a subset of sources and receivers and make it possible to compute them simultaneously using multiple subsets of processors. This divide and conquer approach helps overcome the limited scalability of the parallel direct solver.
tional domain and applying different element sizes to different regions. About 80% of the total elements are used to discretize the relatively small volume of region 5 where the complex distribution of conductivity attributes is mapped from FD to FE. A 200 m long source is placed 100 m above the seafloor at x= ± 4000 m. 10 m long receivers (actual receiver length) are placed with 500 m spacing on the seafloor. For the FE model, the effective receiver size is set to 100 m. The resulting FE and FD solutions are compared in Fig. 13. The two solutions show overall good agreements, demonstrating the accuracy of FEDO. 6.3. SEG salt model The third model is based on the FD-based SEG salt velocity model (Aminzadeh et al., 1997) that has both complex seabed and seafloor structures. The velocity model is translated into a conductivity model (Fig. 14) as done in Um et al. (2014). FEDO generates a corresponding FE model (Fig. 15). For region 3 and 5, the maximum element sizes are set to 500 and 200 m, respectively. For the other regions, the maximum element sizes are determined by δ. The FE seafloor topography is defined by sampling the original coordinates of the FD seafloor topography. We consider hundred sources and two hundred receivers in a grid pattern (Fig. 14b). Ten 10 km long survey lines are apart each other by 1 km from y=1–10 km. Along each line, ten x-oriented electric dipole sources are placed at 1 km interval from x=1–10 km. Twenty xoriented receivers are placed at 500 m interval from x=1–10 km. The lengths of the actual sources and receivers are 100 and 10 m, respectively. In the FE model, the effective receiver size is set to 50 m. FEDO constructs ten independent FE meshes such that each FE mesh can accommodate one survey line (i.e. ten sources and twenty receivers). The ten FE meshes are assigned to ten subsets of processors on a parallel computer. By assigning small problems to subsets of processors and solving them simultaneously, we overcome the limited scalability of the parallel direct solver. On average, one FE mesh consists of 1,206,850 elements and 1,381,135 unknowns. Using 1200 cores (ten groups of ten twelve-core AMD 2.1 GHz processors), the hundred FE solutions are completed in 319 s. The FE and FD solutions along a selected survey line are shown in Fig. 16. Around the sources, the two solutions show large differences because they are computed in the total field mode. Except them, the two solutions are in the right ballpark, but the comparison is not as good as those of the previous two examples. We believe that subtle differences between the two solutions come from the different seafloor topography and different receiver orientations, demonstrating the sensitivity of CSEM to the two factors.
Acknowledgements This work was supported in part by National Natural Science Foundation of China (Grant No. 61303003 and 41374113). S.S.K acknowledges support from Korea National Research Foundation of the Ministry of Science, ICT and Future Planning (Grant No. NRF2013R1A1A1076071 and NRF-2012M2A8A5007440). This work benefited from the constructive comments of Colin Farquharson and three anonymous reviewers. References Amestoy, P.R., Duff, I.S., Koster, J., L’Excellent, J.-Y., 2001. A fully asynchronous multifrontal solver using distributed dynamic scheduling. SIAM J. Matrix Anal. Appl. 23, 15–41. Aminzadeh, F., Brac, J., Kunz, T., 1997. 3-D Salt and Overthrust Models. SEG/EAGE, Tulsa, Oklahoma. Fu, H., Wang, Y., Um, E.S., Fang, J., Wei, T., Huang, X., Yang, G., 2015. A parallel finiteelement time-domain method for transient electromagnetic simulation. Geophysics 80, E213–E224. Jin, J., Riley, D., 2008. Finite Element Analysis of Antennas and Arrays. Wiley-IEEE Press, Hoboken, New Jersey. Key, K., Ovall, J., 2011. A parallel goal-oriented adaptive finite element method for 2.5-D electromagnetic modelling. Geophys. J. Int. 186, 137–154. Li, Y., Key, K., 2007. 2D marine controlled-source electromagnetic modeling, part 1: an adaptive finite element algorithm. Geophysics 72, WA51–WA62. Newman, G.A., Alumbaugh, D.L., 1995. Frequency domain modeling of airborne electromagnetic responses using staggered finite differences. Geophys. Prospect. 43, 1021–1042. Schwarzbach, C., Börner, R.-U., Spitzer, K., 2011. Three-dimensional adaptive higher order finite element simulation for geo-electromagnetics-a marine CSEM example. Geophys. J. Int. 187, 63–74. Si, H., 2006. TetGen: A Quality Tetrahedral Mesh Generator and 3D Delaunay Trianulator. Weierstrass Institute for Applied Analysis and Stochastics, Berlin, Germany. Um, E., Commer, M., Newman, G., 2013. Efficient pre-conditioned iterative solution strategies for the electromagnetic diffusion in the Earth: finite-element frequencydomain approach. Geophys. J. Int. 193, 1460–1473. Um, E., Commer, M., Newman, G., 2014. A strategy for coupled 3D imaging of largescale seismic and electromagnetic data sets: application to subsalt imaging. Geophysics 79, ID1–ID13. Yang, D., Oldenburg, D.W., Haber, E., 2014. 3-D inversion of airborne electromagnetic data parallelized and accelerated by local mesh and adaptive soundings. Geophys. J. Int. 196, (1942-1507). Yang, D., Oldenburg, D.W., 2016. Survey decomposition: a scalable framework for 3D controlled-source electromagnetic inversion. Geophysics 81, E69–E87.
7. Conclusions We have analyzed numerical characteristics of tetrahedral meshes for CSEM modeling. A major characteristic of tetrahedral meshes is that the tetrahedral mesh density required for accurate solutions does not strongly depend on the maximum source-receiver offset because its numerical dispersion errors oscillate. Based on the analysis, a sequence of mesh generation tasks is casted into FEDO and is automated. Benefits of FEDO are summarized as follows. First, unlike FD grids, the
9