GPU-accelerated fast forward method for transient electromagnetic method with large loop source

Abstract

The application discloses a GPU acceleration-based large-loop-source transient electromagnetic fast forward method, which comprises the following steps: firstly, parameter setting and pretreatment are performed; then, the large-loop-source is divided into four edges of north, south, east and west, and each edge is discretized into N electric dipoles by using a Gauss-Legendre quadrature method; a four-level parallel computing architecture of loop edge-frequency sampling point-wave number-time channel is constructed to calculate each edge; wherein, the parallel computing architecture of frequency-wave number-time is used in each edge in the GPU to perform parallel calculation on each electric dipole, so that the frequency-time domain response of each electric dipole is obtained; finally, after the calculation of all electric dipoles is completed, the vectorization calculation and summation of the transient responses of all electric dipoles are performed on the GPU to obtain the total field response of the large-loop-source at the current observation point, and the final time domain electromagnetic field value is output. The above method can efficiently and fully exert the large-scale parallel computing advantage of the GPU, so that the calculation efficiency is effectively improved under the premise of ensuring the calculation precision.

Description

Technical Field

[0001] This invention belongs to the field of detection inversion calculation technology, specifically a GPU-accelerated fast forward modeling method for transient electromagnetic sources with large loops. Background Technology

[0002] Transient electromagnetic method (TEM) is an important geophysical exploration method widely used in mineral exploration, hydrogeology, engineering, and environmental surveys. Among these, large-loop source devices are frequently used due to their advantages such as large detection depth and strong signal. The core of its forward modeling is calculating the electromagnetic response of the subsurface medium under transmitted current excitation. Although one-dimensional forward models can theoretically provide accurate solutions, the efficiency of numerical computation is a significant bottleneck in practical applications. Especially in inversion interpretation, which requires massive forward calculations, the computational speed directly determines the overall efficiency of the work.

[0003] Currently, traditional one-dimensional forward modeling of transient electromagnetic fields in large loop sources is typically executed serially on a central processing unit (CPU). The main technical approach is as follows: first, the large loop source is discretized into N electric dipole units, and the Hankel transform is applied to calculate the electromagnetic field response of each dipole in the frequency domain; then, the frequency domain response is transformed to the time domain using sine and cosine numerical filtering algorithms; finally, the discrete electric dipole responses are superimposed using numerical integration methods to obtain the time-domain electromagnetic field response of the large loop source. However, the above-mentioned traditional method has significant technical bottlenecks:

[0004] 1. Low computational efficiency: Forward modeling of large loop sources typically divides the loop into multiple electric dipoles, calculates the contribution of each segment using the Gauss-Legend quadrature method, and then superimposes them. The calculation of the electric dipole response involves Hankel transform and sine / cosine transform, often using 140-point Hankel transform coefficients and 250-point sine and cosine filter coefficients for frequency-to-time domain conversion, resulting in a huge computational load. While traditional CPU serial computing can achieve forward modeling, its computational speed is slow, especially when the number of time channels increases or the model becomes more complex, making it difficult to meet the needs of large-scale data simulation in actual production or the rapid forward modeling requirements in actual exploration.

[0005] 2. Low memory utilization efficiency: Repetitive data transfer and storage operations during forward modeling reduce overall performance.

[0006] 3. The contradiction between accuracy and speed: In the existing technology, improving the calculation accuracy requires increasing the number of filter coefficients, but this will significantly reduce the calculation speed; while reducing the number of filter coefficients can improve efficiency, it will lead to a loss of accuracy, especially in the calculation of late time channels where the error is large.

[0007] While existing technologies include methods for synthesizing dipole loop sources (such as electric dipole integration along loops), their computational processes are not fully optimized and do not fully utilize the parallel capabilities of GPUs, resulting in insufficient timeliness for large-scale forward or inverse modeling. Therefore, this invention aims to address the complex mathematical processes of recursive calculations, multiple integrations, and transformations in transient electromagnetic forward modeling by providing a new, fast forward modeling framework that can efficiently construct a specific algorithm architecture that fully leverages the advantages of GPU's massively parallel computing, thereby effectively improving computational efficiency while maintaining accuracy. Summary of the Invention

[0008] To address the problems existing in the prior art, this invention provides a GPU-accelerated fast forward modeling method for large loop source transient electromagnetics. By constructing a multi-level parallel computing architecture that can efficiently leverage the advantages of GPU's massive parallel computing, computational efficiency can be effectively improved while ensuring computational accuracy.

[0009] To achieve the above objectives, the technical solution adopted by this invention is: a GPU-accelerated fast forward modeling method for large loop source transient electromagnetic events, comprising the following steps:

[0010] Step 1, Parameter Setting and Preprocessing: First, determine the number of parameters to be set, and assign values ​​to each parameter according to the actual situation of the current observation point. If a parameter is missing in the actual situation, a default value will be automatically assigned. Then, the Hankel transform coefficients and sine and cosine filter coefficients are preloaded using persistent variable technology to avoid repeated loading and transmission, reduce file I / O overhead, and transfer the above parameters and coefficients to the GPU memory, thereby reducing CPU-GPU communication transmission.

[0011] Step 2, Discretization and Parallel Task Planning of the Large Loop Source: The large loop source is divided into four edges: north, south, east, and west. Each edge is discretized into N electric dipoles (integration points) using the Gauss-Legend quadrature method. A four-level parallel computing architecture is used to calculate the four edges, fully exploring the multi-core computing potential of the GPU.

[0012] Step 3: Calculate the frequency-time domain response of each electric dipole: Use a frequency-wavenumber-time parallel computing architecture within the GPU to perform parallel computing on each electric dipole in Step 2, thereby obtaining the frequency-time domain response of each electric dipole.

[0013] Step 4: Synthesize and output the great loop source response: After all electric dipoles have been calculated in Step 3, the transient responses of all electric dipoles are vectorized and summed (Gaussian integral) on the GPU to obtain the total field response of the great loop source at the current observation point.

[0014] Furthermore, the parameters set in step one include the parameters of the layered geodetic model (resistivity ρ and thickness h of each layer), loop dimensions (Lud is the vertical length, Llr is the horizontal width), the center coordinates of the wireframe (X0, Y0), the coordinates of the measuring point (x, y), the emission current I, and the observation time t.

[0015] Furthermore, the parallel task planning in step two specifically includes:

[0016] A. Parallel computation layer of loop edges: The computation tasks of the four edges of the large loop source are distributed to different processing units of the GPU (streaming multiprocessors), realizing parallel computation of each edge at the macro task level; the computation of each loop edge is completely independent, and efficient load balancing is achieved through parallel looping.

[0017] B. Gaussian Integral Parallel Layer: In the calculation of each loop edge, the Gauss-Legend quadrature method is used to discretize each edge into N electric dipoles. This calculation is executed in parallel at the GPU thread block level, which can make full use of the advantages of the GPU's SIMT architecture.

[0018] C. Hankel Transform Parallel Layer: In the frequency domain response calculation of a single electric dipole, the calculation of 47 Hankel transform sampling points is distributed to GPU thread block-level parallel processing; through optimized memory access mode, the efficiency of high-frequency data access is ensured.

[0019] D. Time-channel parallel computing layer: For frequency-to-time domain conversion calculations, the computing tasks of different time channels are parallelized to achieve fine-grained data-level parallel processing; through the parallel processing of time channels, the overall computing time is significantly reduced.

[0020] Furthermore, step three specifically includes:

[0021] I. Frequency Domain and Wavenumber Domain Sampling: Logarithmic uniform sampling is performed in the complex frequency plane and the wavenumber domain respectively to obtain multiple frequency sampling points ωi and wavenumber points λj.

[0022] II. Construction of 3D Parameter Array: The parameters, frequencies, and wavenumbers of the layered geodetic model are expanded into a 3D parameter array (number of strata × number of wavenumber points × number of frequency points) and transferred to the GPU. The purpose of this step is to reconstruct the entire computational task into a highly regularized data parallel problem.

[0023] III. Wave Impedance Recursive Calculation Reconstruction: The traditional layer-by-layer recursive calculation is transformed into matrix multiplication operation, and the dispersion coefficient F(λ,ω) of the surface wave impedance is calculated by combining a three-dimensional parameter array; and the matrix calculation unit of the GPU is used to accelerate the calculation, that is, each thread block or thread of the GPU independently calculates a complete recursive process under a specific (λ,ω) combination, thereby transforming the serial recursion into independent parallel recursive tasks.

[0024] IV. Optimize the Hankel Transform: Using the dispersion coefficients F(λ,ω) calculated in step III, perform the Hankel transform through a linear filtering algorithm and vectorized calculation method to obtain the frequency domain electromagnetic field response H. z (ω);

[0025] V. Improved Sine and Cosine Transform: A non-uniform sampling strategy is employed to maintain high-density sampling in key frequency bands, optimizing the traditional 250-point sine and cosine filter coefficients to 87 points. This is combined with a cubic spline interpolation algorithm to compensate for the accuracy loss caused by the reduced coefficients, thereby improving the frequency domain response H. z (ω) is transformed into a time-domain response h through sine and cosine transforms. z (t). This process involves batch interpolation calculations to reduce GPU-CPU data transfer.

[0026] Furthermore, in step IV, the Hankel transform process uses 47-point Hankel transform coefficients instead of the traditional 140-point Hankel transform coefficients, thereby reducing the amount of computation while ensuring accuracy.

[0027] Furthermore, after obtaining the total field response of the large loop source at the current detection point in step four, the response is transferred from the GPU memory back to the host memory, and then the final time-domain electromagnetic field value is output.

[0028] Compared with the prior art, the present invention has the following advantages:

[0029] 1. Extremely high computational speedup: This invention adopts a four-layer parallel architecture of "loop edge-frequency sampling point-wavenumber-time channel" to match the streaming multiprocessor of the GPU, fully tapping the multi-core computing potential of the GPU and achieving an order-of-magnitude performance improvement.

[0030] 2. Good computational accuracy is maintained: This invention dynamically selects the optimal filter coefficient configuration and compensates for the accuracy loss caused by the reduction of filter coefficients through numerical interpolation methods. This ensures that the accuracy meets engineering requirements while improving computational efficiency, thus guaranteeing the numerical accuracy of forward modeling and making it suitable for complex geological conditions.

[0031] 3. Extremely strong practicality: The efficient computing power of this invention enables the method to be widely applied to practical production scenarios such as rapid processing of transient electromagnetic (TEM) data, inversion of large-scale data, and high-density sampling simulation, which greatly improves work efficiency.

[0032] 4. Significant cost-effectiveness: The hardware of this invention only requires existing GPU hardware to achieve high-performance computing without additional customization, thus reducing the investment cost of high-performance computing equipment. Attached Figure Description

[0033] Figure 1This is an overall flowchart of an embodiment of the present invention.

[0034] Figure 2 This is a schematic diagram of the large loop source being discretized into an electric dipole in an embodiment of the present invention.

[0035] Figure 3 This is a schematic diagram of the frequency-wavenumber-time parallel computing architecture used by each electric dipole in this embodiment of the invention.

[0036] Figure 4 This is a schematic diagram of GPU thread grid planning in an embodiment of the present invention.

[0037] Figure 5 This is a comparison chart of the cosine transform filtering coefficients used in the embodiments of the present invention and the traditional method.

[0038] Figure 6 This is a comparison and relative error diagram of the numerical and analytical solutions of the uniform half-space obtained by the embodiments of the present invention and the traditional methods.

[0039] Figure 7 This is a comparison chart of calculation results under different layered geodetic model parameters in the embodiments of the present invention. Detailed Implementation

[0040] The present invention will be further described below.

[0041] The hardware environment used in this embodiment includes: a workstation or server equipped with an NVIDIA GPU, with a video memory capacity of ≥8GB and a system memory of ≥8GB. The software environment includes: a CentOS 8 operating system (Linux-x86_64), with MATLAB 2025b and CUDA Toolkit 12.2 or later installed, and a MATLAB GPU computing environment configured.

[0042] like Figures 1 to 3 As shown, this embodiment includes the following steps:

[0043] Step 1, Parameter Setting and Preprocessing: First, determine the number of parameters to be set, and assign values ​​to each parameter according to the actual situation of the current observation point, including the parameters of the layered geodetic model (resistivity ρ and thickness h of each layer), loop size (Lud is the vertical length, Llr is the horizontal width), the center coordinates of the wireframe (X0, Y0), the coordinates of the measuring point (x, y), the emission current, and the observation time channel.

[0044] In this embodiment, the parameters of the layered earth model are specifically: uniform half-space geoelectric structure ( ρ 1 = 100 Ω·m);

[0045] Typical three-layer geoelectric structure ( ρ 1 = 100 Ω·m,h 1 = 100 m; ρ 2 = 10 Ω·m, h 2 = 200 m; ρ 3 = 200Ω·m);

[0046] Complex five-layer geoelectric structure ( ρ 1 = 100 Ω·m, h 1 = 100 m; ρ 2 = 10 Ω·m, h 2 = 100 m; ρ 3 = 100 Ω·m, h 3 = 200 m; ρ 4 = 10 Ω·m, h 4 = 300 m; ρ 5 = 200 Ω·m);

[0047] Fine nine-layer geoelectric structure ( ρ 1 = 100 Ω·m, h 1 = 100 m; ρ 2 = 10 Ω·m, h 2 = 100 m; ρ 3 = 100 Ω·m, h 3 = 50 m; ρ 4 = 10 Ω·m, h 4 = 50 m; ρ 5 = 200 Ω·m, h 5 = 50 m; ρ 6 = 500 Ω·m h 6 = 100 m; ρ 7 = 300 Ω·m, h 7 = 100 m; ρ 8 = 100 Ω·m h 8 = 200 m; ρ 9 = 200 Ω·m).

[0048] The loop size is 500m×500m, and the transmission current is 20A; the observation point is located at the center of the loop, that is, the coordinates of the center of the loop and the coordinates of the measurement point are the same; the observation time channel is logspace(log10(5*10^-5,-1,61), that is, 61 log-interval time channels from 5*10^-5 seconds to 0.1 seconds.

[0049] If a parameter is missing in practice, a default value is automatically assigned. Then, persistent variable technology is used to preload the Hankel transform coefficients and sine and cosine filter coefficients to avoid repeated loading and transmission, reduce file I / O overhead, and transfer the above parameters and coefficients to the GPU memory, thereby reducing CPU-GPU communication transmission. At the same time, GPU memory space is pre-allocated, memory layout is optimized, and the gpuArray function is used to transfer the calculated data to the GPU memory.

[0050] Step 2: Discretization and Parallel Task Planning of the Large Loop Source: The large loop source is divided into four edges: north, south, east, and west. Each edge is discretized into N electric dipoles (integration points) using the Gauss-Legend quadrature method. A four-level parallel computing architecture is used to compute the four edges, fully utilizing the multi-core computing potential of the GPU. Specifically:

[0051] A. Parallel computation layer of loop edges: The computation tasks of the four edges (north, south, east, and west) of the large loop source are distributed to different processing units (streaming multiprocessors) of the GPU to achieve parallel computation of each edge at the macro task level; the computation of each loop edge is completely independent, and efficient load balancing is achieved through parallel looping.

[0052] B. Gaussian Integral Parallel Layer: In the calculation of each loop edge, the Gauss-Legend quadrature method is used to discretize each edge into N electric dipoles. This calculation is executed in parallel at the GPU thread block level, which can make full use of the advantages of the GPU's SIMT architecture.

[0053] C. Hankel Transform Parallel Layer: In the frequency domain response calculation of a single electric dipole, the calculation of 47 Hankel transform sampling points is distributed to GPU thread block-level parallel processing; through optimized memory access mode, the efficiency of high-frequency data access is ensured. Specifically, the frequency domain response of a single electric dipole is expressed using the following frequency domain magnetic field expression of a horizontal electric dipole on a layered Earth surface:

[0054] (1)

[0055] in, This represents the current in an electric dipole. The length of the electric dipole. The reflection coefficient, , It is a first-order Bessel function.

[0056] D. Time-channel parallel computing layer: For frequency-to-time domain conversion calculations, the computing tasks of different time channels are parallelized to achieve fine-grained data-level parallel processing; through the parallel processing of time channels, the overall computing time is significantly reduced.

[0057] Step 3: Calculate the frequency-time domain response of each electric dipole: Using a frequency-wavenumber-time parallel computing architecture within the GPU, perform parallel computation on each electric dipole from Step 2 to obtain the frequency-time domain response of each electric dipole, specifically:

[0058] I. Frequency Domain and Wavenumber Domain Sampling: Logarithmic uniform sampling is performed in the complex frequency plane and the wavenumber domain respectively to obtain 87 frequency sampling points ωi and 47 wavenumber points λj.

[0059] II. Construction of 3D Parameter Array: The parameters, frequencies, and wavenumbers of the layered geodetic model are expanded into a 3D parameter array (number of strata × number of wavenumber points × number of frequency points) and transferred to the GPU. The purpose of this step is to reconstruct the entire computational task into a highly regularized data parallel problem.

[0060] III. Wave Impedance Recursive Calculation Reconstruction: The traditional layer-by-layer recursive calculation is transformed into matrix multiplication operations, and the dispersion coefficient F(λ,ω) of the surface wave impedance is calculated by combining a three-dimensional parameter array. The matrix calculation unit of the GPU is used to accelerate the calculation, that is, each thread block or thread of the GPU independently calculates a complete recursive process under a specific (λ,ω) combination, thereby transforming the serial recursion into independent parallel recursive tasks. The bsxfun function is used to replace the repmat function for matrix operations, reducing memory usage.

[0061] IV. Optimization of Hankel Transform: Using the dispersion coefficients F(λ,ω) calculated in step III, a Hankel transform is performed through a linear filtering algorithm and vectorized calculation method. In this process, 47-point Hankel transform coefficients are used instead of the traditional 140-point Hankel transform coefficients, thereby reducing the computational load while maintaining accuracy. The final frequency domain electromagnetic field response H is obtained. z (ω);

[0062] V. Improved Sine and Cosine Transform: A non-uniform sampling strategy is employed to maintain high-density sampling in key frequency bands, optimizing the traditional 250-point sine and cosine filter coefficients to 87 points to reduce GPU-CPU data transfer. Furthermore, a cubic spline interpolation algorithm is used to compensate for the accuracy loss caused by the reduced coefficients, thereby improving the frequency domain response H. z (ω) is transformed into a time-domain response h through sine and cosine transforms. z (t); where the sine and cosine transforms will convert the frequency domain response H. z The specific formulas for transforming (ω) to the time domain are as follows:

[0063]

[0064]

[0065] This can be achieved using digital filtering, transforming the integral into a summation:

[0066] (3)

[0067] (4)

[0068] in and These are the sine and cosine filter coefficients. Where n is the sampling interval, n is the frequency sampling point index, and k is the time channel; in the above formula That is, the time domain response h z (t).

[0069] In the above process, the calculation of wavenumber, frequency, and time channel within each task is performed using a three-layer parallel architecture.

[0070] Step 4: Synthesize and output the big loop source response: After all electric dipoles have been calculated in Step 3, the transient responses of all electric dipoles are vectorized and summed (Gaussian integral) on the GPU to obtain the total field response of the big loop source at the current observation point. Finally, the response is transferred from the GPU memory back to the host memory, and the final time-domain electromagnetic field value is output.

[0071] Effect verification:

[0072] To verify the effectiveness of this embodiment, during the forward modeling of transient electromagnetic data at the current observation point, the existing traditional CPU serial method was simultaneously used to process transient electromagnetic data at the same observation point. The specific results are shown in Table 1:

[0073] Table 1: Comparison of computational efficiency between the GPU-accelerated method of this invention and the traditional CPU serial method under the same conditions:

[0074]

[0075] Summary Table 1 and Figure 5 , 6 Comparative analysis reveals that the GPU-accelerated forward modeling method used in this invention exhibits significant advantages over the traditional CPU serial algorithm in both computational efficiency and accuracy: in a uniform half-space geoelectric model, this method can complete a forward simulation in only about 0.025 seconds, achieving a speedup of 156.9 times compared to the 3.923 seconds of the traditional CPU serial algorithm; and as... Figure 6 As shown, the numerical solution and analytical solution are in high agreement, with the relative errors of the early and late time traces controlled within 4% and 0.5% respectively, fully meeting the requirements for high-precision exploration. Furthermore, Figure 7The comparison results of the four geoelectric models reveal a key pattern: as the complexity of the stratigraphic model increases (from a uniform half-space to a multi-layer geoelectric structure), the parallel computing advantage of this invention becomes increasingly prominent: the speedup ratio for the nine-layer model jumps to 1213.9 times, the GPU computation time only increases to 0.067 seconds, while the traditional CPU time soars to 81.336 seconds. This shows that the method has excellent scalability when dealing with complex geological structures, truly achieving the breakthrough that "the more complex the model, the more significant the efficiency improvement", providing reliable technical support for large-scale high-precision electromagnetic exploration.

[0076] The above description is only a preferred embodiment of the present invention. It should be noted that for those skilled in the art, several improvements and modifications can be made without departing from the principle of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A GPU-accelerated fast forward modeling method for large loop source transient electromagnetic events, characterized in that, Includes the following steps: Step 1, Parameter Setting and Preprocessing: First, determine the number of parameters to be set, and assign values ​​to each parameter according to the actual situation of the current observation point. If a parameter is missing in the actual situation, a default value will be automatically assigned. Then, the Hankel transform coefficients and sine and cosine filter coefficients are preloaded using persistent variable technology, and the above parameters and coefficients are transferred to the GPU memory. Step 2: Discretization and Parallel Task Planning of the Large Loop Source: The large loop source is divided into four edges: north, south, east, and west. Each edge is discretized into N electric dipoles using the Gauss-Legendary quadrature method. A four-level parallel computing architecture is then used to compute the results for each of the four edges, specifically as follows: A. Parallel computation layer of loop edges: Distribute the computation tasks of the four edges of the large loop source to different processing units of the GPU to realize parallel computation of each edge at the macro task level; B. Gaussian Integral Point Parallel Layer: In the calculation of each loop edge, the Gauss-Legend quadrature method is used to discretize each edge into N electric dipoles. This calculation is executed in parallel at the GPU thread block level. C. Hankel Transform Parallel Layer: In the frequency domain response calculation of a single electric dipole, the calculation of 47 Hankel transform sampling points is distributed to GPU thread block-level parallel processing; D. Time-channel parallel computing layer: For frequency-to-time domain conversion calculations, the computing tasks of different time channels are parallelized to achieve fine-grained data-level parallel processing; Step 3: Calculate the frequency-time domain response of each electric dipole: Using a frequency-wavenumber-time parallel computing architecture within the GPU, perform parallel computation on each electric dipole from Step 2 to obtain the frequency-time domain response of each electric dipole, specifically: I. Frequency Domain and Wavenumber Domain Sampling: Logarithmic uniform sampling is performed in the complex frequency plane and the wavenumber domain respectively to obtain multiple frequency sampling points ωi and wavenumber points λj; II. Construction of 3D parameter array: The parameters, frequencies, and wavenumbers of the layered geodetic model are expanded into a 3D parameter array and transferred to the GPU; III. Wave Impedance Recursive Calculation Reconstruction: The traditional layer-by-layer recursive calculation is transformed into matrix multiplication operations, and the dispersion coefficient F(λ,ω) of the surface wave impedance is calculated by combining a three-dimensional parameter array; and the matrix calculation unit of the GPU is used to accelerate the calculation, thereby transforming the serial recursion into independent parallel recursive tasks. IV. Optimize the Hankel Transform: Using the dispersion coefficients F(λ,ω) calculated in step III, perform the Hankel transform through a linear filtering algorithm and vectorized calculation method to obtain the frequency domain electromagnetic field response H. z (ω); V. Improved Sine and Cosine Transform: A non-uniform sampling strategy is employed to maintain high-density sampling in key frequency bands, optimizing the traditional 250-point sine and cosine filter coefficients to 87 points. This is combined with a cubic spline interpolation algorithm to compensate for the accuracy loss caused by the reduced coefficients, thereby improving the frequency domain response H. z (ω) is transformed into a time-domain response h through sine and cosine transforms. z (t); Step 4: Synthesize and output the large loop source response: After all electric dipoles have been calculated in Step 3, the transient responses of all electric dipoles are vectorized and summed on the GPU to obtain the total field response of the large loop source at the current observation point.

2. The GPU-accelerated fast forward modeling method for large loop source transient electromagnetic events according to claim 1, characterized in that, The parameters set in step one include the parameters of the layered geodetic model, loop size, wireframe center coordinates, measuring point coordinates, transmission current, and observation time channel.

3. The GPU-accelerated fast forward modeling method for large loop source transient electromagnetic events according to claim 1, characterized in that, The Hankel transform process in step IV uses 47-point Hankel transform coefficients.

4. The GPU-accelerated fast forward modeling method for large loop source transient electromagnetic events according to claim 1, characterized in that, After obtaining the total field response of the large loop source at the current detection point in step four, the response is transferred from the GPU memory back to the host memory, and then the final time-domain electromagnetic field value is output.

Images