mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-04 08:28:41 +08:00
Added a nodal laplacian
This commit is contained in:
Rowan Cockett
@@ -1,5 +1,8 @@
|
||||
import sys
|
||||
sys.path.append('/Users/haber/dropbox/simpegMaster/')
|
||||
try:
|
||||
import sys
|
||||
sys.path.append('/Users/haber/dropbox/simpegMaster/')
|
||||
except Exception, e:
|
||||
pass
|
||||
|
||||
import SimPEG
|
||||
import scipy.sparse as sp
|
||||
@@ -16,12 +19,12 @@ def HelmholtzSol(model,mesh,w,mbc,q,P):
|
||||
|
||||
k = np.sqrt(w**2 * mbc)
|
||||
|
||||
D1 = utils.sdiag(1./h1) * utils.ddx(mesh.nCx)
|
||||
D2 = utils.sdiag(1./h2) * utils.ddx(mesh.nCy)
|
||||
|
||||
L1 = D1.T*D1; L2 = D2.T*D2;
|
||||
print L1.todense()
|
||||
|
||||
# D1 = utils.sdiag(1./h1) * utils.ddx(mesh.nCx)
|
||||
# D2 = utils.sdiag(1./h2) * utils.ddx(mesh.nCy)
|
||||
|
||||
# L1 = - D1.T*D1
|
||||
# L2 = - D2.T*D2
|
||||
|
||||
Av = mesh.aveN2CC
|
||||
|
||||
B1 = utils.spzeros(n1+1,n1+1);
|
||||
@@ -33,21 +36,23 @@ def HelmholtzSol(model,mesh,w,mbc,q,P):
|
||||
B2[0,0] = 2*1j*k*h2[0]; B2[-1,-1] = 2*1j*k*h2[-1]
|
||||
|
||||
# generate the 2D Laplacian
|
||||
L = sp.kron(sp.identity(n2+1),L1+B1) + sp.kron(L2+B2,sp.identity(n1+1))
|
||||
|
||||
# L = sp.kron(sp.identity(n2+1),L1) + sp.kron(L2,sp.identity(n1+1))
|
||||
L = mesh.nodalLaplacian
|
||||
B = sp.kron(sp.identity(n2+1),B1) + sp.kron(B2,sp.identity(n1+1))
|
||||
L = L+B
|
||||
#plt.spy(L)
|
||||
#plt.show()
|
||||
# Generate the Mass matrix
|
||||
M = utils.sdiag(Av.T*utils.mkvc(model))
|
||||
A = L - w**2 * M
|
||||
A = - L - w**2 * M
|
||||
|
||||
mesh.ForModMat = A
|
||||
u = sp.linalg.spsolve(A,q);
|
||||
d = P*u
|
||||
return u, d
|
||||
|
||||
|
||||
def HelmholtzJmatVec(v,model,u,mesh,w,mbc,P):
|
||||
|
||||
|
||||
Cm = -w**2 * utils.sdiag(u)*mesh.aveN2CC.T
|
||||
Cu = mesh.ForModMat
|
||||
Cmv = Cm*v;
|
||||
@@ -55,20 +60,20 @@ def HelmholtzJmatVec(v,model,u,mesh,w,mbc,P):
|
||||
return -P*lam
|
||||
|
||||
def HelmholtzJTmatVec(v,model,u,mesh,w,mbc,P):
|
||||
|
||||
|
||||
Cm = -w**2 * utils.sdiag(u)*mesh.aveN2CC.T
|
||||
Cu = mesh.ForModMat
|
||||
Pv = -P.T*v
|
||||
z = sp.linalg.spsolve(Cu.T,Pv);
|
||||
return Cm.T*z
|
||||
|
||||
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
# odel,mesh,w,mbc,q
|
||||
n1 = 128; n2 = 128
|
||||
h1 = np.ones(n1); h2 = np.ones(n2);
|
||||
P = sp.identity((n1+1)*(n2+1))
|
||||
|
||||
|
||||
mesh = TensorMesh([h1,h2])
|
||||
model = np.ones(mesh.nC)
|
||||
w = 1
|
||||
@@ -76,21 +81,21 @@ if __name__ == '__main__':
|
||||
q = np.zeros((mesh.nNx,mesh.nNy))
|
||||
q[n1/2,n2/2] = 1.0
|
||||
q = q.reshape(mesh.nN,order = 'F')
|
||||
|
||||
|
||||
u, d = HelmholtzSol(model,mesh,w,mbc,q,P)
|
||||
u = u.reshape((mesh.nCx+1,mesh.nCy+1),order = 'F')
|
||||
|
||||
|
||||
plt.imshow(u.real)
|
||||
plt.show()
|
||||
|
||||
|
||||
dm = np.random.rand(mesh.nC)*1e-1+2
|
||||
u1, d1 = HelmholtzSol(model+dm,mesh,w,mbc,q,P)
|
||||
|
||||
dd = HelmholtzJmatVec(dm,model,u,mesh,w,mbc,P)
|
||||
|
||||
|
||||
print np.linalg.norm(d1-d)
|
||||
print np.linalg.norm(d1-d-dd)
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -97,6 +97,38 @@ class DiffOperators(object):
|
||||
_nodalGrad = None
|
||||
nodalGrad = property(**nodalGrad())
|
||||
|
||||
def nodalLaplacian():
|
||||
doc = "Construct laplacian operator (nodes to edges)."
|
||||
|
||||
def fget(self):
|
||||
if(self._nodalLaplacian is None):
|
||||
print 'Warning: Laplacian has not been tested rigorously.'
|
||||
# The number of cell centers in each direction
|
||||
n = self.n
|
||||
# Compute divergence operator on faces
|
||||
if(self.dim == 1):
|
||||
D1 = sdiag(1./self.hx) * ddx(mesh.nCx)
|
||||
L = - D1.T*D1
|
||||
elif(self.dim == 2):
|
||||
D1 = sdiag(1./self.hx) * ddx(n[0])
|
||||
D2 = sdiag(1./self.hy) * ddx(n[1])
|
||||
L1 = sp.kron(speye(n[1]+1), - D1.T * D1)
|
||||
L2 = sp.kron(- D2.T * D2, speye(n[0]+1))
|
||||
L = L1 + L2
|
||||
elif(self.dim == 3):
|
||||
D1 = sdiag(1./self.hx) * ddx(n[0])
|
||||
D2 = sdiag(1./self.hy) * ddx(n[1])
|
||||
D3 = sdiag(1./self.hz) * ddx(n[2])
|
||||
L1 = kron3(speye(n[2]+1), speye(n[1]+1), - D1.T * D1)
|
||||
L2 = kron3(speye(n[2]+1), - D2.T * D2, speye(n[0]+1))
|
||||
L3 = kron3(- D3.T * D3, speye(n[1]+1), speye(n[0]+1))
|
||||
L = L1 + L2 + L3
|
||||
self._nodalLaplacian = L
|
||||
return self._nodalLaplacian
|
||||
return locals()
|
||||
_nodalLaplacian = None
|
||||
nodalLaplacian = property(**nodalLaplacian())
|
||||
|
||||
def setCellGradBC(self, BC):
|
||||
"""
|
||||
Function that sets the boundary conditions for cell-centred derivative operators.
|
||||
@@ -172,7 +204,6 @@ class DiffOperators(object):
|
||||
return locals()
|
||||
cellGradx = property(**cellGradx())
|
||||
|
||||
|
||||
def cellGrady():
|
||||
doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions."
|
||||
def fget(self):
|
||||
@@ -193,8 +224,6 @@ class DiffOperators(object):
|
||||
return locals()
|
||||
cellGrady = property(**cellGrady())
|
||||
|
||||
|
||||
|
||||
def cellGradz():
|
||||
doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions."
|
||||
def fget(self):
|
||||
@@ -212,7 +241,6 @@ class DiffOperators(object):
|
||||
return locals()
|
||||
cellGradz = property(**cellGradz())
|
||||
|
||||
|
||||
def edgeCurl():
|
||||
doc = "Construct the 3D curl operator."
|
||||
|
||||
@@ -255,6 +283,8 @@ class DiffOperators(object):
|
||||
_edgeCurl = None
|
||||
edgeCurl = property(**edgeCurl())
|
||||
|
||||
# --------------- Averaging ---------------------
|
||||
|
||||
def aveF2CC():
|
||||
doc = "Construct the averaging operator on cell faces to cell centers."
|
||||
|
||||
@@ -356,6 +386,8 @@ class DiffOperators(object):
|
||||
_aveN2F = None
|
||||
aveN2F = property(**aveN2F())
|
||||
|
||||
# --------------- Methods ---------------------
|
||||
|
||||
def getMass(self, materialProp=None, loc='e'):
|
||||
""" Produces mass matricies.
|
||||
|
||||
|
||||
Reference in New Issue
Block a user