wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-06-27 21:08:35 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #89 from simpeg/Tx2Src

Browse Source 
Tx2 src
This commit is contained in:
Rowan Cockett
2015-05-01 11:15:01 -07:00
parent afe7e15d09 39fdcab9b4
commit fbcd3e3691
7 changed files with 333 additions and 332 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/Problem.py
+47 -46
View File
@@ -7,7 +7,7 @@ class Fields(object):
"""Fancy Field Storage
u[:,'phi'] = phi
print u[tx0,'phi']
print u[src0,'phi']
"""
@@ -43,14 +43,14 @@ class Fields(object):
return "%e MB"%sz
def _storageShape(self, loc):
nTx = self.survey.nTx
nSrc = self.survey.nSrc
nP = {'CC': self.mesh.nC,
'N': self.mesh.nN,
'F': self.mesh.nF,
'E': self.mesh.nE}[loc]
return (nP, nTx)
return (nP, nSrc)
def _initStore(self, name):
if name in self._fields:
@@ -70,17 +70,17 @@ class Fields(object):
return field
def _txIndex(self, txTestList):
if type(txTestList) is slice:
ind = txTestList
def _srcIndex(self, srcTestList):
if type(srcTestList) is slice:
ind = srcTestList
else:
if type(txTestList) is not list:
txTestList = [txTestList]
for txTest in txTestList:
if txTest not in self.survey.txList:
raise KeyError('Invalid Transmitter, not in survey list.')
if type(srcTestList) is not list:
srcTestList = [srcTestList]
for srcTest in srcTestList:
if srcTest not in self.survey.srcList:
raise KeyError('Invalid Source, not in survey list.')
ind = np.in1d(self.survey.txList, txTestList)
ind = np.in1d(self.survey.srcList, srcTestList)
return ind
def _nameIndex(self, name, accessType):
@@ -107,11 +107,11 @@ class Fields(object):
if len(key) == 1:
key += (None,)
assert len(key) == 2, 'must be [Tx, fieldName]'
assert len(key) == 2, 'must be [Src, fieldName]'
txTestList, name = key
srcTestList, name = key
name = self._nameIndex(name, accessType)
ind = self._txIndex(txTestList)
ind = self._srcIndex(srcTestList)
return ind, name
def __setitem__(self, key, value):
@@ -150,16 +150,16 @@ class Fields(object):
# Aliased fields
alias, loc, func = self.aliasFields[name]
txII = np.array(self.survey.txList)[ind]
if isinstance(txII, np.ndarray):
txII = txII.tolist()
if len(txII) == 1:
txII = txII[0]
srcII = np.array(self.survey.srcList)[ind]
if isinstance(srcII, np.ndarray):
srcII = srcII.tolist()
if len(srcII) == 1:
srcII = srcII[0]
if type(func) is str:
assert hasattr(self, func), 'The alias field function is a string, but it does not exist in the Fields class.'
func = getattr(self, func)
out = func(self._fields[alias][:,ind], txII)
out = func(self._fields[alias][:,ind], srcII)
if out.shape[1] == 1:
out = Utils.mkvc(out)
@@ -175,7 +175,7 @@ class TimeFields(Fields):
"""Fancy Field Storage for time domain problems
u[:,'phi', timeInd] = phi
print u[tx0,'phi']
print u[src0,'phi']
"""
@@ -184,9 +184,9 @@ class TimeFields(Fields):
'N': self.mesh.nN,
'F': self.mesh.nF,
'E': self.mesh.nE}[loc]
nTx = self.survey.nTx
nSrc = self.survey.nSrc
nT = self.survey.prob.nT + 1
return (nP, nTx, nT)
return (nP, nSrc, nT)
def _indexAndNameFromKey(self, key, accessType):
if type(key) is not tuple:
@@ -196,66 +196,66 @@ class TimeFields(Fields):
if len(key) == 2:
key += (slice(None,None,None),)
assert len(key) == 3, 'must be [Tx, fieldName, times]'
assert len(key) == 3, 'must be [Src, fieldName, times]'
txTestList, name, timeInd = key
srcTestList, name, timeInd = key
name = self._nameIndex(name, accessType)
txInd = self._txIndex(txTestList)
srcInd = self._srcIndex(srcTestList)
return (txInd, timeInd), name
return (srcInd, timeInd), name
def _correctShape(self, name, ind, deflate=False):
txInd, timeInd = ind
srcInd, timeInd = ind
if name in self.knownFields:
loc = self.knownFields[name]
else:
loc = self.aliasFields[name][1]
nP, total_nTx, total_nT = self._storageShape(loc)
nTx = np.ones(total_nTx, dtype=bool)[txInd].sum()
nP, total_nSrc, total_nT = self._storageShape(loc)
nSrc = np.ones(total_nSrc, dtype=bool)[srcInd].sum()
nT = np.ones(total_nT, dtype=bool)[timeInd].sum()
shape = nP, nTx, nT
shape = nP, nSrc, nT
if deflate:
shape = tuple([s for s in shape if s > 1])
return shape
def _setField(self, field, val, name, ind):
txInd, timeInd = ind
srcInd, timeInd = ind
shape = self._correctShape(name, ind)
if Utils.isScalar(val):
field[:,txInd,timeInd] = val
field[:,srcInd,timeInd] = val
return
if val.size != np.array(shape).prod():
raise ValueError('Incorrect size for data.')
correctShape = field[:,txInd,timeInd].shape
field[:,txInd,timeInd] = val.reshape(correctShape, order='F')
correctShape = field[:,srcInd,timeInd].shape
field[:,srcInd,timeInd] = val.reshape(correctShape, order='F')
def _getField(self, name, ind):
txInd, timeInd = ind
srcInd, timeInd = ind
if name in self._fields:
out = self._fields[name][:,txInd,timeInd]
out = self._fields[name][:,srcInd,timeInd]
else:
# Aliased fields
alias, loc, func = self.aliasFields[name]
if type(func) is str:
assert hasattr(self, func), 'The alias field function is a string, but it does not exist in the Fields class.'
func = getattr(self, func)
pointerFields = self._fields[alias][:,txInd,timeInd]
pointerFields = self._fields[alias][:,srcInd,timeInd]
pointerShape = self._correctShape(alias, ind)
pointerFields = pointerFields.reshape(pointerShape, order='F')
timeII = np.arange(self.survey.prob.nT + 1)[timeInd]
txII = np.array(self.survey.txList)[txInd]
if isinstance(txII, np.ndarray):
txII = txII.tolist()
if len(txII) == 1:
txII = txII[0]
srcII = np.array(self.survey.srcList)[srcInd]
if isinstance(srcII, np.ndarray):
srcII = srcII.tolist()
if len(srcII) == 1:
srcII = srcII[0]
if timeII.size == 1:
pointerShapeDeflated = self._correctShape(alias, ind, deflate=True)
pointerFields = pointerFields.reshape(pointerShapeDeflated, order='F')
out = func(pointerFields, txII, timeII)
out = func(pointerFields, srcII, timeII)
else: #loop over the time steps
nT = pointerShape[2]
out = range(nT)
@@ -263,7 +263,7 @@ class TimeFields(Fields):
fieldI = pointerFields[:,:,i]
if fieldI.ndim == 2 and fieldI.shape[1] == 1:
fieldI = Utils.mkvc(fieldI)
out[i] = func(fieldI, txII, TIND_i)
out[i] = func(fieldI, srcII, TIND_i)
if out[i].ndim == 1:
out[i] = out[i][:,np.newaxis,np.newaxis]
elif out[i].ndim == 2:
@@ -321,6 +321,7 @@ class BaseProblem(object):
self.survey._prob = None
self._survey = None
deleteTheseOnModelUpdate = [] # List of strings, e.g. ['_MeSigma', '_MeSigmaI']
@property
SimPEG/Survey.py
+51 -51
View File
@@ -6,7 +6,7 @@ class BaseRx(object):
locs = None #: Locations (nRx x nDim)
knownRxTypes = None #: Set this to a list of strings to ensure that txType is known
knownRxTypes = None #: Set this to a list of strings to ensure that srcType is known
projGLoc = 'CC' #: Projection grid location, default is CC
@@ -111,37 +111,37 @@ class BaseTimeRx(BaseRx):
return P
class BaseTx(object):
"""SimPEG Transmitter Object"""
class BaseSrc(object):
"""SimPEG Source Object"""
loc = None #: Location [x,y,z]
rxList = None #: SimPEG Receiver List
rxPair = BaseRx
knownTxTypes = None #: Set this to a list of strings to ensure that txType is known
knownSrcTypes = None #: Set this to a list of strings to ensure that srcType is known
def __init__(self, loc, txType, rxList, **kwargs):
def __init__(self, loc, srcType, rxList, **kwargs):
assert type(rxList) is list, 'rxList must be a list'
for rx in rxList:
assert isinstance(rx, self.rxPair), 'rxList must be a %s'%self.rxPair.__name__
assert len(set(rxList)) == len(rxList), 'The rxList must be unique'
self.loc = loc
self.txType = txType
self.srcType = srcType
self.rxList = rxList
Utils.setKwargs(self, **kwargs)
@property
def txType(self):
"""Transmitter Type"""
return getattr(self, '_txType', None)
@txType.setter
def txType(self, value):
known = self.knownTxTypes
def srcType(self):
"""Source Type"""
return getattr(self, '_srcType', None)
@srcType.setter
def srcType(self, value):
known = self.knownSrcTypes
if known is not None:
assert value in known, "txType must be in ['%s']" % ("', '".join(known))
self._txType = value
assert value in known, "srcType must be in ['%s']" % ("', '".join(known))
self._srcType = value
@property
def nD(self):
@@ -155,59 +155,59 @@ class BaseTx(object):
class Data(object):
"""Fancy data storage by Tx and Rx"""
"""Fancy data storage by Src and Rx"""
def __init__(self, survey, v=None):
self.survey = survey
self._dataDict = {}
for tx in self.survey.txList:
self._dataDict[tx] = {}
for src in self.survey.srcList:
self._dataDict[src] = {}
if v is not None:
self.fromvec(v)
def _ensureCorrectKey(self, key):
if type(key) is tuple:
if len(key) is not 2:
raise KeyError('Key must be [Tx, Rx]')
if key[0] not in self.survey.txList:
raise KeyError('Tx Key must be a transmitter in the survey.')
raise KeyError('Key must be [Src, Rx]')
if key[0] not in self.survey.srcList:
raise KeyError('Src Key must be a source in the survey.')
if key[1] not in key[0].rxList:
raise KeyError('Rx Key must be a receiver for the transmitter.')
raise KeyError('Rx Key must be a receiver for the source.')
return key
elif isinstance(key, self.survey.txPair):
if key not in self.survey.txList:
raise KeyError('Key must be a transmitter in the survey.')
elif isinstance(key, self.survey.srcPair):
if key not in self.survey.srcList:
raise KeyError('Key must be a source in the survey.')
return key, None
else:
raise KeyError('Key must be [Tx] or [Tx,Rx]')
raise KeyError('Key must be [Src] or [Src,Rx]')
def __setitem__(self, key, value):
tx, rx = self._ensureCorrectKey(key)
assert rx is not None, 'set data using [Tx, Rx]'
src, rx = self._ensureCorrectKey(key)
assert rx is not None, 'set data using [Src, Rx]'
assert isinstance(value, np.ndarray), 'value must by ndarray'
assert value.size == rx.nD, "value must have the same number of data as the transmitter."
self._dataDict[tx][rx] = Utils.mkvc(value)
assert value.size == rx.nD, "value must have the same number of data as the source."
self._dataDict[src][rx] = Utils.mkvc(value)
def __getitem__(self, key):
tx, rx = self._ensureCorrectKey(key)
src, rx = self._ensureCorrectKey(key)
if rx is not None:
if rx not in self._dataDict[tx]:
if rx not in self._dataDict[src]:
raise Exception('Data for receiver has not yet been set.')
return self._dataDict[tx][rx]
return self._dataDict[src][rx]
return np.concatenate([self[tx,rx] for rx in tx.rxList])
return np.concatenate([self[src,rx] for rx in src.rxList])
def tovec(self):
return np.concatenate([self[tx] for tx in self.survey.txList])
return np.concatenate([self[src] for src in self.survey.srcList])
def fromvec(self, v):
v = Utils.mkvc(v)
assert v.size == self.survey.nD, 'v must have the correct number of data.'
indBot, indTop = 0, 0
for tx in self.survey.txList:
for rx in tx.rxList:
for src in self.survey.srcList:
for rx in src.rxList:
indTop += rx.nD
self[tx, rx] = v[indBot:indTop]
self[src, rx] = v[indBot:indTop]
indBot += rx.nD
@@ -226,19 +226,19 @@ class BaseSurvey(object):
def __init__(self, **kwargs):
Utils.setKwargs(self, **kwargs)
txPair = BaseTx #: Transmitter Pair
srcPair = BaseSrc #: Source Pair
@property
def txList(self):
"""Transmitter List"""
return getattr(self, '_txList', None)
def srcList(self):
"""Source List"""
return getattr(self, '_srcList', None)
@txList.setter
def txList(self, value):
assert type(value) is list, 'txList must be a list'
assert np.all([isinstance(tx, self.txPair) for tx in value]), 'All transmitters must be instances of %s' % self.txPair.__name__
assert len(set(value)) == len(value), 'The txList must be unique'
self._txList = value
@srcList.setter
def srcList(self, value):
assert type(value) is list, 'srcList must be a list'
assert np.all([isinstance(src, self.srcPair) for src in value]), 'All sources must be instances of %s' % self.srcPair.__name__
assert len(set(value)) == len(value), 'The srcList must be unique'
self._srcList = value
@property
def prob(self):
@@ -282,12 +282,12 @@ class BaseSurvey(object):
@property
def vnD(self):
"""Vector number of data"""
return np.array([tx.nD for tx in self.txList])
return np.array([src.nD for src in self.srcList])
@property
def nTx(self):
"""Number of Transmitters"""
return len(self.txList)
def nSrc(self):
"""Number of Sources"""
return len(self.srcList)
@Utils.count
@Utils.requires('prob')
SimPEG/Tests/test_Survey.py
+128 -128
View File
@@ -7,22 +7,22 @@ class DataAndFieldsTest(unittest.TestCase):
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.]
srcLoc = np.r_[0,0,0.]
rxList0 = Survey.BaseRx(XYZ, 'exi')
Tx0 = Survey.BaseTx(txLoc, 'VMD', [rxList0])
Src0 = Survey.BaseSrc(srcLoc, 'VMD', [rxList0])
rxList1 = Survey.BaseRx(XYZ, 'bxi')
Tx1 = Survey.BaseTx(txLoc, 'VMD', [rxList1])
Src1 = Survey.BaseSrc(srcLoc, 'VMD', [rxList1])
rxList2 = Survey.BaseRx(XYZ, 'bxi')
Tx2 = Survey.BaseTx(txLoc, 'VMD', [rxList2])
Src2 = Survey.BaseSrc(srcLoc, 'VMD', [rxList2])
rxList3 = Survey.BaseRx(XYZ, 'bxi')
Tx3 = Survey.BaseTx(txLoc, 'VMD', [rxList3])
Tx4 = Survey.BaseTx(txLoc, 'VMD', [rxList0, rxList1, rxList2, rxList3])
txList = [Tx0,Tx1,Tx2,Tx3,Tx4]
survey = Survey.BaseSurvey(txList=txList)
Src3 = Survey.BaseSrc(srcLoc, 'VMD', [rxList3])
Src4 = Survey.BaseSrc(srcLoc, 'VMD', [rxList0, rxList1, rxList2, rxList3])
srcList = [Src0,Src1,Src2,Src3,Src4]
survey = Survey.BaseSurvey(srcList=srcList)
self.D = Survey.Data(survey)
self.F = Problem.Fields(mesh, survey, knownFields={'phi':'CC','e':'E','b':'F'}, dtype={"phi":float,"e":complex,"b":complex})
self.Tx0 = Tx0
self.Tx1 = Tx1
self.Src0 = Src0
self.Src1 = Src1
self.mesh = mesh
self.XYZ = XYZ
@@ -33,12 +33,12 @@ class DataAndFieldsTest(unittest.TestCase):
def test_data(self):
V = []
for tx in self.D.survey.txList:
for rx in tx.rxList:
for src in self.D.survey.srcList:
for rx in src.rxList:
v = np.random.rand(rx.nD)
V += [v]
self.D[tx, rx] = v
self.assertTrue(np.all(v == self.D[tx, rx]))
self.D[src, rx] = v
self.assertTrue(np.all(v == self.D[src, rx]))
V = np.concatenate(V)
self.assertTrue(np.all(V == Utils.mkvc(self.D)))
@@ -47,25 +47,25 @@ class DataAndFieldsTest(unittest.TestCase):
def test_contains(self):
F = self.F
nTx = F.survey.nTx
nSrc = F.survey.nSrc
self.assertTrue('b' not in F)
self.assertTrue('e' not in F)
e = np.random.rand(F.mesh.nE, nTx)
e = np.random.rand(F.mesh.nE, nSrc)
F[:, 'e'] = e
self.assertTrue('b' not in F)
self.assertTrue('e' in F)
def test_uniqueTxs(self):
txs = self.D.survey.txList
txs += [txs[0]]
self.assertRaises(AssertionError, Survey.BaseSurvey, txList=txs)
def test_uniqueSrcs(self):
srcs = self.D.survey.srcList
srcs += [srcs[0]]
self.assertRaises(AssertionError, Survey.BaseSurvey, srcList=srcs)
def test_SetGet(self):
F = self.F
nTx = F.survey.nTx
e = np.random.rand(F.mesh.nE, nTx) + np.random.rand(F.mesh.nE, nTx)*1j
nSrc = F.survey.nSrc
e = np.random.rand(F.mesh.nE, nSrc) + np.random.rand(F.mesh.nE, nSrc)*1j
F[:, 'e'] = e
b = np.random.rand(F.mesh.nF, nTx) + np.random.rand(F.mesh.nF, nTx)*1j
b = np.random.rand(F.mesh.nF, nSrc) + np.random.rand(F.mesh.nF, nSrc)*1j
F[:, 'b'] = b
self.assertTrue(np.all(F[:, 'e'] == e))
@@ -79,31 +79,31 @@ class DataAndFieldsTest(unittest.TestCase):
self.assertTrue(np.all(F[:, 'b'] == b*0))
b = np.random.rand(F.mesh.nF,1)
F[self.Tx0, 'b'] = b
self.assertTrue(np.all(F[self.Tx0, 'b'] == Utils.mkvc(b)))
F[self.Src0, 'b'] = b
self.assertTrue(np.all(F[self.Src0, 'b'] == Utils.mkvc(b)))
b = np.random.rand(F.mesh.nF)
F[self.Tx0, 'b'] = b
self.assertTrue(np.all(F[self.Tx0, 'b'] == b))
F[self.Src0, 'b'] = b
self.assertTrue(np.all(F[self.Src0, 'b'] == b))
phi = np.random.rand(F.mesh.nC,2)
F[[self.Tx0,self.Tx1], 'phi'] = phi
self.assertTrue(np.all(F[[self.Tx0,self.Tx1], 'phi'] == phi))
F[[self.Src0,self.Src1], 'phi'] = phi
self.assertTrue(np.all(F[[self.Src0,self.Src1], 'phi'] == phi))
fdict = F[:,:]
self.assertTrue(type(fdict) is dict)
self.assertTrue(sorted([k for k in fdict]) == ['b','e','phi'])
b = np.random.rand(F.mesh.nF, 2)
F[[self.Tx0, self.Tx1],'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]))
self.assertTrue(np.all(F[self.Tx1,'b'] == b[:,1]))
F[[self.Src0, self.Src1],'b'] = b
self.assertTrue(F[self.Src0]['b'].shape == (F.mesh.nF,))
self.assertTrue(F[self.Src0,'b'].shape == (F.mesh.nF,))
self.assertTrue(np.all(F[self.Src0,'b'] == b[:,0]))
self.assertTrue(np.all(F[self.Src1,'b'] == b[:,1]))
def test_assertions(self):
freq = [self.Tx0, self.Tx1]
bWrongSize = np.random.rand(self.F.mesh.nE, self.F.survey.nTx)
freq = [self.Src0, self.Src1]
bWrongSize = np.random.rand(self.F.mesh.nE, self.F.survey.nSrc)
def fun(): self.F[freq, 'b'] = bWrongSize
self.assertRaises(ValueError, fun)
def fun(): self.F[-999.]
@@ -120,69 +120,69 @@ class FieldsTest_Alias(unittest.TestCase):
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.]
srcLoc = np.r_[0,0,0.]
rxList0 = Survey.BaseRx(XYZ, 'exi')
Tx0 = Survey.BaseTx(txLoc, 'VMD', [rxList0])
Src0 = Survey.BaseSrc(srcLoc, 'VMD', [rxList0])
rxList1 = Survey.BaseRx(XYZ, 'bxi')
Tx1 = Survey.BaseTx(txLoc, 'VMD', [rxList1])
Src1 = Survey.BaseSrc(srcLoc, 'VMD', [rxList1])
rxList2 = Survey.BaseRx(XYZ, 'bxi')
Tx2 = Survey.BaseTx(txLoc, 'VMD', [rxList2])
Src2 = Survey.BaseSrc(srcLoc, 'VMD', [rxList2])
rxList3 = Survey.BaseRx(XYZ, 'bxi')
Tx3 = Survey.BaseTx(txLoc, 'VMD', [rxList3])
Tx4 = Survey.BaseTx(txLoc, 'VMD', [rxList0, rxList1, rxList2, rxList3])
txList = [Tx0,Tx1,Tx2,Tx3,Tx4]
survey = Survey.BaseSurvey(txList=txList)
Src3 = Survey.BaseSrc(srcLoc, 'VMD', [rxList3])
Src4 = Survey.BaseSrc(srcLoc, 'VMD', [rxList0, rxList1, rxList2, rxList3])
srcList = [Src0,Src1,Src2,Src3,Src4]
survey = Survey.BaseSurvey(srcList=srcList)
self.D = Survey.Data(survey)
self.F = Problem.Fields(mesh, survey, knownFields={'e':'E'}, aliasFields={'b':['e','F',(lambda e, ind: self.F.mesh.edgeCurl * e)]})
self.Tx0 = Tx0
self.Tx1 = Tx1
self.Src0 = Src0
self.Src1 = Src1
self.mesh = mesh
self.XYZ = XYZ
def test_contains(self):
F = self.F
nTx = F.survey.nTx
nSrc = F.survey.nSrc
self.assertTrue('b' not in F)
self.assertTrue('e' not in F)
e = np.random.rand(F.mesh.nE, nTx)
e = np.random.rand(F.mesh.nE, nSrc)
F[:, 'e'] = e
self.assertTrue('b' in F)
self.assertTrue('e' in F)
def test_simpleAlias(self):
F = self.F
nTx = F.survey.nTx
e = np.random.rand(F.mesh.nE, nTx)
nSrc = F.survey.nSrc
e = np.random.rand(F.mesh.nE, nSrc)
F[:, 'e'] = e
self.assertTrue(np.all(F[:, 'b'] == F.mesh.edgeCurl * e ))
e = np.random.rand(F.mesh.nE,1)
F[self.Tx0, 'e'] = e
self.assertTrue(np.all(F[self.Tx0, 'b'] == F.mesh.edgeCurl * Utils.mkvc(e)))
F[self.Src0, 'e'] = e
self.assertTrue(np.all(F[self.Src0, 'b'] == F.mesh.edgeCurl * Utils.mkvc(e)))
def f():
F[self.Tx0, 'b'] = F[self.Tx0, 'b']
F[self.Src0, 'b'] = F[self.Src0, 'b']
self.assertRaises(KeyError, f) # can't set a alias attr.
def test_aliasFunction(self):
def alias(e, ind):
self.assertTrue(ind is self.Tx0)
self.assertTrue(ind is self.Src0)
return self.F.mesh.edgeCurl * e
F = Problem.Fields(self.F.mesh, self.F.survey, knownFields={'e':'E'}, aliasFields={'b':['e','F',alias]})
e = np.random.rand(F.mesh.nE,1)
F[self.Tx0, 'e'] = e
F[self.Tx0, 'b']
F[self.Src0, 'e'] = e
F[self.Src0, 'b']
def alias(e, ind):
self.assertTrue(type(ind) is list)
self.assertTrue(ind[0] is self.Tx0)
self.assertTrue(ind[1] is self.Tx1)
self.assertTrue(ind[0] is self.Src0)
self.assertTrue(ind[1] is self.Src1)
return self.F.mesh.edgeCurl * e
F = Problem.Fields(self.F.mesh, self.F.survey, knownFields={'e':'E'}, aliasFields={'b':['e','F',alias]})
e = np.random.rand(F.mesh.nE,2)
F[[self.Tx0, self.Tx1], 'e'] = e
F[[self.Tx0, self.Tx1], 'b']
F[[self.Src0, self.Src1], 'e'] = e
F[[self.Src0, self.Src1], 'b']
class FieldsTest_Time(unittest.TestCase):
@@ -191,34 +191,34 @@ class FieldsTest_Time(unittest.TestCase):
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.]
srcLoc = np.r_[0,0,0.]
rxList0 = Survey.BaseRx(XYZ, 'exi')
Tx0 = Survey.BaseTx(txLoc, 'VMD', [rxList0])
Src0 = Survey.BaseSrc(srcLoc, 'VMD', [rxList0])
rxList1 = Survey.BaseRx(XYZ, 'bxi')
Tx1 = Survey.BaseTx(txLoc, 'VMD', [rxList1])
Src1 = Survey.BaseSrc(srcLoc, 'VMD', [rxList1])
rxList2 = Survey.BaseRx(XYZ, 'bxi')
Tx2 = Survey.BaseTx(txLoc, 'VMD', [rxList2])
Src2 = Survey.BaseSrc(srcLoc, 'VMD', [rxList2])
rxList3 = Survey.BaseRx(XYZ, 'bxi')
Tx3 = Survey.BaseTx(txLoc, 'VMD', [rxList3])
Tx4 = Survey.BaseTx(txLoc, 'VMD', [rxList0, rxList1, rxList2, rxList3])
txList = [Tx0,Tx1,Tx2,Tx3,Tx4]
survey = Survey.BaseSurvey(txList=txList)
Src3 = Survey.BaseSrc(srcLoc, 'VMD', [rxList3])
Src4 = Survey.BaseSrc(srcLoc, 'VMD', [rxList0, rxList1, rxList2, rxList3])
srcList = [Src0,Src1,Src2,Src3,Src4]
survey = Survey.BaseSurvey(srcList=srcList)
prob = Problem.BaseTimeProblem(mesh, timeSteps=[(10.,3), (20.,2)])
survey.pair(prob)
self.F = Problem.TimeFields(mesh, survey, knownFields={'phi':'CC','e':'E','b':'F'})
self.Tx0 = Tx0
self.Tx1 = Tx1
self.Src0 = Src0
self.Src1 = Src1
self.mesh = mesh
self.XYZ = XYZ
def test_contains(self):
F = self.F
nTx = F.survey.nTx
nSrc = F.survey.nSrc
nT = F.survey.prob.nT + 1
self.assertTrue('b' not in F)
self.assertTrue('e' not in F)
self.assertTrue('phi' not in F)
e = np.random.rand(F.mesh.nE, nTx, nT)
e = np.random.rand(F.mesh.nE, nSrc, nT)
F[:, 'e', :] = e
self.assertTrue('e' in F)
self.assertTrue('b' not in F)
@@ -226,11 +226,11 @@ class FieldsTest_Time(unittest.TestCase):
def test_SetGet(self):
F = self.F
nTx = F.survey.nTx
nSrc = F.survey.nSrc
nT = F.survey.prob.nT + 1
e = np.random.rand(F.mesh.nE, nTx, nT)
e = np.random.rand(F.mesh.nE, nSrc, nT)
F[:, 'e'] = e
b = np.random.rand(F.mesh.nF, nTx, nT)
b = np.random.rand(F.mesh.nF, nSrc, nT)
F[:, 'b'] = b
self.assertTrue(np.all(F[:, 'e'] == e))
@@ -244,39 +244,39 @@ class FieldsTest_Time(unittest.TestCase):
self.assertTrue(np.all(F[:, 'b'] == b*0))
b = np.random.rand(F.mesh.nF,1,nT)
F[self.Tx0, 'b'] = b
self.assertTrue(np.all(F[self.Tx0, 'b'] == b[:,0,:]))
F[self.Src0, 'b'] = b
self.assertTrue(np.all(F[self.Src0, 'b'] == b[:,0,:]))
b = np.random.rand(F.mesh.nF,1,nT)
F[self.Tx0, 'b', 0] = b[:,:,0]
self.assertTrue(np.all(F[self.Tx0, 'b', 0] == b[:,0,0]))
F[self.Src0, 'b', 0] = b[:,:,0]
self.assertTrue(np.all(F[self.Src0, 'b', 0] == b[:,0,0]))
phi = np.random.rand(F.mesh.nC,2,nT)
F[[self.Tx0,self.Tx1], 'phi'] = phi
self.assertTrue(np.all(F[[self.Tx0,self.Tx1], 'phi'] == phi))
F[[self.Src0,self.Src1], 'phi'] = phi
self.assertTrue(np.all(F[[self.Src0,self.Src1], 'phi'] == phi))
fdict = F[:]
self.assertTrue(type(fdict) is dict)
self.assertTrue(sorted([k for k in fdict]) == ['b','e','phi'])
b = np.random.rand(F.mesh.nF, 2, nT)
F[[self.Tx0, self.Tx1],'b'] = b
self.assertTrue(F[self.Tx0]['b'].shape == (F.mesh.nF,nT))
self.assertTrue(F[self.Tx0,'b'].shape == (F.mesh.nF,nT))
self.assertTrue(np.all(F[self.Tx0,'b'] == b[:,0,:]))
self.assertTrue(np.all(F[self.Tx1,'b'] == b[:,1,:]))
self.assertTrue(np.all(F[self.Tx0,'b',1] == b[:,0,1]))
self.assertTrue(np.all(F[self.Tx1,'b',1] == b[:,1,1]))
self.assertTrue(np.all(F[self.Tx0,'b',4] == b[:,0,4]))
self.assertTrue(np.all(F[self.Tx1,'b',4] == b[:,1,4]))
F[[self.Src0, self.Src1],'b'] = b
self.assertTrue(F[self.Src0]['b'].shape == (F.mesh.nF,nT))
self.assertTrue(F[self.Src0,'b'].shape == (F.mesh.nF,nT))
self.assertTrue(np.all(F[self.Src0,'b'] == b[:,0,:]))
self.assertTrue(np.all(F[self.Src1,'b'] == b[:,1,:]))
self.assertTrue(np.all(F[self.Src0,'b',1] == b[:,0,1]))
self.assertTrue(np.all(F[self.Src1,'b',1] == b[:,1,1]))
self.assertTrue(np.all(F[self.Src0,'b',4] == b[:,0,4]))
self.assertTrue(np.all(F[self.Src1,'b',4] == b[:,1,4]))
b = np.random.rand(F.mesh.nF, 2, nT)
F[[self.Tx0, self.Tx1],'b', 0] = b[:,:,0]
F[[self.Src0, self.Src1],'b', 0] = b[:,:,0]
def test_assertions(self):
freq = [self.Tx0, self.Tx1]
bWrongSize = np.random.rand(self.F.mesh.nE, self.F.survey.nTx)
freq = [self.Src0, self.Src1]
bWrongSize = np.random.rand(self.F.mesh.nE, self.F.survey.nSrc)
def fun(): self.F[freq, 'b'] = bWrongSize
self.assertRaises(ValueError, fun)
def fun(): self.F[-999.]
@@ -293,35 +293,35 @@ class FieldsTest_Time_Aliased(unittest.TestCase):
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.]
srcLoc = np.r_[0,0,0.]
rxList0 = Survey.BaseRx(XYZ, 'exi')
Tx0 = Survey.BaseTx(txLoc, 'VMD', [rxList0])
Src0 = Survey.BaseSrc(srcLoc, 'VMD', [rxList0])
rxList1 = Survey.BaseRx(XYZ, 'bxi')
Tx1 = Survey.BaseTx(txLoc, 'VMD', [rxList1])
Src1 = Survey.BaseSrc(srcLoc, 'VMD', [rxList1])
rxList2 = Survey.BaseRx(XYZ, 'bxi')
Tx2 = Survey.BaseTx(txLoc, 'VMD', [rxList2])
Src2 = Survey.BaseSrc(srcLoc, 'VMD', [rxList2])
rxList3 = Survey.BaseRx(XYZ, 'bxi')
Tx3 = Survey.BaseTx(txLoc, 'VMD', [rxList3])
Tx4 = Survey.BaseTx(txLoc, 'VMD', [rxList0, rxList1, rxList2, rxList3])
txList = [Tx0,Tx1,Tx2,Tx3,Tx4]
survey = Survey.BaseSurvey(txList=txList)
Src3 = Survey.BaseSrc(srcLoc, 'VMD', [rxList3])
Src4 = Survey.BaseSrc(srcLoc, 'VMD', [rxList0, rxList1, rxList2, rxList3])
srcList = [Src0,Src1,Src2,Src3,Src4]
survey = Survey.BaseSurvey(srcList=srcList)
prob = Problem.BaseTimeProblem(mesh, timeSteps=[(10.,3), (20.,2)])
survey.pair(prob)
def alias(b, txInd, timeInd):
def alias(b, srcInd, timeInd):
return self.F.mesh.edgeCurl.T * b + timeInd
self.F = Problem.TimeFields(mesh, survey, knownFields={'b':'F'}, aliasFields={'e':['b','E',alias]})
self.Tx0 = Tx0
self.Tx1 = Tx1
self.Src0 = Src0
self.Src1 = Src1
self.mesh = mesh
self.XYZ = XYZ
def test_contains(self):
F = self.F
nTx = F.survey.nTx
nSrc = F.survey.nSrc
nT = F.survey.prob.nT + 1
self.assertTrue('b' not in F)
self.assertTrue('e' not in F)
b = np.random.rand(F.mesh.nF, nTx, nT)
b = np.random.rand(F.mesh.nF, nSrc, nT)
F[:, 'b', :] = b
self.assertTrue('e' in F)
self.assertTrue('b' in F)
@@ -329,9 +329,9 @@ class FieldsTest_Time_Aliased(unittest.TestCase):
def test_simpleAlias(self):
F = self.F
nTx = F.survey.nTx
nSrc = F.survey.nSrc
nT = F.survey.prob.nT + 1
b = np.random.rand(F.mesh.nF, nTx, nT)
b = np.random.rand(F.mesh.nF, nSrc, nT)
F[:, 'b', :] = b
self.assertTrue(np.all(F[:, 'e', 0] == F.mesh.edgeCurl.T * b[:,:,0] ))
@@ -341,57 +341,57 @@ class FieldsTest_Time_Aliased(unittest.TestCase):
e[i] = e[i][:,:,np.newaxis]
e = np.concatenate(e, axis=2)
self.assertTrue(np.all(F[:, 'e', :] == e ))
self.assertTrue(np.all(F[self.Tx0, 'e', :] == e[:,0,:] ))
self.assertTrue(np.all(F[self.Tx1, 'e', :] == e[:,1,:] ))
self.assertTrue(np.all(F[self.Src0, 'e', :] == e[:,0,:] ))
self.assertTrue(np.all(F[self.Src1, 'e', :] == e[:,1,:] ))
for t in range(nT):
self.assertTrue(np.all(F[self.Tx1, 'e', t] == e[:,1,t] ))
self.assertTrue(np.all(F[self.Src1, 'e', t] == e[:,1,t] ))
b = np.random.rand(F.mesh.nF,nT)
F[self.Tx0, 'b',:] = b
F[self.Src0, 'b',:] = b
Cb = F.mesh.edgeCurl.T * b
for i in range(Cb.shape[1]):
Cb[:,i] += i
self.assertTrue(np.all(F[self.Tx0, 'e',:] == Cb))
self.assertTrue(np.all(F[self.Src0, 'e',:] == Cb))
def f():
F[self.Tx0, 'e'] = F[self.Tx0, 'e']
F[self.Src0, 'e'] = F[self.Src0, 'e']
self.assertRaises(KeyError, f) # can't set a alias attr.
def test_aliasFunction(self):
nT = self.F.survey.prob.nT + 1
count = [0]
def alias(e, txInd, timeInd):
def alias(e, srcInd, timeInd):
count[0] += 1
self.assertTrue(txInd is self.Tx0)
self.assertTrue(srcInd is self.Src0)
return self.F.mesh.edgeCurl * e
F = Problem.TimeFields(self.F.mesh, self.F.survey, knownFields={'e':'E'}, aliasFields={'b':['e','F',alias]})
e = np.random.rand(F.mesh.nE,1,nT)
F[self.Tx0, 'e', :] = e
F[self.Tx0, 'b', :]
F[self.Src0, 'e', :] = e
F[self.Src0, 'b', :]
self.assertTrue(count[0] == nT) # ensure that this is called for every time separately.
e = np.random.rand(F.mesh.nE,1,1)
F[self.Tx0, 'e', 1] = e
F[self.Src0, 'e', 1] = e
count[0] = 0
F[self.Tx0, 'b', 1]
F[self.Src0, 'b', 1]
self.assertTrue(count[0] == 1) # ensure that this is called only once.
def alias(e, txInd, timeInd):
def alias(e, srcInd, timeInd):
count[0] += 1
self.assertTrue(type(txInd) is list)
self.assertTrue(txInd[0] is self.Tx0)
self.assertTrue(txInd[1] is self.Tx1)
self.assertTrue(type(srcInd) is list)
self.assertTrue(srcInd[0] is self.Src0)
self.assertTrue(srcInd[1] is self.Src1)
return self.F.mesh.edgeCurl * e
F = Problem.TimeFields(self.F.mesh, self.F.survey, knownFields={'e':'E'}, aliasFields={'b':['e','F',alias]})
e = np.random.rand(F.mesh.nE,2, nT)
F[[self.Tx0, self.Tx1], 'e', :] = e
F[[self.Src0, self.Src1], 'e', :] = e
count[0] = 0
F[[self.Tx0, self.Tx1], 'b', :]
F[[self.Src0, self.Src1], 'b', :]
self.assertTrue(count[0] == nT) # ensure that this is called for every time separately.
e = np.random.rand(F.mesh.nE,2, 1)
F[[self.Tx0, self.Tx1], 'e', 1] = e
F[[self.Src0, self.Src1], 'e', 1] = e
count[0] = 0
F[[self.Tx0, self.Tx1], 'b', 1]
F[[self.Src0, self.Src1], 'b', 1]
self.assertTrue(count[0] == 1) # ensure that this is called only once.
SimPEG/Tests/test_boundaryPoisson.py
+63 -63
View File
@@ -20,9 +20,9 @@ class Test1D_InhomogeneousDirichlet(OrderTest):
j_fun = lambda x: -np.pi*np.sin(np.pi*x)
q_fun = lambda x: -(np.pi**2)*np.cos(np.pi*x)
xc_anal = phi(self.M.gridCC)
q_anal = q_fun(self.M.gridCC)
j_anal = j_fun(self.M.gridFx)
xc_ana = phi(self.M.gridCC)
q_ana = q_fun(self.M.gridCC)
j_ana = j_fun(self.M.gridFx)
#TODO: Check where our boundary conditions are CCx or Nx
# vec = self.M.vectorNx
@@ -38,32 +38,32 @@ class Test1D_InhomogeneousDirichlet(OrderTest):
V = Utils.sdiag(self.M.vol)
G = -Pin.T*Pin*self.M.faceDiv.T * V
D = self.M.faceDiv
j = McI*(G*xc_anal + P*phi_bc)
j = McI*(G*xc_ana + P*phi_bc)
q = V*D*Pin.T*Pin*j + V*D*Pout.T*j_bc
# Rearrange if we know q to solve for x
A = V*D*Pin.T*Pin*McI*G
rhs = V*q_anal - V*D*Pin.T*Pin*McI*P*phi_bc - V*D*Pout.T*j_bc
rhs = V*q_ana - V*D*Pin.T*Pin*McI*P*phi_bc - V*D*Pout.T*j_bc
# A = D*McI*G
# rhs = q_anal - D*McI*P*phi_bc
# rhs = q_ana - D*McI*P*phi_bc
if self.myTest == 'j':
err = np.linalg.norm((j-j_anal), np.inf)
err = np.linalg.norm((j-j_ana), np.inf)
elif self.myTest == 'q':
err = np.linalg.norm((q-V*q_anal), np.inf)
err = np.linalg.norm((q-V*q_ana), np.inf)
elif self.myTest == 'xc':
#TODO: fix the null space
solver = SolverCG(A, maxiter=1000)
xc = solver * (rhs)
print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
err = np.linalg.norm((xc-xc_anal), np.inf)
err = np.linalg.norm((xc-xc_ana), np.inf)
elif self.myTest == 'xcJ':
#TODO: fix the null space
xc = Solver(A) * (rhs)
print np.linalg.norm(Utils.mkvc(A*xc) - rhs)
j = McI*(G*xc + P*phi_bc)
err = np.linalg.norm((j-j_anal), np.inf)
err = np.linalg.norm((j-j_ana), np.inf)
return err
@@ -102,11 +102,11 @@ class Test2D_InhomogeneousDirichlet(OrderTest):
j_funY = lambda x: -np.pi*np.cos(np.pi*x[:,0])*np.sin(np.pi*x[:,1])
q_fun = lambda x: -2*(np.pi**2)*phi(x)
xc_anal = phi(self.M.gridCC)
q_anal = q_fun(self.M.gridCC)
jX_anal = j_funX(self.M.gridFx)
jY_anal = j_funY(self.M.gridFy)
j_anal = np.r_[jX_anal,jY_anal]
xc_ana = phi(self.M.gridCC)
q_ana = q_fun(self.M.gridCC)
jX_ana = j_funX(self.M.gridFx)
jY_ana = j_funY(self.M.gridFy)
j_ana = np.r_[jX_ana,jY_ana]
#TODO: Check where our boundary conditions are CCx or Nx
# fxm,fxp,fym,fyp = self.M.faceBoundaryInd
@@ -126,26 +126,26 @@ class Test2D_InhomogeneousDirichlet(OrderTest):
McI = Utils.sdInv(self.M.getFaceInnerProduct())
G = -self.M.faceDiv.T * Utils.sdiag(self.M.vol)
D = self.M.faceDiv
j = McI*(G*xc_anal + P*bc)
j = McI*(G*xc_ana + P*bc)
q = D*j
# self.M.plotImage(j, 'FxFy', showIt=True)
# Rearrange if we know q to solve for x
A = D*McI*G
rhs = q_anal - D*McI*P*bc
rhs = q_ana - D*McI*P*bc
if self.myTest == 'j':
err = np.linalg.norm((j-j_anal), np.inf)
err = np.linalg.norm((j-j_ana), np.inf)
elif self.myTest == 'q':
err = np.linalg.norm((q-q_anal), np.inf)
err = np.linalg.norm((q-q_ana), np.inf)
elif self.myTest == 'xc':
xc = Solver(A) * (rhs)
err = np.linalg.norm((xc-xc_anal), np.inf)
err = np.linalg.norm((xc-xc_ana), np.inf)
elif self.myTest == 'xcJ':
xc = Solver(A) * (rhs)
j = McI*(G*xc + P*bc)
err = np.linalg.norm((j-j_anal), np.inf)
err = np.linalg.norm((j-j_ana), np.inf)
return err
@@ -182,9 +182,9 @@ class Test1D_InhomogeneousNeumann(OrderTest):
j_fun = lambda x: np.pi*np.cos(np.pi*x)
q_fun = lambda x: -(np.pi**2)*np.sin(np.pi*x)
xc_anal = phi(self.M.gridCC)
q_anal = q_fun(self.M.gridCC)
j_anal = j_fun(self.M.gridFx)
xc_ana = phi(self.M.gridCC)
q_ana = q_fun(self.M.gridCC)
j_ana = j_fun(self.M.gridFx)
#TODO: Check where our boundary conditions are CCx or Nx
vecN = self.M.vectorNx
@@ -200,24 +200,24 @@ class Test1D_InhomogeneousNeumann(OrderTest):
V = Utils.sdiag(self.M.vol)
G = -Pin.T*Pin*self.M.faceDiv.T * V
D = self.M.faceDiv
j = McI*(G*xc_anal + P*phi_bc)
j = McI*(G*xc_ana + P*phi_bc)
q = V*D*Pin.T*Pin*j + V*D*Pout.T*j_bc
# Rearrange if we know q to solve for x
A = V*D*Pin.T*Pin*McI*G
rhs = V*q_anal - V*D*Pin.T*Pin*McI*P*phi_bc - V*D*Pout.T*j_bc
rhs = V*q_ana - V*D*Pin.T*Pin*McI*P*phi_bc - V*D*Pout.T*j_bc
# A = D*McI*G
# rhs = q_anal - D*McI*P*phi_bc
# rhs = q_ana - D*McI*P*phi_bc
if self.myTest == 'j':
err = np.linalg.norm((Pin*j-Pin*j_anal), np.inf)
err = np.linalg.norm((Pin*j-Pin*j_ana), np.inf)
elif self.myTest == 'q':
err = np.linalg.norm((q-V*q_anal), np.inf)
err = np.linalg.norm((q-V*q_ana), np.inf)
elif self.myTest == 'xc':
#TODO: fix the null space
xc, info = sp.linalg.minres(A, rhs, tol = 1e-6)
err = np.linalg.norm((xc-xc_anal), np.inf)
err = np.linalg.norm((xc-xc_ana), np.inf)
if info > 0:
print 'Solve does not work well'
print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
@@ -225,7 +225,7 @@ class Test1D_InhomogeneousNeumann(OrderTest):
#TODO: fix the null space
xc, info = sp.linalg.minres(A, rhs, tol = 1e-6)
j = McI*(G*xc + P*phi_bc)
err = np.linalg.norm((Pin*j-Pin*j_anal), np.inf)
err = np.linalg.norm((Pin*j-Pin*j_ana), np.inf)
if info > 0:
print 'Solve does not work well'
print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
@@ -261,11 +261,11 @@ class Test2D_InhomogeneousNeumann(OrderTest):
j_funY = lambda x: np.pi*np.sin(np.pi*x[:,0])*np.cos(np.pi*x[:,1])
q_fun = lambda x: -2*(np.pi**2)*phi(x)
xc_anal = phi(self.M.gridCC)
q_anal = q_fun(self.M.gridCC)
jX_anal = j_funX(self.M.gridFx)
jY_anal = j_funY(self.M.gridFy)
j_anal = np.r_[jX_anal,jY_anal]
xc_ana = phi(self.M.gridCC)
q_ana = q_fun(self.M.gridCC)
jX_ana = j_funX(self.M.gridFx)
jY_ana = j_funY(self.M.gridFy)
j_ana = np.r_[jX_ana,jY_ana]
#TODO: Check where our boundary conditions are CCx or Nx
@@ -290,21 +290,21 @@ class Test2D_InhomogeneousNeumann(OrderTest):
V = Utils.sdiag(self.M.vol)
G = -Pin.T*Pin*self.M.faceDiv.T * V
D = self.M.faceDiv
j = McI*(G*xc_anal + P*phi_bc)
j = McI*(G*xc_ana + P*phi_bc)
q = V*D*Pin.T*Pin*j + V*D*Pout.T*j_bc
# Rearrange if we know q to solve for x
A = V*D*Pin.T*Pin*McI*G
rhs = V*q_anal - V*D*Pin.T*Pin*McI*P*phi_bc - V*D*Pout.T*j_bc
rhs = V*q_ana - V*D*Pin.T*Pin*McI*P*phi_bc - V*D*Pout.T*j_bc
if self.myTest == 'j':
err = np.linalg.norm((Pin*j-Pin*j_anal), np.inf)
err = np.linalg.norm((Pin*j-Pin*j_ana), np.inf)
elif self.myTest == 'q':
err = np.linalg.norm((q-V*q_anal), np.inf)
err = np.linalg.norm((q-V*q_ana), np.inf)
elif self.myTest == 'xc':
#TODO: fix the null space
xc, info = sp.linalg.minres(A, rhs, tol = 1e-6)
err = np.linalg.norm((xc-xc_anal), np.inf)
err = np.linalg.norm((xc-xc_ana), np.inf)
if info > 0:
print 'Solve does not work well'
print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
@@ -312,7 +312,7 @@ class Test2D_InhomogeneousNeumann(OrderTest):
#TODO: fix the null space
xc, info = sp.linalg.minres(A, rhs, tol = 1e-6)
j = McI*(G*xc + P*phi_bc)
err = np.linalg.norm((Pin*j-Pin*j_anal), np.inf)
err = np.linalg.norm((Pin*j-Pin*j_ana), np.inf)
if info > 0:
print 'Solve does not work well'
print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
@@ -346,9 +346,9 @@ class Test1D_InhomogeneousMixed(OrderTest):
j_fun = lambda x: -0.5*np.pi*np.sin(0.5*np.pi*x)
q_fun = lambda x: -0.25*(np.pi**2)*np.cos(0.5*np.pi*x)
xc_anal = phi(self.M.gridCC)
q_anal = q_fun(self.M.gridCC)
j_anal = j_fun(self.M.gridFx)
xc_ana = phi(self.M.gridCC)
q_ana = q_fun(self.M.gridCC)
j_ana = j_fun(self.M.gridFx)
#TODO: Check where our boundary conditions are CCx or Nx
vecN = self.M.vectorNx
@@ -364,24 +364,24 @@ class Test1D_InhomogeneousMixed(OrderTest):
V = Utils.sdiag(self.M.vol)
G = -Pin.T*Pin*self.M.faceDiv.T * V
D = self.M.faceDiv
j = McI*(G*xc_anal + P*phi_bc)
j = McI*(G*xc_ana + P*phi_bc)
q = V*D*Pin.T*Pin*j + V*D*Pout.T*j_bc
# Rearrange if we know q to solve for x
A = V*D*Pin.T*Pin*McI*G
rhs = V*q_anal - V*D*Pin.T*Pin*McI*P*phi_bc - V*D*Pout.T*j_bc
rhs = V*q_ana - V*D*Pin.T*Pin*McI*P*phi_bc - V*D*Pout.T*j_bc
# A = D*McI*G
# rhs = q_anal - D*McI*P*phi_bc
# rhs = q_ana - D*McI*P*phi_bc
if self.myTest == 'j':
err = np.linalg.norm((Pin*j-Pin*j_anal), np.inf)
err = np.linalg.norm((Pin*j-Pin*j_ana), np.inf)
elif self.myTest == 'q':
err = np.linalg.norm((q-V*q_anal), np.inf)
err = np.linalg.norm((q-V*q_ana), np.inf)
elif self.myTest == 'xc':
#TODO: fix the null space
xc, info = sp.linalg.minres(A, rhs, tol = 1e-6)
err = np.linalg.norm((xc-xc_anal), np.inf)
err = np.linalg.norm((xc-xc_ana), np.inf)
if info > 0:
print 'Solve does not work well'
print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
@@ -389,7 +389,7 @@ class Test1D_InhomogeneousMixed(OrderTest):
#TODO: fix the null space
xc, info = sp.linalg.minres(A, rhs, tol = 1e-6)
j = McI*(G*xc + P*phi_bc)
err = np.linalg.norm((Pin*j-Pin*j_anal), np.inf)
err = np.linalg.norm((Pin*j-Pin*j_ana), np.inf)
if info > 0:
print 'Solve does not work well'
print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
@@ -425,11 +425,11 @@ class Test2D_InhomogeneousMixed(OrderTest):
j_funY = lambda x: -0.5*np.pi*np.cos(0.5*np.pi*x[:,0])*np.sin(0.5*np.pi*x[:,1])
q_fun = lambda x: -2*((0.5*np.pi)**2)*phi(x)
xc_anal = phi(self.M.gridCC)
q_anal = q_fun(self.M.gridCC)
jX_anal = j_funX(self.M.gridFx)
jY_anal = j_funY(self.M.gridFy)
j_anal = np.r_[jX_anal,jY_anal]
xc_ana = phi(self.M.gridCC)
q_ana = q_fun(self.M.gridCC)
jX_ana = j_funX(self.M.gridFx)
jY_ana = j_funY(self.M.gridFy)
j_ana = np.r_[jX_ana,jY_ana]
#TODO: Check where our boundary conditions are CCx or Nx
@@ -454,21 +454,21 @@ class Test2D_InhomogeneousMixed(OrderTest):
V = Utils.sdiag(self.M.vol)
G = -Pin.T*Pin*self.M.faceDiv.T * V
D = self.M.faceDiv
j = McI*(G*xc_anal + P*phi_bc)
j = McI*(G*xc_ana + P*phi_bc)
q = V*D*Pin.T*Pin*j + V*D*Pout.T*j_bc
# Rearrange if we know q to solve for x
A = V*D*Pin.T*Pin*McI*G
rhs = V*q_anal - V*D*Pin.T*Pin*McI*P*phi_bc - V*D*Pout.T*j_bc
rhs = V*q_ana - V*D*Pin.T*Pin*McI*P*phi_bc - V*D*Pout.T*j_bc
if self.myTest == 'j':
err = np.linalg.norm((Pin*j-Pin*j_anal), np.inf)
err = np.linalg.norm((Pin*j-Pin*j_ana), np.inf)
elif self.myTest == 'q':
err = np.linalg.norm((q-V*q_anal), np.inf)
err = np.linalg.norm((q-V*q_ana), np.inf)
elif self.myTest == 'xc':
#TODO: fix the null space
xc, info = sp.linalg.minres(A, rhs, tol = 1e-6)
err = np.linalg.norm((xc-xc_anal), np.inf)
err = np.linalg.norm((xc-xc_ana), np.inf)
if info > 0:
print 'Solve does not work well'
print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
@@ -476,7 +476,7 @@ class Test2D_InhomogeneousMixed(OrderTest):
#TODO: fix the null space
xc, info = sp.linalg.minres(A, rhs, tol = 1e-6)
j = McI*(G*xc + P*phi_bc)
err = np.linalg.norm((Pin*j-Pin*j_anal), np.inf)
err = np.linalg.norm((Pin*j-Pin*j_ana), np.inf)
if info > 0:
print 'Solve does not work well'
print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
SimPEG/Tests/test_cylMesh.py
+4 -4
View File
@@ -234,9 +234,9 @@ class TestFaceDiv2D(OrderTest):
F = self.M.projectFaceVector(Fc)
divF = self.M.faceDiv.dot(F)
divF_anal = call3(sol, self.M.gridCC)
divF_ana = call3(sol, self.M.gridCC)
err = np.linalg.norm((divF-divF_anal), np.inf)
err = np.linalg.norm((divF-divF_ana), np.inf)
return err
def test_order(self):
@@ -272,9 +272,9 @@ class TestEdgeCurl2D(OrderTest):
Fc = cylF2(self.M, solR, solZ)
Fc = np.c_[Fc[:,0],np.zeros(self.M.nF),Fc[:,1]]
curlE_anal = self.M.projectFaceVector(Fc)
curlE_ana = self.M.projectFaceVector(Fc)
err = np.linalg.norm((curlE-curlE_anal), np.inf)
err = np.linalg.norm((curlE-curlE_ana), np.inf)
return err
def test_order(self):
SimPEG/Tests/test_interpolation.py
+16 -16
View File
@@ -34,7 +34,7 @@ class TestInterpolation1D(OrderTest):
def getError(self):
funX = lambda x: np.cos(2*np.pi*x)
anal = call1(funX, self.LOCS)
ana = call1(funX, self.LOCS)
if 'CC' == self.type:
grid = call1(funX, self.M.gridCC)
@@ -43,7 +43,7 @@ class TestInterpolation1D(OrderTest):
comp = self.M.getInterpolationMat(self.LOCS, self.type)*grid
err = np.linalg.norm((comp - anal), 2)
err = np.linalg.norm((comp - ana), 2)
return err
def test_orderCC(self):
@@ -82,11 +82,11 @@ class TestInterpolation2d(OrderTest):
funY = lambda x, y: np.cos(2*np.pi*x)
if 'x' in self.type:
anal = call2(funX, self.LOCS)
ana = call2(funX, self.LOCS)
elif 'y' in self.type:
anal = call2(funY, self.LOCS)
ana = call2(funY, self.LOCS)
else:
anal = call2(funX, self.LOCS)
ana = call2(funX, self.LOCS)
if 'F' in self.type:
Fc = cartF2(self.M, funX, funY)
@@ -101,7 +101,7 @@ class TestInterpolation2d(OrderTest):
comp = self.M.getInterpolationMat(self.LOCS, self.type)*grid
err = np.linalg.norm((comp - anal), np.inf)
err = np.linalg.norm((comp - ana), np.inf)
return err
def test_orderCC(self):
@@ -165,13 +165,13 @@ class TestInterpolation2dCyl(OrderTest):
funY = lambda x, y: np.cos(2*np.pi*x)
if 'x' in self.type:
anal = call2(funX, self.LOCS)
ana = call2(funX, self.LOCS)
elif 'y' in self.type:
anal = call2(funY, self.LOCS)
ana = call2(funY, self.LOCS)
elif 'z' in self.type:
anal = call2(funY, self.LOCS)
ana = call2(funY, self.LOCS)
else:
anal = call2(funX, self.LOCS)
ana = call2(funX, self.LOCS)
if 'Fx' == self.type:
Fc = cartF2Cyl(self.M, funX, funY)
@@ -192,7 +192,7 @@ class TestInterpolation2dCyl(OrderTest):
comp = self.M.getInterpolationMat(self.LOCS, self.type)*grid
err = np.linalg.norm((comp - anal), np.inf)
err = np.linalg.norm((comp - ana), np.inf)
return err
def test_orderCC(self):
@@ -234,13 +234,13 @@ class TestInterpolation3D(OrderTest):
funZ = lambda x, y, z: np.cos(2*np.pi*x)
if 'x' in self.type:
anal = call3(funX, self.LOCS)
ana = call3(funX, self.LOCS)
elif 'y' in self.type:
anal = call3(funY, self.LOCS)
ana = call3(funY, self.LOCS)
elif 'z' in self.type:
anal = call3(funZ, self.LOCS)
ana = call3(funZ, self.LOCS)
else:
anal = call3(funX, self.LOCS)
ana = call3(funX, self.LOCS)
if 'F' in self.type:
Fc = cartF3(self.M, funX, funY, funZ)
@@ -255,7 +255,7 @@ class TestInterpolation3D(OrderTest):
comp = self.M.getInterpolationMat(self.LOCS, self.type)*grid
err = np.linalg.norm((comp - anal), np.inf)
err = np.linalg.norm((comp - ana), np.inf)
return err
def test_orderCC(self):
SimPEG/Tests/test_operators.py
+24 -24
View File
@@ -35,15 +35,15 @@ class TestCurl(OrderTest):
E = self.M.projectEdgeVector(Ec)
Fc = cartF3(self.M, solX, solY, solZ)
curlE_anal = self.M.projectFaceVector(Fc)
curlE_ana = self.M.projectFaceVector(Fc)
curlE = self.M.edgeCurl.dot(E)
if self._meshType == 'rotateLRM':
# Really it is the integration we should be caring about:
# So, let us look at the l2 norm.
err = np.linalg.norm(self.M.area*(curlE - curlE_anal), 2)
err = np.linalg.norm(self.M.area*(curlE - curlE_ana), 2)
else:
err = np.linalg.norm((curlE - curlE_anal), np.inf)
err = np.linalg.norm((curlE - curlE_ana), np.inf)
return err
def test_order(self):
@@ -63,8 +63,8 @@ class TestCurl2D(OrderTest):
sol_curl2d = call2(sol, self.M.gridCC)
Ec = cartE2(self.M, ex, ey)
sol_anal = self.M.edgeCurl*self.M.projectFaceVector(Ec)
err = np.linalg.norm((sol_curl2d-sol_anal), np.inf)
sol_ana = self.M.edgeCurl*self.M.projectFaceVector(Ec)
err = np.linalg.norm((sol_curl2d-sol_ana), np.inf)
return err
@@ -86,13 +86,13 @@ class TestCellGrad1D_InhomogeneousDirichlet(OrderTest):
xc = sol(self.M.gridCC)
gradX_anal = fx(self.M.gridFx)
gradX_ana = fx(self.M.gridFx)
bc = np.array([1,1])
self.M.setCellGradBC('dirichlet')
gradX = self.M.cellGrad.dot(xc) + self.M.cellGradBC*bc
err = np.linalg.norm((gradX-gradX_anal), np.inf)
err = np.linalg.norm((gradX-gradX_ana), np.inf)
return err
@@ -114,12 +114,12 @@ class TestCellGrad2D_Dirichlet(OrderTest):
xc = call2(sol, self.M.gridCC)
Fc = cartF2(self.M, fx, fy)
gradX_anal = self.M.projectFaceVector(Fc)
gradX_ana = self.M.projectFaceVector(Fc)
self.M.setCellGradBC('dirichlet')
gradX = self.M.cellGrad.dot(xc)
err = np.linalg.norm((gradX-gradX_anal), np.inf)
err = np.linalg.norm((gradX-gradX_ana), np.inf)
return err
@@ -143,12 +143,12 @@ class TestCellGrad3D_Dirichlet(OrderTest):
xc = call3(sol, self.M.gridCC)
Fc = cartF3(self.M, fx, fy, fz)
gradX_anal = self.M.projectFaceVector(Fc)
gradX_ana = self.M.projectFaceVector(Fc)
self.M.setCellGradBC('dirichlet')
gradX = self.M.cellGrad.dot(xc)
err = np.linalg.norm((gradX-gradX_anal), np.inf)
err = np.linalg.norm((gradX-gradX_ana), np.inf)
return err
@@ -170,12 +170,12 @@ class TestCellGrad2D_Neumann(OrderTest):
xc = call2(sol, self.M.gridCC)
Fc = cartF2(self.M, fx, fy)
gradX_anal = self.M.projectFaceVector(Fc)
gradX_ana = self.M.projectFaceVector(Fc)
self.M.setCellGradBC('neumann')
gradX = self.M.cellGrad.dot(xc)
err = np.linalg.norm((gradX-gradX_anal), np.inf)
err = np.linalg.norm((gradX-gradX_ana), np.inf)
return err
@@ -199,12 +199,12 @@ class TestCellGrad3D_Neumann(OrderTest):
xc = call3(sol, self.M.gridCC)
Fc = cartF3(self.M, fx, fy, fz)
gradX_anal = self.M.projectFaceVector(Fc)
gradX_ana = self.M.projectFaceVector(Fc)
self.M.setCellGradBC('neumann')
gradX = self.M.cellGrad.dot(xc)
err = np.linalg.norm((gradX-gradX_anal), np.inf)
err = np.linalg.norm((gradX-gradX_ana), np.inf)
return err
@@ -227,14 +227,14 @@ class TestFaceDiv3D(OrderTest):
F = self.M.projectFaceVector(Fc)
divF = self.M.faceDiv.dot(F)
divF_anal = call3(sol, self.M.gridCC)
divF_ana = call3(sol, self.M.gridCC)
if self._meshType == 'rotateLRM':
# Really it is the integration we should be caring about:
# So, let us look at the l2 norm.
err = np.linalg.norm(self.M.vol*(divF-divF_anal), 2)
err = np.linalg.norm(self.M.vol*(divF-divF_ana), 2)
else:
err = np.linalg.norm((divF-divF_anal), np.inf)
err = np.linalg.norm((divF-divF_ana), np.inf)
return err
def test_order(self):
@@ -257,9 +257,9 @@ class TestFaceDiv2D(OrderTest):
F = self.M.projectFaceVector(Fc)
divF = self.M.faceDiv.dot(F)
divF_anal = call2(sol, self.M.gridCC)
divF_ana = call2(sol, self.M.gridCC)
err = np.linalg.norm((divF-divF_anal), np.inf)
err = np.linalg.norm((divF-divF_ana), np.inf)
return err
@@ -283,9 +283,9 @@ class TestNodalGrad(OrderTest):
gradE = self.M.nodalGrad.dot(phi)
Ec = cartE3(self.M, solX, solY, solZ)
gradE_anal = self.M.projectEdgeVector(Ec)
gradE_ana = self.M.projectEdgeVector(Ec)
err = np.linalg.norm((gradE-gradE_anal), np.inf)
err = np.linalg.norm((gradE-gradE_ana), np.inf)
return err
@@ -309,9 +309,9 @@ class TestNodalGrad2D(OrderTest):
gradE = self.M.nodalGrad.dot(phi)
Ec = cartE2(self.M, solX, solY)
gradE_anal = self.M.projectEdgeVector(Ec)
gradE_ana = self.M.projectEdgeVector(Ec)
err = np.linalg.norm((gradE-gradE_anal), np.inf)
err = np.linalg.norm((gradE-gradE_ana), np.inf)
return err