diff --git a/simpegPF/BaseMag.py b/simpegPF/BaseMag.py index fdd4aea1..23e91a13 100644 --- a/simpegPF/BaseMag.py +++ b/simpegPF/BaseMag.py @@ -9,20 +9,69 @@ class BaseMagData(Data.BaseData): rxType = None #: receiver type def __init__(self, **kwargs): - Data.BaseData.__init__(**kwargs) - + Data.BaseData.__init__(self, **kwargs) #TODO: change to inc, dec, intensity - def setBackgroundField(self, x, y, z): + def setBackgroundField(self, x=1., y=0., z=0.): # Primary field in x-direction (background) - self.B0 = np.r_[1.,0,0] + self.B0 = np.r_[x,y,z] + + @property + def Qfx(self): + if getattr(self, '_Qfx', None) is None: + self._Qfx = self.prob.mesh.getInterpolationMat(self.rxLoc,'Fx') + return self._Qfx + + @property + def Qfy(self): + if getattr(self, '_Qfy', None) is None: + self._Qfy = self.prob.mesh.getInterpolationMat(self.rxLoc,'Fy') + return self._Qfy + + @property + def Qfz(self): + if getattr(self, '_Qfz', None) is None: + self._Qfz = self.prob.mesh.getInterpolationMat(self.rxLoc,'Fz') + return self._Qfz + + def projectFields(self, B): + """ + This function projects the fields onto the data space. + + + .. math:: + d_\\text{pred} = \mathbf{P} u(m) + """ + + bfx = self.Qfx*B + bfy = self.Qfy*B + bfz = self.Qfz*B + + return np.sqrt(bfx**2 + bfy**2 + bfz**2) + + def projectFieldsAsVector(self, B): + + bfx = self.Qfx*B + bfy = self.Qfy*B + bfz = self.Qfz*B + + return np.c_[bfx, bfy, bfz] + +class MagDataBx(object): + """docstring for MagDataBx""" + def __init__(self, **kwargs): + Data.BaseData.__init__(self, **kwargs) + + def projectFields(self, B): + bfx = self.Qfx*B + return bfx class BaseMagModel(Model.BaseModel): """BaseMagModel""" def __init__(self, mesh, **kwargs): - Model.BaseModel.__init__(mesh, **kwargs) + Model.BaseModel.__init__(self, mesh) def transform(self, m, asMu=True): if asMu: diff --git a/simpegPF/MagAnalytics.py b/simpegPF/MagAnalytics.py index 08fad64e..41f2ad48 100644 --- a/simpegPF/MagAnalytics.py +++ b/simpegPF/MagAnalytics.py @@ -1,5 +1,5 @@ from scipy.constants import mu_0 -from SimPEG.Utils.sputils import kron3, speye, sdiag +from SimPEG.Utils.matutils import kron3, speye, sdiag import matplotlib.pyplot as plt from SimPEG import * import numpy as np diff --git a/simpegPF/Magnetics.py b/simpegPF/Magnetics.py index 3447e02f..c81af597 100644 --- a/simpegPF/Magnetics.py +++ b/simpegPF/Magnetics.py @@ -1,6 +1,7 @@ -from SimPEG import Problem, Utils +from SimPEG import Mesh, Problem, Utils, np, sp import BaseMag from scipy.constants import mu_0 +from MagAnalytics import spheremodel, CongruousMagBC @@ -11,27 +12,27 @@ class MagneticsDiffSecondary(Problem.BaseProblem): modelPair = BaseMag.BaseMagModel def __init__(self, mesh, model, **kwargs): - Problem.BaseProblem.__init__(mesh, model, **kwargs) + Problem.BaseProblem.__init__(self, mesh, model, **kwargs) - self._Pbc, self._Pin, self._Pout = / + Pbc, Pin, self._Pout = \ self.mesh.getBCProjWF('neumann', discretization='CC') Dface = self.mesh.faceDiv Mc = Utils.sdiag(self.mesh.vol) - self._Div = Mc*Dface*self._Pin.T*self._Pin + self._Div = Mc*Dface*Pin.T*Pin @property def MfMuI(self): return self._MfMuI @property - def Mfmu0(self): return self._Mfmu0 + def MfMu0(self): return self._MfMu0 def makeMassMatrices(self, m): mu = self.model.transform(m) - Mfmui = self.mesh.getFaceInnerProduct(1./mu) + MfMui = self.mesh.getFaceInnerProduct(1./mu) #TODO: this will break if tensor mu - self._MfmuI = Utils.sdiag(1./Mfmui.diagonal()) - self._Mfmu0 = self.mesh.getFaceInnerProduct(1/mu_0) + self._MfMuI = Utils.sdiag(1./MfMui.diagonal()) + self._MfMu0 = self.mesh.getFaceInnerProduct(1/mu_0) def getRHS(self, m): b0 = self.data.B0 @@ -45,7 +46,7 @@ class MagneticsDiffSecondary(Problem.BaseProblem): chi = self.model.transform(m, asMu=False) Bbc = CongruousMagBC(self.mesh, self.data.B0, chi) - return -self._Div*self.MfmuI*self.Mfmu0*B0 + self._Div*B0 - Mc*Dface*self._Pout.T*Bbc + return -self._Div*self.MfMuI*self.MfMu0*B0 + self._Div*B0 - Mc*Dface*self._Pout.T*Bbc def getA(self, m): """ @@ -60,15 +61,74 @@ class MagneticsDiffSecondary(Problem.BaseProblem): """ - return -self._Div*self.MfmuI*self._Div.T + return -self._Div*self.MfMuI*self._Div.T def fields(self, m): self.makeMassMatrices(m) + #TODO: change to pos def A A = self.getA(m) rhs = self.getRHS(m) + m1 = sp.linalg.interface.aslinearoperator(Utils.sdiag(-1/A.diagonal())) + phi, info = sp.linalg.bicgstab(A, rhs, tol=1e-6, maxiter=1000, M=m1) + + #TODO: make onPair function call + b0 = self.data.B0 + B0 = np.r_[b0[0]*np.ones(self.mesh.nFx), + b0[1]*np.ones(self.mesh.nFy), + b0[2]*np.ones(self.mesh.nFz)] + + B = self.MfMuI*self.MfMu0*B0-B0-self.MfMuI*self._Div.T*phi + + return B + # F = self.getInitialFields() # return self.forward(m, self.getRHS, self.calcFields, F=F) +if __name__ == '__main__': + import matplotlib.pyplot as plt + hxind = ((5,25,1.3),(41, 12.5),(5,25,1.3)) + hyind = ((5,25,1.3),(41, 12.5),(5,25,1.3)) + hzind = ((5,25,1.3),(40, 12.5),(1,25,1.3)) + hx, hy, hz = Utils.meshTensors(hxind, hyind, hzind) + mesh = Mesh.TensorMesh([hx, hy, hz], [-hx.sum()/2,-hy.sum()/2,-hz.sum()/2]) + + chibkg = 0. + chiblk = 0.01 + chi = np.ones(mesh.nC)*chibkg + sph_ind = spheremodel(mesh, 0., 0., 0., 100) + chi[sph_ind] = chiblk + model = BaseMag.BaseMagModel(mesh) + # mu = (1.+chi)*mu_0 + + data = BaseMag.BaseMagData() + data.setBackgroundField(x=1., y=1., z=0.) + xr = np.linspace(-300, 300, 41) + yr = np.linspace(-300, 300, 41) + X, Y = np.meshgrid(xr, yr) + Z = np.ones((xr.size, yr.size))*150 + rxLoc = np.c_[Utils.mkvc(X), Utils.mkvc(Y), Utils.mkvc(Z)] + data.rxLoc = rxLoc + + prob = MagneticsDiffSecondary(mesh, model) + + prob.pair(data) + + B = prob.fields(chi) + mesh.plotSlice(B, 'F', view='vec', showIt=True) + + dpred = data.dpred(chi, u=B) + + plt.pcolor(X, Y, dpred.reshape(X.shape, order='F')) + + plt.show() + + + + + + + +