diff --git a/SimPEG/EMforward.py b/SimPEG/EMforward.py new file mode 100644 index 00000000..9c332c9b --- /dev/null +++ b/SimPEG/EMforward.py @@ -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) diff --git a/SimPEG/InnerProducts.py b/SimPEG/InnerProducts.py index bc482de4..5e53961f 100644 --- a/SimPEG/InnerProducts.py +++ b/SimPEG/InnerProducts.py @@ -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) + diff --git a/SimPEG/interpmat.py b/SimPEG/interpmat.py index 68b6c4c1..bfb80291 100644 --- a/SimPEG/interpmat.py +++ b/SimPEG/interpmat.py @@ -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]) diff --git a/SimPEG/ops.py b/SimPEG/ops.py new file mode 100644 index 00000000..4997cde9 --- /dev/null +++ b/SimPEG/ops.py @@ -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) + + + + + + +