Merge branch 'develop' of https://github.com/simpeg/simpeg into cylClean

Conflicts:
	SimPEG/Mesh/View.py
	SimPEG/Tests/test_utils.py
This commit is contained in:
rowanc1
2014-02-26 10:10:29 -08:00
23 changed files with 894 additions and 131 deletions
+7 -13
View File
@@ -52,6 +52,7 @@ class BaseData(object):
instead of recalculating the fields (which may be expensive!).
.. math::
d_\\text{pred} = P(u(m))
Where P is a projection of the fields onto the data space.
@@ -67,29 +68,22 @@ class BaseData(object):
.. math::
d_\\text{pred} = \mathbf{P} u(m)
"""
return u
@Utils.count
def projectFieldsAdjoint(self, d):
def projectFieldsDeriv(self, u):
"""
This function is the adjoint of the projection.
**projectFieldsAdjoint** is used in the
calculation of the sensitivities.
This function projects the fields onto the data space.
.. math::
u = \mathbf{P}^\\top d
:param numpy.array d: data
:param numpy.array u: fields (ish)
:rtype: fields like object
:return: data
\\frac{\partial d_\\text{pred}}{\partial u} = \mathbf{P}
"""
return d
#TODO: def projectFieldDeriv(self, u): Does this need to be made??!
return sp.identity(u.size)
@Utils.count
def residual(self, m, u=None):
+104 -2
View File
@@ -291,8 +291,8 @@ class DiffOperators(object):
"""
if(type(BC) is str):
BC = [BC for _ in self.vnC] # Repeat the str self.dim times
elif(type(BC) is list):
BC = [BC]*self.dim
if(type(BC) is list):
assert len(BC) == self.dim, 'BC list must be the size of your mesh'
else:
raise Exception("BC must be a str or a list.")
@@ -460,6 +460,108 @@ class DiffOperators(object):
_edgeCurl = None
edgeCurl = property(**edgeCurl())
def getBCProjWF(self, BC, discretization='CC'):
"""
The weak form boundary condition projection matrices.
Examples::
BC = 'neumann' # Neumann in all directions
BC = ['neumann', 'dirichlet', 'neumann'] # 3D, Dirichlet in y Neumann else
BC = [['neumann', 'dirichlet'], 'dirichlet', 'dirichlet'] # 3D, Neumann in x on bottom of domain,
# Dirichlet else
"""
if discretization is not 'CC':
raise NotImplementedError('Boundary conditions only implemented for CC discretization.')
if(type(BC) is str):
BC = [BC for _ in self.vnC] # Repeat the str self.dim times
elif(type(BC) is list):
assert len(BC) == self.dim, 'BC list must be the size of your mesh'
else:
raise Exception("BC must be a str or a list.")
for i, bc_i in enumerate(BC):
BC[i] = checkBC(bc_i)
def projDirichlet(n, bc):
bc = checkBC(bc)
ij = ([0,n], [0,1])
vals = [0,0]
if(bc[0] == 'dirichlet'):
vals[0] = -1
if(bc[1] == 'dirichlet'):
vals[1] = 1
return sp.csr_matrix((vals, ij), shape=(n+1,2))
def projNeumannIn(n, bc):
bc = checkBC(bc)
P = sp.identity(n+1).tocsr()
if(bc[0] == 'neumann'):
P = P[1:,:]
if(bc[1] == 'neumann'):
P = P[:-1,:]
return P
def projNeumannOut(n, bc):
bc = checkBC(bc)
ij = ([0, 1],[0, n])
vals = [0,0]
if(bc[0] == 'neumann'):
vals[0] = 1
if(bc[1] == 'neumann'):
vals[1] = 1
return sp.csr_matrix((vals, ij), shape=(2,n+1))
n = self.vnC
indF = self.faceBoundaryInd
if(self.dim == 1):
Pbc = projDirichlet(n[0], BC[0])
indF = indF[0] | indF[1]
Pbc = Pbc*sdiag(self.area[indF])
Pin = projNeumannIn(n[0], BC[0])
Pout = projNeumannOut(n[0], BC[0])
elif(self.dim == 2):
Pbc1 = sp.kron(speye(n[1]), projDirichlet(n[0], BC[0]))
Pbc2 = sp.kron(projDirichlet(n[1], BC[1]), speye(n[0]))
Pbc = sp.block_diag((Pbc1, Pbc2), format="csr")
indF = np.r_[(indF[0] | indF[1]), (indF[2] | indF[3])]
Pbc = Pbc*sdiag(self.area[indF])
P1 = sp.kron(speye(n[1]), projNeumannIn(n[0], BC[0]))
P2 = sp.kron(projNeumannIn(n[1], BC[1]), speye(n[0]))
Pin = sp.block_diag((P1, P2), format="csr")
P1 = sp.kron(speye(n[1]), projNeumannOut(n[0], BC[0]))
P2 = sp.kron(projNeumannOut(n[1], BC[1]), speye(n[0]))
Pout = sp.block_diag((P1, P2), format="csr")
elif(self.dim == 3):
Pbc1 = kron3(speye(n[2]), speye(n[1]), projDirichlet(n[0], BC[0]))
Pbc2 = kron3(speye(n[2]), projDirichlet(n[1], BC[1]), speye(n[0]))
Pbc3 = kron3(projDirichlet(n[2], BC[2]), speye(n[1]), speye(n[0]))
Pbc = sp.block_diag((Pbc1, Pbc2, Pbc3), format="csr")
indF = np.r_[(indF[0] | indF[1]), (indF[2] | indF[3]), (indF[4] | indF[5])]
Pbc = Pbc*sdiag(self.area[indF])
P1 = kron3(speye(n[2]), speye(n[1]), projNeumannIn(n[0], BC[0]))
P2 = kron3(speye(n[2]), projNeumannIn(n[1], BC[1]), speye(n[0]))
P3 = kron3(projNeumannIn(n[2], BC[2]), speye(n[1]), speye(n[0]))
Pin = sp.block_diag((P1, P2, P3), format="csr")
P1 = kron3(speye(n[2]), speye(n[1]), projNeumannOut(n[0], BC[0]))
P2 = kron3(speye(n[2]), projNeumannOut(n[1], BC[1]), speye(n[0]))
P3 = kron3(projNeumannOut(n[2], BC[2]), speye(n[1]), speye(n[0]))
Pout = sp.block_diag((P1, P2, P3), format="csr")
return Pbc, Pin, Pout
# --------------- Averaging ---------------------
@property
+50
View File
@@ -396,6 +396,56 @@ class TensorMesh(BaseTensorMesh, TensorView, DiffOperators, InnerProducts):
Q[indZeros, :] = 0
return Q.tocsr()
@property
def faceBoundaryInd(self):
"""
Find indices of boundary faces in each direction
"""
if self.dim==1:
indxd = (self.gridFx==min(self.gridFx))
indxu = (self.gridFx==max(self.gridFx))
return indxd, indxu
elif self.dim==2:
indxd = (self.gridFx[:,0]==min(self.gridFx[:,0]))
indxu = (self.gridFx[:,0]==max(self.gridFx[:,0]))
indyd = (self.gridFy[:,1]==min(self.gridFy[:,1]))
indyu = (self.gridFy[:,1]==max(self.gridFy[:,1]))
return indxd, indxu, indyd, indyu
elif self.dim==3:
indxd = (self.gridFx[:,0]==min(self.gridFx[:,0]))
indxu = (self.gridFx[:,0]==max(self.gridFx[:,0]))
indyd = (self.gridFy[:,1]==min(self.gridFy[:,1]))
indyu = (self.gridFy[:,1]==max(self.gridFy[:,1]))
indzd = (self.gridFz[:,2]==min(self.gridFz[:,2]))
indzu = (self.gridFz[:,2]==max(self.gridFz[:,2]))
return indxd, indxu, indyd, indyu, indzd, indzu
@property
def cellBoundaryInd(self):
"""
Find indices of boundary faces in each direction
"""
if self.dim==1:
indxd = (self.gridCC==min(self.gridCC))
indxu = (self.gridCC==max(self.gridCC))
return indxd, indxu
elif self.dim==2:
indxd = (self.gridCC[:,0]==min(self.gridCC[:,0]))
indxu = (self.gridCC[:,0]==max(self.gridCC[:,0]))
indyd = (self.gridCC[:,1]==min(self.gridCC[:,1]))
indyu = (self.gridCC[:,1]==max(self.gridCC[:,1]))
return indxd, indxu, indyd, indyu
elif self.dim==3:
indxd = (self.gridCC[:,0]==min(self.gridCC[:,0]))
indxu = (self.gridCC[:,0]==max(self.gridCC[:,0]))
indyd = (self.gridCC[:,1]==min(self.gridCC[:,1]))
indyu = (self.gridCC[:,1]==max(self.gridCC[:,1]))
indzd = (self.gridCC[:,2]==min(self.gridCC[:,2]))
indzu = (self.gridCC[:,2]==max(self.gridCC[:,2]))
return indxd, indxu, indyd, indyu, indzd, indzu
if __name__ == '__main__':
print('Welcome to tensor mesh!')
+69 -11
View File
@@ -74,7 +74,7 @@ class TensorView(object):
if I.size != np.prod(self.vnEz): ex, ey, I = self.r(I,'E','E','M')
elif imageType[0] == 'E':
plotAll = len(imageType) == 1
options = {"direction":direction,"numbering":numbering,"annotationColor":annotationColor,"showIt":showIt}
options = {"direction":direction,"numbering":numbering,"annotationColor":annotationColor,"showIt":False}
fig = plt.figure(figNum)
# Determine the subplot number: 131, 121
numPlots = 130 if plotAll else len(imageType)/2*10+100
@@ -92,10 +92,11 @@ class TensorView(object):
ax_z = plt.subplot(numPlots+pltNum)
self.plotImage(ez, imageType='Ez', ax=ax_z, **options)
pltNum +=1
if showIt: plt.show()
return
elif imageType[0] == 'F':
plotAll = len(imageType) == 1
options = {"direction":direction,"numbering":numbering,"annotationColor":annotationColor,"showIt":showIt}
options = {"direction":direction,"numbering":numbering,"annotationColor":annotationColor,"showIt":False}
fig = plt.figure(figNum)
# Determine the subplot number: 131, 121
numPlots = 130 if plotAll else len(imageType)/2*10+100
@@ -113,6 +114,7 @@ class TensorView(object):
ax_z = plt.subplot(numPlots+pltNum)
self.plotImage(fxyz[2], imageType='Fz', ax=ax_z, **options)
pltNum +=1
if showIt: plt.show()
return
else:
raise Exception("imageType must be 'CC', 'N','Fx','Fy','Fz','Ex','Ey','Ez'")
@@ -241,11 +243,30 @@ class TensorView(object):
def plotSlice(self, v, vType='CC',
normal='Z', ind=None, grid=False, view='real',
ax=None, clim=None, showIt=False,
pcolorOpts={},
streamOpts={'color':'k'},
gridOpts={'color':'k'}
):
"""
Plots a slice of a 3D mesh.
.. plot::
from SimPEG import *
mT = Utils.meshTensors(((2,5),(4,2),(2,5)),((2,2),(6,2),(2,2)),((2,2),(6,2),(2,2)))
M = Mesh.TensorMesh(mT)
q = np.zeros(M.vnC)
q[[4,4],[4,4],[2,6]]=[-1,1]
q = Utils.mkvc(q)
A = M.faceDiv*M.cellGrad
b = Solver(A).solve(q)
M.plotSlice(M.cellGrad*b, 'F', view='vec', grid=True, showIt=True, pcolorOpts={'alpha':0.8})
"""
viewOpts = ['real','imag','abs','vec']
normalOpts = ['X', 'Y', 'Z']
vTypeOpts = ['CC','F','E']
vTypeOpts = ['CC', 'CCv','F','E']
# Some user error checking
assert vType in vTypeOpts, "vType must be in ['%s']" % "','".join(vTypeOpts)
@@ -277,11 +298,15 @@ class TensorView(object):
def doSlice(v):
if vType == 'CC':
return getIndSlice(self.r(v,'CC','CC','M'))
# Now just deal with 'F' and 'E'
aveOp = 'ave' + vType + ('2CCV' if view == 'vec' else '2CC')
v = getattr(self,aveOp)*v # average to cell centers (might be a vector)
if view == 'vec':
elif vType == 'CCv':
v = self.r(v.reshape((self.nC,3),order='F'),'CC','CC','M')
assert view == 'vec', 'Other types for CCv not yet supported'
else:
# Now just deal with 'F' and 'E'
aveOp = 'ave' + vType + ('2CCV' if view == 'vec' else '2CC')
v = getattr(self,aveOp)*v # average to cell centers (might be a vector)
v = self.r(v.reshape((self.nC,3),order='F'),'CC','CC','M')
if view == 'vec':
outSlice = []
if 'X' not in normal: outSlice.append(getIndSlice(v[0]))
if 'Y' not in normal: outSlice.append(getIndSlice(v[1]))
@@ -302,13 +327,31 @@ class TensorView(object):
v = doSlice(v)
if clim is None:
clim = [v.min(),v.max()]
out += (ax.pcolormesh(tM.vectorNx, tM.vectorNy, v.T, vmin=clim[0], vmax=clim[1]),)
out += (ax.pcolormesh(tM.vectorNx, tM.vectorNy, v.T, vmin=clim[0], vmax=clim[1], **pcolorOpts),)
elif view in ['vec']:
U, V = doSlice(v)
if clim is None:
clim = [v.min(),v.max()]
out += (ax.pcolormesh(tM.vectorNx, tM.vectorNy, 0.5*(U+V).T, vmin=clim[0], vmax=clim[1]),)
out += (plt.streamplot(tM.vectorCCx, tM.vectorCCy, U.T, V.T),)
uv = np.r_[mkvc(U), mkvc(V)]
uv = np.sqrt(uv**2)
clim = [uv.min(),uv.max()]
# Matplotlib seems to not support irregular
# spaced vectors at the moment. So we will
# Interpolate down to a regular mesh at the
# smallest mesh size in this 2D slice.
nxi = int(tM.hx.sum()/tM.hx.min())
nyi = int(tM.hy.sum()/tM.hy.min())
tMi = self.__class__([np.ones(nxi)*tM.hx.sum()/nxi,
np.ones(nyi)*tM.hy.sum()/nyi])
P = tM.getInterpolationMat(tMi.gridCC,'CC',zerosOutside=True)
Ui = P*mkvc(U)
Vi = P*mkvc(V)
Ui = tMi.r(Ui, 'CC', 'CC', 'M')
Vi = tMi.r(Vi, 'CC', 'CC', 'M')
# End Interpolation
out += (ax.pcolormesh(tM.vectorNx, tM.vectorNy, np.sqrt(U**2+V**2).T, vmin=clim[0], vmax=clim[1], **pcolorOpts),)
out += (ax.streamplot(tMi.vectorCCx, tMi.vectorCCy, Ui.T, Vi.T, **streamOpts),)
if grid:
xXGrid = np.c_[tM.vectorNx,tM.vectorNx,np.nan*np.ones(tM.nNx)].flatten()
@@ -320,6 +363,8 @@ class TensorView(object):
ax.set_xlabel('y' if normal == 'X' else 'x')
ax.set_ylabel('y' if normal == 'Z' else 'z')
ax.set_title('Slice %d' % ind)
ax.set_xlim(*tM.vectorNx[[0,-1]])
ax.set_ylim(*tM.vectorNy[[0,-1]])
if showIt: plt.show()
return out
@@ -610,3 +655,16 @@ class LomView(object):
ax.set_ylabel('x2')
if showIt: plt.show()
if __name__ == '__main__':
from SimPEG import *
mT = Utils.meshTensors(((2,5),(4,2),(2,5)),((2,2),(6,2),(2,2)),((2,2),(6,2),(2,2)))
M = Mesh.TensorMesh(mT)
q = np.zeros(M.vnC)
q[[4,4],[4,4],[2,6]]=[-1,1]
q = Utils.mkvc(q)
A = M.faceDiv*M.cellGrad
b = Solver(A).solve(q)
M.plotSlice(M.cellGrad*b, 'F', view='vec', grid=True, showIt=True, pcolorOpts={'alpha':0.8})
+62
View File
@@ -64,6 +64,68 @@ class BaseModel(object):
m = self.example()
return checkDerivative(lambda m : [self.transform(m), self.transformDeriv(m)], m, plotIt=False)
class BaseNonLinearModel(object):
"""
SimPEG BaseNonLinearModel
"""
__metaclass__ = Utils.SimPEGMetaClass
counter = None #: A SimPEG.Utils.Counter object
mesh = None #: A SimPEG Mesh
def __init__(self, mesh):
self.mesh = mesh
def transform(self, u, m):
"""
:param numpy.array u: fields
:param numpy.array m: model
:rtype: numpy.array
:return: transformed model
The *transform* changes the model into the physical property.
"""
return m
def transformDerivU(self, u, m):
"""
:param numpy.array u: fields
:param numpy.array m: model
:rtype: scipy.csr_matrix
:return: derivative of transformed model
The *transform* changes the model into the physical property.
The *transformDerivU* provides the derivative of the *transform* with respect to the fields.
"""
raise NotImplementedError('The transformDerivU is not implemented.')
def transformDerivM(self, u, m):
"""
:param numpy.array u: fields
:param numpy.array m: model
:rtype: scipy.csr_matrix
:return: derivative of transformed model
The *transform* changes the model into the physical property.
The *transformDerivU* provides the derivative of the *transform* with respect to the model.
"""
raise NotImplementedError('The transformDerivM is not implemented.')
@property
def nP(self):
"""Number of parameters in the model."""
return self.mesh.nC
def example(self):
raise NotImplementedError('The example is not implemented.')
def test(self, m=None):
raise NotImplementedError('The test is not implemented.')
class LogModel(BaseModel):
"""SimPEG LogModel"""
+4 -2
View File
@@ -1,5 +1,5 @@
import Utils, Data, numpy as np, scipy.sparse as sp
import Model
class BaseProblem(object):
"""
@@ -39,10 +39,12 @@ class BaseProblem(object):
counter = None #: A SimPEG.Utils.Counter object
dataPair = Data.BaseData
modelPair = Model.BaseModel
def __init__(self, mesh, model, *args, **kwargs):
def __init__(self, mesh, model, **kwargs):
Utils.setKwargs(self, **kwargs)
self.mesh = mesh
assert isinstance(model, self.modelPair), "Model object must be an instance of a %s class."%(self.modelPair.__name__)
self.model = model
@property
+7 -5
View File
@@ -90,9 +90,9 @@ class OrderTest(unittest.TestCase):
if 'uniform' in self._meshType:
h = [nc, nc, nc]
elif 'random' in self._meshType:
h1 = np.random.rand(nc)
h2 = np.random.rand(nc)
h3 = np.random.rand(nc)
h1 = np.random.rand(nc)*nc*0.5 + nc*0.5
h2 = np.random.rand(nc)*nc*0.5 + nc*0.5
h3 = np.random.rand(nc)*nc*0.5 + nc*0.5
h = [hi/np.sum(hi) for hi in [h1, h2, h3]] # normalize
else:
raise Exception('Unexpected meshType')
@@ -122,10 +122,12 @@ class OrderTest(unittest.TestCase):
kwrd = 'rotate'
else:
raise Exception('Unexpected meshType')
if self.meshDimension == 2:
if self.meshDimension == 1:
raise Exception('Lom not supported for 1D')
elif self.meshDimension == 2:
X, Y = Utils.exampleLomGird([nc, nc], kwrd)
self.M = LogicallyOrthogonalMesh([X, Y])
if self.meshDimension == 3:
elif self.meshDimension == 3:
X, Y, Z = Utils.exampleLomGird([nc, nc, nc], kwrd)
self.M = LogicallyOrthogonalMesh([X, Y, Z])
return 1./nc
+501
View File
@@ -0,0 +1,501 @@
import numpy as np
import scipy.sparse as sp
import unittest
from TestUtils import OrderTest
import matplotlib.pyplot as plt
from SimPEG import Utils, Solver
MESHTYPES = ['uniformTensorMesh']
class Test1D_InhomogeneousDirichlet(OrderTest):
name = "1D - Dirichlet"
meshTypes = MESHTYPES
meshDimension = 1
expectedOrders = 2
meshSizes = [4, 8, 16, 32, 64, 128]
def getError(self):
#Test function
phi = lambda x: np.cos(np.pi*x)
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)
#TODO: Check where our boundary conditions are CCx or Nx
# vec = self.M.vectorNx
vec = self.M.vectorCCx
phi_bc = phi(vec[[0,-1]])
j_bc = j_fun(vec[[0,-1]])
P, Pin, Pout = self.M.getBCProjWF([['dirichlet', 'dirichlet']])
Mc = self.M.getFaceInnerProduct()
McI = Utils.sdInv(self.M.getFaceInnerProduct())
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)
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
# A = D*McI*G
# rhs = q_anal - D*McI*P*phi_bc
if self.myTest == 'j':
err = np.linalg.norm((j-j_anal), np.inf)
elif self.myTest == 'q':
err = np.linalg.norm((q-V*q_anal), np.inf)
elif self.myTest == 'xc':
#TODO: fix the null space
solver = Solver(A, doDirect=False, options={'M':'J','iterSolver':'CG','backend':'scipy','maxIter':1000})
xc = solver.solve(rhs)
print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
err = np.linalg.norm((xc-xc_anal), np.inf)
elif self.myTest == 'xcJ':
#TODO: fix the null space
xc = Solver(A).solve(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)
return err
def test_orderJ(self):
self.name = "1D - InhomogeneousDirichlet_Forward j"
self.myTest = 'j'
self.orderTest()
def test_orderQ(self):
self.name = "1D - InhomogeneousDirichlet_Forward q"
self.myTest = 'q'
self.orderTest()
def test_orderX(self):
self.name = "1D - InhomogeneousDirichlet_Inverse"
self.myTest = 'xc'
self.orderTest()
def test_orderXJ(self):
self.name = "1D - InhomogeneousDirichlet_Inverse J"
self.myTest = 'xcJ'
self.orderTest()
class Test2D_InhomogeneousDirichlet(OrderTest):
name = "2D - Dirichlet"
meshTypes = MESHTYPES
meshDimension = 2
expectedOrders = 2
meshSizes = [4, 8, 16, 32]
def getError(self):
#Test function
phi = lambda x: np.cos(np.pi*x[:,0])*np.cos(np.pi*x[:,1])
j_funX = lambda x: -np.pi*np.sin(np.pi*x[:,0])*np.cos(np.pi*x[:,1])
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]
#TODO: Check where our boundary conditions are CCx or Nx
# fxm,fxp,fym,fyp = self.M.faceBoundaryInd
# gBFx = self.M.gridFx[(fxm|fxp),:]
# gBFy = self.M.gridFy[(fym|fyp),:]
fxm,fxp,fym,fyp = self.M.cellBoundaryInd
gBFx = self.M.gridCC[(fxm|fxp),:]
gBFy = self.M.gridCC[(fym|fyp),:]
bc = phi(np.r_[gBFx,gBFy])
# P = sp.csr_matrix(([-1,1],([0,self.M.nF-1],[0,1])), shape=(self.M.nF, 2))
P, Pin, Pout = self.M.getBCProjWF('dirichlet')
Mc = self.M.getFaceInnerProduct()
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)
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
if self.myTest == 'j':
err = np.linalg.norm((j-j_anal), np.inf)
elif self.myTest == 'q':
err = np.linalg.norm((q-q_anal), np.inf)
elif self.myTest == 'xc':
xc = Solver(A).solve(rhs)
err = np.linalg.norm((xc-xc_anal), np.inf)
elif self.myTest == 'xcJ':
xc = Solver(A).solve(rhs)
j = McI*(G*xc + P*bc)
err = np.linalg.norm((j-j_anal), np.inf)
return err
def test_orderJ(self):
self.name = "2D - InhomogeneousDirichlet_Forward j"
self.myTest = 'j'
self.orderTest()
def test_orderQ(self):
self.name = "2D - InhomogeneousDirichlet_Forward q"
self.myTest = 'q'
self.orderTest()
def test_orderX(self):
self.name = "2D - InhomogeneousDirichlet_Inverse"
self.myTest = 'xc'
self.orderTest()
def test_orderXJ(self):
self.name = "2D - InhomogeneousDirichlet_Inverse J"
self.myTest = 'xcJ'
self.orderTest()
class Test1D_InhomogeneousNeumann(OrderTest):
name = "1D - Neumann"
meshTypes = MESHTYPES
meshDimension = 1
expectedOrders = 2
meshSizes = [4, 8, 16, 32, 64, 128]
def getError(self):
#Test function
phi = lambda x: np.sin(np.pi*x)
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)
#TODO: Check where our boundary conditions are CCx or Nx
vecN = self.M.vectorNx
vecC = self.M.vectorCCx
phi_bc = phi(vecC[[0,-1]])
j_bc = j_fun(vecN[[0,-1]])
P, Pin, Pout = self.M.getBCProjWF([['neumann', 'neumann']])
Mc = self.M.getFaceInnerProduct()
McI = Utils.sdInv(self.M.getFaceInnerProduct())
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)
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
# A = D*McI*G
# rhs = q_anal - D*McI*P*phi_bc
if self.myTest == 'j':
err = np.linalg.norm((Pin*j-Pin*j_anal), np.inf)
elif self.myTest == 'q':
err = np.linalg.norm((q-V*q_anal), 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)
if info > 0:
print 'Solve does not work well'
print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
elif self.myTest == 'xcJ':
#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)
if info > 0:
print 'Solve does not work well'
print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
return err
def test_orderJ(self):
self.name = "1D - InhomogeneousNeumann_Forward j"
self.myTest = 'j'
self.orderTest()
def test_orderQ(self):
self.name = "1D - InhomogeneousNeumann_Forward q"
self.myTest = 'q'
self.orderTest()
def test_orderXJ(self):
self.name = "1D - InhomogeneousNeumann_Inverse J"
self.myTest = 'xcJ'
self.orderTest()
class Test2D_InhomogeneousNeumann(OrderTest):
name = "2D - Neumann"
meshTypes = MESHTYPES
meshDimension = 2
expectedOrders = 2
meshSizes = [4, 8, 16, 32]
# meshSizes = [4]
def getError(self):
#Test function
phi = lambda x: np.sin(np.pi*x[:,0])*np.sin(np.pi*x[:,1])
j_funX = lambda x: np.pi*np.cos(np.pi*x[:,0])*np.sin(np.pi*x[:,1])
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]
#TODO: Check where our boundary conditions are CCx or Nx
cxm,cxp,cym,cyp = self.M.cellBoundaryInd
fxm,fxp,fym,fyp = self.M.faceBoundaryInd
gBFx = self.M.gridFx[(fxm|fxp),:]
gBFy = self.M.gridFy[(fym|fyp),:]
gBCx = self.M.gridCC[(cxm|cxp),:]
gBCy = self.M.gridCC[(cym|cyp),:]
phi_bc = phi(np.r_[gBFx,gBFy])
j_bc = np.r_[j_funX(gBFx), j_funY(gBFy)]
# P = sp.csr_matrix(([-1,1],([0,self.M.nF-1],[0,1])), shape=(self.M.nF, 2))
P, Pin, Pout = self.M.getBCProjWF('neumann')
Mc = self.M.getFaceInnerProduct()
McI = Utils.sdInv(self.M.getFaceInnerProduct())
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)
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
if self.myTest == 'j':
err = np.linalg.norm((Pin*j-Pin*j_anal), np.inf)
elif self.myTest == 'q':
err = np.linalg.norm((q-V*q_anal), 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)
if info > 0:
print 'Solve does not work well'
print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
elif self.myTest == 'xcJ':
#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)
if info > 0:
print 'Solve does not work well'
print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
return err
def test_orderJ(self):
self.name = "2D - InhomogeneousNeumann_Forward j"
self.myTest = 'j'
self.orderTest()
def test_orderQ(self):
self.name = "2D - InhomogeneousNeumann_Forward q"
self.myTest = 'q'
self.orderTest()
def test_orderXJ(self):
self.name = "2D - InhomogeneousNeumann_Inverse J"
self.myTest = 'xcJ'
self.orderTest()
class Test1D_InhomogeneousMixed(OrderTest):
name = "1D - Mixed"
meshTypes = MESHTYPES
meshDimension = 1
expectedOrders = 2
meshSizes = [4, 8, 16, 32, 64, 128]
def getError(self):
#Test function
phi = lambda x: np.cos(0.5*np.pi*x)
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)
#TODO: Check where our boundary conditions are CCx or Nx
vecN = self.M.vectorNx
vecC = self.M.vectorCCx
phi_bc = phi(vecC[[0,-1]])
j_bc = j_fun(vecN[[0,-1]])
P, Pin, Pout = self.M.getBCProjWF([['dirichlet', 'neumann']])
Mc = self.M.getFaceInnerProduct()
McI = Utils.sdInv(self.M.getFaceInnerProduct())
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)
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
# A = D*McI*G
# rhs = q_anal - D*McI*P*phi_bc
if self.myTest == 'j':
err = np.linalg.norm((Pin*j-Pin*j_anal), np.inf)
elif self.myTest == 'q':
err = np.linalg.norm((q-V*q_anal), 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)
if info > 0:
print 'Solve does not work well'
print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
elif self.myTest == 'xcJ':
#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)
if info > 0:
print 'Solve does not work well'
print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
return err
def test_orderJ(self):
self.name = "1D - InhomogeneousMixed_Forward j"
self.myTest = 'j'
self.orderTest()
def test_orderQ(self):
self.name = "1D - InhomogeneousMixed_Forward q"
self.myTest = 'q'
self.orderTest()
def test_orderXJ(self):
self.name = "1D - InhomogeneousMixed_Inverse J"
self.myTest = 'xcJ'
self.orderTest()
class Test2D_InhomogeneousMixed(OrderTest):
name = "2D - Mixed"
meshTypes = MESHTYPES
meshDimension = 2
expectedOrders = 2
meshSizes = [2, 4, 8, 16]
# meshSizes = [4]
def getError(self):
#Test function
phi = lambda x: np.cos(0.5*np.pi*x[:,0])*np.cos(0.5*np.pi*x[:,1])
j_funX = lambda x: -0.5*np.pi*np.sin(0.5*np.pi*x[:,0])*np.cos(0.5*np.pi*x[:,1])
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]
#TODO: Check where our boundary conditions are CCx or Nx
cxm,cxp,cym,cyp = self.M.cellBoundaryInd
fxm,fxp,fym,fyp = self.M.faceBoundaryInd
gBFx = self.M.gridFx[(fxm|fxp),:]
gBFy = self.M.gridFy[(fym|fyp),:]
gBCx = self.M.gridCC[(cxm|cxp),:]
gBCy = self.M.gridCC[(cym|cyp),:]
phi_bc = phi(np.r_[gBCx,gBCy])
j_bc = np.r_[j_funX(gBFx), j_funY(gBFy)]
# P = sp.csr_matrix(([-1,1],([0,self.M.nF-1],[0,1])), shape=(self.M.nF, 2))
P, Pin, Pout = self.M.getBCProjWF([['dirichlet', 'neumann'], ['dirichlet', 'neumann']])
Mc = self.M.getFaceInnerProduct()
McI = Utils.sdInv(self.M.getFaceInnerProduct())
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)
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
if self.myTest == 'j':
err = np.linalg.norm((Pin*j-Pin*j_anal), np.inf)
elif self.myTest == 'q':
err = np.linalg.norm((q-V*q_anal), 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)
if info > 0:
print 'Solve does not work well'
print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
elif self.myTest == 'xcJ':
#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)
if info > 0:
print 'Solve does not work well'
print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
return err
def test_orderJ(self):
self.name = "2D - InhomogeneousMixed_Forward j"
self.myTest = 'j'
self.orderTest()
def test_orderQ(self):
self.name = "2D - InhomogeneousMixed_Forward q"
self.myTest = 'q'
self.orderTest()
def test_orderXJ(self):
self.name = "2D - InhomogeneousMixed_Inverse J"
self.myTest = 'xcJ'
self.orderTest()
if __name__ == '__main__':
unittest.main()
+1 -1
View File
@@ -178,7 +178,7 @@ class TestEdgeCurl2D(OrderTest):
meshDimension = 2
def getError(self):
pass
#TODO!
# def test_order(self):
+7 -6
View File
@@ -3,6 +3,7 @@ from SimPEG.Utils import *
from SimPEG import Mesh, np, sp
from SimPEG.Tests import checkDerivative
TOL = 1e-8
class TestCheckDerivative(unittest.TestCase):
@@ -104,7 +105,7 @@ class TestSequenceFunctions(unittest.TestCase):
sp.hstack((sdiag(a[2]), sdiag(a[3])))))
Z2 = B*A - sp.eye(10, 10)
self.assertTrue(np.linalg.norm(Z2.todense().ravel(), 2) < 1e-10)
self.assertTrue(np.linalg.norm(Z2.todense().ravel(), 2) < TOL)
a = [np.random.rand(5, 1) for i in range(9)]
B = inv3X3BlockDiagonal(*a)
@@ -115,7 +116,7 @@ class TestSequenceFunctions(unittest.TestCase):
Z3 = B*A - sp.eye(15, 15)
self.assertTrue(np.linalg.norm(Z3.todense().ravel(), 2) < 1e-10)
self.assertTrue(np.linalg.norm(Z3.todense().ravel(), 2) < TOL)
def test_invPropertyTensor2D(self):
@@ -134,9 +135,9 @@ class TestSequenceFunctions(unittest.TestCase):
B2 = invPropertyTensor(M, prop, returnMatrix=True)
Z = B1*A - sp.identity(M.nC*2)
self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < 1e-12)
self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < TOL)
Z = B2*A - sp.identity(M.nC*2)
self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < 1e-12)
self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < TOL)
def test_invPropertyTensor3D(self):
@@ -158,9 +159,9 @@ class TestSequenceFunctions(unittest.TestCase):
B2 = invPropertyTensor(M, prop, returnMatrix=True)
Z = B1*A - sp.identity(M.nC*3)
self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < 1e-12)
self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < TOL)
Z = B2*A - sp.identity(M.nC*3)
self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < 1e-12)
self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < TOL)
if __name__ == '__main__':
-8
View File
@@ -1,8 +0,0 @@
.. _api_BaseMesh:
Base Mesh
*********
.. automodule:: SimPEG.Mesh.BaseMesh
:members:
:undoc-members:
-8
View File
@@ -1,8 +0,0 @@
.. _api_Cyl1DMesh:
Cylindrical 1D Mesh
*******************
.. automodule:: SimPEG.Mesh.Cyl1DMesh
:members:
:undoc-members:
-8
View File
@@ -1,8 +0,0 @@
.. _api_DiffOperators:
Differential Operators
**********************
.. automodule:: SimPEG.Mesh.DiffOperators
:members:
:undoc-members:
+3 -3
View File
@@ -2,7 +2,7 @@
Model
*****
=====
.. automodule:: SimPEG.Model
:show-inheritance:
@@ -11,7 +11,7 @@ Model
:inherited-members:
Data
****
====
.. automodule:: SimPEG.Data
:show-inheritance:
@@ -20,7 +20,7 @@ Data
:inherited-members:
Problem
*******
=======
.. automodule:: SimPEG.Problem
:show-inheritance:
-8
View File
@@ -1,8 +0,0 @@
.. _api_InnerProducts:
Inner Products
**************
.. automodule:: SimPEG.Mesh.InnerProducts
:members:
:undoc-members:
+4 -4
View File
@@ -2,7 +2,7 @@
Regularization
**************
==============
.. automodule:: SimPEG.Regularization
:show-inheritance:
@@ -11,7 +11,7 @@ Regularization
Objective Function
******************
==================
.. automodule:: SimPEG.ObjFunction
:members:
@@ -19,7 +19,7 @@ Objective Function
Optimize
********
========
.. automodule:: SimPEG.Optimization
:show-inheritance:
@@ -28,7 +28,7 @@ Optimize
:undoc-members:
Inversion
*********
=========
.. automodule:: SimPEG.Inversion
:show-inheritance:
-10
View File
@@ -1,10 +0,0 @@
.. _api_LogicallyOrthogonalMesh:
Logically Orthogonal Mesh
*************************
.. automodule:: SimPEG.Mesh.LogicallyOrthogonalMesh
:show-inheritance:
:members:
:undoc-members:
:inherited-members:
+53
View File
@@ -0,0 +1,53 @@
.. _api_Mesh:
Tensor Mesh
===========
.. automodule:: SimPEG.Mesh.TensorMesh
:show-inheritance:
:members:
:undoc-members:
:inherited-members:
Cylindrical 1D Mesh
===================
.. automodule:: SimPEG.Mesh.Cyl1DMesh
:members:
:undoc-members:
Logically Orthogonal Mesh
=========================
.. automodule:: SimPEG.Mesh.LogicallyOrthogonalMesh
:show-inheritance:
:members:
:undoc-members:
:inherited-members:
Base Mesh
=========
.. automodule:: SimPEG.Mesh.BaseMesh
:members:
:undoc-members:
Inner Products
==============
.. automodule:: SimPEG.Mesh.InnerProducts
:members:
:undoc-members:
Differential Operators
======================
.. automodule:: SimPEG.Mesh.DiffOperators
:members:
:undoc-members:
+1 -1
View File
@@ -2,7 +2,7 @@
Parameters
**********
==========
.. automodule:: SimPEG.Parameters
:show-inheritance:
-10
View File
@@ -1,10 +0,0 @@
.. _api_TensorMesh:
Tensor Mesh
***********
.. automodule:: SimPEG.Mesh.TensorMesh
:show-inheritance:
:members:
:undoc-members:
:inherited-members:
+1 -1
View File
@@ -1,7 +1,7 @@
.. _api_Tests:
Testing SimPEG
**************
==============
.. automodule:: SimPEG.Tests.TestUtils
:members:
+5 -12
View File
@@ -17,42 +17,35 @@ Utilities
:undoc-members:
Matrix Utilities
****************
================
.. automodule:: SimPEG.Utils.matutils
:members:
:undoc-members:
Sparse Utilities
****************
.. automodule:: SimPEG.Utils.sputils
:members:
:undoc-members:
LOM Utilities
*************
=============
.. automodule:: SimPEG.Utils.lomutils
:members:
:undoc-members:
Mesh Utilities
**************
==============
.. automodule:: SimPEG.Utils.meshutils
:members:
:undoc-members:
Model Builder Utilities
***********************
=======================
.. automodule:: SimPEG.Utils.ModelBuilder
:members:
:undoc-members:
Interpolation Utilities
***********************
=======================
.. automodule:: SimPEG.Utils.interputils
:members:
+15 -18
View File
@@ -11,10 +11,12 @@ The vision is to create a package for finite volume simulation and parameter est
applications to geophysical imaging and subsurface flow. To enable
these goals, this package has the following features:
* is modular with respect to discretization, physics, optimization, and regularization
* is built with the (large-scale) inverse problem in mind
* provides a framework for geophysical and hydrogeologic problems
* supports 1D, 2D and 3D problems
* is modular with respect to discretization, physics, optimization, and regularization
* is built with the (large-scale) inverse problem in mind
* provides a framework for geophysical and hydrogeologic problems
* supports 1D, 2D and 3D problems
.. image:: simpeg-framework.png
:width: 400 px
@@ -22,20 +24,15 @@ these goals, this package has the following features:
:align: center
Meshing & Operators
===================
*******************
.. toctree::
:maxdepth: 2
api_BaseMesh
api_TensorMesh
api_LogicallyOrthogonalMesh
api_Cyl1DMesh
api_DiffOperators
api_InnerProducts
api_Mesh
Forward Problems
================
****************
.. toctree::
:maxdepth: 2
@@ -43,7 +40,7 @@ Forward Problems
api_Forward
Inversion
=========
*********
.. toctree::
:maxdepth: 2
@@ -52,7 +49,7 @@ Inversion
api_Parameters
Testing SimPEG
==============
**************
.. toctree::
:maxdepth: 2
@@ -60,20 +57,20 @@ Testing SimPEG
api_Tests
* Master Branch
.. image:: https://travis-ci.org/simpeg/simpeg.png?branch=master
.. image:: https://travis-ci.org/simpeg/simpeg.png?branch*master
:target: https://travis-ci.org/simpeg/simpeg
:alt: Master Branch
:align: center
* Develop Branch
.. image:: https://travis-ci.org/simpeg/simpeg.png?branch=develop
.. image:: https://travis-ci.org/simpeg/simpeg.png?branch*develop
:target: https://travis-ci.org/simpeg/simpeg
:alt: Develop Branch
:align: center
Utility Codes
=============
*************
.. toctree::
:maxdepth: 2
@@ -82,7 +79,7 @@ Utility Codes
Project Index & Search
======================
**********************
* :ref:`genindex`
* :ref:`modindex`