From dcae030d811b994b1d2786f5ab849f05d9bf7c7e Mon Sep 17 00:00:00 2001 From: ehaber99 Date: Fri, 2 Aug 2013 00:18:46 -0700 Subject: [PATCH] 1. Added output for the edgeInnerProduct so we can use it in the derivative phase 2. The getMisfit may be working but needs a little bit more work We need to go through our code to make sure that all vectors that can be multiplied by a matrix are set to vec = mkvc(vec,2) --- SimPEG/EMforward.py | 47 ++++++++++++++++++++++++++++++++--------- SimPEG/InnerProducts.py | 12 ++++++----- 2 files changed, 44 insertions(+), 15 deletions(-) diff --git a/SimPEG/EMforward.py b/SimPEG/EMforward.py index fc143841..0f65cb21 100644 --- a/SimPEG/EMforward.py +++ b/SimPEG/EMforward.py @@ -1,13 +1,15 @@ 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): +def getMisfit(m,mesh,forward,param): mu0 = 4*np.pi*1e-7 omega = forward['omega'] #[param['indomega']] rhs = forward['rhs'] #[:,param['indrhs']] - misfit = 0 + mis = 0 + dmis = m*0 # Maxwell's system for E for i in range(len(omega)): @@ -15,8 +17,8 @@ def getMisfit(m,mesh,forward): Curl = mesh.edgeCurl #Grad = mesh.nodalGrad sigma = np.exp(m) - Me,PP = mesh.getEdgeMass(sigma) - Mf = 1/mu0 * mesh.getFaceMass() # assume mu = mu0 + 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])) @@ -30,20 +32,45 @@ def getMisfit(m,mesh,forward): mis = mis + 0.5*(r.T*r) # get derivatives lam = dsl.spsolve(A.T,P.T*r) - Gij = PP.T*diag((PP*e)*mesh.vol) + 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) - ne = np.sum(mesh.nE) - Q = np.matrix(np.random.randn(ne,5)) - P = np.matrix(Q.T) - forward = {'omega':[1,2,3], 'rhs':Q,'projection':P} + 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] + Px = interpmat(x,y,z,xs,ys,zs) + xyz = mesh.gridEy + x = xyz[:,0]; y = xyz[:,1]; z = xyz[:,2] + Py = interpmat(x,y,z,xs,ys,zs) + xyz = mesh.gridEz + x = xyz[:,0]; y = xyz[:,1]; z = xyz[:,2] + 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),np.shape(Q,2),np.size(omega)) + param = {'dobs':dobs} + + m = np.ones(mesh.nC) - getMisfit(m,mesh,forward) + getMisfit(m,mesh,forward,param) diff --git a/SimPEG/InnerProducts.py b/SimPEG/InnerProducts.py index bc482de4..d4b2e9ad 100644 --- a/SimPEG/InnerProducts.py +++ b/SimPEG/InnerProducts.py @@ -169,13 +169,15 @@ 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 + #V = sdiag(v3)*Sigma*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 = 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 + #A = 0.125*A + 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*Sigma*P) + return A, P