Large-Scale Parallel Forward Modeling Method for Electromagnetic Anisotropy of the Base Solved by PDE

Abstract

The present application relates to an electromagnetic anisotropic large-scale parallel forward modeling method based on a PDE solving base. The method includes: constructing an anisotropic electromagnetic model according to the constitutive relationship between the electric and magnetic fields when facing anisotropic media; the constitutive parameters of the constitutive relationship include the dielectric constant, magnetic permeability, and conductivity of the anisotropic media; inputting the configuration file information of the electromagnetic problem to be solved into the PDE solving base, and the PDE solving base loads the electromagnetic model plug-in according to the configuration file information, constructs an electromagnetic field solving module, divides the calculation area into multiple sub-domains, allocates processes from the process set to each sub-domain, and in the processes corresponding to each sub-domain, parallelly solves the anisotropic electromagnetic model of the sub-domain according to the multi-level multi-granularity parallel algorithm and the configuration information in the electromagnetic field solving module, and outputs the electromagnetic field results of the sub-domain; obtaining the forward modeling result according to the electromagnetic field results of each sub-domain. Using this method can improve the efficiency and reliability of large-scale electromagnetic forward modeling.

Description

Technical Field

[0001] The present application relates to the technical field of numerical simulation, and particularly to an electromagnetic anisotropic large-scale parallel forward modeling method based on a PDE solving base. Background Art

[0002] Based on the electromagnetic properties of media, especially the differences in conductivity, permittivity, permeability, etc., a wide range of applications have emerged, especially in the fields of geophysics, communication technology, and imaging technology. In the field of geophysics, electromagnetic exploration is a technique that uses the conductivity differences of underground media to detect geological structures. This method searches for ore deposits, identifies geological structures, and solves geological problems by observing the spatial distribution laws and time characteristics of artificial or natural electric fields and electromagnetic fields. In communication technology, artificial media based on layer structures can achieve unidirectional electromagnetic wave transmission, which performs excellently in anti-backscattering and is applicable to fields such as communication and imaging. Electromagnetic properties are also widely used in imaging technology, such as magnetic resonance imaging (MRI) in medical imaging. MRI uses the different responses of human tissues to magnetic fields to generate detailed internal images. In addition, electromagnetic stealth technology relies on specially designed electromagnetic media to reduce radar wave reflection, thereby improving stealth performance.

[0003] The above research and applications all focus on the electromagnetic property differences between the target body and the background medium, and the electromagnetic properties of the target body are defaulted to be isotropic. Taking electromagnetic exploration as an example: In recent decades, electromagnetic exploration has mainly focused on the research of isotropic media, that is, it is assumed that the earth's media are isotropic. Based on this assumption of isotropic media, there have been many studies and applications in electromagnetic method exploration. However, anisotropy objectively exists, and many studies have proven that underground media are anisotropic. Rocks in nature all have varying degrees of anisotropy, and factors such as the earth's tectonic stress field, the earth's medium deformation zone, rock fractures, pore water, and geological sedimentation are all causes of anisotropy.

[0004] In recent years, in the field of geophysical exploration and other fields, with the deepening of exploration depth, the increase in detection area, and the improvement of accuracy requirements, the computational scale of numerical simulation has increased sharply, and the computational time has also increased exponentially. In order to ensure the timeliness of applications and be able to quickly solve problems, the research on parallel algorithms has begun to attract the attention of scholars. High-performance computing can accelerate the solution speed or solve large-scale problems that ordinary computers cannot solve.

[0005] Partial differential equations (PDEs) play a crucial role in various fields such as physics, engineering, financial mathematics, and biology, and are used to describe various natural phenomena and processes. For example, the Navier-Stokes equations in fluid mechanics, the Maxwell equations in electromagnetism, and the heat equation in heat conduction processes are all typical application examples of partial differential equations. Common PDE solving methods include the finite difference method (FDM), the finite element method (FEM), the finite volume method (FVM), etc. With the development of computer technology, numerical simulation has become one of the important means to solve complex PDE problems.

[0006] Most of the existing research and applications of electromagnetic methods are based on the isotropy of the target body, and there is little research on the differences in the electromagnetic characteristics of the target body itself. Taking electromagnetic exploration in geophysical exploration as an example, most are based on the assumption of conductivity isotropy. However, in actual exploration, affected by pore water, sedimentation, plate tectonic movements, etc., the underground medium shows a certain degree of conductivity anisotropy. It is obviously not feasible to analyze and interpret anisotropic data based on isotropic theory. Although there have been some studies on anisotropy in recent years, most are forward modeling for conductivity anisotropy, and there are very few electromagnetic forward models considering permeability anisotropy. At the same time, the electromagnetic forward modeling considering the anisotropy of conductivity, permeability, and permittivity is still blank. And the existing electromagnetic parallel algorithms face problems such as single parallel granularity, single parallel means, and poor parallel scalability. Most can only be parallelized on a small scale, and very few large-scale parallelizations avoid the low-frequency problems with limited parallel efficiency. However, small-scale parallelization obviously cannot meet the timeliness requirements.

[0007] Although significant progress has been made in existing PDE solving software frameworks, there are still some limitations and challenges. To support a wide range of PDE types, many frameworks are designed to be more general, but this often comes at the cost of computational efficiency. As the problem scale grows, the memory and computing power on a single machine gradually become bottlenecks. Although some frameworks have begun to utilize multi-core processors and GPUs to accelerate calculations, the support in a distributed computing environment is still insufficient, limiting the effective solution of large-scale problems. Summary of the Invention

[0008] Based on this, it is necessary to provide a large-scale parallel electromagnetic anisotropic forward modeling method based on a PDE solving base for the above technical problems.

[0009] A large-scale parallel electromagnetic anisotropic forward modeling method based on a PDE solving base, the method includes:

[0010] Construct an anisotropic electromagnetic model according to the constitutive relationship between the electric field and the magnetic field in an anisotropic medium; the constitutive parameters of the constitutive relationship include the permittivity, permeability, and conductivity of the anisotropic medium;

[0011] Input the configuration file information of the electromagnetic problem to be solved into the PDE solving base. The PDE solving base loads the electromagnetic model plug-in according to the configuration file information, constructs an electromagnetic field solving module, divides the calculation region into multiple subdomains, allocates processes from the process set to each subdomain, and in the processes corresponding to each subdomain, parallelly solve the anisotropic electromagnetic model corresponding to the subdomain according to the preset multi-level and multi-granularity parallel algorithm and the configuration information in the electromagnetic field solving module, and output the electromagnetic field results of the subdomain;

[0012] Obtain the forward modeling result according to the electromagnetic field results of each subdomain.

[0013] The above electromagnetic anisotropic large-scale parallel forward modeling method based on the PDE solving base constructs an electromagnetic model by considering the permittivity, permeability and conductivity of anisotropic media, which can make the model more in line with the actual underground medium situation and improve the accuracy of the forward modeling result. Then, input the configuration file information of the electromagnetic problem to be solved into the PDE solving base, load the electromagnetic model plug-in and construct an electromagnetic field solving module, which can utilize the advantages of the PDE solving base and integrate various numerical algorithms, such as finite element, finite difference algorithms, etc. This enables the program to flexibly select appropriate algorithms according to different calculation scenarios, enhancing the applicability and flexibility of the program. Divide the calculation region into multiple subdomains and allocate processes, and combine the multi-level and multi-granularity parallel algorithm to parallelly solve the anisotropic electromagnetic model corresponding to the subdomain, which can flexibly allocate computing resources at different levels, greatly improve the computing efficiency, break through the bottleneck of single-machine memory and computing power, and solve the problem of long large-scale computing time. Finally, integrate the electromagnetic field results of each subdomain to obtain the forward modeling result, which can summarize the calculation results of each part and output a comprehensive and high-precision forward modeling result. The embodiment of the present invention can improve the efficiency and reliability of large-scale electromagnetic forward modeling. Description of the Drawings

[0014] Figure 1 It is a schematic flow chart of the electromagnetic anisotropic large-scale parallel forward modeling method based on the PDE solving base in one embodiment;

[0015] Figure 2 It is a schematic diagram of the corresponding relationship between the edges and nodes of the Nédélec vector shape function of tetrahedral elements in one embodiment;

[0016] Figure 3 It is a schematic diagram of the ParMETIS division of the anomaly model in one embodiment;

[0017] Figure 4 It is a schematic flow chart of the multi-level and multi-granularity algorithm based on MPI and OpenMP in one embodiment. Detailed Embodiment

[0018] In order to make the objectives, technical solutions, and advantages of this application clearer, the following further details this application in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely used to explain this application and are not used to limit this application.

[0019] In one embodiment, as Figure 1 shown, an electromagnetic anisotropic large-scale parallel forward modeling method based on a PDE solving base is provided, including the following steps:

[0020] Step 102: Construct an anisotropic electromagnetic model according to the constitutive relationship between the electric field and the magnetic field in the case of an anisotropic medium.

[0021] The constitutive parameters of the constitutive relationship include the permittivity, permeability, and conductivity of the anisotropic medium. It can be understood that using the permittivity, permeability, and conductivity of the anisotropic medium as constitutive parameters to construct the model fully considers the anisotropic changes of the target body, and at the same time takes into account the anisotropy of conductivity, permeability, and permittivity, making the electromagnetic model closer to the real earth medium, capable of improving the accuracy of electromagnetic forward modeling results, and being able to more comprehensively consider the characteristics of underground media in actual exploration, providing a basis for subsequent accurate analysis and interpretation of electromagnetic data.

[0022] Step 104: Input the configuration file information of the electromagnetic problem to be solved into the PDE solving base. The PDE solving base loads the electromagnetic model plug-in according to the configuration file information, constructs an electromagnetic field solving module, divides the calculation area into multiple sub-domains, allocates processes from the process set to each sub-domain, and in the processes corresponding to each sub-domain, parallelly solve the anisotropic electromagnetic model corresponding to the sub-domain according to the pre-set multi-level multi-granularity parallel algorithm and the configuration information in the electromagnetic field solving module, and output the sub-domain electromagnetic field results.

[0023] The PDE solving base includes a core module, a plug-in management module, and a numerical calculation module. The specific component structures and functions of each module are described in the Chinese invention patent "Multi-physics field coupling analysis system and method based on a general PDE solving base" with the application number 2025100793285 in the prior application, and will not be elaborated here. Based on the PDE programming framework, the finite element and finite difference algorithms are encapsulated, and appropriate numerical methods can be selected according to specific calculation scenarios, increasing the flexibility of the program. The PDE solving base uses a multi-level and multi-granularity parallel algorithm, combining various parallel strategies. For example, different parallel programming techniques such as MPI and OpenMP can be adopted at different calculation levels, making full use of the advantages of distributed memory systems and shared memory systems to achieve multi-node and multi-thread parallel computing. At the same time, the distributed storage technology can disperse large-scale data and store it on different nodes, avoiding the limitation of the memory of a single machine, being able to handle ultra-large-scale electromagnetic forward problems, and meeting the requirements of large-scale electromagnetic forward for timeliness.

[0024] Step 106: Obtain the forward result based on the electromagnetic field results of each subdomain.

[0025] Using the obtained high-precision electromagnetic forward result can provide accurate data support for electromagnetic exploration in fields such as mineral resource exploration and marine geophysics, helping to more accurately identify targets and improve the success rate and accuracy of exploration.

[0026] In the above electromagnetic anisotropic large-scale parallel forward method based on the PDE solving base, by constructing an electromagnetic model considering the permittivity, permeability, and conductivity of anisotropic media, the model can better fit the actual underground medium situation, improving the accuracy of the forward result. Then, inputting the configuration file information of the electromagnetic problem to be solved into the PDE solving base, loading the electromagnetic model plug-in and constructing the electromagnetic field solving module can utilize the advantages of the PDE solving base to integrate various numerical algorithms, such as finite element and finite difference algorithms. This enables the program to flexibly select appropriate algorithms according to different calculation scenarios, enhancing the applicability and flexibility of the program. Dividing the calculation area into multiple subdomains and allocating processes, and combining the multi-level and multi-granularity parallel algorithm to parallelly solve the anisotropic electromagnetic models corresponding to the subdomains can flexibly allocate computing resources at different levels, greatly improving the computing efficiency, breaking through the bottleneck of single-machine memory and computing power, and solving the problem of long calculation time for large-scale calculations. Finally, integrating the electromagnetic field results of each subdomain to obtain the forward result can summarize the calculation results of each part and output a comprehensive and high-precision forward result. The embodiments of the present invention can improve the efficiency and reliability of large-scale electromagnetic forward.

[0027] The present invention is implemented based on the PDE programming framework. When building the PDE solution base, first, the two basic algorithms of the finite element method and the finite difference method are modularly encapsulated. The finite element method has advantages in dealing with complex geometric shapes and multi-physical boundary conditions, while the finite difference method has higher computational efficiency for problems with regular geometric shapes. Through modular encapsulation, these two algorithms can be independently developed, tested, and optimized, and can work together when needed. For electromagnetic problems with regular geometric shapes, we use structured grid meshing and select the finite difference method for solution, which can reduce memory occupancy and improve computational efficiency. For electromagnetic problems involving complex geometric shapes and various physical property materials, we use unstructured grid meshing and select the finite element method for solution, and improve the solution accuracy of specific regions by locally refining the grid in the regions of interest or where the field quantities change violently.

[0028] Electromagnetic problems are usually described by Maxwell's equations. For the solution of electromagnetic problems, first, prepare the model grid file of the electromagnetic problem to be solved externally, and set the physical property parameters and solution parameters of the electromagnetic model through the configuration file; after the program starts, the core module (Base module) in the PDE solution base starts to work. First, it calls the command-line parameter parsing module to read the path of the configuration file input from the command line, and then the PDE solution base completes the reading of the configuration file and the model grid file by calling the parallel communication and file I / O module. According to the settings in the configuration file, the PDE solution base dynamically loads the electromagnetic model plug-in and its related components, completes the integration of the plug-in and the core module, and constructs the electromagnetic field solution module. The electromagnetic model in the electromagnetic field solution module is a mathematical model established for electromagnetic problems, and this mathematical model is established based on Maxwell's equations. The electromagnetic field solution module extracts the curl operator ( ) and the vector dot product operator (·) from the discrete operator library already discretized in the PDE base according to the Helmholtz equation in the mathematical model, and combines the operators through operator overloading to form a complete discrete equation, completes the finite element discretization of the electromagnetic mathematical model, obtains a linear equation set describing the electromagnetic field distribution, and with the help of the linear system solver in the PDE solution base, completes the solution of the linear equation set to obtain the electromagnetic field strength on each unit.

[0029] In one embodiment, loading the electromagnetic model plug-in according to the configuration file information and constructing the electromagnetic field solution module includes: loading the electromagnetic model plug-in according to the input configuration file information, and completing the integration of the plug-in and the core module through the standardized interface to construct the electromagnetic field solution module.

[0030] In one embodiment, the configuration file information includes the electromagnetic model subclass and related components; the electromagnetic model subclass includes the anisotropic electromagnetic model and boundary conditions; the related components include the grid file processing sub-plug-in, the solution algorithm sub-plug-in, and the PDE discretization sub-plug-in.

[0031] In one embodiment, according to the constitutive relationship between the electric field and the magnetic field in an anisotropic medium, constructing an anisotropic electromagnetic model includes: obtaining the definite solution form of Maxwell's equations according to the constitutive relationship between the electric field and the magnetic field in an anisotropic medium, converting the definite solution form of Maxwell's equations into the control differential equation of the electric field, and obtaining the anisotropic electromagnetic model as:

[0032]

[0033] Wherein, is the electric field strength, represents the permittivity of the anisotropic medium, represents the permeability of the anisotropic medium, represents the conductivity of the anisotropic medium, is the curl operator, is the angular frequency, is the imaginary unit.

[0034] Specifically, Maxwell's equations are a set of fundamental equations that govern all macroscopic electromagnetic phenomena. For a general time-varying field, the differential form of Maxwell's equations can be written as

[0035] (1)

[0036] (2)

[0037] (3)

[0038] (4)

[0039] In the formula:

[0040]

[0041] Another fundamental equation is the current continuity equation, which can be written as

[0042] (5)

[0043] When the field quantities in Maxwell's equations are single-frequency harmonic functions, we obtain a time-harmonic field. Using the complex phase factor representation method and the time convention , is the angular frequency. Equations (1), (2), and (5) can be written as

[0044] (6)

[0045] (7)

[0046] (8)

[0047] In this case, the electric field and the magnetic field must obviously exist simultaneously and interact with each other. The static field is the limiting case when the angular frequency tends to zero for the time-harmonic field. Among the first five equations, only three are independent. Taking the divergence of both sides of equation (1) gives equation (4), and taking the divergence of both sides of equation (2) gives equation (5). Since the number of equations is less than the number of unknowns, the three independent equations are in an underdetermined form. When the constitutive relations between the field quantities are determined, Maxwell's equations become in a determined form. The constitutive relations describe the macroscopic properties of the medium under consideration. For a simple medium, there are

[0048] (9)

[0049] (10)

[0050] (11)

[0051] In the above equations, the constitutive parameters 、 and represent the permittivity (F / m), permeability (H / m), and conductivity (S / m) of the medium, respectively. In the algorithm of the present invention, when calculating according to the anisotropic medium, the corresponding permittivity, permeability, and conductivity parameters are all tensors, which are represented by 、 and respectively. When directly dealing with the time-harmonic field situation using the electric field and the magnetic field, starting from Maxwell's equations containing the electric field and the magnetic field, a control differential equation containing only any one field quantity is derived. Using the constitutive relations (9)-(11), taking the curl of both sides of equation (6), and then substituting equation (7) into it, we can get

[0052] (12)

[0053] Since the in factor (12) is the conduction current formed under the action of the electric field, combining with equation (11), equation (12) can be further written as

[0054] (13)

[0055] The above equation is the double curl equation of the electric field including the displacement current, and is also called the electric field wave equation.

[0056] In one embodiment, the anisotropic electromagnetic model corresponding to the subdomain is solved in parallel according to the preset multi-level and multi-granularity parallel algorithm and the configuration information in the electromagnetic field solving module, and the output of the subdomain electromagnetic field results includes: each subdomain performs data preprocessing in parallel according to the subdomain grid data and subdomain matrix data stored in the corresponding process; each subdomain performs unit analysis on the anisotropic electromagnetic model corresponding to the subdomain in parallel according to the configuration information in the electromagnetic field solving module and the preprocessed data, and performs unit analysis on the processed subdomain grid in parallel to obtain a system of linear equations; each subdomain solves the system of linear equations in parallel to obtain the subdomain electromagnetic field distribution; each subdomain performs post-processing on the subdomain electromagnetic field distribution in parallel to obtain the subdomain electromagnetic field results. In this embodiment, grid distributed storage is adopted to effectively reduce the memory consumption on a single node and achieve load balancing.

[0057] It should be noted that currently, the mainstream parallel programming environments are MPI, OpenMP, and CUDA / GPU. Among them, MPI is a message-passing programming model suitable for distributed memory systems, and its ultimate goal is to serve the goal of inter-process communication. It mainly uses coarse-grained parallelism, allocating tasks to multiple nodes of the distributed memory system to complete parallel computing. However, in the large-scale parallel forward modeling of electromagnetic anisotropy, relying solely on coarse-grained parallelism cannot meet complex computing requirements, and there are efficiency issues in the face of situations requiring finer-grained parallelism. Moreover, it only focuses on inter-process communication and has relatively single functions when dealing with complex electromagnetic forward problems. The OpenMP parallel programming technology is suitable for shared memory systems and adopts a fork-join parallel execution mode. When encountering large amounts of computation, it generates a thread parallel region to complete parallel computing tasks. However, for large-scale parallel computing in electromagnetic methods, especially in scenarios involving multiple nodes and distribution, its adaptability is limited and it cannot be well applied to large-scale computing tasks of distributed memory systems. GPUs focus more on data parallel computing and are not suitable for logical processing. Moreover, they need to be used through the CUDA architecture. In electromagnetic forward modeling, they may not be able to handle the complex electromagnetic model logic well, and their application scenarios are relatively limited and cannot independently solve all problems in the large-scale parallel forward modeling of electromagnetic methods.

[0058] In this embodiment, a multi-level hybrid parallel method is adopted, combining the advantages of various parallel modes. It not only utilizes the parallel advantage of a distributed memory system similar to MPI for task allocation, but also combines the thread parallel advantage of OpenMP in a shared memory system in some links to achieve more flexible parallel computing, so as to adapt to the different requirements of parallel granularity for different computing tasks in large-scale parallel forward modeling of the electromagnetic method. The method of the present invention stores grid data and equation system matrices distributively. First, the global grid is divided into multiple subdomains, and each process only stores the grid data of one subdomain, and the grid data at the junction of adjacent subdomains is transmitted through MPI. In order to obtain better parallel efficiency, the distributed grid storage needs to meet load balancing. Considering that the number of grids is usually proportional to the computational complexity, the number of grids assigned to each subdomain is close. The method of the present invention adopts ParMETIS that can be executed in parallel to implement partitioning. As Figure 3 shown, the entire three-dimensional grid area of the anomaly body model is divided into multiple subdomains by using ParMETIS. Finally, a new grid data structure is reconstructed in each subdomain, and the cells or edges in the same subdomain are numbered continuously.

[0059] The multi-level and multi-granularity parallel algorithm involves seven parts, namely mesh generation, input, subdomain division, preprocessing, finite element calculation, solution of equations, postprocessing, and output: ① Use open-source software Tetgen or Gmsh to generate mesh files; ② In the input part, read the mesh data and calculation parameters and distribute them to all processes; ③ Use ParMETIS to divide the calculation domain into multiple subdomains. The processes in the same subdomain only store the mesh information of that subdomain to prepare for subdomain parallelization; ④ Perform preprocessing on the data, continuously number the elements or edges in the same subdomain. The data structures involved include elements, nodes, edges, edges on the boundary, and some other data; ⑤ The finite element calculation part includes element analysis, matrix assembly, imposition of the first kind of boundary conditions, and matrix format conversion; ⑥ Use Mumps or SuperLU to solve large sparse linear equations to obtain the electric field; ⑦ Calculate the magnetic field, tensor impedance, apparent resistivity, and phase through the electric field calculation and output the results. To fully improve the calculation efficiency, parallel computing is adopted for a total of six parts, namely subdomain division, preprocessing, finite element calculation, solution of equations, and postprocessing. In the part of solving the linear system, first evenly divide the processes into N process groups and evenly distribute the frequencies to each process group. Due to the natural parallelism of frequencies, the calculations between process groups are independent of each other. Within a process group, all processes multiply the stored matrix elements by the frequency coefficient to generate a frequency-dependent sparse matrix, and parallelly solve the sparse linear system through a direct solver. In the postprocessing part, the calculations between process groups are independent of each other. Each process group only processes the frequencies assigned to it. The main process collects the solutions of other processes in the group, calculates the impedance, apparent resistivity, and phase response, and finally outputs the results. If a process group is assigned multiple frequencies, the frequencies need to be calculated sequentially until all frequencies are calculated.

[0060] As Figure 4 shown, a flow schematic diagram of an MPI and OpenMP multi-level and multi-granularity hybrid parallel algorithm is provided. Among them, the parallel modes include thread-level parallelism and process-level parallelism. OpenMP thread-level parallelism: Use OpenMP to implement parts such as element analysis, matrix assembly, and imposition of the first kind of boundary conditions, mainly involving the parallelism of loops without data dependencies: ① First, enable OpenMP; ② Read the data required for the loop; ③ Specify or dynamically allocate the number of threads for the loop; ④ Determine whether variables are private or shared; ⑤ The allocated number of threads calculates together for the DO / for loop; ⑥ Finally, the OpenMP parallel calculation ends. MPI process-level parallelism: ① First, use ParMetis to parallelly divide the calculation area into multiple subdomains; ② Perform data preprocessing in parallel for each subdomain. Each subdomain corresponds to at least one process, and the processes in each subdomain only store the mesh and matrix data of that subdomain; ③ All processes in each subdomain parallelly solve the equations of that subdomain; ④ Finally, the main process collects the data of other processes in the same subdomain.

[0061] In one embodiment, unit analysis is performed in parallel on the anisotropic electromagnetic model corresponding to the subdomain, and unit analysis is performed in parallel on the subdomain grid obtained after processing to obtain a system of linear equations, including: in a single subdomain, multiple units to be processed are assigned to different processes, each process processes the corresponding unit to be processed in parallel, and different processes communicate through a message passing interface; in each process, different iterations in the unit loop are assigned to different threads, and multi-threaded parallel execution is used to achieve parallel unit analysis of the anisotropic electromagnetic model corresponding to the subdomain, and unit analysis is performed in parallel on the subdomain grid obtained after processing to obtain a system of linear equations.

[0062] In one embodiment, the system of linear equations is solved in parallel for each subdomain to obtain the electromagnetic field distribution of the subdomain, including: the processes corresponding to each subdomain are divided into multiple process groups, frequencies are assigned to each process group, and all processes within the process group multiply the stored matrix elements by the corresponding frequencies to obtain a frequency-dependent sparse matrix; each process group uses a direct solver to solve the system of linear equations composed of the frequency-dependent sparse matrix in parallel. In each process group, the master process collects the solutions of other processes within the process group to obtain the electromagnetic field distribution of the subdomain. In this embodiment, frequency-point parallelism is used. For problems that need to be calculated at multiple frequency points, the calculation tasks at different frequency points can be allocated among different processing units. By calling a third-party solver, the system of linear equations is solved in parallel. Inside the third-party solver, MPI is used for communication and coordination between processes, while OpenMP is used for multi-threaded parallelization within the same node, thus realizing hybrid parallel computing.

[0063] In one embodiment, post-processing is performed in parallel on the electromagnetic field distribution of each subdomain to obtain the electromagnetic field results of the subdomain, including: each subdomain calculates impedance, apparent resistivity, and phase response based on the electromagnetic field distribution of the subdomain to obtain the electromagnetic field results of the subdomain. In this embodiment, for the computationally intensive process of constructing an impedance matrix for computational electromagnetics, the calculation tasks are divided, and each process calculates the matrix elements of its responsible part respectively, and then the results are merged through communication. This not only accelerates the matrix filling process but also allows for the handling of larger-scale problems.

[0064] In one embodiment, unit analysis is performed in parallel on the anisotropic electromagnetic model corresponding to the subdomain, and unit analysis is performed in parallel on the subdomain grid obtained after processing to obtain a system of linear equations, including: the space is discretized using tetrahedral meshes, the integral of the computational region is decomposed into the accumulation of discrete element integrals and processed using the weighted residual method, and unit analysis is performed in parallel on the subdomain grid obtained after processing to obtain a system of linear equations.

[0065] Specifically, the weighted residual method is used to derive the finite element equations from the governing equations. The weighted residual method substitutes the approximate solution into the governing equations. Since the approximate solution is not exactly the same as the true solution, a residual will be generated. The weighted residual method makes the residual minimum through mathematical methods to obtain an approximate solution that is infinitely close to the true solution. Equation (13) is converted into a linear form, and the unknown quantity to be solved in the equation is the electric field vector, and the corresponding residual can be expressed as

[0066] (14)

[0067] By taking the inner product of the residual and the weighting function and setting the integral within each unit region composing the model equal to 0

[0068] (15)

[0069] Substituting equation (14) into equation (15) gives the weighted residual equation for the electric field

[0070] (16)

[0071] According to the vector identity (17)

[0072] (17)

[0073] and the divergence theorem (Gauss's theorem)

[0074] (18)

[0075] Equation (16) can be transformed into

[0076] (19)

[0077] According to Stokes' theorem: the integral of the normal vector of the curl of a vector along a surface is equal to the line integral of the tangential component of the vector along the boundary of the surface, i.e.,

[0078] (20)

[0079] Applying Stokes' theorem to the first term in equation (19), since the calculation of this term is for the closed surface of the element, the closed surface can be divided into upper and lower parts by a closed curve, and then two integrations in opposite directions are performed along the closed curve ( ), so the result of this term's integration is 0

[0080] (21)

[0081] Thus, equation (19) can be simplified to

[0082] (22)

[0083] The above equation is the weighted residual equation satisfied within the element, and equation (22) is also satisfied for the entire study area.

[0084] The space is discretized using tetrahedral meshes, and the Nédélec vector shape functions shown in equation (23) are selected

[0085] (23)

[0086] denotes the edge number of the tetrahedron ( ), and denote and are the nodal shape functions of denotes the th edge length, Figure 2 describes the correspondence between the edge and the nodes. According to Figure 1 the correspondence between the edge and the nodes, combined with the Nédélec vector shape function expression (23), it can be known that

[0087] (24)

[0088] the expression of the nodal shape function (node number )

[0089] (25)

[0090] where denotes the element volume,[[]] denotes an arbitrary point within the tetrahedron,[[]] ( ) denote the coefficients related to the element respectively, and the specific numerical calculations are as follows:

[0091]

[0092]

[0093]

[0094]

[0095]

[0096]

[0097]

[0098] The above equation gives of expression. Expressions for other values can be deduced by analogy. Note the plus and minus signs in front of the coefficients. where represents the row, represents the column, and this is the cofactor in the determinant.

[0099] Taking the gradient of the nodal shape function in Equation (25) gives:

[0100] (26)

[0101] Substituting Equation (26) and Equation (25) into Equation (24), the specific expression of the Nédélec vector shape function can be obtained as follows:

[0102]

[0103]

[0104]

[0105]

[0106]

[0107] Decomposing the integral over the entire computational domain into the form of discrete element integral summation, Equation (22) becomes

[0108] (27)

[0109] In the method of weighted residuals, different weight functions are adopted, corresponding to different calculation methods, such as the collocation method, the least squares method, and the Galerkin method, etc. Since the Galerkin method can be well used for the finite element method, the Galerkin method is selected. In the Galerkin method, the weight function is the same as the basis function, that is:

[0110] (28)

[0111] Therefore, Equation (27) becomes

[0112] (29)

[0113] The electric field at any point within a tetrahedral element can be expressed as

[0114] (30)

[0115] In one embodiment, by solving the linear equations, the output of the electromagnetic field results in the subdomain includes: Dividing each term of the linear equations by the magnetic permeability simultaneously to obtain a new set of equations:

[0116]

[0117] where is the curl operator, is the shape function, is the solution region, is the electric field strength, represents the permittivity of the anisotropic medium, represents the magnetic permeability of the anisotropic medium, represents the conductivity of the anisotropic medium, is the angular frequency, is the imaginary unit; performing element analysis on each term in the new set of equations and accumulating them to obtain the matrix elements corresponding to the subdomain grid. After each subdomain completes the element analysis and format conversion for the solver in parallel, adding boundary conditions, and calling the direct solver to solve in parallel, each subdomain obtains the corresponding electromagnetic field distribution in the subdomain.

[0118] Specifically, the six basic parameters of conductivity anisotropy are: ① The three principal conductivities along the three anisotropic principal axes ; ② The three Euler rotation angles and , representing the anisotropic strike angle, anisotropic dip angle, and anisotropic dip direction angle respectively. The conductivity tensor is jointly determined by the three principal conductivities and the three Euler rotation angles and can represent any anisotropy.

[0119] (31)

[0120] where each element is respectively:

[0121]

[0122]

[0123]

[0124]

[0125] .

[0126] For anisotropic media, the magnetic susceptibility It is a second-order tensor related to spatial position, has 9 components, and is available in the rectangular coordinate system The matrix is ​​represented as

[0127]

[0128] The magnetic induction intensity vector at this time and magnetic field strength The relationship between

[0129]

[0130] The six basic parameters of magnetic permeability anisotropy are: ① The three main axis magnetic susceptibilities along the three main anisotropy axes: ; ②Three Euler rotation angles and , represent the anisotropic strike angle, anisotropic dip angle, and anisotropic inclination angle, respectively. The expressions of the nine components in are consistent with the calculation method of the conductivity tensor, which are all determined by the principal axis magnetic susceptibility and the Euler angle.

[0131] For anisotropic media, the polarizability It is a second-order tensor related to spatial position, has 9 components, and is available in the rectangular coordinate system The matrix is ​​represented as

[0132]

[0133] The electric displacement vector and electric field strength The relationship between

[0134]

[0135] The six basic parameters of dielectric anisotropy are: ① The three main axis polarizabilities along the three main anisotropy axes: ; ②Three Euler rotation angles and , represent the anisotropic strike angle, anisotropic dip angle, and anisotropic inclination angle, respectively. Polarizability The expressions of the nine components in are consistent with the calculation method of the conductivity tensor, which are all determined by the principal axis polarizability and the Euler angle.

[0136] Perform unit analysis on equation (29) and divide each term in the equation by the magnetic permeability

[0137] (32)

[0138] For the convenience of derivation, let

[0139] (33)

[0140] Substitute equation (33) into equation (32), and perform element analysis on the first term therein. The element matrix thus calculated is denoted as

[0141]

[0142]

[0143]

[0144]

[0145] After simplification, we get:

[0146] Due to the symmetry of the matrix elements, here we only give the expressions of the upper diagonal matrix elements

[0147]

[0148]

[0149]

[0150]

[0151]

[0152]

[0153] Thus, the element analysis of the first term in Equation (32) is completed. Next, the element analysis of the second term in Equation (32) is carried out, and the element matrix calculated therefrom is denoted as .

[0154]

[0155] According to the formula (Zienkiewicz and Taylor, 1989)

[0156]

[0157]

[0158] Let

[0159]

[0160] Due to the symmetry of the dielectric constant, it can be seen from the above formula that , so there is:

[0161]

[0162] The other terms are similar and will not be expanded for proof. Since the matrix is symmetric, the values of the elements in the upper triangular position of the tetrahedral element are given here.

[0163]

[0164]

[0165]

[0166]

[0167]

[0168]

[0169]

[0170]

[0171]

[0172]

[0173]

[0174]

[0175]

[0176]

[0177]

[0178]

[0179]

[0180]

[0181]

[0182]

[0183]

[0184] The element analysis of the second term in Equation (32) is thus completed. Next, the element analysis of the third term in Equation (32) is carried out, and the element matrix calculated therefrom is denoted as .

[0185]

[0186]

[0187]

[0188] According to the formula (Zienkiewicz and Taylor, 1989)

[0189]

[0190]

[0191] Let

[0192]

[0193] Due to the symmetry of the conductivity, it can be seen from the above formula that , so there is:

[0194]

[0195] The other terms are similar and will not be expanded for proof. Since the matrix is symmetric, the values of the elements in the upper triangular position of the tetrahedral element are given here.

[0196]

[0197]

[0198]

[0199]

[0200]

[0201]

[0202]

[0203]

[0204]

[0205]

[0206]

[0207]

[0208]

[0209]

[0210]

[0211]

[0212]

[0213]

[0214]

[0215]

[0216]

[0217] The above is based on the anisotropic algorithm of the full-wave electromagnetic vector finite element for tetrahedral unstructured grids, and the element calculation expressions for the analysis of each tetrahedral element are given above.

[0218] It should be understood that although Figure 1The steps in the flowchart are shown in sequence according to the arrows, but these steps are not necessarily executed in the order indicated by the arrows. Unless otherwise specified in this document, there is no strict order restriction for the execution of these steps, and these steps can be executed in other orders. Moreover, Figure 1 At least some of the steps 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 alternately or in turn with at least some of the sub-steps or stages of other steps or other steps.

[0219] The technical features of the above embodiments can be combined arbitrarily. For the sake of brevity of description, 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, it should be considered as the scope recorded in this specification.

[0220] The above-described embodiments merely represent several implementation manners of the present application. The description is relatively specific and detailed, but it should not be construed as a limitation on the scope of the invention. It should be noted that for those of ordinary skill in the art, without departing from the concept of the present application, several modifications and improvements can still be made, and these all belong to the protection scope of the present application. Therefore, the protection scope of the present application should be subject to the appended claims.

Claims

1. A large-scale parallel forward modeling method for solving electromagnetic anisotropy of a base based on PDE, characterized in that: The method comprises: An anisotropic electromagnetic model is constructed according to a constitutive relationship between an electric field and a magnetic field facing an anisotropic medium; the constitutive parameters of the constitutive relationship include the dielectric constant, magnetic permeability and electrical conductivity of the anisotropic medium; Inputting configuration file information of the electromagnetic problem to be solved into the PDE solving base, the PDE solving base loads the electromagnetic model plug-in according to the configuration file information, constructs an electromagnetic field solving module, divides the calculation area into multiple subdomains, and assigns a process from a process set to each subdomain. In the process corresponding to each subdomain, according to the pre-set multi-level and multi-granularity parallel algorithm and the configuration information in the electromagnetic field solving module, the anisotropic electromagnetic model corresponding to the subdomain is solved in parallel, and the electromagnetic field result of the subdomain is output; According to the electromagnetic field results of each subdomain, the forward modeling results are obtained; The anisotropic electromagnetic model corresponding to the subdomain is solved in parallel according to the pre-set multi-level and multi-granularity parallel algorithm and the configuration information in the electromagnetic field solution module. The output subdomain electromagnetic field results include: Each subdomain performs data pre-processing in parallel according to the subdomain grid data and subdomain matrix data stored in the corresponding process; Each subdomain performs unit analysis on the anisotropic electromagnetic model corresponding to the subdomain in parallel according to the configuration information in the electromagnetic field solution module and the pre-processed data, and performs unit analysis on the processed subdomain grid in parallel to obtain a linear equation system; Solving the linear equations in parallel in each subdomain to obtain the subdomain electromagnetic field distribution; Each subdomain performs post-processing on the subdomain electromagnetic field distribution in parallel to obtain a subdomain electromagnetic field result.

2. The method according to claim 1, characterized in that: According to the constitutive relationship between electric and magnetic fields facing anisotropic media, the anisotropic electromagnetic model is constructed including: According to the constitutive relationship between the electric field and the magnetic field facing the anisotropic medium, the definite form of the Maxwell equations is obtained. The definite form of the Maxwell equations is converted into the governing differential equations of the electric field, and the anisotropic electromagnetic model is obtained as follows: in, is the electric field strength, represents the dielectric constant of anisotropic medium, represents the magnetic permeability of anisotropic medium, represents the conductivity of anisotropic media, is the curl operator, is the angular frequency, Is an imaginary unit.

3. The method according to claim 1, characterized in that The anisotropic electromagnetic model corresponding to the subdomain is analyzed in parallel, and the linear equations obtained by performing unit analysis on the processed subdomain grid are as follows: In a single subdomain, multiple processing units are assigned to different processes. Each process processes the corresponding processing units in parallel. Different processes communicate with each other through a message passing interface. In each process, different iterations in the unit loop are assigned to different threads, and multiple threads are executed in parallel to realize parallel unit analysis of the anisotropic electromagnetic model corresponding to the subdomain, and parallel unit analysis of the processed subdomain grid to obtain a linear equation system.

4. The method according to claim 1, characterized in that: The anisotropic electromagnetic model corresponding to the subdomain is analyzed in parallel, and the linear equations obtained by performing unit analysis on the processed subdomain grid are as follows: The space is discretized using tetrahedral meshes, and the integral of the calculation area is decomposed into the accumulation of discrete unit integrals and processed using the weighted residual method. In parallel, the subdomain meshes obtained after processing are subjected to unit analysis to obtain a set of linear equations.

5. The method according to claim 1, characterized in that: Each subdomain solves the linear equations in parallel, and obtains the subdomain electromagnetic field distribution including: The processes corresponding to each subdomain are divided into multiple process groups, and the frequencies are allocated to each process group. All processes in the process group multiply the stored matrix elements by the corresponding frequencies to obtain a frequency-related sparse matrix; Each process group solves the linear equations consisting of frequency-dependent sparse matrices in parallel through a direct solver. In each process group, the main process collects the solutions of other processes in the process group to obtain the subdomain electromagnetic field distribution.

6. The method according to claim 1, characterized in that Solve the linear equations and output the subdomain electromagnetic field results including: Dividing each term of the linear equations by the magnetic permeability, we get a new set of equations: in, is the curl operator, is the shape function, To solve the area, is the electric field strength, represents the dielectric constant of anisotropic medium, represents the magnetic permeability of anisotropic medium, represents the conductivity of anisotropic media, is the angular frequency, is an imaginary unit; Each item in the new set of equations is analyzed and accumulated separately to obtain the matrix elements corresponding to the subdomain grid. After each subdomain completes the unit analysis in parallel and converts the format of the solver, boundary conditions are added and the direct solver is called for parallel solution. Each subdomain obtains the corresponding subdomain electromagnetic field distribution.

7. The method according to claim 5, characterized in that Each subdomain performs post-processing on the subdomain electromagnetic field distribution in parallel, and obtains the subdomain electromagnetic field results including: The impedance, apparent resistivity and phase response of each subdomain are calculated according to the subdomain electromagnetic field distribution to obtain the subdomain electromagnetic field results.

8. The method according to claim 1, characterized in that: Loading the electromagnetic model plug-in according to the configuration file information and constructing the electromagnetic field solution module includes: The electromagnetic model plug-in is loaded according to the input configuration file information, and the plug-in is integrated with the core module through the standardized interface to build the electromagnetic field solution module.

9. The method according to claim 1, characterized in that: The configuration file information includes electromagnetic model subclasses and related components; the electromagnetic model subclasses include anisotropic electromagnetic models and boundary conditions; the related components include a grid file processing sub-plug-in, a solution algorithm sub-plug-in and a PDE discretization sub-plug-in.

Images