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

Move Fields object to base simpeg.

Browse Source
This commit is contained in:
rowanc1
2014-04-26 17:08:24 -07:00
parent 06de6a346f
commit f3f571edbd
10 changed files with 126 additions and 214 deletions
Download Patch File Download Diff File Expand all files Collapse all files
simpegEM/FDEM/FDEM.py
+5 -3
View File
@@ -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
simpegEM/FDEM/SurveyFDEM.py
+7 -110
View File
@@ -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."""
simpegEM/TDEM/BaseTDEM.py
-1
View File
@@ -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
simpegEM/TDEM/FieldsTDEM.py
-73
View File
@@ -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
simpegEM/TDEM/SurveyTDEM.py
+76 -1
View File
@@ -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
simpegEM/TDEM/TDEM_b.py
-1
View File
@@ -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
simpegEM/TDEM/__init__.py
-1
View File
@@ -1,4 +1,3 @@
from BaseTDEM import ProblemBaseTDEM
from SurveyTDEM import SurveyTDEM1D
from FieldsTDEM import FieldsTDEM
from TDEM_b import ProblemTDEM_b
simpegEM/Tests/test_FDEM_analytics.py
-1
View File
@@ -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)
simpegEM/Tests/test_FieldsObject.py
+25 -18
View File
@@ -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):
simpegEM/Tests/test_TDEM_b_DerivAdjoint.py
+13 -5
View File
@@ -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):