Initial fields and data implementation.

This commit is contained in:
rowanc1
2014-03-10 22:47:43 -07:00
parent f10c029161
commit e1010f9eca
4 changed files with 35 additions and 48 deletions
+3 -7
View File
@@ -93,8 +93,6 @@ class ProblemFDEM_e(Problem.BaseProblem):
def Jvec(self, m, v, u=None):
# TODO: only 1 transmitter for now
# TODO: all P's the same
if u is None:
u = self.fields(m)
@@ -108,17 +106,15 @@ class ProblemFDEM_e(Problem.BaseProblem):
for tx in self.survey.getTransmitters(freq):
dMe_dsig = self.mesh.getEdgeInnerProductDeriv(m, v=e)
dsig_dm = self.model.transformDeriv(m)
b = 1j*omega(freq) * ( dMe_dsig * ( dsig_dm * v ) )
Ab = solver.solve(b)
#TODO: look at Rx for this...
P = self.survey.projectFieldsDeriv(u)
P = tx.projectFieldsDeriv(self.mesh, u)
Jv[tx] = -P*Ab
Jv = np.concatenate(Jvs)
return Jv
return Utils.mkvc(Jv)
def Jtvec(self, m, v, u=None):
+18 -32
View File
@@ -35,15 +35,18 @@ class TxFDEM(Survey.BaseTx):
def projectFields(self, mesh, u):
if self.rxList.rxType in ['Ex', 'Ey', 'Ez']:
u_part = u[self, 'e']
else:
raise NotImplemented('Unknown receiver type.')
P = self.rxList.getP(mesh)
u_part = u[self]
Pu = P*u_part
return Pu
def projectFieldsDeriv(self, mesh, u):
pass
P = self.rxList.getP(mesh)
return P
class FieldsFDEM(object):
@@ -119,8 +122,10 @@ class FieldsFDEM(object):
elif isinstance(key, TxFDEM):
raise Exception('Cannot set one transmitter at a time.')
for field in newFields:
for name in newFields:
field = self._initStore(name)
if field[freq].shape[1] == 1:
newFields[name] = Utils.mkvc(newFields[name],2)
assert field[freq].shape == newFields[name].shape, 'Must be correct shape (n%s x nTx[freq])' % self.knownFields[name]
field[freq] = newFields[name]
@@ -173,13 +178,13 @@ class DataFDEM(object):
self._ensureCorrectKey(key)
return self._dataDict[key]
def toarray(self):
def tovec(self):
D = self._dataDict
return np.concatenate([D[k] for k in D])
class SurveyFDEM(Survey.MixinFancyProjection, Survey.BaseSurvey):
class SurveyFDEM(Survey.BaseSurvey):
"""
docstring for SurveyFDEM
"""
@@ -228,35 +233,16 @@ class SurveyFDEM(Survey.MixinFancyProjection, Survey.BaseSurvey):
self._nTx[freq] = len(self.getTransmitters(freq))
return self._nTx
def getTransmitters(self, freq):
"""Returns the transmitters associated with a specific frequency."""
assert freq in self._freqDict, "The requested frequency is not in this survey."
return self._freqDict[freq]
def projectFields(self, u):
P = sp.identity(self.prob.mesh.nE)
Pes = range(self.nFreq)
#TODO: this is hardcoded to 1Tx
for i, freqInd in enumerate(range(self.nFreq)):
e = u.get_e(freqInd)
Pes[i] = P*e
Pe = np.concatenate(Pes)
return Pe
data = DataFDEM(self)
for tx in self.txList:
data[tx] = tx.projectFields(self.mesh, u)
return data
def projectFieldsDerivVec(self, u, v=None):
raise NotImplemented('projectFieldsDerivVec is not yet implemented')
def projectAdjointFieldsDerivVec(self, u, v=None):
raise NotImplemented('projectAdjointFieldsDerivVec is not yet implemented')
####################################################
# Interpolation Matrices
####################################################
@Utils.requires('prob')
def getP(self, Tx):
# TODO: store these in a mesh lookup table by transmitter?
pass
def projectFieldsDeriv(self, u):
raise Exception('Use Transmitters to project fields deriv.')
+6 -7
View File
@@ -21,13 +21,12 @@ class FDEM_bDerivTests(unittest.TestCase):
model = Model.LogModel(mesh)
opts = {'txLoc':0.,
'txType':'VMD_MVP',
'rxLoc': rxLoc,
'rxType':'bz',
'freq': np.logspace(0,3,4),
}
dat = EM.FDEM.SurveyFDEM(**opts)
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)
dat = EM.FDEM.SurveyFDEM([Tx0])
prb = EM.FDEM.ProblemFDEM_e(model)
prb.pair(dat)
+8 -2
View File
@@ -22,10 +22,16 @@ class FieldsTest(unittest.TestCase):
def test_SetGet(self):
F = self.F
for freq in F.survey.freqs:
e = np.random.rand(F.mesh.nE, F.survey.nTx[freq])
nFreq = F.survey.nTx[freq]
e = np.random.rand(F.mesh.nE, nFreq)
F[freq, 'e'] = e
b = np.random.rand(F.mesh.nF, F.survey.nTx[freq])
b = np.random.rand(F.mesh.nF, nFreq)
F[freq, 'b'] = b
if nFreq == 1:
F[freq, 'b'] = Utils.mkvc(b)
self.assertTrue(np.all(F[freq, 'e'] == e))
self.assertTrue(np.all(F[freq, 'b'] == b))
F[freq] = {'b':b,'e':e}
self.assertTrue(np.all(F[freq, 'e'] == e))
self.assertTrue(np.all(F[freq, 'b'] == b))