wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-07-08 19:29:00 +08:00
Code Issues Packages Projects Releases Wiki Activity

forward problem changes in data projection. have broken the derivative test.

Browse Source
This commit is contained in:
rowanc1
2014-03-21 17:09:39 -07:00
parent 837a114ad7
commit a8f09ee869
5 changed files with 130 additions and 47 deletions
Download Patch File Download Diff File Expand all files Collapse all files
simpegEM/FDEM/FDEM.py
+9 -9
View File
@@ -110,14 +110,15 @@ class BaseProblemFDEM(Problem.BaseProblem):
u_tx = u[tx, self.solType]
w = self.getADeriv(freq, u_tx, v)
Ainvw = solver.solve(w)
fAinvw = self.calcFields(Ainvw, freq, tx.rxList.fieldType)
P = tx.projectFieldsDeriv(self.mesh, u)
fAinvw = self.calcFieldsList(Ainvw, freq, tx.rxList.fieldTypes)
P = lambda v: tx.projectFieldsDeriv(self.mesh, u, v)
df_dm = self.calcFieldsDeriv(u_tx, freq, tx.rxList.fieldType, v)
df_dm = self.calcFieldsDerivList(u_tx, freq, tx.rxList.fieldTypes, v)
#TODO: this is now a list?
if df_dm is None:
Jv[tx] = - P*fAinvw
Jv[tx] = - P(fAinvw)
else:
Jv[tx] = - P*fAinvw + P*df_dm
Jv[tx] = - P(fAinvw) + P(df_dm)
return Utils.mkvc(Jv)
@@ -140,14 +141,13 @@ class BaseProblemFDEM(Problem.BaseProblem):
for tx in self.survey.getTransmitters(freq):
u_tx = u[tx, self.solType]
P = tx.projectFieldsDeriv(self.mesh, u)
PTv = P.T * v[tx]
fPTv = self.calcFields(PTv, freq, tx.rxList.fieldType, adjoint=True)
PTv = tx.projectFieldsDeriv(self.mesh, u, v[tx], adjoint=True)
fPTv = self.calcFields(PTv, freq, tx.rxList.fieldTypes, adjoint=True)
w = solver.solve( fPTv )
Jtv_tx = self.getADeriv(freq, u_tx, w, adjoint=True)
df_dm = self.calcFieldsDeriv(u_tx, freq, tx.rxList.fieldType, PTv, adjoint=True)
df_dm = self.calcFieldsDeriv(u_tx, freq, tx.rxList.fieldTypes, PTv, adjoint=True)
if df_dm is None:
Jtv += - Jtv_tx
simpegEM/FDEM/SurveyFDEM.py
+81 -25
View File
@@ -3,29 +3,71 @@ from SimPEG import Survey, Utils, np, sp
class RxListFDEM(Survey.BaseRxList):
knownRxTypes = {
'ex':'Ex',
'ey':'Ey',
'ez':'Ez',
'exr':['e', 'Ex', 'real'],
'eyr':['e', 'Ey', 'real'],
'ezr':['e', 'Ez', 'real'],
'exi':['e', 'Ex', 'imag'],
'eyi':['e', 'Ey', 'imag'],
'ezi':['e', 'Ez', 'imag'],
'bx':'Fx',
'by':'Fy',
'bz':'Fz',
'bxr':['b', 'Fx', 'real'],
'byr':['b', 'Fy', 'real'],
'bzr':['b', 'Fz', 'real'],
'bxi':['b', 'Fx', 'imag'],
'byi':['b', 'Fy', 'imag'],
'bzi':['b', 'Fz', 'imag'],
}
def __init__(self, locs, rxType):
Survey.BaseRxList.__init__(self, locs, rxType)
self._Ps = {}
for rx in self.rxTypes:
self._Ps[self._projGLoc(rx)] = {}
def _projField(self, rx):
"""Field Type projection (e.g. e b ...)"""
if type(rx) is int: rx = self.rxTypes[rx]
return self.knownRxTypes[rx][0]
def _projGLoc(self, rx):
"""Grid Location projection (e.g. Ex Fy ...)"""
if type(rx) is int: rx = self.rxTypes[rx]
return self.knownRxTypes[rx][1]
def _projComp(self, rx):
"""Component projection (real/imag)"""
if type(rx) is int: rx = self.rxTypes[rx]
return self.knownRxTypes[rx][2]
@property
def fieldType(self):
def rxTypes(self):
"""A list of the recieve types that are collected at this rxList locations."""
return self.rxType.split(',')
@property
def fieldTypes(self):
#TODO: This assumes that it has the structure ex, by ...
return self.rxType[0]
return [self._projField(rx) for rx in self.rxTypes]
def getP(self, mesh):
if mesh not in self._Ps:
self._Ps[mesh] = mesh.getInterpolationMat(self.locs, self.knownRxTypes[self.rxType])
return self._Ps[mesh]
"""
Returns the projection matrices as a
list for all components collected by
the receivers.
.. note::
Projection matrices are stored as a nested dict,
First gridLocation, then mesh.
"""
P = []
for rx in self.rxTypes:
gloc = self._projGLoc(rx)
if mesh not in self._Ps[gloc]:
self._Ps[gloc][mesh] = mesh.getInterpolationMat(self.locs, gloc)
P += [self._Ps[gloc][mesh]]
return P
class TxFDEM(Survey.BaseTx):
@@ -43,25 +85,39 @@ class TxFDEM(Survey.BaseTx):
@property
def nD(self):
"""Number of data"""
return self.rxList.locs.shape[0]
return self.rxList.locs.shape[0]*len(self.rxList.rxTypes)
def projectFields(self, mesh, u):
if self.rxList.rxType in ['ex', 'ey', 'ez']:
u_part = u[self, 'e']
elif self.rxList.rxType in ['bx', 'by', 'bz']:
u_part = u[self, 'b']
else:
raise NotImplemented('Unknown receiver type.')
nRt = len(self.rxList.rxTypes)
Pu = range(nRt)
Ps = self.rxList.getP(mesh)
P = self.rxList.getP(mesh)
Pu = P*u_part
return Pu
for ii, rx in enumerate(self.rxList.rxTypes):
fieldType = self.rxList._projField(rx)
u_part_complex = u[self, fieldType]
# get the real or imag component
real_or_imag = self.rxList._projComp(rx)
u_part = getattr(u_part_complex, real_or_imag)
Pu[ii] = Ps[ii]*u_part
def projectFieldsDeriv(self, mesh, u):
P = self.rxList.getP(mesh)
return P
return np.concatenate(Pu)
def projectFieldsDeriv(self, mesh, u, v, adjoint=False):
Ps = self.rxList.getP(mesh)
V = v.reshape((-1, len(Ps)), order='F')
Pvs = range(len(Ps))
for ii, rx in enumerate(self.rxList.rxTypes):
Pv = Ps[ii] * V[:, ii]
real_or_imag = self.rxList._projComp(rx)
if real_or_imag == 'imag':
Pv = 1j*Pv
elif real_or_imag == 'real':
Pv = Pv.astype(complex)
else:
raise NotImplementedError('must be real or imag')
return np.concatenate(Pvs)
class FieldsFDEM(object):
simpegEM/Tests/test_FDEM.py
+4 -4
View File
@@ -18,9 +18,9 @@ def getProblem(fdemType):
x = np.linspace(5,10,3)
XYZ = Utils.ndgrid(x,x,np.r_[0])
if fdemType == 'e':
rxList = EM.FDEM.RxListFDEM(XYZ, 'bx')
rxList = EM.FDEM.RxListFDEM(XYZ, 'bxr,exi')
elif fdemType == 'b':
rxList = EM.FDEM.RxListFDEM(XYZ, 'ex')
rxList = EM.FDEM.RxListFDEM(XYZ, 'exi')
else:
raise NotImplementedError()
@@ -29,9 +29,9 @@ def getProblem(fdemType):
x = np.linspace(5,10,3)
XYZ = Utils.ndgrid(x,x,np.r_[0])
if fdemType == 'e':
rxList = EM.FDEM.RxListFDEM(XYZ, 'ey')
rxList = EM.FDEM.RxListFDEM(XYZ, 'eyi')
elif fdemType == 'b':
rxList = EM.FDEM.RxListFDEM(XYZ, 'ey')
rxList = EM.FDEM.RxListFDEM(XYZ, 'eyr')
else:
raise NotImplementedError()
simpegEM/Tests/test_FDEM_analytics.py
+1 -1
View File
@@ -21,7 +21,7 @@ class FDEM_analyticTests(unittest.TestCase):
x = np.linspace(-10,10,5)
XYZ = Utils.ndgrid(x,np.r_[0],np.r_[0])
rxList = EM.FDEM.RxListFDEM(XYZ, 'ex')
rxList = EM.FDEM.RxListFDEM(XYZ, 'exi')
Tx0 = EM.FDEM.TxFDEM(np.r_[0.,0.,0.], 'VMD', 1e2, rxList)
survey = EM.FDEM.SurveyFDEM([Tx0])
simpegEM/Tests/test_FieldsObject.py
+35 -8
View File
@@ -5,20 +5,27 @@ import simpegEM as EM
class FieldsTest(unittest.TestCase):
def setUp(self):
x = np.linspace(5,10,3)
XYZ = Utils.ndgrid(x,x,np.r_[0])
rxList = EM.FDEM.RxListFDEM(XYZ, 'ex')
Tx0 = EM.FDEM.TxFDEM(None, 'VMD', 3, rxList)
Tx1 = EM.FDEM.TxFDEM(None, 'VMD', 3, rxList)
Tx2 = EM.FDEM.TxFDEM(None, 'VMD', 2, rxList)
Tx3 = EM.FDEM.TxFDEM(None, 'VMD', 1, rxList)
txList = [Tx0,Tx1,Tx2,Tx3]
mesh = Mesh.TensorMesh([np.ones(n)*5 for n in [10,11,12]],[0,0,-30])
x = np.linspace(5,10,3)
XYZ = Utils.ndgrid(x,x,np.r_[0.])
txLoc = np.r_[0,0,0.]
rxList0 = EM.FDEM.RxListFDEM(XYZ, 'exi,exr,eyi,eyr,ezi,ezr')
Tx0 = EM.FDEM.TxFDEM(txLoc, 'VMD', 3., rxList0)
rxList1 = EM.FDEM.RxListFDEM(XYZ, 'bxi,bxr,byi,byr,bzi,bzr')
Tx1 = EM.FDEM.TxFDEM(txLoc, 'VMD', 3., rxList1)
rxList2 = EM.FDEM.RxListFDEM(XYZ, 'bxi,eyr')
Tx2 = EM.FDEM.TxFDEM(txLoc, 'VMD', 2., rxList2)
rxList3 = EM.FDEM.RxListFDEM(XYZ, 'bxi')
Tx3 = EM.FDEM.TxFDEM(txLoc, 'VMD', 2., rxList3)
Tx4 = EM.FDEM.TxFDEM(txLoc, 'VMD', 1., rxList0)
txList = [Tx0,Tx1,Tx2,Tx3,Tx4]
survey = EM.FDEM.SurveyFDEM(txList)
self.F = EM.FDEM.FieldsFDEM(mesh, survey)
self.D = EM.FDEM.DataFDEM(survey)
self.Tx0 = Tx0
self.Tx1 = Tx1
self.mesh = mesh
self.XYZ = XYZ
def test_SetGet(self):
F = self.F
@@ -82,7 +89,27 @@ class FieldsTest(unittest.TestCase):
D2 = EM.FDEM.DataFDEM(self.D.survey, V)
self.assertTrue(np.all(Utils.mkvc(D2) == Utils.mkvc(self.D)))
def test_FieldProjections(self):
F = self.F
for freq in F.survey.freqs:
nFreq = F.survey.nTx[freq]
e = np.random.rand(F.mesh.nE, nFreq)
b = np.random.rand(F.mesh.nF, nFreq)
F[freq] = {'b':b,'e':e}
Txs = F.survey.getTransmitters(freq)
for ii, tx in enumerate(Txs):
dat = tx.projectFields(self.mesh, F)
self.assertTrue(dat.dtype == float)
dat = dat.reshape((self.XYZ.shape[0], len(tx.rxList.rxTypes)), order='F')
for jj, rx in enumerate(tx.rxList.rxTypes):
fieldType = tx.rxList._projField(rx)
u = {'b':b[:,ii], 'e': e[:,ii]}[fieldType]
real_or_imag = tx.rxList._projComp(rx)
u = getattr(u, real_or_imag)
gloc = tx.rxList._projGLoc(rx)
d = self.mesh.getInterpolationMat(self.XYZ, gloc)*u
self.assertTrue(np.all(dat[:, jj] == d))