mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-31 12:40:46 +08:00
Some analytics and set up for forward problem
This commit is contained in:
SEOGI KANG
@@ -0,0 +1,221 @@
|
||||
from scipy.constants import mu_0
|
||||
from SimPEG.Utils.sputils import kron3, speye, sdiag
|
||||
import matplotlib.pyplot as plt
|
||||
from SimPEG import *
|
||||
import numpy as np
|
||||
import scipy.sparse as sp
|
||||
|
||||
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 MagSphereAnalFun(x, y, z, R, x0, y0, z0, mu1, mu2, H0, flag):
|
||||
"""
|
||||
Analytic function for Magnetics problem. The set up here is
|
||||
magnetic sphere in whole-space.
|
||||
|
||||
* (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 == 'total'):
|
||||
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]**2-y[~ind]**2-z[~ind]**2))
|
||||
elif (flag == 'secondary'):
|
||||
Bx[~ind] = mu1*(H0/r[~ind]**5*(R**3)*rf1*(2*x[~ind]**2-y[~ind]**2-z[~ind]**2))
|
||||
|
||||
By[~ind] = mu1*(H0/r[~ind]**5*(R**3)*rf1*(3*x[~ind]*y[~ind]))
|
||||
Bz[~ind] = mu1*(H0/r[~ind]**5*(R**3)*rf1*(3*x[~ind]*z[~ind]))
|
||||
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()
|
||||
|
||||
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]
|
||||
|
||||
|
||||
def MagSphereAnalFunA(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
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
|
||||
hxind = ((0,25,1.3),(21, 12.5),(0,25,1.3))
|
||||
hyind = ((0,25,1.3),(21, 12.5),(0,25,1.3))
|
||||
hzind = ((0,25,1.3),(20, 12.5),(0,25,1.3))
|
||||
hx, hy, hz = Utils.meshTensors(hxind, hyind, hzind)
|
||||
M3 = Mesh.TensorMesh([hx, hy, hz], [-sum(hx)/2,-sum(hy)/2,-sum(hz)/2])
|
||||
indxd, indxu, indyd, indyu, indzd, indzu = M3.faceBoundaryInd
|
||||
mu0 = 4*np.pi*1e-7
|
||||
chibkg = 0.
|
||||
chiblk = 0.01
|
||||
chi = np.ones(M3.nC)*chibkg
|
||||
sph_ind = spheremodel(M3, 0, 0, 0, 100)
|
||||
chi[sph_ind] = chiblk
|
||||
mu = (1.+chi)*mu0
|
||||
Bbc = CongruousMagBC(M3, np.array([1., 0., 0.]), chi)
|
||||
|
||||
flag = 'secondary'
|
||||
Box = 1.
|
||||
H0 = Box/mu_0
|
||||
Bbcxx, Bbcxy, Bbcxz = MagSphereAnalFun(M3.gridFx[(indxd|indxu),0], M3.gridFx[(indxd|indxu),1], M3.gridFx[(indxd|indxu),2], 100, 0., 0., 0., mu_0, mu_0*(1+chiblk), H0, flag)
|
||||
Bbcyx, Bbcyy, Bbcyz = MagSphereAnalFun(M3.gridFy[(indyd|indyu),0], M3.gridFy[(indyd|indyu),1], M3.gridFy[(indyd|indyu),2], 100, 0., 0., 0., mu_0, mu_0*(1+chiblk), H0, flag)
|
||||
Bbczx, Bbczy, Bbczz = MagSphereAnalFun(M3.gridFz[(indzd|indzu),0], M3.gridFz[(indzd|indzu),1], M3.gridFz[(indzd|indzu),2], 100, 0., 0., 0., mu_0, mu_0*(1+chiblk), H0, flag)
|
||||
Bbc_anal = np.r_[Bbcxx, Bbcyy, Bbczz]
|
||||
|
||||
fig, ax = plt.subplots(1,1, figsize = (10, 10))
|
||||
ax.plot(Bbc_anal)
|
||||
ax.plot(Bbc)
|
||||
plt.show()
|
||||
err = np.linalg.norm(Bbc-Bbc_anal)/np.linalg.norm(Bbc_anal)
|
||||
|
||||
if err < 0.1:
|
||||
print 'Mag Boundary computation is valid, err = ', err
|
||||
else:
|
||||
print 'Mag Boundary computation is wrong!!, err = ', err
|
||||
|
||||
pass
|
||||
@@ -0,0 +1,26 @@
|
||||
from SimPEG import *
|
||||
|
||||
|
||||
class Magnetics(object):
|
||||
"""docstring for Magnetics"""
|
||||
def __init__(self, arg):
|
||||
super(Magnetics, self).__init__()
|
||||
self.arg = arg
|
||||
|
||||
def getA(self, m):
|
||||
"""
|
||||
GetA creates and returns the A matrix for the Magnetics problem
|
||||
|
||||
The A matrix has the form:
|
||||
|
||||
.. math::
|
||||
|
||||
\mathbf{A} = \mathbf{D}\mu\mathbf{G}u
|
||||
|
||||
|
||||
|
||||
"""
|
||||
return self.mesh.faceDiv
|
||||
|
||||
|
||||
|
||||
@@ -11,12 +11,12 @@ class MagProblemTests(unittest.TestCase):
|
||||
prob = PF.Mag.MagProblem(M, mod, None)
|
||||
|
||||
self.prob = prob
|
||||
self.M = M
|
||||
self.M = M
|
||||
|
||||
def test_forward(self):
|
||||
|
||||
|
||||
passed = True
|
||||
self.assertTrue(passed)
|
||||
self.assertTrue(passed)
|
||||
|
||||
|
||||
def test_DirchletBC(self):
|
||||
|
||||
Reference in New Issue
Block a user