diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index a457c84b..f88af098 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -56,9 +56,10 @@ class BaseProblemFDEM(Problem.BaseProblem): @property def MeSigmaI(self): - # TODO: this will not work if tensor conductivity + #TODO: hardcoded to sigma as the model if getattr(self, '_MeSigmaI', None) is None: - self._MeSigmaI = Utils.sdiag(1/self.MeSigma.diagonal()) + sigma = self.curTModel + self._MeSigmaI = self.mesh.getEdgeInnerProduct(sigma, invMat=True) return self._MeSigmaI curModel = Utils.dependentProperty('_curModel', None, ['_MeSigma', '_MeSigmaI', '_curTModel', '_curTModelDeriv'], 'Sets the current model, and removes dependent mass matrices.') @@ -90,7 +91,8 @@ class BaseProblemFDEM(Problem.BaseProblem): solver = self.Solver(A, **self.solverOpts) sol = solver.solve(rhs) for fieldType in self.storeTheseFields: - F[freq, fieldType] = CalcFields(sol, freq, fieldType) + Txs = self.survey.getTransmitters(freq) + F[Txs, fieldType] = CalcFields(sol, freq, fieldType) return F diff --git a/simpegEM/FDEM/SurveyFDEM.py b/simpegEM/FDEM/SurveyFDEM.py index d9c460f5..48fa1f4d 100644 --- a/simpegEM/FDEM/SurveyFDEM.py +++ b/simpegEM/FDEM/SurveyFDEM.py @@ -79,113 +79,10 @@ class TxFDEM(Survey.BaseTx): -class FieldsFDEM(object): +class FieldsFDEM(Survey.Fields): """Fancy Field Storage for a FDEM survey.""" - knownFields = {'b': 'F', 'e': 'E'} - - def __init__(self, mesh, survey): - self.survey = survey - self.mesh = mesh - self._fields = {} - - def _initStore(self, name): - if name in self._fields: - return self._fields[name] - - assert name in self.knownFields, 'field name is not known.' - - loc = self.knownFields[name] - - nP = {'CC': self.mesh.nC, - 'F': self.mesh.nF, - 'E': self.mesh.nE}[loc] - - field = {} - for freq in self.survey.freqs: - nTx_f = len(self.survey.getTransmitters(freq)) - field[freq] = np.empty((nP, nTx_f)) - - self._fields[name] = field - - return field - - def _ensureCorrectKey(self, key): - if type(key) is tuple: - assert len(key) == 2, 'must be [freq, fieldName]' - freqTest, name = key - - if name not in self.knownFields: - raise KeyError('Invalid field name') - - if type(freqTest) is float: - freq = freqTest - elif isinstance(freqTest, TxFDEM): - freq = freqTest.freq - if freqTest not in self.survey.txList: - raise KeyError('Invalid Transmitter') - else: - raise KeyError('Invalid Frequency Key') - - elif type(key) is float: - freq = key - elif isinstance(key, TxFDEM): - freq = key.freq - if key not in self.survey.txList: - raise KeyError('Invalid Transmitter') - else: - raise KeyError('Unexpected key use [freq, fieldName]') - if freq not in self.survey.freqs: - raise KeyError('Invalid frequency') - - def __setitem__(self, key, value): - self._ensureCorrectKey(key) - if type(key) is tuple: - freq, name = key - assert type(freq) is float, 'Frequency must be a float for setter.' - assert type(value) is np.ndarray, 'Must be set to a numpy array' - newFields = {name: value} - elif type(key) is float: - freq = key - assert type(value) is dict, 'New fields must be a dictionary' - newFields = value - elif isinstance(key, TxFDEM): - raise Exception('Cannot set one transmitter at a time.') - - 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] - - def __getitem__(self, key): - self._ensureCorrectKey(key) - if type(key) is tuple: - freqTest, name = key - if type(freqTest) is float: - return self._fields[name][freqTest] - elif isinstance(freqTest, TxFDEM): - key = freqTest - ind = np.array([tx is key for tx in self.survey.getTransmitters(key.freq)]) - return Utils.mkvc(self._fields[name][key.freq][:,ind]) - elif type(key) is float: - freq = key - out = {} - for name in self._fields: - out[name] = self._fields[name][freq] - return out - elif isinstance(key, TxFDEM): - freq = key.freq - ind = np.array([tx is key for tx in self.survey.getTransmitters(freq)]) - out = {} - for name in self._fields: - out[name] = Utils.mkvc(self._fields[name][freq][:,ind]) - return out - - - def __contains__(self, key): - return key in self.children + dtype = complex class SurveyFDEM(Survey.BaseSurvey): @@ -220,12 +117,12 @@ class SurveyFDEM(Survey.BaseSurvey): return len(self._freqDict) @property - def nTx(self): - if getattr(self, '_nTx', None) is None: - self._nTx = {} + def nTxByFreq(self): + if getattr(self, '_nTxByFreq', None) is None: + self._nTxByFreq = {} for freq in self.freqs: - self._nTx[freq] = len(self.getTransmitters(freq)) - return self._nTx + self._nTxByFreq[freq] = len(self.getTransmitters(freq)) + return self._nTxByFreq def getTransmitters(self, freq): """Returns the transmitters associated with a specific frequency.""" diff --git a/simpegEM/TDEM/BaseTDEM.py b/simpegEM/TDEM/BaseTDEM.py index dcee4705..72a6e6fa 100644 --- a/simpegEM/TDEM/BaseTDEM.py +++ b/simpegEM/TDEM/BaseTDEM.py @@ -1,7 +1,6 @@ from SimPEG import Solver from SimPEG.Problem import BaseTimeProblem from simpegEM.Utils import Sources -from FieldsTDEM import FieldsTDEM from scipy.constants import mu_0 from SimPEG.Utils import sdiag, mkvc from SimPEG import Utils, Mesh diff --git a/simpegEM/TDEM/FieldsTDEM.py b/simpegEM/TDEM/FieldsTDEM.py deleted file mode 100644 index 79e54ed2..00000000 --- a/simpegEM/TDEM/FieldsTDEM.py +++ /dev/null @@ -1,73 +0,0 @@ -import numpy as np - - -class FieldsTDEM(object): - """docstring for FieldsTDEM""" - - phi0 = None #: Initial electric potential - A0 = None #: Initial magnetic vector potential - e0 = None #: Initial electric field - b0 = None #: Initial magnetic flux density - j0 = None #: Initial current density - h0 = None #: Initial magnetic field - - phi = None #: Electric potential - A = None #: Magnetic vector potential - e = None #: Electric field - b = None #: Magnetic flux density - j = None #: Current density - h = None #: Magnetic field - - def __init__(self, mesh, nTx, nT, store='b'): - - self.nT = nT #: Number of times - self.nTx = nTx #: Number of transmitters - self.mesh = mesh - - def update(self, newFields, tInd): - self.set_b(newFields['b'], tInd) - self.set_e(newFields['e'], tInd) - - def fieldVec(self): - u = np.ndarray((0, self.nTx)) - for i in range(self.nT): - u = np.r_[u, self.get_b(i), self.get_e(i)] - if self.nTx == 1: - u = u.flatten() - return u - - #################################################### - # Get Methods - #################################################### - - def get_b(self, ind): - if ind == -1: - return self.b0 - else: - return self.b[ind,:,:] - - def get_e(self, ind): - if ind == -1: - return self.e0 - else: - return self.e[ind,:,:] - - #################################################### - # Set Methods - #################################################### - - def set_b(self, b, ind): - if self.b is None: - self.b = np.zeros((self.nT, np.sum(self.mesh.nF), self.nTx)) - self.b[:] = np.nan - if len(b.shape) == 1: - b = b[:, np.newaxis] - self.b[ind,:,:] = b - - def set_e(self, e, ind): - if self.e is None: - self.e = np.zeros((self.nT, np.sum(self.mesh.nE), self.nTx)) - self.e[:] = np.nan - if len(e.shape) == 1: - e = e[:, np.newaxis] - self.e[ind,:,:] = e diff --git a/simpegEM/TDEM/SurveyTDEM.py b/simpegEM/TDEM/SurveyTDEM.py index 4e712fd4..81c2421e 100644 --- a/simpegEM/TDEM/SurveyTDEM.py +++ b/simpegEM/TDEM/SurveyTDEM.py @@ -1,6 +1,5 @@ from SimPEG import Utils, np from SimPEG.Survey import BaseSurvey -from FieldsTDEM import FieldsTDEM class SurveyTDEM1D(BaseSurvey): """ @@ -51,3 +50,79 @@ class SurveyTDEM1D(BaseSurvey): self._Qrx = self.prob.mesh.getInterpolationMat(self.rxLoc, locType=locType) return self._Qrx _Qrx = None + + +class FieldsTDEM(object): + """docstring for FieldsTDEM""" + + phi0 = None #: Initial electric potential + A0 = None #: Initial magnetic vector potential + e0 = None #: Initial electric field + b0 = None #: Initial magnetic flux density + j0 = None #: Initial current density + h0 = None #: Initial magnetic field + + phi = None #: Electric potential + A = None #: Magnetic vector potential + e = None #: Electric field + b = None #: Magnetic flux density + j = None #: Current density + h = None #: Magnetic field + + def __init__(self, mesh, nTx, nT, store='b'): + + self.nT = nT #: Number of times + self.nTx = nTx #: Number of transmitters + self.mesh = mesh + + def update(self, newFields, tInd): + self.set_b(newFields['b'], tInd) + self.set_e(newFields['e'], tInd) + + def fieldVec(self): + u = np.ndarray((0, self.nTx)) + for i in range(self.nT): + u = np.r_[u, self.get_b(i), self.get_e(i)] + if self.nTx == 1: + u = u.flatten() + return u + + #################################################### + # Get Methods + #################################################### + + def get_b(self, ind): + if ind == -1: + return self.b0 + else: + return self.b[ind,:,:] + + def get_e(self, ind): + if ind == -1: + return self.e0 + else: + return self.e[ind,:,:] + + #################################################### + # Set Methods + #################################################### + + def set_b(self, b, ind): + if self.b is None: + self.b = np.zeros((self.nT, np.sum(self.mesh.nF), self.nTx)) + self.b[:] = np.nan + if len(b.shape) == 1: + b = b[:, np.newaxis] + self.b[ind,:,:] = b + + def set_e(self, e, ind): + if self.e is None: + self.e = np.zeros((self.nT, np.sum(self.mesh.nE), self.nTx)) + self.e[:] = np.nan + if len(e.shape) == 1: + e = e[:, np.newaxis] + self.e[ind,:,:] = e + + + def __contains__(self, key): + return key in self.children diff --git a/simpegEM/TDEM/TDEM_b.py b/simpegEM/TDEM/TDEM_b.py index 9e52bdf7..df05f9e6 100644 --- a/simpegEM/TDEM/TDEM_b.py +++ b/simpegEM/TDEM/TDEM_b.py @@ -1,5 +1,4 @@ from BaseTDEM import ProblemBaseTDEM -from FieldsTDEM import FieldsTDEM from SimPEG.Utils import mkvc import numpy as np from SurveyTDEM import SurveyTDEM1D diff --git a/simpegEM/TDEM/__init__.py b/simpegEM/TDEM/__init__.py index d1a7d0ac..d489d934 100644 --- a/simpegEM/TDEM/__init__.py +++ b/simpegEM/TDEM/__init__.py @@ -1,4 +1,3 @@ from BaseTDEM import ProblemBaseTDEM from SurveyTDEM import SurveyTDEM1D -from FieldsTDEM import FieldsTDEM from TDEM_b import ProblemTDEM_b diff --git a/simpegEM/Tests/test_FDEM_analytics.py b/simpegEM/Tests/test_FDEM_analytics.py index e129d2f7..d9d7408d 100644 --- a/simpegEM/Tests/test_FDEM_analytics.py +++ b/simpegEM/Tests/test_FDEM_analytics.py @@ -16,7 +16,6 @@ class FDEM_analyticTests(unittest.TestCase): hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)] hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)] mesh = Mesh.TensorMesh([hx,hy,hz], 'CCC') - print mesh.vectorCCx mapping = Maps.ExpMap(mesh) diff --git a/simpegEM/Tests/test_FieldsObject.py b/simpegEM/Tests/test_FieldsObject.py index f371c36d..f28d1ad5 100644 --- a/simpegEM/Tests/test_FieldsObject.py +++ b/simpegEM/Tests/test_FieldsObject.py @@ -29,29 +29,34 @@ class FieldsTest(unittest.TestCase): def test_SetGet(self): F = self.F for freq in F.survey.freqs: - nFreq = F.survey.nTx[freq] + nFreq = F.survey.nTxByFreq[freq] + Txs = F.survey.getTransmitters(freq) e = np.random.rand(F.mesh.nE, nFreq) - F[freq, 'e'] = e + F[Txs, 'e'] = e b = np.random.rand(F.mesh.nF, nFreq) - F[freq, 'b'] = b + F[Txs, '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)) + F[Txs, 'b'] = Utils.mkvc(b) + if e.shape[1] == 1: + e, b = Utils.mkvc(e), Utils.mkvc(b) + self.assertTrue(np.all(F[Txs, 'e'] == e)) + self.assertTrue(np.all(F[Txs, 'b'] == b)) + F[Txs] = {'b':b,'e':e} + self.assertTrue(np.all(F[Txs, 'e'] == e)) + self.assertTrue(np.all(F[Txs, 'b'] == b)) - lastFreq = F[freq] + lastFreq = F[Txs] self.assertTrue(type(lastFreq) is dict) self.assertTrue(sorted([k for k in lastFreq]) == ['b','e']) self.assertTrue(np.all(lastFreq['b'] == b)) self.assertTrue(np.all(lastFreq['e'] == e)) - self.assertTrue(F[3.,'b'].shape == (F.mesh.nF, 2)) + Tx_f3 = F.survey.getTransmitters(3.) + self.assertTrue(F[Tx_f3,'b'].shape == (F.mesh.nF, 2)) b = np.random.rand(F.mesh.nF, 2) - F[self.Tx0.freq,'b'] = b + Tx_f0 = F.survey.getTransmitters(self.Tx0.freq) + F[Tx_f0,'b'] = b self.assertTrue(F[self.Tx0]['b'].shape == (F.mesh.nF,)) self.assertTrue(F[self.Tx0,'b'].shape == (F.mesh.nF,)) self.assertTrue(np.all(F[self.Tx0,'b'] == b[:,0])) @@ -59,23 +64,25 @@ class FieldsTest(unittest.TestCase): def test_assertions(self): freq = self.F.survey.freqs[0] - bWrongSize = np.random.rand(self.F.mesh.nE, self.F.survey.nTx[freq]) - def fun(): self.F[freq, 'b'] = bWrongSize - self.assertRaises(AssertionError, fun) + Txs = self.F.survey.getTransmitters(freq) + bWrongSize = np.random.rand(self.F.mesh.nE, self.F.survey.nTxByFreq[freq]) + def fun(): self.F[Txs, 'b'] = bWrongSize + self.assertRaises(ValueError, fun) def fun(): self.F[-999.] self.assertRaises(KeyError, fun) def fun(): self.F['notRight'] self.assertRaises(KeyError, fun) - def fun(): self.F[freq,'notThere'] + def fun(): self.F[Txs,'notThere'] self.assertRaises(KeyError, fun) def test_FieldProjections(self): F = self.F for freq in F.survey.freqs: - nFreq = F.survey.nTx[freq] + nFreq = F.survey.nTxByFreq[freq] + Txs = F.survey.getTransmitters(freq) e = np.random.rand(F.mesh.nE, nFreq) b = np.random.rand(F.mesh.nF, nFreq) - F[freq] = {'b':b,'e':e} + F[Txs] = {'b':b,'e':e} Txs = F.survey.getTransmitters(freq) for ii, tx in enumerate(Txs): diff --git a/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py b/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py index 993991b9..a1665324 100644 --- a/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py +++ b/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py @@ -2,6 +2,8 @@ import unittest from SimPEG import * import simpegEM as EM +plotIt = False + class TDEM_bDerivTests(unittest.TestCase): def setUp(self): @@ -29,7 +31,7 @@ class TDEM_bDerivTests(unittest.TestCase): self.dat = EM.TDEM.SurveyTDEM1D(**opts) self.prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping) - self.prb.timeSteps = [(1e-05, 10), (5e-05, 10), (0.00025, 10)] + self.prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)] self.sigma = np.ones(mesh.nCz)*1e-8 self.sigma[mesh.vectorCCz<0] = 1e-1 @@ -146,7 +148,7 @@ class TDEM_bDerivTests(unittest.TestCase): derChk = lambda m: [self.prb.AhVec(m, f).fieldVec(), lambda mx: self.prb.Gvec(sigma, mx, u=f).fieldVec()] print '\ntest_DerivG' - passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=6, eps=1e-20) + passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20) self.assertTrue(passed) def test_Deriv_dUdM(self): @@ -162,7 +164,7 @@ class TDEM_bDerivTests(unittest.TestCase): derChk = lambda m: [self.prb.fields(m).fieldVec(), lambda mx: -prb.solveAh(sigma, prb.Gvec(sigma, mx, u=f)).fieldVec()] print '\n' print 'test_Deriv_dUdM' - passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=6, eps=1e-20) + passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20) self.assertTrue(passed) def test_Deriv_J(self): @@ -179,7 +181,7 @@ class TDEM_bDerivTests(unittest.TestCase): derChk = lambda m: [prb.survey.dpred(m), lambda mx: -prb.Jvec(sigma, mx)] print '\n' print 'test_Deriv_J' - passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=6, eps=1e-20) + passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=4, eps=1e-20) self.assertTrue(passed) def test_projectAdjoint(self): @@ -225,15 +227,21 @@ class TDEM_bDerivTests(unittest.TestCase): sigma = np.random.rand(prb.mapping.nP) f1 = EM.TDEM.FieldsTDEM(mesh, 1, prb.nT, 'b') - for i in range(f1.nT): + for i in range(prb.nT): f1.set_b(np.random.rand(mesh.nF, 1), i) f1.set_e(np.random.rand(mesh.nE, 1), i) f2 = prb.solveAht(sigma, f1) f3 = prb.AhtVec(sigma, f2) + if plotIt: + import matplotlib.pyplot as plt + plt.plot(f3.fieldVec()) + plt.plot(f1.fieldVec()) + plt.show() V1 = np.linalg.norm(f3.fieldVec()-f1.fieldVec()) V2 = np.linalg.norm(f1.fieldVec()) + print V1, V2 self.assertLess(V1/V2, 1e-6) def test_adjointsolveAhVssolveAht(self):