Three-dimensional magnetotelluric anisotropic parallel finite element forward method coupled with pml

Abstract

The application relates to a three-dimensional magnetotelluric anisotropy parallel finite element forward method coupled with PML. The method comprises the following steps: determining a non-reflection condition required by PML according to the continuity condition of an electromagnetic field and the dispersion relation of a plane wave in TE and TM modes; determining the parameter form of the stable absorption performance of PML by using the relationship between the incident field and the transmission field in a lossy medium; constructing an improved wide-band PML parameter tensor according to the non-reflection condition, the parameter form of the stable absorption performance and preset characteristic parameters; and realizing finite element forward by using the improved wide-band PML parameter tensor under a preset multi-level parallel architecture to obtain a finite element forward result. The method can be accurate, efficient and widely applicable.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of three-dimensional electromagnetic numerical simulation, in particular to a three-dimensional magnetotelluric anisotropy parallel finite element forward method coupled with PML. BACKGROUND

[0002] In resource exploration, mainly by measuring the change rule of natural or artificial electromagnetic field to explore the electrical structure of the underground or seabed by using the electromagnetic difference of the medium. Electromagnetic forward simulation provides technical support for inversion, and fast, efficient and accurate forward simulation can provide reference for the interpretation of actual data, help to understand the geological / sea bottom structure (such as identifying oil and gas based on the large resistivity difference between the oil-bearing reservoir and the surrounding water-saturated stratum), and reduce the risk and cost of land exploration and seabed drilling.

[0003] At present, the two-dimensional electromagnetic forward and inverse problems have been relatively mature, while the three-dimensional forward and inverse problems are not perfect, and the response law of complex geological bodies needs further understanding. Therefore, efficient and high-precision three-dimensional forward is a research hotspot and difficulty. In three-dimensional magnetotelluric forward, the methods widely used at present mainly include integral equation method, finite difference method and finite element method, etc. The integral equation method only discretizes the abnormal body, and is generally used for the calculation and analysis of electromagnetic field of simple model. The finite difference method uses difference to approximate differential, which has certain difficulty in dealing with arbitrary undulating terrain. The finite element method has the advantages of flexible grid division and accurate modeling of complex boundary, and can be used for models with complex physical property distribution or geometric characteristics, and has become one of the mainstream methods of magnetotelluric forward. The finite element method can be divided into node finite element and vector finite element. The vector finite element can automatically satisfy the continuity of tangential electric field or magnetic field and the condition of divergence being zero, thereby avoiding the problem of pseudo-solution that may be produced by node finite element. When using finite element for electromagnetic numerical simulation, whether using total field method, secondary field method or potential method, the open space needs to be truncated to reduce the calculation requirements. In order to avoid the influence of truncation, the outer boundary usually needs to be extended for a sufficient distance (usually several times the skin depth) and then Dirichlet, Neumann or mixed boundary conditions are used, which will greatly increase the calculation resources and time. Therefore, under the premise of ensuring calculation accuracy, improving the boundary conditions of magnetotelluric forward and compressing the numerical simulation solving space are very helpful to improve the efficiency of forward. Perfectly Matched Layer (PML) is a free space simulation method proposed by Berenger in 1994 based on field splitting theory. It is a specially designed virtual medium that can absorb electromagnetic waves without reflection. Berenger's PML is not perfect, and the control equation used in PML is a non-Maxwell equation, and the absorption effect of evanescent wave is not ideal. Chew and Weedon systematically analyzed Berenger's PML from the perspective of complex coordinate stretching, thereby providing a theoretical basis. Sacks et al. and Gedney proposed Uniaxial Anisotropic Perfectly Matched Layer (UPML) based on Maxwell's equations, which does not need to split the field and can better absorb evanescent waves. Kuzuoglu and Mittra introduced a complex frequency shift stretching operator and proposed a complex frequency shifted perfectly matched layer (CFS-PML), which can effectively absorb grazing waves and low-frequency waves.Roden and Gedney proposed the Convolutional Perfectly Matched Layer (CPML) by calculating the convolution terms in the time-domain PML equations using recursive formulas, thus improving computational efficiency. After continuous refinement, PML has been widely applied in recent years. Some researchers have used PML based on complex frequency shift scaling tensors to significantly improve the performance of absorbing boundaries. Others have applied CFS-PML to time-domain finite element simulations of Debye dispersive medium ground-penetrating radar (GPR), improving simulation accuracy. Some researchers have used staggered-grid finite-difference methods coupled with PML boundary conditions to improve the computational accuracy and efficiency of 3D magnetotelluric forward modeling. Some researchers have applied CFS-PML to high-order finite element methods to solve 3D GPR numerical simulation problems. Some researchers have applied PML to model order reduction algorithms, reducing the dimension of the Krylov subspace and improving convergence. Some researchers have proposed GMI-PML, avoiding the problem of improper time step synchronization and achieving higher absorption performance. UPML, derived from high-frequency wave equations, exhibits strong absorption in low-frequency diffuse fields, leading to significant numerical reflection. CFS-PML is also applicable in low-frequency diffuse fields, but its absorption effect is influenced by the model's conductivity and frequency, requiring adjustment of PML parameters based on the model. To address these issues, researchers have developed a broadband PML metric that improves the cutoff effect in diffuse fields, further expanding the application scenarios of PML.

[0004] However, the above techniques or methods have two drawbacks: PML is limited to scenarios involving only a single conductivity parameter variation and has low computational efficiency. Current PML is limited to scenarios involving only a single conductivity parameter variation, failing to consider magnetic permeability parameters and more complex anisotropic scenarios. However, in real-world scenarios, the magnetotelluric response is usually influenced by both the conductivity and magnetic permeability of the medium. Numerous studies have shown that in highly magnetic regions, the influence of magnetic permeability on the magnetotelluric response cannot be ignored. As magnetic susceptibility increases, its influence on the electromagnetic forward modeling response increases significantly. The magnetic susceptibility parameter has a stronger impact on the TM model than on the TE model, and high-resistivity anomalies are more susceptible to the influence of magnetic susceptibility than low-resistivity anomalies. Some scholars have pointed out that electromagnetic data contains inherent information about the depth distribution of magnetic susceptibility. Others have used three-dimensional models with conductivity and magnetic permeability anisotropy to compare the impact of different magnetic susceptibilities on apparent resistivity. Traditional three-dimensional electromagnetic forward modeling computation is inefficient. In electromagnetic forward modeling of anisotropic media based on vector finite element method, numerical simulation methods often require truncation of the infinite space to reduce computational demands. While inner boundary conditions are automatically satisfied in electromagnetic forward modeling, outer boundary conditions are typically set after extending the space sufficiently to avoid the effects of truncation. Boundary conditions can be categorized into three types: Type I, Type II, and Type III. Type I boundary conditions, also known as Dirichlet boundary conditions, directly impose known values ​​on the boundary. Type II boundary conditions, also known as Neumann boundary conditions, provide the normal derivative values ​​of the physical quantities on the boundary. Type III boundary conditions, also known as Robin boundary conditions, provide a combination of physical quantities and their normal derivative values ​​on the boundary, hence they are also called hybrid boundary conditions. Even with the computationally efficient Type III boundary conditions, the boundary still needs to be set far from the anomaly, resulting in a large computational scale. Furthermore, three-dimensional electromagnetic forward modeling is a large problem with a long computation time. Direct solution methods can achieve higher solution accuracy, but this further increases the consumption of computer memory and other system resources. A common approach is to use multi-frequency parallel solution to speed up the solution process, but the acceleration potential of intra-frequency parallelism and distributed storage has not yet been fully explored. Summary of the Invention

[0005] Therefore, it is necessary to provide an accurate, efficient, and widely applicable three-dimensional magnetotelluric anisotropic parallel finite element forward modeling method coupled with PML to address the aforementioned technical problems.

[0006] A three-dimensional magnetotelluric anisotropic parallel finite element forward modeling method coupled with PML, the method comprising:

[0007] Based on the continuity condition of the electromagnetic field and the dispersion relation of plane waves in TE and TM modes, the non-reflection condition that PML needs to satisfy is determined; the parameter form that PML can satisfy stable absorption performance is determined by using the relationship between the incident field and the transmitted field in the lossy medium.

[0008] An improved broadband PML parameter tensor is constructed based on the non-reflection condition, the parameter form that satisfies stable absorption performance, and the pre-set characteristic parameters.

[0009] Finite element forward modeling was performed using the improved wideband PML parameter tensor under a pre-set multi-level parallel architecture, and the finite element forward modeling results were obtained.

[0010] The aforementioned parallel finite element forward modeling method for three-dimensional magnetotelluric anisotropy coupled with PML, as described in this application, determines the reflection-free condition based on the electromagnetic field continuity condition and the dispersion relation of TE and TM modes, and determines the parameter form by combining the relationship between the incident and transmitted fields in the lossy medium. This allows the improved PML to simultaneously consider the changes in both conductivity and permeability. In the high magnetic region, permeability has a significant impact on the magnetotelluric response, and the improved PML model is no longer limited to a single conductivity scenario, making it suitable for more complex anisotropic models. Furthermore, the parameter form determined in this application, which satisfies stable absorption performance, enables the improved broadband PML to maintain stable absorption performance, overcoming the shortcomings of traditional PML at low frequencies and improving overall absorption performance. Finite element forward modeling is achieved using the improved PML model under a multi-level parallel architecture, fully exploring the acceleration potential of intra-frequency parallelism and distributed storage. Compared with traditional methods, this significantly reduces the number of grids, computation time, and memory usage. An improved PML parameter tensor is constructed by integrating the non-reflection condition, improved absorption performance, and pre-set characteristic parameters, making its performance more stable and ensuring the accuracy and reliability of the calculation results under different parameter variations and complex models. Attached Figure Description

[0011] Figure 1 This is a flowchart illustrating a three-dimensional magnetotelluric anisotropic parallel finite element forward modeling method coupled with PML in one embodiment.

[0012] Figure 2 This is a schematic diagram of plane wave incident in one embodiment;

[0013] Figure 3 This is a schematic diagram of the division of various regions of the PML in a three-dimensional model in one embodiment;

[0014] Figure 4 This is a schematic diagram of multi-level parallelism in another embodiment;

[0015] Figure 5 This is a diagram illustrating the vector finite element solution process for coupled PML boundary conditions in one embodiment.

[0016] Figure 6 This is a geometric schematic diagram of the experimental model in one embodiment;

[0017] Figure 7 This is a schematic diagram of the apparent resistivity, phase, and error calculated by Model 1 using the improved PML and Dirichlet boundary conditions of this application in one embodiment.

[0018] Figure 8 This is a schematic diagram of the apparent resistivity, phase, and error calculated by Model 2 using the PML and Dirichlet boundary conditions of this application in one embodiment.

[0019] Figure 9 XOY cross-sectional views of apparent resistivity and phase in XY and YX modes in one embodiment;

[0020] Figure 10 This is a schematic diagram illustrating the parallel acceleration results of PML and the conventional method at a single frequency point in one embodiment.

[0021] Figure 11 This is a schematic diagram illustrating the parallel acceleration of PML and conventional methods at 8 frequency points in one embodiment. Detailed Implementation

[0022] To make the objectives, technical solutions, and advantages of this application clearer, the following detailed description is provided in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the scope of this application.

[0023] In one embodiment, such as Figure 1 As shown, a parallel finite element forward modeling method for three-dimensional magnetotelluric anisotropy coupled with PML is provided. This method can be used for scenarios with varying conductivity and permeability as well as three-dimensional electromagnetic forward modeling containing arbitrary anisotropic anomalies, and can significantly improve the truncation effect. The method includes the following steps:

[0024] Step 102: Based on the continuity condition of the electromagnetic field and the dispersion relation of the plane wave in TE and TM modes, determine the non-reflection condition that the PML needs to satisfy; use the relationship between the incident field and the transmission field in the lossy medium to determine the parameter form of the PML to satisfy the stable absorption performance.

[0025] like Figure 2 The diagram shows a plane wave incident from the simulated region onto the PML. To simplify the derivation, displacement current is neglected in the diffused field, and the time harmonic factor is taken as . And assuming that PML is a uniaxial anisotropic medium, with anisotropic parameters and about Rotationally symmetric, using tensor form

[0026] (1)

[0027] The incident surface is The interface is , The region is the simulated region. The region is defined as PML. The conductivity of the simulated region is... Magnetic permeability is The electrical conductivity and magnetic permeability of PML are respectively and . Let be the angle of incidence relative to the normal direction. Represents the incident wave vector. Represents the reflected wave vector. This represents the transmitted wave vector. For ease of derivation, plane waves are represented using wave vectors and position vectors, i.e. Thus, the phase matching principle and operators can be utilized. The substitution relation represents the double curl equation. To make... and If a non-zero solution exists, then the coefficient matrix must be of partial rank. Based on this condition, we obtain... and The dispersion relation of the mode. Combining the continuity condition of the tangential components of the electric and magnetic fields at the interface, the non-reflection condition that the PML needs to satisfy is comprehensively solved. The purpose of using the PML is to truncate the computational domain as efficiently as possible; therefore, based on the first step, it is necessary to further determine the form of the PML parameters so that they satisfy both the non-reflection condition and have good absorption performance. Using the parameter relationships determined by the non-reflection condition and the electric and magnetic field relationships in the lossy medium, the parameter form that satisfies stable absorption performance is further solved as follows:

[0028] (2)

[0029] Should This form is applicable to both high-frequency wave fields and low-frequency diffused fields, and can be used across a wide frequency range. Ultimately, the three-dimensional parameter matching matrix of a broadband PML that meets the conditions can be obtained as follows:

[0030] (3)

[0031] in, , , These are characteristic parameters of the PML. The discretization of the simulation space using the finite element method leads to discontinuities on both sides of the interface between the simulation region and the PML, causing stray wave reflections. Therefore, in the direction perpendicular to the interface, the absorption of the PML should gradually increase from zero. This application... , , The parameters adopt a gradually changing exponential spatial distribution, and as they move away from the boundary, and The value gradually increases. The value gradually decreases. Typically, the total thickness of a PML is finite, and the number of layers is not excessive. When the grid cell size of the PML is uniform, characteristic parameters are given. , , The spatial distribution is as follows

[0032] (4)

[0033] in, The dielectric constant in vacuum. It is the thickness of the PML mesh cell. This represents the distance from the center of each PML element to the simulation region-PML interface. It is the total thickness of PML. It is a constant. At the outer boundary of the PML, the entire computational domain is terminated using a perfect electrical conductor boundary condition. In equation (4), , , They are respectively , , The maximum value. In order to control the absorption capacity of PML and balance the reflection error and numerical discretization error of PEC boundary, the maximum value is... , , The optimal value is denoted as , , When the outer boundary of PML is truncated by the PEC condition, the parameters... PML depth The function of PML. The reflection coefficient of PML is defined as...

[0034] (5)

[0035] in, It is a gradually changing distribution, which can be written as the coefficient is The function, i.e. ,in This is the thickness of each PML layer. The integral in equation (5) can be written in discrete form.

[0036] (6)

[0037] in It is the total number of PML layers on each side. This is the number of each PML layer. Equation (6) shows that the reflectance of the PML is related to its thickness. It is irrelevant to the number of layers in PML. This indicates that the absorption performance of the PML is insensitive to its thickness; as long as the number of PML layers remains constant, the thickness can be set as needed, thus allowing its application in various scenarios. This application considers both electrical conductivity and magnetic permeability parameters, resulting in a PML suitable for scenarios with varying dual parameters. Since the non-reflection and attenuation conditions of the PML are theoretically only related to the medium of adjacent simulation regions, it can be used for truncating simulation regions containing arbitrary anisotropic anomalous bodies.

[0038] Step 104: Construct the improved broadband PML parameter tensor based on the non-reflection condition, the parameter form that satisfies the stable absorption performance, and the pre-set spatial distribution of characteristic parameters.

[0039] The improved broadband PML parameter tensor is no longer limited to a single conductivity scenario, making it suitable for more complex anisotropic models. Furthermore, the parameter form determined in this application, which satisfies stable absorption performance, enables the improved PML to maintain stable absorption performance in low-frequency diffusion fields, overcoming the shortcomings of traditional PML at low frequencies and improving overall absorption performance.

[0040] Step 106: finite element forward modeling is performed using the improved wideband PML parameter tensor under a pre-set multi-level parallel architecture to obtain the finite element forward modeling results.

[0041] In a 3D model, such as Figure 3 As shown, the PML region consists of 6 planar regions, 12 edge regions, and 8 corner regions. The parameter forms of each region of the PML are shown in Table 1.

[0042] Table 1

[0043]

[0044] In calculation When a frequency point (e.g.) Figure 4 As shown), the program runs within the overall communication domain, using... Each node and Each process. Since the computation at each frequency point is independent, and considering sufficient parallelism within each frequency point, this application divides the total communication domain according to the number of frequency points. Each sub-communication domain contains [number] sub-communication domains. There are 10 frequency points to be calculated, and each frequency point is assigned a frequency point. Multiple processes run concurrently when calculating each frequency point to accelerate matrix assembly, right-hand side construction, and matrix solving. The matrix is ​​distributed and stored on different nodes responsible for calculating that frequency point to reduce the memory requirements of individual nodes. The petsc distributed toolkit is used for distributed processing in the matrix assembly and right-hand side construction parts, while the superlu-dist solver combined with petsc is used for distributed parallel solving in the matrix solving part.

[0045] The process of implementing finite element forward modeling using the improved PML model under a pre-set multi-level parallel architecture is as follows: Figure 5 As shown, when truncating the boundary in an open-domain problem, applying a perfectly matched layer to absorb outward traveling waves yields better results. Therefore, the calculation of the total electric field is decomposed into a primary field here. and secondary field Independent calculations. The relationship is...

[0046] (7)

[0047] In the total field, Maxwell's equations and the double curl equations of the electric field are respectively...

[0048] (8)

[0049] (9)

[0050] in, The relative permeability of the simulated region. It is the conductivity of the simulated region. The dielectric constant of the simulated region, For displacement current, The identity matrix is ​​denoted as . In a single field, the Maxwell's equations (governing equations) in the simulation region and the PML region are:

[0051] (10)

[0052] in, For background conductivity, Given the background relative permeability, the double curl equation of the electric field can be expressed as:

[0053] (11)

[0054] From equations (7), (9), and (11), the double curl equation of the secondary electric field can be obtained as follows:

[0055] (12)

[0056] In the simulation area, Equal to the identity matrix; in the PML region, The result is presented in the form shown in Table 1. Using the Galerkin weighted residual method, the entire region is integrated. Utilizing vector identities and the divergence theorem, the final result can be obtained.

[0057] (13)

[0058] in, for Variation, The unit coefficient matrix, The corresponding overall matrix equation is:

[0059] (14)

[0060] in, This is the overall coefficient matrix. After obtaining the electric field from equation (14), the magnetic field can be further obtained.

[0061] The field source used in the magnetotelluric method of this application is: and Two orthogonal sources. for , , ;source for , , .source The calculated electric and magnetic fields are as follows , , , ;source The calculated electric and magnetic fields are as follows , , , The formula for calculating tensor impedance is:

[0062] (15)

[0063] From equation (15), the apparent resistivity and phase under different modes can be obtained as follows:

[0064] (16)

[0065] The vector finite element method with coupled PML boundary conditions only requires adding PML constraints to the governing equations and multiplying the element stiffness matrix by the PML parameter matching matrix when establishing the element stiffness matrix. The rest is the same as the ordinary vector finite element method.

[0066] The aforementioned parallel finite element forward modeling method for three-dimensional magnetotelluric anisotropy coupled with PML, as described in this application, determines the reflection-free condition based on the electromagnetic field continuity condition and the dispersion relation of TE and TM modes, and determines the parameter form by combining the relationship between the incident and transmitted fields in the lossy medium. This allows the improved PML to simultaneously consider the changes in both conductivity and permeability. In the high magnetic region, permeability has a significant impact on the magnetotelluric response, and the improved PML model is no longer limited to a single conductivity scenario, making it suitable for more complex anisotropic models. Furthermore, the parameter form determined in this application, which satisfies stable absorption performance, enables the improved PML to maintain stable absorption performance in low-frequency diffusion fields, overcoming the shortcomings of traditional PML at low frequencies and improving overall absorption performance. Finite element forward modeling is achieved using the improved wideband PML model under a multi-level parallel architecture, fully exploring the acceleration potential of intra-frequency parallelism and distributed storage. Compared with traditional methods, this significantly reduces the number of grids, computation time, and memory usage. An improved PML was constructed by integrating non-reflection conditions, improved absorption performance, and pre-set characteristic parameters, making its performance more stable and ensuring the accuracy and reliability of calculation results under different parameter variations and complex models.

[0067] In one embodiment, based on the continuity condition of the electromagnetic field and the dispersion relation of the plane wave in TE and TM modes, the reflection-free conditions that the PML needs to satisfy are determined, including:

[0068] Plane waves are represented using wave vectors and position vectors, i.e. Using the phase matching principle and operators The substitution relation represents the double curl equation, and to make and If there is a non-zero solution, then the coefficient matrix must be of incomplete rank. Based on the condition that the coefficient matrix is ​​of incomplete rank, the dispersion relation in the TE and TM modes is obtained. Then, combined with the continuity condition of the tangential components of the electric and magnetic fields at the interface, the reflection-free condition that the PML needs to satisfy is solved.

[0069] In one embodiment, the reflection-free condition that needs to be satisfied to solve the PML is:

[0070]

[0071] in, Represents the anisotropic parameter tensor. Indicates parameters, express , and direction.

[0072] In one embodiment, the parameter form for determining the stable absorption performance of the PML in a low-frequency diffusion field is determined by utilizing the relationship between the incident field and the transmission field in the lossy medium, including:

[0073] The parameter form for determining the stable absorption performance of PML is as follows, based on the relationship between the incident field and the transmitted field in a lossy medium:

[0074]

[0075] In a low-frequency diffusion field, it simplifies to

[0076]

[0077] in, , , These represent the different characteristic parameters of PML. For electrical conductivity, Permeability, It represents angular frequency.

[0078] In one embodiment, the pre-set feature parameters are:

[0079]

[0080] in, The dielectric constant in vacuum. It is the thickness of the PML mesh cell. This represents the distance from the center of each PML element to the simulation region-PML interface. It is the total thickness of PML. It is a constant. , , They are respectively , , The maximum value, It represents the magnetic permeability.

[0081] In one embodiment, the pre-configured multi-level parallel architecture includes:

[0082] In calculation When using a single frequency point, within the total communication domain, using Each node and The process divides the total communication domain into several parts based on the number of frequency points. Each sub-communication domain contains [number] sub-communication domains. There are 10 frequency points to be calculated, and each frequency point is assigned a frequency point. Multiple processes run simultaneously to accelerate matrix assembly, right-hand side construction, and matrix solving when calculating each frequency point. The matrix is ​​distributed and stored on different nodes responsible for calculating that frequency point. The petsc distributed toolkit is used for distributed processing in the matrix assembly and right-hand side construction parts, and the superlu-dist solver is used in combination with petsc for distributed parallel solving in the matrix solving part.

[0083] In specific embodiments, air-to-ground model 1, air-to-sea-to-ground model 2, and different magnetic permeability model groups 3 were used, and all numerical experiments were completed on a supercomputer. Figure 6 As shown, the geometric dimensions of the simulation region in the model are all 2000 m × 2000 m × 2000 m, and the background resistivity is 100. The background relative permeability is 1. The anomaly is located at the center of the simulation region, with a top burial depth of 600 m and dimensions of 800 m × 400 m × 400 m. The simulation region and the anomaly are uniformly divided into hexahedral units of 100 m × 100 m × 100 m. When using PML, six layers of PML are set outside the simulation region, each with a thickness of 0.05 m. The outermost layer is truncated under perfect electrical conductor conditions. All parameters are set as follows. , , When using the traditional mesh extension method, extending the mesh about 5 times the skin depth outside the simulation region (at this extension distance, the secondary field attenuates to less than 1% at the outer boundary), the size growth rate of the extended mesh is 1.25. Finding the optimal parameters for the PML is a current technique and will not be elaborated here. Figure 6 shows the model geometry; the blue box represents the computational region, the purple box represents the simulation region, and the region in between is the PML or extended region. The green part in the middle represents the anomaly.

[0084] Model 1 is a simplified geological model. The top four layers of the background medium in the simulated area are air layers with a resistivity of [missing value]. The relative permeability is 1. The anomalous body is anisotropic, with principal axis resistivity of... , , Euler rotation angles are respectively , , The principal spindle susceptibility (principal spindle susceptibility + 1 = principal spindle permeability) are 0.3, 0.2, and 0.1, respectively, and the Euler rotation angles are respectively... , , The test frequency was 0.125 Hz. Figure 7The apparent resistivity, phase, and error calculated using the improved PML and Dirichlet boundary conditions of this application are shown. The test frequency was 0.125 Hz. ResXY-PML represents the PML calculation result of this application, ResXY-Diri represents the Dirichlet boundary condition calculation result (reference value), and RE represents the relative error (‰). Indicates absolute error ( The first line shows the calculation results of the survey line along the X direction at Y=50 m on the ground. The second to fourth lines are schematic diagrams of the XOY section.

[0085] Model 2 is a simplified layered ocean model. Unlike Model 1, the simulated region contains two layers of seawater beneath the air layer of the background medium, with a resistivity of [missing information]. The relative permeability is 1. The anomalous body is anisotropic, with principal axis resistivity of... , , Euler rotation angles are respectively , , The principal spindle susceptibility is 0.2, 0.1, and 0.3, respectively, and the Euler rotation angles are respectively... , , The test frequency was 0.1 Hz. Figure 8 The apparent resistivity, phase, and error calculated using the PML and Dirichlet boundary conditions of this application are shown. The test frequency was 0.1 Hz. ResXY-PML represents the PML calculation result of this application, ResXY-Diri represents the Dirichlet boundary condition calculation result (reference value), and RE represents the relative error (‰). Indicates absolute error ( The first line shows the calculation results of the seabed survey line along the X direction at Y=50 m, and the second to fourth lines are schematic diagrams of the XOY section.

[0086] Model group 3 is a simplified geological model used to test the effects of different magnetic susceptibility, illustrating the significance of the dual-parameter PML considering both electrical conductivity and magnetic permeability in this application, compared to previous PMLs that only considered electrical conductivity. These models differ from Model 1 only in the magnetic susceptibility of the anomaly; all other settings are identical, and they all use the same solution method with the PML of this application as boundary conditions. Based on the actual conditions of the natural world, four settings—negative magnetic susceptibility, zero magnetic susceptibility, low magnetic susceptibility, and high magnetic susceptibility—are considered for comparison. The corresponding magnetic susceptibility of the principal axes of the anomaly are: -0.1, -0.2, -0.3; 0, 0, 0; 0.3, 0.2, 0.1; 0.7, 0.6, 0.5. Figure 9The XOY cross-sectional plots of apparent resistivity and phase in XY and YX modes are shown. Figure 9 shows the XOY cross-sectional plots of apparent resistivity and phase for model group 4 under different permeabilities. The test frequency was 0.125 Hz. ResXY represents the apparent resistivity in the XY mode, PhaseXY represents the phase in the XY mode, and ms = -0.1, -0.2, and -0.3 represent the principal axis susceptibility of the anomalous body as -0.1, -0.2, and -0.3, respectively.

[0087] To fully test the acceleration effect brought by multi-level parallelism and PML boundary conditions, two sets of tests were conducted: single-frequency and 8-frequency tests. The computational model is as follows: Figure 6 As shown, the number of unknowns in the traditional method is 888,822, and the test frequency is 0.125Hz. Figure 10 shows the parallel acceleration results of PML and the traditional method at a single frequency point. The test frequency is 0.125Hz. The blue solid line with circles represents the parallel speedup ratio of the PML method of this application, the black solid line with triangles represents the parallel speedup ratio of the traditional method, and the red dashed line with squares represents the total speedup ratio of the parallel PML method compared to the serial traditional method. Figure 11 shows the parallel acceleration of PML and the traditional method at 8 frequency points. The test frequency is 0.125Hz. The blue solid line with circles represents the parallel speedup ratio of the PML method of this application, the black solid line with triangles represents the parallel speedup ratio of the traditional method, and the red dashed line with squares represents the total speedup ratio of the parallel PML method compared to the serial traditional method. Table 2 shows the performance comparison between the PML method of this application and the traditional method. Table 3 shows the comparison of parallel computing time at a single frequency point. Table 4 shows the comparison of parallel computing time at 8 frequency points.

[0088] Table 2

[0089]

[0090] Table 3

[0091]

[0092] Table 4

[0093]

[0094] exist Figure 7 and 8In Figure 7, (a) and (b) show the apparent resistivity in the XY and YX modes, respectively, while (c) and (d) show the phase in the XY and YX modes, respectively. The first row of the image shows the apparent resistivity and phase curves for the survey line at Y=50m on the ground and seabed. The second and third rows show the XOY cross-sectional diagrams of the apparent resistivity and phase on the ground. The fourth row shows the relative error of the apparent resistivity and the absolute error of the phase. The results show that the calculated results of the two methods are in excellent agreement. In Figure 7, the maximum relative error of the apparent resistivity in both modes is less than 0.2‰, and the absolute error of the phase is... In Figure 8, the maximum relative error of apparent resistivity under both modes is less than 0.03‰, and the absolute error of phase is... This indicates that the PML method of this application is applicable to the forward modeling of homogeneous media and anisotropic anomalous bodies, as well as layered media and anisotropic anomalous bodies. Moreover, it only requires 6 layers of PML absorption, totaling 0.3m, to achieve calculation results that are almost identical to those of the first type of boundary conditions using the large-distance continuation method, demonstrating very high computational efficiency and accuracy.

[0095] exist Figure 9 In the figure, (a) to (d) show the calculated results from negative to high magnetic susceptibility, respectively. The first to fourth rows show the apparent resistivity and phase in the XY and YX modes, respectively. It can be observed that as the magnetic susceptibility increases, the dark areas of low resistivity gradually fade, while the yellow areas of high resistivity become brighter. The phase cross-section plots show the same trend. This indicates that it is essential to improve upon previous PMLs that only considered conductivity parameters and propose a PML that considers both conductivity and permeability variations.

[0096] Parallel test results of the two methods at single and eight frequencies are as follows: Figure 10 and Figure 11 As shown, the PML speedup curve represents the ratio of single-process computation time to parallel computation time of the PML method in this application, the traditional speedup curve represents the ratio of single-process computation time to parallel computation time of the traditional method, and the overall speedup curve represents the ratio of single-process computation time of the traditional method to parallel computation time of the PML method in this application. It is the acceleration effect brought about by the combined application of hierarchical parallelization design and the PML boundary conditions in this application. Figure 10The results show that the PML speedup reaches its maximum at 32 processes, approximately 5.46; the traditional speedup is approximately 6.73 at 128 processes; and the overall speedup reaches its maximum at 32 processes, approximately 85.15. For 1-4 processes, the hierarchical parallel design provides almost the same speedup effect for both methods. From 4-32 processes, the speedup of both methods increases rapidly with the number of processes, but the speedup curve of the PML method in this application is above that of the traditional method. This may be because the advantage of the PML method in saving computational resources is converted into additional parallel speedup effect in the additional system overhead brought by the parallel design. From 32-168 processes, the traditional method still shows some improvement in speedup, but the improvement effect from increasing the number of processes becomes weaker. The speedup of the PML method in this application no longer increases but instead begins to decrease slightly. At this point, the computation time of the PML method in this application has decreased from 106.2 seconds for a single process to 19.44 seconds for 32 processes. This is because the additional system overhead of dividing the computation into more processes is too great, resulting in a slight decrease in speedup.

[0097] The overall speedup can be divided into three parts: the boundary condition acceleration brought by PML in this application, the parallel acceleration brought by the hierarchical parallelization design, and the overhead acceleration brought by the smaller parallel overhead of the PML method. For 1-4 processes, the overall speedup curve has almost the same slope as the other two curves, indicating that the parallel acceleration effect is almost the same at this point, with the boundary condition acceleration causing the overall speedup curve to be above the other two curves. For 4-32 processes, apart from the boundary condition acceleration, the overall speedup and PML speedup have almost the same parallel acceleration and overhead acceleration. However, further observation of the three curves reveals that the overhead acceleration effect of both the overall speedup and PML speedup is best for 4-8 processes, gradually decreasing for 8-32 processes. For 32-128 processes, for the same reason, the overall speedup curve and the PML speedup curve show the same trend.

[0098] Figure 11The results show that the PML speedup reaches its maximum value of approximately 37.07 with 512 processes; the traditional speedup is approximately 33.62 with 512 processes; and the total speedup also reaches its maximum value of approximately 649.63 with 512 processes. From 1 to 32 processes, the PML speedup increases rapidly with the number of processes, while the traditional speedup increases slowly. This may be because the traditional method requires more system resources, and the current number of processes cannot fully realize its acceleration potential compared to the additional overhead of parallelism. From 32 to 256 processes, all three speedup curves rise rapidly with the number of processes. From 256 to 512 processes, the PML and total speedup curves gradually slow down, which may be due to the same reason as the speedup bottleneck that occurs in single-frequency parallel testing. In the figure, the PML speedup curve of this application is consistently above the traditional method speedup curve. This may be because the increased computational frequency requires more computational resources, leading to a more significant acceleration effect. As shown in Table 2, in the three test models, the difference in mesh count between the two methods exceeded 8 times, the time difference exceeded 13 times, and the average memory usage exceeded 12 times. This indicates that, for the same truncation effect, the PML of this application has significant advantages in terms of mesh count, computation time, and memory usage compared to the traditional boundary conditions. Tables 3 and 4 show that the hierarchical parallel design of this application has a significant speedup effect for both single-frequency and multi-frequency points. Compared to the traditional boundary conditions, the maximum speedup of the parallelized PML design of this application is approximately 85.24 for single-frequency points and approximately 649.63 for 8-frequency points. Further observation reveals that the PML method of this application achieves near-maximum speedup with a relatively small number of processes (e.g., the speedup is almost the same with 32 and 128 processes at a single frequency point, and with 128 and 512 processes at an 8-frequency point), while traditional boundary conditions require more processes to achieve maximum speedup (e.g., more than 128 processes are needed at a single frequency point, and more than 512 processes are needed at an 8-frequency point). This further demonstrates the advantage of the multi-level parallelized PML method of this application in saving system resources.

[0099] Previous PML methods were derived based on conductivity parameters, without considering the influence of magnetic permeability. This application proposes a parallel finite element forward modeling method for 3D magnetotelluric conductivity and magnetic permeability anisotropy, coupled with PML boundary conditions. Unlike previous methods, this application considers the magnetic permeability parameter, proposing a PML applicable to variations in both conductivity and magnetic permeability, and improving it for low-frequency bands, applying it to anisotropic models. Furthermore, to further improve 3D forward modeling efficiency, a multi-level parallel design was implemented. Numerical experiments yielded the following conclusions: 1. The proposed PML method is correct, efficient, and accurate, suitable for homogeneous and layered models with anisotropic conductivity and magnetic permeability anomalies; 2. Compared with traditional boundary conditions, the proposed PML reduces the number of meshes by more than 85%, computation time by more than 90%, and memory usage by more than 90%, improving the efficiency of 3D electromagnetic forward modeling; 3. Compared to traditional methods with 888,822 unknowns, this application's PML achieves a speedup of approximately 85.24 when solving a single frequency point using 32 processes and approximately 649.63 when solving 8 frequency points using 512 processes, and can achieve maximum speedup with fewer processes; 4. The parameters of this application's PML have a wide range of applications, and the same set of optimal parameters has a good truncation effect in models of different sizes, mesh divisions, and test frequencies; 5. Changes in magnetic susceptibility will affect apparent resistivity and phase. As magnetic susceptibility increases, apparent resistivity and phase will also be affected and increase, even producing drastically different characteristics.

[0100] Compared to traditional boundary conditions, the Perfectly Matched Layer (PML) is a more efficient and accurate truncation method. However, current PMLs are limited to scenarios with only a single conductivity parameter variation and cannot be applied to dual-parameter variations of conductivity and permeability, or complex anisotropic models. Therefore, this application considers both conductivity and permeability parameters, as well as anisotropy, and proposes a PML applicable to complex anisotropic models. Furthermore, it combines MPI design with a multi-level parallel scheme to achieve parallel vector finite element forward modeling of three-dimensional magnetotelluric conductivity and permeability anisotropy with coupled PML boundary conditions. Comparison with previous results verifies that the PML boundary conditions proposed in this application have the advantages of high efficiency, high accuracy, and stable performance. Numerical experimental results of several models in this application show that, compared with traditional boundary conditions, the PML mesh number in this application is reduced by more than 85%, and the computation time and memory usage are reduced by more than 90%. Compared with the traditional method with 888,822 unknowns, by applying the PML and multi-level parallel design in this application, the problem of low efficiency in three-dimensional electromagnetic forward modeling is alleviated. When solving a single frequency point with 32 processes, the speedup ratio is approximately 85.24, and when solving 8 frequency points with 512 processes, the speedup ratio is approximately 649.63. Moreover, the maximum speedup effect can be achieved with fewer processes. The PML proposed in this application has a wider range of applications and better performance, and has a broader application prospect.

[0101] It should be understood that, although Figure 1 The steps in the flowchart are shown sequentially as indicated by the arrows, but these steps are not necessarily executed in the order indicated by the arrows. Unless explicitly stated in this application, there is no strict order restriction on the execution of these steps, and they can be executed in other orders. Furthermore, Figure 1 At least some of the steps in the process may include multiple sub-steps or multiple stages. These sub-steps or stages are not necessarily completed at the same time, but can be executed at different times. The execution order of these sub-steps or stages is not necessarily sequential, but can be executed in turn or alternately with other steps or at least some of the sub-steps or stages of other steps.

[0102] The technical features of the above embodiments can be combined in any way. For the sake of brevity, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0103] The embodiments described above are merely illustrative of several implementation methods of this application, and while the descriptions are specific and detailed, they should not be construed as limiting the scope of the invention. It should be noted that those skilled in the art can make various modifications and improvements without departing from the concept of this application, and these all fall within the protection scope of this application. Therefore, the protection scope of this application should be determined by the appended claims.

Claims

1. A three-dimensional magnetotelluric anisotropic parallel finite element forward modeling method coupled with PML, characterized in that, The method includes: Based on the continuity condition of the electromagnetic field and the dispersion relation of plane waves in TE and TM modes, the reflection-free condition that the PML needs to satisfy is determined; the parametric form for the stable absorption performance of the PML is determined using the relationship between the incident field and the transmitted field in the lossy medium; the reflection-free condition that the PML needs to satisfy is then solved. in, Represents the anisotropic parameter tensor. Indicates parameters, express , and direction; The parametric form for ensuring stable absorption performance of a PML is determined by utilizing the relationship between the incident and transmitted fields in a lossy medium, including: The parameter form for determining the stable absorption performance of PML is as follows, based on the relationship between the incident field and the transmitted field in a lossy medium: in, , , These represent the different characteristic parameters of PML. For electrical conductivity, Permeability, Indicates angular frequency; The pre-set feature parameters are: in, The dielectric constant in vacuum. It is the thickness of the PML mesh cell. This represents the distance from the center of each PML element to the simulation region-PML interface. It is the total thickness of PML. It is a constant. , , They are respectively , , The maximum value, Indicates magnetic permeability; The improved broadband PML parameter tensor is determined based on the non-reflection condition, the parameter form that satisfies stable absorption performance, and the pre-set spatial distribution of characteristic parameters. Finite element forward modeling is performed using the improved wideband PML parameter tensor under a pre-set multi-level parallel architecture to obtain finite element forward modeling results.

2. The method according to claim 1, characterized in that, Based on the continuity condition of the electromagnetic field and the dispersion relation of plane waves in TE and TM modes, the reflection-free conditions that PML needs to satisfy are determined, including: Plane waves are represented using wave vectors and position vectors, i.e. Using the phase matching principle and operators The substitution relation represents the double curl equation, and to make and If there is a non-zero solution, then the coefficient matrix must be of incomplete rank. Based on the condition that the coefficient matrix is ​​of incomplete rank, the dispersion relation in the TE and TM modes is obtained. Then, combined with the continuity condition of the tangential components of the electric and magnetic fields at the interface, the reflection-free condition that the PML needs to satisfy is solved.

3. The method according to claim 1, characterized in that, The pre-configured multi-level parallel architecture includes: In calculation When using a single frequency point, within the total communication domain, using Each node and The process divides the total communication domain into several parts based on the number of frequency points. Each sub-communication domain contains [number] sub-communication domains. There are 10 frequency points to be calculated, and each frequency point is assigned a frequency point. Multiple processes run simultaneously to accelerate matrix assembly, right-hand side construction, and matrix solving when calculating each frequency point. The matrix is ​​distributed and stored on different nodes responsible for calculating that frequency point. The petsc distributed toolkit is used for distributed processing in the matrix assembly and right-hand side construction parts, and the superlu-dist solver combined with petsc is used for distributed parallel solving in the matrix solving part.

Images