mirror of
https://github.com/wassname/simpeg.git
synced 2026-06-28 09:13:21 +08:00
Removed EMforward code for merging into master. (This is still being developed in branch: eldadsWork)
This commit is contained in:
Rowan Cockett
@@ -1,79 +0,0 @@
|
||||
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)
|
||||
@@ -1,94 +0,0 @@
|
||||
from scipy import sparse as sp
|
||||
import numpy as np
|
||||
|
||||
|
||||
def interpmat(x,y,z,xr,yr,zr):
|
||||
#
|
||||
# This function does local linear interpolation
|
||||
# computed for each receiver point in turn
|
||||
#
|
||||
# [Q] = linint(x,y,z,xr,yr,zr)
|
||||
# Interpolation matrix
|
||||
#
|
||||
|
||||
nx = np.size(x)
|
||||
ny = np.size(y)
|
||||
nz = np.size(z)
|
||||
|
||||
nps = np.size(xr)
|
||||
|
||||
#Q = spalloc(np,nx*ny*nz,8*np);
|
||||
Q = sp.lil_matrix((nps,nx*ny*nz))
|
||||
ind_x = np.array([0,0])
|
||||
ind_y = np.array([0,0])
|
||||
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.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
|
||||
ind_x[0] = im-1; ind_x[1] = im
|
||||
|
||||
|
||||
dx[0] = xr[i] - x[ind_x[0]]
|
||||
dx[1] = x[ind_x[1]] - xr[i]
|
||||
|
||||
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
|
||||
ind_y[0] = im-1; ind_y[1] = im
|
||||
|
||||
|
||||
dy[0] = yr[i] - y[ind_y[0]]
|
||||
dy[1] = y[ind_y[1]] - yr[i];
|
||||
|
||||
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
|
||||
ind_z[0] = im-1; ind_z[1] = im;
|
||||
|
||||
dz[0] = zr[i] - z[ind_z[0]]; dz[1] = z[ind_z[1]] - zr[i]
|
||||
|
||||
Dx = x[ind_x[1]] - x[ind_x[0]]
|
||||
Dy = y[ind_y[1]] - y[ind_y[0]]
|
||||
Dz = z[ind_z[1]] - z[ind_z[0]]
|
||||
#dv = Dx*Dy*Dz
|
||||
|
||||
# Get the row in the matrix
|
||||
v = np.zeros([nx, ny,nz])
|
||||
|
||||
v[ ind_x[0], ind_y[0], ind_z[0]] = (1-dx[0]/Dx)*(1-dy[0]/Dy)*(1-dz[0]/Dz)
|
||||
v[ ind_x[0], ind_y[1], ind_z[0]] = (1-dx[0]/Dx)*(1-dy[1]/Dy)*(1-dz[0]/Dz)
|
||||
v[ ind_x[1], ind_y[0], ind_z[0]] = (1-dx[1]/Dx)*(1-dy[0]/Dy)*(1-dz[0]/Dz)
|
||||
v[ ind_x[1], ind_y[1], ind_z[0]] = (1-dx[1]/Dx)*(1-dy[1]/Dy)*(1-dz[0]/Dz)
|
||||
v[ ind_x[0], ind_y[0], ind_z[1]] = (1-dx[0]/Dx)*(1-dy[0]/Dy)*(1-dz[1]/Dz)
|
||||
v[ ind_x[0], ind_y[1], ind_z[1]] = (1-dx[0]/Dx)*(1-dy[1]/Dy)*(1-dz[1]/Dz)
|
||||
v[ ind_x[1], ind_y[0], ind_z[1]] = (1-dx[1]/Dx)*(1-dy[0]/Dy)*(1-dz[1]/Dz)
|
||||
v[ ind_x[1], ind_y[1], ind_z[1]] = (1-dx[1]/Dx)*(1-dy[1]/Dy)*(1-dz[1]/Dz)
|
||||
|
||||
|
||||
print(np.shape(v.flatten('F')))
|
||||
print(np.shape(Q))
|
||||
|
||||
Q[i,:] = v.flatten('F')
|
||||
|
||||
|
||||
return Q.tocsr()
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
|
||||
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])
|
||||
zr = np.array([2.5,3.9])
|
||||
|
||||
A = interpmat(x,y,z,xr,yr,zr)
|
||||
|
||||
@@ -1,62 +0,0 @@
|
||||
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)
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
Reference in New Issue
Block a user