forked from simpeg/simpeg
-
Notifications
You must be signed in to change notification settings - Fork 1
/
MagAnalytics.py
239 lines (178 loc) · 7.46 KB
/
MagAnalytics.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
from __future__ import print_function
from scipy.constants import mu_0
from SimPEG import *
from SimPEG import Utils
from SimPEG import Mesh
import numpy as np
from SimPEG.Utils import kron3, speye, sdiag
import matplotlib.pyplot as plt
def spheremodel(mesh, x0, y0, z0, r):
"""
Generate model indicies for sphere
- (x0, y0, z0 ): is the center location of sphere
- r: is the radius of the sphere
- it returns logical indicies of cell-center model
"""
ind = np.sqrt( (mesh.gridCC[:, 0]-x0)**2+(mesh.gridCC[:, 1]-y0)**2+(mesh.gridCC[:, 2]-z0)**2 ) < r
return ind
def MagSphereAnaFun(x, y, z, R, x0, y0, z0, mu1, mu2, H0, flag='total'):
"""
test
Analytic function for Magnetics problem. The set up here is
magnetic sphere in whole-space assuming that the inducing field is oriented in the x-direction.
* (x0, y0, z0)
* (x0, y0, z0 ): is the center location of sphere
* r: is the radius of the sphere
.. math::
\mathbf{H}_0 = H_0\hat{x}
"""
if (~np.size(x)==np.size(y)==np.size(z)):
print("Specify same size of x, y, z")
return
dim = x.shape
x = Utils.mkvc(x)
y = Utils.mkvc(y)
z = Utils.mkvc(z)
ind = np.sqrt((x-x0)**2+(y-y0)**2+(z-z0)**2) < R
r = Utils.mkvc(np.sqrt((x-x0)**2+(y-y0)**2+(z-z0)**2))
Bx = np.zeros(x.size)
By = np.zeros(x.size)
Bz = np.zeros(x.size)
# Inside of the sphere
rf2 = 3*mu1/(mu2+2*mu1)
if flag is 'total' and any(ind):
Bx[ind] = mu2*H0*(rf2)
elif (flag == 'secondary'):
Bx[ind] = mu2*H0*(rf2)-mu1*H0
By[ind] = 0.
Bz[ind] = 0.
# Outside of the sphere
rf1 = (mu2-mu1)/(mu2+2*mu1)
if (flag == 'total'):
Bx[~ind] = mu1*(H0+H0/r[~ind]**5*(R**3)*rf1*(2*(x[~ind]-x0)**2-(y[~ind]-y0)**2-(z[~ind]-z0)**2))
elif (flag == 'secondary'):
Bx[~ind] = mu1*(H0/r[~ind]**5*(R**3)*rf1*(2*(x[~ind]-x0)**2-(y[~ind]-y0)**2-(z[~ind]-z0)**2))
By[~ind] = mu1*(H0/r[~ind]**5*(R**3)*rf1*(3*(x[~ind]-x0)*(y[~ind]-y0)))
Bz[~ind] = mu1*(H0/r[~ind]**5*(R**3)*rf1*(3*(x[~ind]-x0)*(z[~ind]-z0)))
return np.reshape(Bx, x.shape, order='F'), np.reshape(By, x.shape, order='F'), np.reshape(Bz, x.shape, order='F')
def CongruousMagBC(mesh, Bo, chi):
"""Computing boundary condition using Congrous sphere method.
This is designed for secondary field formulation.
>> Input
* mesh: Mesh class
* Bo: np.array([Box, Boy, Boz]): Primary magnetic flux
* chi: susceptibility at cell volume
.. math::
\\vec{B}(r) = \\frac{\mu_0}{4\pi} \\frac{m}{ \| \\vec{r} - \\vec{r}_0\|^3}[3\hat{m}\cdot\hat{r}-\hat{m}]
"""
ind = chi > 0.
V = mesh.vol[ind].sum()
gamma = 1/V*(chi*mesh.vol).sum() # like a mass!
Bot = np.sqrt(sum(Bo**2))
mx = Bo[0]/Bot
my = Bo[1]/Bot
mz = Bo[2]/Bot
mom = 1/mu_0*Bot*gamma*V/(1+gamma/3)
xc = sum(chi[ind]*mesh.gridCC[:, 0][ind])/sum(chi[ind])
yc = sum(chi[ind]*mesh.gridCC[:, 1][ind])/sum(chi[ind])
zc = sum(chi[ind]*mesh.gridCC[:, 2][ind])/sum(chi[ind])
indxd, indxu, indyd, indyu, indzd, indzu = mesh.faceBoundaryInd
const = mu_0/(4*np.pi)*mom
rfun = lambda x: np.sqrt((x[:, 0]-xc)**2 + (x[:, 1]-yc)**2 + (x[:, 2]-zc)**2)
mdotrx = (mx*(mesh.gridFx[(indxd|indxu), 0]-xc)/rfun(mesh.gridFx[(indxd|indxu), :]) +
my*(mesh.gridFx[(indxd|indxu), 1]-yc)/rfun(mesh.gridFx[(indxd|indxu), :]) +
mz*(mesh.gridFx[(indxd|indxu), 2]-zc)/rfun(mesh.gridFx[(indxd|indxu), :]))
Bbcx = const/(rfun(mesh.gridFx[(indxd|indxu), :])**3)*(3*mdotrx*(mesh.gridFx[(indxd|indxu), 0]-xc)/rfun(mesh.gridFx[(indxd|indxu), :])-mx)
mdotry = (mx*(mesh.gridFy[(indyd|indyu), 0]-xc)/rfun(mesh.gridFy[(indyd|indyu), :]) +
my*(mesh.gridFy[(indyd|indyu), 1]-yc)/rfun(mesh.gridFy[(indyd|indyu), :]) +
mz*(mesh.gridFy[(indyd|indyu), 2]-zc)/rfun(mesh.gridFy[(indyd|indyu), :]))
Bbcy = const/(rfun(mesh.gridFy[(indyd|indyu), :])**3)*(3*mdotry*(mesh.gridFy[(indyd|indyu), 1]-yc)/rfun(mesh.gridFy[(indyd|indyu), :])-my)
mdotrz = (mx*(mesh.gridFz[(indzd|indzu), 0]-xc)/rfun(mesh.gridFz[(indzd|indzu), :]) +
my*(mesh.gridFz[(indzd|indzu), 1]-yc)/rfun(mesh.gridFz[(indzd|indzu), :]) +
mz*(mesh.gridFz[(indzd|indzu), 2]-zc)/rfun(mesh.gridFz[(indzd|indzu), :]))
Bbcz = const/(rfun(mesh.gridFz[(indzd|indzu), :])**3)*(3*mdotrz*(mesh.gridFz[(indzd|indzu), 2]-zc)/rfun(mesh.gridFz[(indzd|indzu), :])-mz)
return np.r_[Bbcx, Bbcy, Bbcz], (1/gamma-1/(3+gamma))*1/V
def MagSphereAnaFunA(x, y, z, R, xc, yc, zc, chi, Bo, flag):
"""Computing boundary condition using Congrous sphere method.
This is designed for secondary field formulation.
>> Input
mesh: Mesh class
Bo: np.array([Box, Boy, Boz]): Primary magnetic flux
Chi: susceptibility at cell volume
.. math::
\\vec{B}(r) = \\frac{\mu_0}{4\pi}\\frac{m}{\| \\vec{r}-\\vec{r}_0\|^3}[3\hat{m}\cdot\hat{r}-\hat{m}]
"""
if (~np.size(x)==np.size(y)==np.size(z)):
print("Specify same size of x, y, z")
return
dim = x.shape
x = Utils.mkvc(x)
y = Utils.mkvc(y)
z = Utils.mkvc(z)
Bot = np.sqrt(sum(Bo**2))
mx = Bo[0]/Bot
my = Bo[1]/Bot
mz = Bo[2]/Bot
ind = np.sqrt((x-xc)**2+(y-yc)**2+(z-zc)**2 ) < R
Bx = np.zeros(x.size)
By = np.zeros(x.size)
Bz = np.zeros(x.size)
# Inside of the sphere
rf2 = 3/(chi+3)*(1+chi)
if (flag == 'total'):
Bx[ind] = Bo[0]*(rf2)
By[ind] = Bo[1]*(rf2)
Bz[ind] = Bo[2]*(rf2)
elif (flag == 'secondary'):
Bx[ind] = Bo[0]*(rf2)-Bo[0]
By[ind] = Bo[1]*(rf2)-Bo[1]
Bz[ind] = Bo[2]*(rf2)-Bo[2]
r = Utils.mkvc(np.sqrt((x-xc)**2+(y-yc)**2+(z-zc)**2 ))
V = 4*np.pi*R**3/3
mom = Bot/mu_0*chi/(1+chi/3)*V
const = mu_0/(4*np.pi)*mom
mdotr = (mx*(x[~ind]-xc)/r[~ind] + my*(y[~ind]-yc)/r[~ind] + mz*(z[~ind]-zc)/r[~ind])
Bx[~ind] = const/(r[~ind]**3)*(3*mdotr*(x[~ind]-xc)/r[~ind]-mx)
By[~ind] = const/(r[~ind]**3)*(3*mdotr*(y[~ind]-yc)/r[~ind]-my)
Bz[~ind] = const/(r[~ind]**3)*(3*mdotr*(z[~ind]-zc)/r[~ind]-mz)
return Bx, By, Bz
def IDTtoxyz(Inc, Dec, Btot):
"""Convert from Inclination, Declination,
Total intensity of earth field to x, y, z
"""
Bx = Btot*np.cos(Inc/180.*np.pi)*np.sin(Dec/180.*np.pi)
By = Btot*np.cos(Inc/180.*np.pi)*np.cos(Dec/180.*np.pi)
Bz = -Btot*np.sin(Inc/180.*np.pi)
return np.r_[Bx, By, Bz]
def MagSphereFreeSpace(x, y, z, R, xc, yc, zc, chi, Bo):
"""Computing the induced response of magnetic sphere in free-space.
>> Input
x, y, z: Observation locations
R: radius of the sphere
xc, yc, zc: Location of the sphere
chi: Susceptibility of sphere
Bo: Inducing field components [bx, by, bz]*|H0|
"""
if (~np.size(x) == np.size(y) == np.size(z)):
print("Specify same size of x, y, z")
return
x = Utils.mkvc(x)
y = Utils.mkvc(y)
z = Utils.mkvc(z)
nobs = len(x)
Bot = np.sqrt(sum(Bo**2))
mx = np.ones([nobs]) * Bo[0] * R**3 / 3. * chi
my = np.ones([nobs]) * Bo[1] * R**3 / 3. * chi
mz = np.ones([nobs]) * Bo[2] * R**3 / 3. * chi
M = np.c_[mx, my, mz]
rx = (x - xc)
ry = (y - yc)
rz = (zc - z)
rvec = np.c_[rx, ry, rz]
r = np.sqrt((rx)**2+(ry)**2+(rz)**2)
B = -Utils.sdiag(1./r**3)*M + \
Utils.sdiag((3 * np.sum(M*rvec, axis=1))/r**5)*rvec
Bx = B[:, 0]
By = B[:, 1]
Bz = B[:, 2]
return Bx, By, Bz