Hybrid precision domain decomposition solution method based on electrical source electromagnetic method exploration

Abstract

The application discloses a mixed precision area decomposition solving method based on an electric source electromagnetic method exploration, and comprises the following steps: establishing a forward model and performing tetrahedral mesh division; assigning resistivity to any tetrahedral mesh; dividing the tetrahedral mesh into a plurality of non-overlapping sub-regions; establishing a discrete equation of an electric field control equation for any sub-region; coupling the sub-regions according to a boundary condition, and constructing a single-precision interface equation and a double-precision interface equation; performing single-precision iterative solving on the single-precision interface equation by using a minimum residual method; when Arnoldi residual norm stagnation or a preset condition is reached, a double-precision iterative residual calculated by using the double-precision interface equation is used to correct the single-precision iteration, and a corrected parameter is obtained; performing single-precision iterative solving on the corrected parameter by using the single-precision interface equation; obtaining a solution of the interface equation; and substituting the solution of the interface equation into a corresponding receiving point in the sub-region, so as to obtain an electromagnetic response of the receiving point.

Description

Technical Field

[0001] This invention relates to the field of geophysical exploration technology, and in particular to a mixed-precision domain decomposition solution method based on electromagnetic exploration using electrical sources. Background Technology

[0002] As is well known, controlled-source electromagnetic methods have advantages such as high transmission power and strong anti-interference ability, and are widely used in mineral exploration and engineering geology. Compared with magnetic source emission, electric source emission has a greater exploration depth, is more sensitive to high-resistivity media, and has better vertical resolution, and is widely used in electromagnetic exploration fields such as semi-airborne electromagnetic methods and marine electromagnetic methods. With increasingly complex exploration environments and higher requirements for exploration accuracy, conventional one- and two-dimensional forward and inverse methods can no longer meet exploration requirements. Currently, three-dimensional forward and inverse methods have become a research hotspot and challenge in this field. Three-dimensional forward modeling is the foundation of three-dimensional inverse modeling, and it is essential to study accurate and efficient three-dimensional forward modeling methods. Compared with forward modeling methods such as finite volume and finite difference methods, the finite element method has higher computational accuracy.

[0003] Currently, the finite element method (FEM) based on unstructured tetrahedral meshes can accurately simulate complex models and is used for fine-grained exploration of complex terrain areas, making it one of the most commonly used methods in electromagnetic forward modeling. However, to obtain an accurate solution, the FEM requires fine meshing of the computational domain, resulting in millions or even tens of millions of meshes, consuming enormous computational resources and exhibiting low computational efficiency. The FETI-DP method (Dual-Primal Finite Element Tearing and Interconnecting method), which combines the FEM and domain decomposition, can decompose the solution domain into several sub-regions and couple all sub-regions through internal boundary conditions to form reduced-order interface equations. This method reduces the solution dimensionality and is easy to parallelize, improving computational efficiency. The reduced-order interface equations solve for internal boundary conditions, and the requirements for computational accuracy are not high.

[0004] Currently, existing research uses double-precision calculations, while single-precision calculations are twice as fast as double-precision calculations. This process wastes computing power and reduces computational efficiency.

[0005] Therefore, there is an urgent need to propose a hybrid accuracy solution method for interface equations that can both guarantee computational accuracy and utilize low computational power and speed. Summary of the Invention

[0006] To address the aforementioned problems, the present invention aims to provide a mixed-precision domain decomposition solution method based on electromagnetic exploration using electrical sources. The technical solution adopted by the present invention is as follows:

[0007] The mixed-accuracy domain decomposition solution method based on electromagnetic exploration using electrical sources includes the following steps:

[0008] Based on the underground geoelectric structure, a forward model of the area to be tested is established and tetrahedral mesh is generated to obtain the node and mesh information of the tetrahedral mesh.

[0009] Assign resistivity values ​​to any tetrahedral mesh;

[0010] Divide the tetrahedral mesh into several non-overlapping sub-regions;

[0011] Discrete equations for the electric field control equations are established for any sub-region;

[0012] The sub-regions are coupled according to the boundary conditions, and single-precision interface equations and double-precision interface equations are constructed.

[0013] The minimum residual method is used to solve the single-precision interface equation using single-precision iterative solution.

[0014] In single-precision iteration, when the Arnoldi residual norm stagnates or reaches a preset condition, the single-precision iteration residual calculated by the double-precision interface equation is used to correct the single-precision iteration, and the corrected parameters are obtained.

[0015] The corrected parameters are solved using single-precision iterative methods based on the single-precision interface equation until the convergence condition is met; the solution to the interface equation is then obtained.

[0016] Substituting the solution of the interface equation into the corresponding receiving point within the subregion yields the electromagnetic response of the receiving point.

[0017] Furthermore, in the single-precision iteration, when the Arnoldi residual norm stagnates or reaches a preset condition, the single-precision iteration is corrected using the double-precision iterative residual calculated by the double-precision interface equation, resulting in corrected parameters, including:

[0018] The convergence threshold k for single-precision iteration and the relative change threshold m for the residual norm of adjacent iterations are preset;

[0019] The single-precision interface equation is solved by single-precision iterative solution.

[0020] If the single-precision iteration residual norm is less than the convergence threshold k or the relative change in the residual norm between adjacent iterations is less than the relative change threshold m, then the single-precision iteration is corrected using the double-precision iteration residual calculated by the double-precision interface equation, and the corrected parameters are obtained.

[0021] Furthermore, the correction of the single-precision iteration using the double-precision residual calculated by the double-precision interface equation also includes:

[0022] The convergence threshold p for double-precision iteration is preset;

[0023] If the relative iterative residual norm calculated by the double-precision interface equation is less than the convergence threshold p, then the iteration is stopped.

[0024] If the relative iterative residual norm calculated by the double-precision interface equation is greater than or equal to the convergence threshold p, then the corrected parameters are solved iteratively using the single-precision interface equation until the convergence condition is met.

[0025] Furthermore, the expression for the electric field control equation is as follows:

[0026]

[0027] in, Represents the differential operator; E represents the electric field strength; ω represents the angular frequency; μ represents the permeability; σ represents the conductivity; J s This indicates the intensity of the external electrical source.

[0028] Furthermore, the expression for the discrete equation is:

[0029] (A p +iωB p E p =-iωS p -λ p

[0030] Where p represents the p-th sub-region; A p B represents the stiffness matrix of the p-th subregion; p S represents the quality matrix of the p-th sub-region; p Let λ represent the discrete emission source term of the p-th sub-region. p This represents the boundary term of the p-th sub-region.

[0031] Furthermore, the expression for the boundary condition is:

[0032]

[0033] Where n represents the unit normal vector at the interface; μ represents the magnetic permeability; Λ p Γ represents an unknown value; p Indicates the internal boundary.

[0034] Furthermore, the expression for the single-precision interface equation is:

[0035] F s λ s =b s

[0036] Among them, F sλ represents the coefficient matrix formed by coupling the discrete equations of a sub-region under single-precision conditions. s In single-precision, based on the boundary term λ p The resulting unknown terms to be solved; b s This represents the known terms formed by the coupling of source terms in the discrete equations of a sub-region under single-precision conditions.

[0037] The expression for the double-precision interface equation is:

[0038] F d λ d =b d

[0039] Among them, F d λ represents the coefficient matrix formed by the coupling of discrete equations in a sub-region under double precision. d In double precision, based on the boundary term λ p The resulting unknown terms to be solved; b d This represents the known terms formed by the coupling of source terms in the discrete equations of a subregion under double precision.

[0040] Furthermore, the single-precision iterative residual norm is the Arnoldi residual norm used in the minimum residual method calculation process.

[0041] Furthermore, the double-precision iterative residual calculated by the double-precision interface equation is the residual ε of the double-precision equation system to be solved, and its expression is:

[0042] ε=F d λ n -b d

[0043] Where, λ n F represents the iteration update amount of the nth iteration of double precision; d b represents the coefficient matrix formed by coupling the discrete equations of the sub-region under double precision; d This represents the known terms formed by the coupling of source terms in the discrete equations of a sub-region under double precision.

[0044] The double-precision iterative residual calculated using the double-precision interface equation is used to correct the single-precision iteration, resulting in corrected parameters. These corrected parameters include the iteration update amount λ for the nth double-precision iteration. n The iterative error of calculating the double-precision interface equation.

[0045] Furthermore, the relative iterative residual norm L calculated by the double-precision interface equation is expressed as follows:

[0046]

[0047] Compared with the prior art, the present invention has the following beneficial effects:

[0048] (1) This invention combines the domain decomposition method and mixed precision solution technology, which can make full use of the low precision computing power of the computer while ensuring the calculation accuracy and improving the calculation efficiency.

[0049] (2) This invention utilizes the single-precision computing power of a computer to quickly solve the interface equations stored in single precision, thereby fully leveraging the computational speed of single precision to improve the computational efficiency of equation solving. In addition, this invention uses the high-precision stored interface equations to adaptively correct computational parameters by utilizing the iteration error of single precision, thus ensuring computational accuracy while utilizing single-precision computing power.

[0050] In summary, this invention has the advantages of simple logic and high calculation accuracy. It is applicable to different electromagnetic exploration fields, such as marine electromagnetic exploration and semi-airborne electromagnetic exploration, and has high practical and promotional value in the field of geophysical exploration technology. Attached Figure Description

[0051] To more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings used in the embodiments will be briefly introduced below. It should be understood that the following drawings only show some embodiments of the present invention and should not be regarded as a limitation on the scope of protection. For those skilled in the art, other related drawings can be obtained based on these drawings without creative effort.

[0052] Figure 1 This is a logic flowchart of the present invention.

[0053] Figure 2 This is a comparison of the computational accuracy and efficiency of the method of the present invention with that of traditional methods. Detailed Implementation

[0054] To make the objectives, technical solutions, and advantages of this application clearer, the present invention will be further described below with reference to the accompanying drawings and embodiments. The embodiments of the present invention include, but are not limited to, the following embodiments. All other embodiments obtained by those skilled in the art based on the embodiments in this application without inventive effort are within the scope of protection of this application.

[0055] In this embodiment, the term "and / or" is merely a description of the relationship between related objects, indicating that there can be three relationships. For example, A and / or B can represent three situations: A exists alone, A and B exist simultaneously, and B exists alone.

[0056] The terms "first" and "second," etc., used in the specification and claims of this embodiment are used to distinguish different objects, not to describe a specific order of objects. For example, "first target object" and "second target object," etc., are used to distinguish different target objects, not to describe a specific order of target objects.

[0057] In the embodiments of this application, the terms "exemplary" or "for example" are used to indicate that something is an example, illustration, or description. Any embodiment or design that is described as "exemplary" or "for example" in the embodiments of this application should not be construed as being more preferred or advantageous than other embodiments or design. Specifically, the use of the terms "exemplary" or "for example" is intended to present the relevant concepts in a specific manner.

[0058] In the description of the embodiments in this application, unless otherwise stated, "multiple" means two or more. For example, multiple processing units means two or more processing units; multiple systems means two or more systems.

[0059] like Figures 1 to 2 As shown, this embodiment provides a mixed-accuracy domain decomposition solution method based on electromagnetic exploration using electrical sources, which includes the following steps:

[0060] The first step is to establish a three-dimensional forward model according to the exploration requirements, and to perform tetrahedral meshing based on the characteristics of the finite element method.

[0061] The second step is to divide the tetrahedral mesh into sub-regions to obtain the partitioned mesh information.

[0062] The third step is to independently establish discrete equations for the electric field control equations for each sub-region, and store them in both single-precision and double-precision formats.

[0063] In this embodiment, the expression for the electric field control equation is:

[0064]

[0065] in, Represents the differential operator; E represents the electric field strength; ω represents the angular frequency; μ represents the permeability; σ represents the conductivity; J s This indicates the intensity of the external electrical source.

[0066] Furthermore, the expression for the discrete equation in this embodiment is:

[0067] (A p +iωB p E p =-iωS p -λ p

[0068] Where p represents the p-th sub-region; A p B represents the stiffness matrix of the p-th subregion; p S represents the quality matrix of the p-th sub-region; p Let λ represent the discrete emission source term of the p-th sub-region. p This represents the boundary term of the p-th sub-region.

[0069] The fourth step is to couple all sub-regions according to the boundary conditions and construct single-precision and double-precision interface equations.

[0070] In this embodiment, the expression for the boundary condition is:

[0071]

[0072] Where n represents the unit normal vector at the interface; μ represents the magnetic permeability; Λ p Γ represents an unknown value; p Indicates the internal boundary.

[0073] In this embodiment, the expression for the single-precision interface equation is:

[0074] F s λ s =b s

[0075] Among them, F s λ represents the coefficient matrix formed by coupling the discrete equations of a sub-region under single-precision conditions. s In single-precision, based on the boundary term λ p The resulting unknown terms to be solved; b s This represents the known terms formed by the coupling of source terms in the discrete equations of a sub-region under single-precision conditions.

[0076] In this embodiment, the expression for the double-precision interface equation is:

[0077] F d λ d =b d

[0078] Among them, F d λ represents the coefficient matrix formed by the coupling of discrete equations in a sub-region under double precision. d In double precision, based on the boundary term λ p The resulting unknown terms to be solved; b d This represents the known terms formed by the coupling of source terms in the discrete equations of a subregion under double precision.

[0079] The fifth step involves setting a single-precision iteration convergence threshold k and a relative change threshold m for the residual norm between adjacent iterations. The single-precision interface equation is then iteratively solved. If the single-precision iteration residual norm is less than k or the relative change in residual norm between adjacent iterations is less than m, the iteration is corrected using the double-precision interface equation. In this embodiment, the single-precision iteration residual norm is the Arnoldi residual norm used in the generalized minimum residual method calculation.

[0080] Step 6: Preset the double-precision iteration convergence threshold p. If the relative iteration residual norm calculated by the double-precision interface equation is less than p, then stop the iteration. If the relative iteration residual norm calculated by the double-precision interface equation is greater than or equal to p, use the corrected parameters to perform single-precision interface equation iteration again until the convergence condition is met.

[0081] The expression for the relative iterative residual norm of double precision is:

[0082]

[0083] Where, λ n This represents the update amount in the nth iteration of double precision.

[0084] Step 7: Substitute the solution of the interface equation into the discrete equation of the sub-region containing the receiving point to calculate the electromagnetic response at the receiving point. The double-precision correction parameter is: the double-precision updated λ... n The iterative error of calculating the double-precision interface equation.

[0085] The following is a real-world example:

[0086] This embodiment provides a domain decomposition solution method for calculating the semi-airborne electromagnetic response. A 1000-meter-long conductor source is used, with a transmitting current of 1A. The effective area of ​​the receiving coil is 1 square meter, the aircraft's flight altitude is 30 meters, and three measuring lines are parallel to the long conductor source, at distances of 200 meters, 600 meters, and 1000 meters from the source, respectively.

[0087] A geological model of a low-resistivity aquifer is designed, and it is subjected to unstructured mesh generation and sub-region division. Discrete equations of the electric field control equations are independently established for each sub-region, and stored in both single-precision and double-precision formats. All sub-regions are coupled according to boundary conditions to construct single-precision and double-precision interface equations. A preset single-precision iteration convergence threshold of 10 is used. -6 The single-precision interface equation is iteratively solved by taking a threshold of 1% relative change in error between adjacent iterations. If the single-precision iteration error is less than 10%, the solution is considered successful. -6 If the relative change in error between adjacent iterations is less than one percent, then the iteration is corrected using the double-precision interface equation; the preset double-precision iteration convergence threshold is 10. -6If the iteration error norm of the double-precision interface equation calculation is less than 10 -6 Then stop iterating; when the iteration error norm calculated by the double-precision interface equation is greater than or equal to 10... -6 Then, the corrected parameters are used to iteratively solve the single-precision interface equation until the convergence condition is met; the solution of the interface equation is substituted into the discrete equation of the sub-region containing the receiving point to calculate the electromagnetic response at the receiving point. Figure 2 (a) Comparison of the computational accuracy and efficiency of the method of the present invention with that of conventional methods. Figure 2 (b) Comparison of computational efficiency between the method of the present invention and the conventional method.

[0088] Through the above examples of this invention, compared with the traditional domain decomposition method in the prior art, the improved domain decomposition method of this invention can effectively improve computational efficiency while ensuring computational accuracy. The practical application results of this invention demonstrate that the improved domain decomposition method proposed in this invention has significant application value and practical significance in electromagnetic exploration.

[0089] The above embodiments are merely preferred embodiments of the present invention and are not intended to limit the scope of protection of the present invention. Any changes made based on the design principles of the present invention, or any non-creative modifications made thereon, shall fall within the scope of protection of the present invention.

Claims

1. A mixed-precision domain decomposition solution method based on electromagnetic exploration using electrical sources, characterized in that, Includes the following steps: Based on the underground geoelectric structure, a forward model of the area to be tested is established and tetrahedral mesh is generated to obtain the node and mesh information of the tetrahedral mesh. Assign resistivity values ​​to any tetrahedral mesh; Divide the tetrahedral mesh into several non-overlapping sub-regions; Discrete equations for the electric field control equations are established for any sub-region; The sub-regions are coupled according to the boundary conditions, and single-precision interface equations and double-precision interface equations are constructed. The minimum residual method is used to solve the single-precision interface equation using single-precision iterative solution. In single-precision iteration, when the Arnoldi residual norm stagnates or reaches a preset condition, the single-precision iteration residual calculated by the double-precision interface equation is used to correct the single-precision iteration, and the corrected parameters are obtained. The corrected parameters are solved iteratively using a single-precision interface equation until the convergence condition is met. The solution to the interface equation is obtained; Substituting the solution of the interface equation into the corresponding receiving point within the subregion yields the electromagnetic response of the receiving point.

2. The hybrid-accuracy domain decomposition solution method based on electromagnetic exploration using electrical sources according to claim 1, characterized in that, In single-precision iteration, when the Arnoldi residual norm stagnates or reaches a preset condition, the single-precision iteration is corrected using the double-precision iterative residual calculated by the double-precision interface equation, resulting in corrected parameters, including: The convergence threshold k for single-precision iteration and the relative change threshold m for the residual norm of adjacent iterations are preset; The single-precision interface equation is solved by single-precision iterative solution. If the single-precision iteration residual norm is less than the convergence threshold k or the relative change in the residual norm between adjacent iterations is less than the relative change threshold m, then the single-precision iteration is corrected using the double-precision iteration residual calculated by the double-precision interface equation, and the corrected parameters are obtained.

3. The hybrid-accuracy domain decomposition solution method based on electromagnetic exploration using electrical sources according to claim 1 or 2, characterized in that, The correction of single-precision iterations using the double-precision residuals calculated with double-precision interface equations also includes: The convergence threshold p for double-precision iteration is preset; If the relative iterative residual norm calculated by the double-precision interface equation is less than the convergence threshold p, then the iteration is stopped. If the relative iterative residual norm calculated by the double-precision interface equation is greater than or equal to the convergence threshold p, then the corrected parameters are solved iteratively using the single-precision interface equation until the convergence condition is met.

4. The hybrid-accuracy domain decomposition solution method based on electromagnetic exploration using electrical sources according to claim 3, characterized in that, The expression for the electric field control equation is as follows: in, Represents the differential operator; E represents the electric field strength; ω represents the angular frequency; μ represents the permeability; σ represents the conductivity; J s This indicates the intensity of the external electrical source.

5. The hybrid-accuracy domain decomposition solution method based on electromagnetic exploration using electrical sources according to claim 4, characterized in that, The expression for the discrete equation is: (A p +iωB p )E p =-iωS p -l p Where p represents the p-th sub-region; A p B represents the stiffness matrix of the p-th subregion; p S represents the quality matrix of the p-th sub-region; p Let λ represent the discrete emission source term of the p-th sub-region. p This represents the boundary term of the p-th sub-region.

6. The hybrid-accuracy domain decomposition solution method based on electromagnetic exploration using electrical sources according to claim 5, characterized in that, The expression for the boundary condition is: Where n represents the unit normal vector at the interface; μ represents the magnetic permeability; Λ p Γ represents an unknown value; p Indicates the internal boundary.

7. The hybrid-accuracy domain decomposition solution method based on electromagnetic exploration using electrical sources according to claim 1, characterized in that, The expression for the single-precision interface equation is: F s l s =b s Among them, F s λ represents the coefficient matrix formed by coupling the discrete equations of a sub-region under single-precision conditions. s In single-precision, based on the boundary term λ p The resulting unknown terms to be solved; b s This represents the known terms formed by the coupling of source terms in the discrete equations of a sub-region under single precision. The expression for the double-precision interface equation is: F d l d =b d Among them, F d λ represents the coefficient matrix formed by the coupling of discrete equations in a sub-region under double precision. d In double precision, based on the boundary term λ p The resulting unknown terms to be solved; b d This represents the known terms formed by the coupling of source terms in the discrete equations of a subregion under double precision.

8. The hybrid-accuracy domain decomposition solution method based on electromagnetic exploration using electrical sources according to claim 2, characterized in that, The single-precision iterative residual norm is the Arnoldi residual norm used in the minimum residual method calculation process.

9. The hybrid-accuracy domain decomposition solution method based on electromagnetic exploration using electrical sources according to claim 1, characterized in that, The double-precision iterative residual calculated from the double-precision interface equation is the residual ε of the double-precision equation system to be solved, and its expression is: e=F d l n -b d Where, λ n F represents the iteration update amount of the nth iteration of double precision; d b represents the coefficient matrix formed by coupling the discrete equations of the sub-region under double precision; d This represents the known terms formed by the coupling of source terms in the discrete equations of a sub-region under double precision. The double-precision iterative residual calculated using the double-precision interface equation is used to correct the single-precision iteration, resulting in corrected parameters. These corrected parameters include the iteration update amount λ for the nth double-precision iteration. n The iterative error of calculating the double-precision interface equation.

10. The hybrid-accuracy domain decomposition solution method based on electromagnetic exploration using electrical sources according to claim 9, characterized in that, The relative iterative residual norm L calculated by the double-precision interface equation is expressed as follows:

Images