傳統大地電磁法(Magnetotellurics)，是使用中央處理器(CPU)演算為主的有限差分法(Finite Different)，其將馬克士威偏微分方程式離散化(Discretize)產生的線性系統(Linear System)，再使用迭代解法(Iterative Methods)估計準靜態電磁場(Quasi-Static Electromagnetic Field)的分布，並利用地層結構逆推演算法(Stratum Structure Inversion Algorithm)找出估計的電性地層(Electric Stratum)電導結構。但是，其迭代解法(Iterative Methods)卻需要花費極高的運算時間成本，且會有相依性，無法切割演算以達到平行加速。為了解決運算時間的問題，本研究提出適用於圖形處理器(GPU)為主的有限元素法(Finite Element)，在不依賴無法平行化迭代解法(Iterative Methods)條件下，估計出準靜態電磁場(Quasi-Static Electromagnetic Field)的分布，同時切割演算以達到平行加速的功能。

本研究主要分為兩個區塊的加速：模擬電磁波擴散分佈演算法(Simulation of Electromagnetic Diffusion Distribution Algorithm)和地層結構逆推演算法(Stratum Structure Inversion Algorithm)。模擬電磁波擴散分佈演算法(Simulation of Electromagnetic Diffusion Distribution Algorithm)主要的工作是將所有的電性地層網格視為獨立的電磁場源，再利用此特性估計地面觀測站所偵測到的電磁場值。地層結構逆推演算法(Stratum Structure Inversion Algorithm)則是基於模擬電磁波擴散分佈演算法(Simulation of Electromagnetic Diffusion Distribution Algorithm)求出網格的雅可比矩陣(Jacobian Matrix)，使得本研究可以使用牛頓法找出下一個的遞推搜尋方向。在這兩個演算法中因為每個電磁場源都是獨立的所以適用於圖形處理器(GPU)平行運算。實驗結果如下，一、在模擬電磁波擴散分佈演算法(Simulation of Electromagnetic Diffusion Distribution Algorithm)中，圖形處理器(GPU)相較於中央處理器(CPU)可加速達50~71倍。二、簡化雅可比矩陣(Jacobian Matrix)的地層結構逆推演算法(Stratum Structure Inversion Algorithm)在圖形處理器(GPU)相較於傳統的地層結構逆推演算法(Stratum Structure Inversion Algorithm)在中央處理器(CPU)上的運算，於一次遞迴條件下可加速達11,000~55,000倍。三、簡化雅可比矩陣(Jacobian Matrix)的地層結構逆推演算法(Stratum Structure Inversion Algorithm)，圖形處理器(GPU)相較於中央處理器(CPU)在估算整個電性地層條件時，則加速可達到68~82倍。

The 3D Magnetotelluric Method is a CPU-based finite difference method which uses Maxwell partial differential equations discretization to produce a linear system used in the iterative solution of quasi-static electromagnetic field distribution estimation and a Stratum Structure Inversion Algorithm to estimate the stratum's electric conductivity structure. But one of the drawbacks of the 3D Magnetotelluric Method is the solving process of linear systems iterative methods has a very high computation time cost. In order to solve the problem of computation time, this study proposes a graphics processor (GPU) based finite element method that estimates the distribution of the quasi-static electromagnetic field without relying on the iterative solving of linear systems for achieving a faster estimation process.

This study can be separated into two parts: Simulation of Electromagnetic Diffusion Distribution Algorithm and Stratum Structure Inversion Algorithm. Simulation of Electromagnetic Diffusion Distribution Algorithm considers each cell in the grid of electromagnetic structure as an independent source of an electromagnetic field and uses this concept to estimate the electromagnetic field values detected at the observation points on the ground. Stratum Structure Inversion Algorithm calculates the grid’s Jacobian matrix based on the result of Simulation of Electromagnetic Diffusion Distribution Algorithm and find the next recursive search direction with Newton's Method. Because all electromagnetic field sources is independent in both algorithms, so our method is suitable for GPU's parallel computing capability.
The results of our experiment are as follows: 1. For the Simulation of Electromagnetic Diffusion Distribution Algorithm, the GPU version is 50 ~ 71 times faster than the CPU version. 2. For a single iteration in the improved Stratum Structure Inversion Algorithm, the GPU version with simplified Jacobian matrix can speed 11,000 ~ 55,000 times faster than the CPU version with normal Jacobian matrix. 3. The improved Stratum Structure Inversion Algorithm on GPU completes the stratum condition estimation 68 ~ 82 times faster than the CPU version.

[1] M.ˇCuma and M. S. Zhdanov, “Massively parallel regularized 3d inversion ofpotential fields on cpus and gpus,” Computers & Geosciences, vol. 62, pp. 80–87, 2014.
[2]U. S. Inan and A. Inan, Engineering Electromagnetics. Prentice Hall; 1 edition,1998.
[3] F. Simpson and K. Bahr, Practical Magnetotellurics. Cambridge University
Press, 2005.
[4] P. E. Wannamaker, J. A. Stodt, and L. Rijo, “A stable finite element solution for two-dimensional magnetotelluric modelling,” Geophysical Journal International, vol. 88, no. 1, pp. 277–296, 1987.
[5] A. Kelbert, N. Meqbel, G. D. Egbert, and K. Tandon, “Modem: A modular system for inversion of elecgtromagnetic geophysical data,” Computers &
Geosciences, vol. 66, pp. 40–53, 2014.
[6] W. Siripunvaraporn and G. Egbert, “Wsinv3dmt: vertical magnetic field transfer function inversion and parallel implementation,” Physics of the Earth and Planetary Interiors, vol. 173, no. 3, pp. 317–329, 2009.
[7] R. L. Mackie, T. R. Madden, and P. E. Wannamaker, “Three-dimensional magnetotelluric modeling using difference equations-theory and comparisons to integral equation solutions,” Geophysics, vol. 58, no. 2, pp. 215–226, 1993.
[8] V. Puzyrev, J. Koldan, J. de la Puente, G. Houzeaux, M. V´azquez, and
J. M. Cela, “A parallel finite-element method for three-dimensional controlled-source electromagnetic forward modelling,” Geophysical Journal International,p. ggt027, 2013.
[9] P. N. Glaskowsky, “Nvidia’s fermi: the first complete gpu computing architecture,” Nvidia, September, p. 23, 2009.
[10] N. Goodnight, C. Woolley, G. Lewin, D. Luebke, and G. Humphreys, “A multigrid solver for boundary value problems using programmable graphics hardware,” ACM SIGGRAPH 2005 Courses, p. 193, 2005.
[11] Y. Wang, M. Baboulin, K. Rupp, O. Le Maふıtre, and Y. Fraigneau, “Solving 3d incompressible navier-stokes equations on hybrid cpu/gpu systems,” Society for Computer Simulation International, p. 12, 2014.
[12] W. Kahan, “Ieee standard 754 for binary floating-point arithmetic,” Lecture Notes on the Status of IEEE, vol. 754, no. 94720-1776, p. 11, 1996.
[13] J. T. Smith and J. R. Booker, “Magnetotelluric inversion for minimum structure,” Geophysics, vol. 53, no. 12, pp. 1565–1576, 1988.
[14] C. Brewitt-Taylor and J. Weaver, “On the finite difference solution of twodimensional induction problems,” Geophysical Journal International, vol. 47,no. 2, pp. 375–396, 1976.
[15] M. Sadiku, Numerical Techniques in Electromagnetics, Second Edition. Taylor & Francis, 2000.