mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-05 10:34:26 +08:00
Changed Data --> Survey
This commit is contained in:
+9
-9
@@ -1,15 +1,15 @@
|
||||
from SimPEG import Model, Data, Utils, np, sp
|
||||
from SimPEG import Model, Survey, Utils, np, sp
|
||||
from scipy.constants import mu_0
|
||||
|
||||
|
||||
class BaseMagData(Data.BaseData):
|
||||
"""Base Magnetics Data"""
|
||||
class BaseMagSurvey(Survey.BaseSurvey):
|
||||
"""Base Magnetics Survey"""
|
||||
|
||||
rxLoc = None #: receiver locations
|
||||
rxType = None #: receiver type
|
||||
|
||||
def __init__(self, **kwargs):
|
||||
Data.BaseData.__init__(self, **kwargs)
|
||||
Survey.BaseSurvey.__init__(self, **kwargs)
|
||||
|
||||
|
||||
def setBackgroundField(self, Inc, Dec, Btot):
|
||||
@@ -65,7 +65,7 @@ class BaseMagData(Data.BaseData):
|
||||
bfz = self.Qfz*u['B']
|
||||
|
||||
# Generate unit vector
|
||||
B0 = self.prob.data.B0
|
||||
B0 = self.prob.survey.B0
|
||||
Bot = np.sqrt(B0[0]**2+B0[1]**2+B0[2]**2)
|
||||
box = B0[0]/Bot
|
||||
boy = B0[1]/Bot
|
||||
@@ -88,7 +88,7 @@ class BaseMagData(Data.BaseData):
|
||||
|
||||
"""
|
||||
# Generate unit vector
|
||||
B0 = self.prob.data.B0
|
||||
B0 = self.prob.survey.B0
|
||||
Bot = np.sqrt(B0[0]**2+B0[1]**2+B0[2]**2)
|
||||
box = B0[0]/Bot
|
||||
boy = B0[1]/Bot
|
||||
@@ -105,10 +105,10 @@ class BaseMagData(Data.BaseData):
|
||||
|
||||
return np.r_[bfx, bfy, bfz]
|
||||
|
||||
class MagDataBx(object):
|
||||
"""docstring for MagDataBx"""
|
||||
class MagSurveyBx(object):
|
||||
"""docstring for MagSurveyBx"""
|
||||
def __init__(self, **kwargs):
|
||||
Data.BaseData.__init__(self, **kwargs)
|
||||
Survey.BaseData.__init__(self, **kwargs)
|
||||
|
||||
def projectFields(self, B):
|
||||
bfx = self.Qfx*B
|
||||
|
||||
@@ -10,11 +10,11 @@ class MagneticsDiffSecondary(Problem.BaseProblem):
|
||||
Secondary field approach using differential equations!
|
||||
"""
|
||||
|
||||
dataPair = BaseMag.BaseMagData
|
||||
surveyPair = BaseMag.BaseMagSurvey
|
||||
modelPair = BaseMag.BaseMagModel
|
||||
|
||||
def __init__(self, mesh, model, **kwargs):
|
||||
Problem.BaseProblem.__init__(self, mesh, model, **kwargs)
|
||||
def __init__(self, model, **kwargs):
|
||||
Problem.BaseProblem.__init__(self, model, **kwargs)
|
||||
|
||||
Pbc, Pin, self._Pout = \
|
||||
self.mesh.getBCProjWF('neumann', discretization='CC')
|
||||
@@ -41,9 +41,9 @@ class MagneticsDiffSecondary(Problem.BaseProblem):
|
||||
self._MfMu0 = self.mesh.getFaceMass(1/mu_0)
|
||||
# self._MfMu0 = self.mesh.getFaceInnerProduct(1/mu_0)
|
||||
|
||||
@Utils.requires('data')
|
||||
@Utils.requires('survey')
|
||||
def getB0(self):
|
||||
b0 = self.data.B0
|
||||
b0 = self.survey.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)]
|
||||
@@ -64,7 +64,7 @@ class MagneticsDiffSecondary(Problem.BaseProblem):
|
||||
mu = self.model.transform(m)
|
||||
chi = mu/mu_0-1
|
||||
|
||||
Bbc, Bbc_const = CongruousMagBC(self.mesh, self.data.B0, chi)
|
||||
Bbc, Bbc_const = CongruousMagBC(self.mesh, self.survey.B0, chi)
|
||||
self.Bbc = Bbc
|
||||
self.Bbc_const = Bbc_const
|
||||
return self._Div*self.MfMuI*self.MfMu0*B0 - self._Div*B0 + Mc*Dface*self._Pout.T*Bbc
|
||||
@@ -191,7 +191,7 @@ class MagneticsDiffSecondary(Problem.BaseProblem):
|
||||
vol = self.mesh.vol
|
||||
Div = self._Div
|
||||
Dface = self.mesh.faceDiv
|
||||
P = self.data.projectFieldsDeriv(B) # Projection matrix
|
||||
P = self.survey.projectFieldsDeriv(B) # Projection matrix
|
||||
B0 = self.getB0()
|
||||
|
||||
MfMuIvec = 1/self.MfMui.diagonal()
|
||||
@@ -268,7 +268,7 @@ class MagneticsDiffSecondary(Problem.BaseProblem):
|
||||
vol = self.mesh.vol
|
||||
Div = self._Div
|
||||
Dface = self.mesh.faceDiv
|
||||
P = self.data.projectFieldsDeriv(B) # Projection matrix
|
||||
P = self.survey.projectFieldsDeriv(B) # Projection matrix
|
||||
B0 = self.getB0()
|
||||
|
||||
MfMuIvec = 1/self.MfMui.diagonal()
|
||||
@@ -389,7 +389,7 @@ if __name__ == '__main__':
|
||||
prob.pair(data)
|
||||
|
||||
dpred = data.dpred(chi)
|
||||
|
||||
|
||||
# fig = plt.figure( figsize = (8,5) )
|
||||
# ax = plt.subplot(111)
|
||||
# dat = plt.imshow(np.reshape(dpred, (xr.size, yr.size), order='F'), extent=[min(xr), max(xr), min(yr), max(yr)])
|
||||
|
||||
@@ -20,7 +20,7 @@ class MagFwdProblemTests(unittest.TestCase):
|
||||
sph_ind = PF.MagAnalytics.spheremodel(M, 0., 0., 0., 100)
|
||||
chi[sph_ind] = chiblk
|
||||
model = PF.BaseMag.BaseMagModel(M)
|
||||
prob = PF.Magnetics.MagneticsDiffSecondary(M, model)
|
||||
prob = PF.Magnetics.MagneticsDiffSecondary(model)
|
||||
self.prob = prob
|
||||
self.M = M
|
||||
self.chi = chi
|
||||
@@ -28,39 +28,31 @@ class MagFwdProblemTests(unittest.TestCase):
|
||||
|
||||
def test_anal_forward(self):
|
||||
|
||||
data = PF.BaseMag.BaseMagData()
|
||||
survey = PF.BaseMag.BaseMagSurvey()
|
||||
|
||||
Inc = 45.
|
||||
Dec = 45.
|
||||
Btot = 51000
|
||||
|
||||
b0 = PF.MagAnalytics.IDTtoxyz(Inc, Dec, Btot)
|
||||
data.setBackgroundField(Inc, Dec, Btot)
|
||||
survey.setBackgroundField(Inc, Dec, Btot)
|
||||
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
|
||||
survey.rxLoc = rxLoc
|
||||
|
||||
self.prob.pair(data)
|
||||
self.prob.pair(survey)
|
||||
u = self.prob.fields(self.chi)
|
||||
B = u['B']
|
||||
|
||||
bxa,bya,bza = PF.MagAnalytics.MagSphereAnalFunA(rxLoc[:,0],rxLoc[:,1],rxLoc[:,2],100.,0.,0.,0.,0.01, b0,'secondary')
|
||||
|
||||
dpred = data.projectFieldsAsVector(B)
|
||||
dpred = survey.projectFieldsAsVector(B)
|
||||
err = np.linalg.norm(dpred-np.r_[bxa, bya, bza])/np.linalg.norm(np.r_[bxa, bya, bza])
|
||||
|
||||
# plt.plot(dpred)
|
||||
# plt.plot(np.r_[bxa, bya, bza])
|
||||
# plt.show()
|
||||
|
||||
if err > 0.05:
|
||||
raise Exception('Anaytic test is failed T.T')
|
||||
else:
|
||||
print "Anaytic test is passed"
|
||||
pass
|
||||
self.assertTrue(err < 0.05)
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
|
||||
@@ -28,27 +28,27 @@ class MagSensProblemTests(unittest.TestCase):
|
||||
chi[sph_ind] = chiblk
|
||||
model = PF.BaseMag.BaseMagModel(M)
|
||||
|
||||
data = BaseMag.BaseMagData()
|
||||
survey = BaseMag.BaseMagSurvey()
|
||||
|
||||
data.setBackgroundField(Inc, Dec, Btot)
|
||||
survey.setBackgroundField(Inc, Dec, Btot)
|
||||
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
|
||||
survey.rxLoc = rxLoc
|
||||
|
||||
|
||||
|
||||
prob = PF.Magnetics.MagneticsDiffSecondary(M, model)
|
||||
prob.pair(data)
|
||||
dpre = data.dpred(chi)
|
||||
prob = PF.Magnetics.MagneticsDiffSecondary(model)
|
||||
prob.pair(survey)
|
||||
dpre = survey.dpred(chi)
|
||||
|
||||
fields = prob.fields(chi)
|
||||
self.u = fields['u']
|
||||
self.B = fields['B']
|
||||
|
||||
self.data = data
|
||||
self.survey = survey
|
||||
self.model = model
|
||||
self.prob = prob
|
||||
self.M = M
|
||||
@@ -145,7 +145,7 @@ class MagSensProblemTests(unittest.TestCase):
|
||||
self.prob.makeMassMatrices(chi)
|
||||
dmudm = self.model.transformDeriv(chi)
|
||||
dchidmu = Utils.sdiag(1/(dmudm.diagonal()))
|
||||
Bbc, Bbc_const = PF.MagAnalytics.CongruousMagBC(self.prob.mesh, self.data.B0, chi)
|
||||
Bbc, Bbc_const = PF.MagAnalytics.CongruousMagBC(self.prob.mesh, self.survey.B0, chi)
|
||||
MfMuIvec = 1/self.prob.MfMui.diagonal()
|
||||
dMfMuI = Utils.sdiag(MfMuIvec**2)*aveF2CC.T*Utils.sdiag(vol*1./mu**2)
|
||||
RHS1 = Div*self.prob.MfMuI*self.prob.MfMu0*B0
|
||||
@@ -161,7 +161,7 @@ class MagSensProblemTests(unittest.TestCase):
|
||||
self.prob.makeMassMatrices(chi)
|
||||
dmudm = self.model.transformDeriv(chi)
|
||||
dmdmu = Utils.sdiag(1/(dmudm.diagonal()))
|
||||
Bbc, Bbc_const = PF.MagAnalytics.CongruousMagBC(self.prob.mesh, self.data.B0, chi)
|
||||
Bbc, Bbc_const = PF.MagAnalytics.CongruousMagBC(self.prob.mesh, self.survey.B0, chi)
|
||||
MfMuIvec = 1/self.prob.MfMui.diagonal()
|
||||
dMfMuI = Utils.sdiag(MfMuIvec**2)*aveF2CC.T*Utils.sdiag(vol*1./mu**2)
|
||||
dCdm_RHS1 = Div * (Utils.sdiag( self.prob.MfMu0*B0 ) * dMfMuI)
|
||||
@@ -199,7 +199,7 @@ class MagSensProblemTests(unittest.TestCase):
|
||||
u = self.u
|
||||
dmudm = self.model.transformDeriv(chi)
|
||||
dmdmu = Utils.sdiag(1/(dmudm.diagonal()))
|
||||
Bbc, Bbc_const = PF.MagAnalytics.CongruousMagBC(self.prob.mesh, self.data.B0, chi)
|
||||
Bbc, Bbc_const = PF.MagAnalytics.CongruousMagBC(self.prob.mesh, self.survey.B0, chi)
|
||||
MfMuIvec = 1/self.prob.MfMui.diagonal()
|
||||
dMfMuI = Utils.sdiag(MfMuIvec**2)*aveF2CC.T*Utils.sdiag(vol*1./mu**2)
|
||||
dCdu = self.prob.getA(chi)
|
||||
@@ -248,7 +248,7 @@ class MagSensProblemTests(unittest.TestCase):
|
||||
u = self.u
|
||||
dmudm = self.model.transformDeriv(chi)
|
||||
dmdmu = Utils.sdiag(1/(dmudm.diagonal()))
|
||||
Bbc, Bbc_const = PF.MagAnalytics.CongruousMagBC(self.prob.mesh, self.data.B0, chi)
|
||||
Bbc, Bbc_const = PF.MagAnalytics.CongruousMagBC(self.prob.mesh, self.survey.B0, chi)
|
||||
MfMuIvec = 1/self.prob.MfMui.diagonal()
|
||||
dMfMuI = Utils.sdiag(MfMuIvec**2)*aveF2CC.T*Utils.sdiag(vol*1./mu**2)
|
||||
dCdu = self.prob.getA(chi)
|
||||
@@ -289,7 +289,7 @@ class MagSensProblemTests(unittest.TestCase):
|
||||
|
||||
a = self.prob.Jvec(self.chi, d_chi)
|
||||
|
||||
derChk = lambda m: [self.data.dpred(m), lambda mx: self.prob.Jvec(self.chi, mx)]
|
||||
derChk = lambda m: [self.survey.dpred(m), lambda mx: self.prob.Jvec(self.chi, mx)]
|
||||
# TODO: I am not sure why the order get worse as step decreases .. --;
|
||||
passed = Tests.checkDerivative(derChk, self.chi, num=2, dx = d_chi, plotIt=False)
|
||||
self.assertTrue(passed)
|
||||
@@ -297,10 +297,10 @@ class MagSensProblemTests(unittest.TestCase):
|
||||
def test_Jtvec(self):
|
||||
print ">> Derivative test for Jtvec"
|
||||
mu = self.model.transform(self.chi)
|
||||
dobs = self.data.dpred(self.chi)
|
||||
dobs = self.survey.dpred(self.chi)
|
||||
|
||||
def misfit (m, dobs):
|
||||
dpre = self.data.dpred(m)
|
||||
dpre = self.survey.dpred(m)
|
||||
misfit = 0.5*np.linalg.norm(dpre-dobs)**2
|
||||
residual = dpre-dobs
|
||||
dmisfit = self.prob.Jtvec(self.chi, residual)
|
||||
|
||||
Reference in New Issue
Block a user