High-performance parallel solver for 3D electromagnetic integral equations Mikhail Kruglyakov
[email protected] Lomonsov Moscow State University
Electromagnetic sounding
COMMEMI3d3: High conductivity contrast modeling I I I
σb(z)
I
Coefficients calculation (continuation) αβ Calculation of Wl
7 blocks in 3-layered media Layers resistivity: 103, 104 and 10 Ω·m Layers thickness: 1 and 6.5 km The maximum conductivity contrast is more than 3 · 104
F =
a
Top view
σa(M )
6
T
c
b
F˜ (p, α, β, ψ)ep =
1 2
8
I
Large size of modeling domain
hx α= , R
6 4 2 0
High conductivity contrast, more than 104
-2
Y,km
-3
1
0
-1
-2
2
3
4
5
6
7
300 Ω·m I
H — total magnetic field
I
E0 — normal electrical field for the layered media H0 — normal electrical field for the layered media
I
σb(z) — conductivity of layered media
I
T — bounded domain with anomalous conductivity σa(M ) ∆σ (M0) = σa(M0) − σb(M0) ˆE, G ˆ H — electrical and magnetic Green’s tensors G
I
I I I
I I
Tn — nonoverlapping rectangular cells
X,km
-2
-1
-2
0
1
30 Ω·m
100 Ω·m
ρxy at y = 3830m
103
3
0.1 Ω·m
5
4
6
7
a−b ψ = arctan , c−d
hy β= , R
R = ep ,
λ = e−q .
0.3 Ω·m
Matrix storage
Profiles
ρyx at y = 3830m
103
102
Block Toeplitz matrix n,k y y n,k x x z z b =B b G In − Ik , In − Ik , In, Ik E
(9)
b n,k induced by 4NxNy blocks Q(I, J), I = −(Ny − 1) . . . Ny − 1, B J = −(Nx − 1) . . . Nx − 1 Qxx Qxy Qxz Q(I, J) = Qyx Qyy Qyz (10) Qzx Qzy Qzz Here Qαβ are the matrices of the order Nz , α, β = x, y, z
102
Symmetric and antisymmetric properties Using Lorentz reciprocity we obtain
N S
101
Qzx = −QTxz , Qzy = −QTyz , Qxx = QTxx, Qyy = QTyy , Qzz = QTzz ,
101
(11)
100
Integral Equations
Qxy = QTxy = Qyx = QTyx,
For any point M ∈ R3 Z bE (M, M0)∆σ (M0)E(M0)dTM + E0(M ) E(M ) = G 0 T
H(M ) =
300 Ω·m
30 Ω·m
2
X,km
The solid lines are the results of modeling with different cell diameters d by presented IE solver. The green crosses are the result of modeling by the 2nd order finite elements solver from [2].
Tn = T n=1 h i y y Tn = [hxInx, hx (Inx + 1)] × hy In , hy In + 1 × zInz , zInz +1 N = NxNy Nz , Ikτ = 1 . . . Nτ , τ = x, y, z
Z
(8)
f˜(q) = f (λ),
αβ
Notations E — total electrical field
p−q Q(p − q, α, β, ψ)e f˜(q)dq,
Wl are the coefficients of digital filter approximating convolution integral (8). To compute them we need only one pair of analytically αβ computed functions F and f satisfied (7). Then Wl computed as ratio between discrete Fourier spectrums of F˜ and f˜.
0
-3
I
F˜ (p, α, β, ψ) = F (a, b, c, d, hx, hy ), 2
Complex conductivity distribution
I
Z∞
100 −3000 −2000 −1000
(1)
10
0
1000
2000
m
3000
4000
5000
6000
7000
10−1 −3000 −2000 −1000
ρxy at x = 1900m
3
10
0
1000
2000
m
3000
4000
5000
6000
7000
ρyx at x = 1900m
3
bH (M, M0)∆σ (M0)E(M0)dTM + H0(M ) G 0
Parallelization
T
If M ∈ T we obtain 3D singular vector equation Z bE (M, M0)∆σ (M0)E(M0)dTM = E0(M ) E(M ) − G 0
where T indicates a matrix transpose. I Only Qxz , Qyz and upper diagonal parts of Qxx , Qxy , Qyy , Qzz are stored. I We need 8 · Nx Ny Nz · (2Nz + 1) · 16 bytes to store matricies, which is practically a half to [3].
10
Distributed storage Node 1 Node 2 Q(0, 0) Q(1, 0) Q(0, 1) Q(1, 1) .. .. Q(0, 2Nx − 1) Q(1, 2Nx − 1)
2
(2)
T
102
Projective method
Efficiency of parallelization Matrix calculation
101
Idea Approximate integral averages of E over Tn instead of E in the centres of Tn as in collocation method
10
Main stages I Construct the system of linear equations I Solve this system by FGRMES approach I Compute EM field in observation points by (1)
4
10
... Node 2Ny ... Q(2Ny − 1, 0) ... Q(2Ny − 1, 1) .. .. . . . Q(2Ny − 1, 2Nx − 1) Equation solving
3
103 101
100 −2000
0
2000
8000
0
−2000
2000
Φxy at y = 3830m
84
Challenges I Fast and accurate coefficients calculations I Storage of dense matrix
m
6000
4000
m
4000
6000
102
8000
Φyx at y = 3830m
Time, s
I
0
4
3
Main challenges of 3D electromagnetic modeling in the modern geophysics:
d
(7)
−∞ Y,km
Z,km
0
Z
J0(λρ)f (λ)dλ dxdx0dydy0
q Let R = (c − d)2 + (a − b))2 and transform RHS of (7) to convolution integral:
8
3D view
a+h Z y b+h Z y c+h Z x d+h Z x Z∞
102
−96
101 −98 101
82 −100
Coefficients calculation
−102
80
To approximate integral averages we need double volumetric integrals Z Z bE (M, M0)dTM dTM bn,k = G G E 0 (3)
78
76
(4)
0
q ρ = (x − x0)2 + (y − y0)2 ν = 1, σ
82
Uν (z, z0, λ) = Uν0(z, z0, λ) is a fundamental function of layered media from [1].
0
1000
2000
m
3000
4000
5000
6000
=
1000
2000
m
3000
4000
5000
6000
7000
Nλ X
γ,ν
103
104
Number of processes Lomonosov Nx=128 Ny =128 Nz =50
Linear speed up
Gnu Integral Equation Modeling in ElectroMagnetic Geophysics I
Φyx at x = 1900m −99
I
81.5 I
80.5
−101 I
(5)
−102
requires only half of an amount of RAM compared with the other existing 3D IE solvers with matrix storage has high parallelization degree to perform computations for large scale models with a complex conductivity distribution demonstrate the good fit with the high-order finite elements method, including the solving of the high conductivity contrast problems is an open source software distributed under the GPLv2 license and it is available at https://github.com/DarthLaran/GIEM2G
79.5
α,β
Wl
y
y
Contacts
γ,ν
(Inx − Ikx, In − Ik )κl (Inz , Ikz ) γ,ν
κl (Inz , Ikz ) = fl
79
Ikz
(Inz ) ×
Y
−103
[email protected]
78.5
plm(ν) × glν (Ikz )
m=Inz
(6)
−104 78 77.5
I
0
Φxy at x = 1900m
l=1
I
−116 −3000 −2000 −1000
80
Tn Tk
I
7000
81
∂ α+β+γ ν I (M, M0)dTM dTM0 = α β γ ∂x ∂y ∂z
102
The presented solver
−100
Z Z
Bluegene/P Nx=128 Ny =1024 Nz =50
101
Conclusion
−114 72 −3000 −2000 −1000
104
−112
Z∞ J0(λρ)Uν (z, z0, λ)λdλ
103
100 100
Time of matricies calculation and IE solving for different models at different HPC from MSU depends on number of processes
−110
74
I ν (x − x0, y − y0, z, z0) =
102
Bluegene/P Nx=1024 Ny =1024 Nz =10
−106 −108
ˆ E (M, M0), G ˆ H (M, M0) are expressed as followThe components of G ing integral and its partial derivatives [1]:
101
Number of processes
−104
Tk Tn
I
100 100
α,β −2000 0 2000 Wl (I, J) depends on hx, hy , I, J, α, β, l only m d = 86.6m γ,ν ν l The factors fl (s), gl (s), pm(ν), s, m = 1 . . . Nz depends on References z0, . . . zNz , normal section and frequency
This factors can be stably computed analytically
4000
6000
8000
d = 43.3m
−105 −2000
0
d = 28.5m
2000
m
4000
6000
8000
2nd order FE
[1] V. Dmitriev, A. Silkin, and R. Farzan, “Tensor green function for the system of maxwell´s equations in a layered medium,” Computational Mathematics and Modeling, vol. 13, no. 2, pp. 107–118, 2002.
γ,ν z z A. Grayver and T. Kolev, “Large-scale 3d geo-electromagnetic modeling using parallel adaptive high-order finite element method,” Geophysics, To compute all factors for all κl (In, Ik ) for all Inz , Ikz = 1 . . . Nz [2] in review. [3] D. B. Avdeev, A. V. Kuvshinov, O. V. Pankratov, and G. A. Newman, “High-performance three-dimensional electromagnetic modelling using we need the order of O(Nz ) operations for fixed l, γ, ν modified neumann series. wide-band numerical solution and examples,” J. Geomagn. Geoelectr., vol. 49, pp. 1519–1539, 1997.
Feel free to take my card and email me
