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

Merge branch 'eldadswork' of https://bitbucket.org/rcockett/simpeg into LOM

Browse Source
This commit is contained in:
Rowan Cockett
2013-08-03 13:55:04 -07:00
parent 306e31ef04 caaa7dcb41
commit 333b484cfb
4 changed files with 183 additions and 19 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/EMforward.py
+79
View File
@@ -0,0 +1,79 @@
import numpy as np
from utils import mkvc
import scipy.sparse.linalg.dsolve as dsl
from InnerProducts import getFaceInnerProduct, getEdgeInnerProduct
def getMisfit(m,mesh,forward,param):
mu0 = 4*np.pi*1e-7
omega = forward['omega'] #[param['indomega']]
rhs = forward['rhs'] #[:,param['indrhs']]
mis = 0
dmis = m*0
# Maxwell's system for E
for i in range(len(omega)):
for j in range(rhs.shape[1]):
Curl = mesh.edgeCurl
#Grad = mesh.nodalGrad
sigma = np.exp(m)
Me,PP = getEdgeInnerProduct(mesh,sigma)
Mf = 1/mu0 * getFaceInnerProduct(mesh) # assume mu = mu0
A = Curl.T * Mf * Curl - 1j * omega[i] * Me
b = mkvc(np.array(rhs[:,j]))
e = dsl.spsolve(A,b)
e = mkvc(e,2)
#print np.linalg.norm(A*e-b)/np.linalg.norm(b)
P = forward['projection']
d = P*e
r = mkvc(d - param.dobs[i,j,:],2)
mis = mis + 0.5*(r.T*r)
# get derivatives
lam = dsl.spsolve(A.T,P.T*r)
lam = mkvc(lam,2)
Gij = - 1j * omega[i] * PP.T*sp.diag((PP*e)*mesh.vol)
dmis = dmis - Gij.T*lam
return mis, dmis, d
if __name__ == '__main__':
from TensorMesh import TensorMesh
from interpmat import interpmat
from scipy import sparse as sp
h = [np.ones(7),np.ones(8),np.ones(9)]
mesh = TensorMesh(h)
xs = np.array([3.1,4.3,5.4,6.5])
ys = np.array([3.2,4.1,5.4,6.2])
zs = np.array([4.3,4.2,4.1,4.1]);
xyz = mesh.gridEx
x = xyz[:,0]; y = xyz[:,1]; z = xyz[:,2]
x = list(set(x)); y = list(set(y)); z = list(set(z))
Px = interpmat(x,y,z,xs,ys,zs)
xyz = mesh.gridEy
x = xyz[:,0]; y = xyz[:,1]; z = xyz[:,2]
x = list(set(x)); y = list(set(y)); z = list(set(z))
Py = interpmat(x,y,z,xs,ys,zs)
xyz = mesh.gridEz
x = xyz[:,0]; y = xyz[:,1]; z = xyz[:,2]
x = list(set(x)); y = list(set(y)); z = list(set(z))
Pz = interpmat(x,y,z,xs,ys,zs)
P = sp.hstack((Px,Py,Pz))
ne = np.sum(mesh.nE)
Q = np.matrix(np.random.randn(ne,5))
omega = [1,2,3]
forward = {'omega':omega, 'rhs':Q,'projection':P}
dobs = np.ones([np.size(xs),5,np.size(omega)])
param = {'dobs':dobs}
m = np.ones(mesh.nC)
getMisfit(m,mesh,forward,param)
SimPEG/InnerProducts.py
+34 -13
View File
@@ -81,23 +81,34 @@ def getFaceInnerProduct(mesh, mu):
if mu.size == mesh.nC: # Isotropic!
mu = mkvc(mu) # ensure it is a vector.
mu = sdiag(np.r_[mu, mu, mu])
Mu = sdiag(np.r_[mu, mu, mu])
elif mu.shape[1] == 3: # Diagonal tensor
mu = sdiag(np.r_[mu[:, 0], mu[:, 1], mu[:, 2]])
Mu = sdiag(np.r_[mu[:, 0], mu[:, 1], mu[:, 2]])
elif mu.shape[1] == 6: # Fully anisotropic
row1 = sp.hstack((sdiag(mu[:, 0]), sdiag(mu[:, 3]), sdiag(mu[:, 4])))
row2 = sp.hstack((sdiag(mu[:, 3]), sdiag(mu[:, 1]), sdiag(mu[:, 5])))
row3 = sp.hstack((sdiag(mu[:, 4]), sdiag(mu[:, 5]), sdiag(mu[:, 2])))
mu = sp.vstack((row1, row2, row3))
Mu = sp.vstack((row1, row2, row3))
# Cell volume
v = np.sqrt(mesh.vol)
v3 = np.r_[v, v, v]
V = sdiag(v3)*mu*sdiag(v3) # to keep symmetry
v3 = (0.125)**(0.5)*np.r_[v, v, v]
#V = sdiag(v3)*mu*sdiag(v3) # to keep symmetry
#A = P000.T*V*P000 + P001.T*V*P001 + P010.T*V*P010 + P011.T*V*P011 + P100.T*V*P100 + P101.T*V*P101 + P110.T*V*P110 + P111.T*V*P111
#A = 0.125*A
P000 = sdiag(v3)*P000; P001 = sdiag(v3)*P001; P010 = sdiag(v3)*P010; P011 = sdiag(v3)*P011
P100 = sdiag(v3)*P100; P101 = sdiag(v3)*P101; P110 = sdiag(v3)*P110; P111 = sdiag(v3)*P111
A = P000.T*Mu*P000 + P001.T*Mu*P001 + P010.T*Mu*P010 + P011.T*Mu*P011 + P100.T*Mu*P100 + P101.T*Mu*P101 + P110.T*Mu*P110 + P111.T*Mu*P111
#P = sp.vstack((sdiag(v3)*P000,sdiag(v3)*P001,sdiag(v3)*P010,sdiag(v3)*P011,
# sdiag(v3)*P100,sdiag(v3)*P101,sdiag(v3)*P110,sdiag(v3)*P111))
#A = 0.125*(P.T * sp.kron(sp.eye(8),Sigma) * P)
P = [P000,P001,P010,P011,P100,P101,P110,P111]
return A, P
A = P000.T*V*P000 + P001.T*V*P001 + P010.T*V*P010 + P011.T*V*P011 + P100.T*V*P100 + P101.T*V*P101 + P110.T*V*P110 + P111.T*V*P111
A = 0.125*A
return A
@@ -168,14 +179,22 @@ def getEdgeInnerProduct(mesh, sigma):
# Cell volume
v = np.sqrt(mesh.vol)
v3 = np.r_[v, v, v]
V = sdiag(v3)*Sigma*sdiag(v3) # to keep symmetry
v3 = (0.125)**(0.5)*np.r_[v, v, v]
P000 = sdiag(v3)*P000; P001 = sdiag(v3)*P001; P010 = sdiag(v3)*P010; P011 = sdiag(v3)*P011
P100 = sdiag(v3)*P100; P101 = sdiag(v3)*P101; P110 = sdiag(v3)*P110; P111 = sdiag(v3)*P111
A = P000.T*V*P000 + P001.T*V*P001 + P010.T*V*P010 + P011.T*V*P011 + P100.T*V*P100 + P101.T*V*P101 + P110.T*V*P110 + P111.T*V*P111
A = 0.125*A
#V = sdiag(v3)*Sigma*sdiag(v3) # to keep symmetry
return A
A = P000.T*Sigma*P000 + P001.T*Sigma*P001 + P010.T*Sigma*P010 + P011.T*Sigma*P011 + P100.T*Sigma*P100 + P101.T*Sigma*P101 + P110.T*Sigma*P110 + P111.T*Sigma*P111
#P = sp.vstack((sdiag(v3)*P000,sdiag(v3)*P001,sdiag(v3)*P010,sdiag(v3)*P011,
# sdiag(v3)*P100,sdiag(v3)*P101,sdiag(v3)*P110,sdiag(v3)*P111))
#A = 0.125*(P.T * sp.kron(sp.eye(8),Sigma) * P)
P = [P000,P001,P010,P011,P100,P101,P110,P111]
return A, P
@@ -184,4 +203,6 @@ if __name__ == '__main__':
h = [np.array([1, 2, 3, 4]), np.array([1, 2, 1, 4, 2]), np.array([1, 1, 4, 1])]
mesh = TensorMesh(h)
mu = np.ones((mesh.nC, 6))
A = mesh.getFaceInnerProduct(mu)
A = getFaceInnerProduct(mesh,mu)
B, P = getEdgeInnerProduct(mesh,mu)
SimPEG/interpmat.py
+8 -6
View File
@@ -1,6 +1,7 @@
from scipy import sparse as sp
import numpy as np
def interpmat(x,y,z,xr,yr,zr):
#
# This function does local linear interpolation
@@ -23,7 +24,8 @@ def interpmat(x,y,z,xr,yr,zr):
ind_z = np.array([0,0])
dx, dy, dz = np.zeros(2), np.zeros(2), np.zeros(2)
for i in range(0, nps):
im = np.amin(abs(xr[i]-x))
im = np.argmin(abs(xr[i]-x))
print i,im
if xr[i] - x[im] >= 0: # Point on the left
ind_x[0] = im; ind_x[1] = im+1
else: # Point on the right
@@ -33,7 +35,7 @@ def interpmat(x,y,z,xr,yr,zr):
dx[0] = xr[i] - x[ind_x[0]]
dx[1] = x[ind_x[1]] - xr[i]
im = np.amin(abs(yr[i] - y))
im = np.argmin(abs(yr[i] - y))
if yr[i] - y[im] >= 0: # Point on the left
ind_y[0] = im; ind_y[1] = im+1
else: # Point on the right
@@ -43,7 +45,7 @@ def interpmat(x,y,z,xr,yr,zr):
dy[0] = yr[i] - y[ind_y[0]]
dy[1] = y[ind_y[1]] - yr[i];
im = np.amin(abs(zr[i] - z));
im = np.argmin(abs(zr[i] - z));
if zr[i] -z[im] >= 0: # Point on the left
ind_z[0] = im; ind_z[1] = im+1
else: # Point on the right
@@ -80,9 +82,9 @@ def interpmat(x,y,z,xr,yr,zr):
if __name__ == '__main__':
x = np.array([1, 2, 3, 4])
y = np.array([1, 2, 3, 4, 5])
z = np.array([0, 1, 4, 6])
x = np.array([1.1, 2.1, 3.6, 4.9])
y = np.array([1.2, 2.2, 3.3, 4.9, 5.6])
z = np.array([0.8, 1.7, 4.9, 6.5])
xr = np.array([2.5,3.2])
yr = np.array([2.4,3.6])
SimPEG/ops.py
+62
View File
@@ -0,0 +1,62 @@
import numpy as np
from utils import mkvc
import scipy.sparse.linalg as spla
import scipy.sparse as sp
def matmul(A,B):
# first check shape
if np.shape(A)[1] != np.shape(B)[0]:
print 'error in sizes'
return
# Check types
sA = sp.issparse(A)
sB = sp.issparse(B)
if ((sA == False) & (sB == True)): # doesno't work unless we trick it
return (B.T.dot(A.T)).T
else:
return A.dot(B)
def dot(A,B):
A = mkvc(A,1)
B = mkvc(B,1)
return np.dot(A,B)
def inner(A,B):
A = mkvc(A,1)
B = mkvc(B,1)
return np.dot(A,B)
if __name__ == '__main__':
import numpy as np
from utils import mkvc
import scipy.sparse as sp
# generate sparse and dense matrices
A = sp.rand(100, 200, density=0.05, format='csr', dtype=None)
B = sp.rand(200, 150, density=0.05, format='csr', dtype=None)
C = np.random.rand(200,150)
D = np.random.rand(150,100)
b = mkvc(np.arange(200),1)
c = np.reshape(b,(1,200))
matmul(A,B)
matmul(A,C)
matmul(C,D)
matmul(D,A)
matmul(A,b)
dot(c,b)
dot(C,C)
print np.shape(c), np.shape(b)[0]
print matmul(c,b),dot(c,b)