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

Testing differential operators (Div, Grad, Curl)

Browse Source
This commit is contained in:
SEOGI KANG
2013-07-17 15:03:56 -07:00
parent 4e1e56b15a
commit 082df8f110
5 changed files with 319 additions and 16 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/getDiffop.py
+198
View File
@@ -0,0 +1,198 @@
import numpy as np
from scipy import sparse
from utils import mkvc
from sputils import ddx, sdiag, speye, kron3, spzeros, av
def getvol(h):
"""Construct cell volumes of the 3D model as 1d array."""
# Cell sizes in each direction
h1 = h[0]
h2 = h[1]
h3 = h[2]
# Compute cell volumes
v12 = h1.T*h2
V = mkvc(v12.reshape(-1,1)*h3)
return V
def getarea(h):
"""Construct face areas of the 3D model as 1d array."""
# Cell sizes in each direction
h1 = h[0]
h2 = h[1]
h3 = h[2]
# The number of cell centers in each direction
n1 = np.size(h1)
n2 = np.size(h2)
n3 = np.size(h3)
# Compute areas of cell faces
area1 = np.ones((n1+1,1))*mkvc(h2.T*h3)
area2 = h1.T*mkvc(np.ones((n2+1,1))*h3)
area3 = h1.T*mkvc(h2.T*np.ones(n3+1))
area = np.concatenate((mkvc(area1), mkvc(area2), mkvc(area3)), axis=0)
return area
def getlength_e(h):
"""Construct edge legnths of the 3D model as 1d array."""
# Cell sizes in each direction
h1 = h[0]
h2 = h[1]
h3 = h[2]
# The number of cell centers in each direction
n1 = np.size(h1)
n2 = np.size(h2)
n3 = np.size(h3)
# Compute areas of cell faces
l1 = h1.T*mkvc(np.ones((n2+1,1))*np.ones(n3+1))
l2 = np.ones((n1+1,1))*mkvc(h2.T*np.ones(n3+1))
l3 = np.ones((n1+1,1))*mkvc(np.ones((n2+1,1))*h3)
#l = np.hstack((np.hstack((mkvc(area1), mkvc(area2))), mkvc(area3)))
l = np.concatenate((mkvc(l1), mkvc(l2), mkvc(l3)), axis=0)
return l
def getDivMatrix(h):
"""Construct the 3D divergence operator on Faces."""
# Cell sizes in each direction
h1 = h[0]
h2 = h[1]
h3 = h[2]
# The number of cell centers in each direction
n1 = np.size(h1)
n2 = np.size(h2)
n3 = np.size(h3)
# Compute areas of cell faces
S = getarea(h)
# Compute cell volumes
V = getvol(h)
# Compute divergence operator on faces
d1 = ddx(n1)
d2 = ddx(n2)
d3 = ddx(n3)
D1 = kron3(speye(n3), speye(n2), d1)
D2 = kron3(speye(n3), d2, speye(n1))
D3 = kron3(d3, speye(n2), speye(n1))
D = sparse.hstack((D1, D2, D3), format="csr")
return sdiag(1/V)*D*sdiag(S)
def getGradMatrix(h):
"""Construct the 3D nodal gradient operator."""
# Cell sizes in each direction
h1 = h[0]
h2 = h[1]
h3 = h[2]
# The number of cell centers in each direction
n1 = np.size(h1)
n2 = np.size(h2)
n3 = np.size(h3)
# Compute lengths of cell edges
L = getlength_e(h)
# Compute divergence operator on faces
d1 = ddx(n1)
d2 = ddx(n2)
d3 = ddx(n3)
D1 = kron3(speye(n3+1), speye(n2+1), d1)
D2 = kron3(speye(n3+1), d2, speye(n1+1))
D3 = kron3(d3, speye(n2+1), speye(n1+1))
G = sparse.vstack((D1, D2, D3), format="csr")
return sdiag(1/L)*G
def getCurlMatrix(h):
"""Construct the 3D curl operator."""
# Cell sizes in each direction
h1 = h[0]
h2 = h[1]
h3 = h[2]
# The number of cell centers in each direction
n1 = np.size(h1)
n2 = np.size(h2)
n3 = np.size(h3)
# Compute lengths of cell edges
L = getlength_e(h)
# Compute areas of cell faces
S = getarea(h)
# Compute divergence operator on faces
d1 = ddx(n1)
d2 = ddx(n2)
d3 = ddx(n3)
D32 = kron3(d3, speye(n2), speye(n1+1))
D23 = kron3(speye(n3), d2, speye(n1+1))
D31 = kron3(d3, speye(n2+1), speye(n1))
D13 = kron3(speye(n3), speye(n2+1), d1)
D21 = kron3(speye(n3+1), d2, speye(n1))
D12 = kron3(speye(n3+1), speye(n2), d1)
O1 = spzeros(np.shape(D32)[0], np.shape(D31)[1])
O2 = spzeros(np.shape(D31)[0], np.shape(D32)[1])
O3 = spzeros(np.shape(D21)[0], np.shape(D13)[1])
C = sparse.vstack((sparse.hstack((O1,-D32, D23)),
sparse.hstack((D31,O2, -D13)),
sparse.hstack((-D21,D12, O3))), format="csr")
return sdiag(1/S)*(C*sdiag(L))
def getAverageMatrixF(h):
"""Construct the 3D averaging operator on cell faces."""
# Cell sizes in each direction
h1 = h[0]
h2 = h[1]
h3 = h[2]
# The number of cell centers in each direction
n1 = np.size(h1)
n2 = np.size(h2)
n3 = np.size(h3)
av1 = av(n1)
av2 = av(n2)
av3 = av(n3)
AvF = sparse.hstack(kron3(speye(n3), speye(n2), av1),
kron3(speye(n3), av2, speye(n3)),
kron3(av3, speye(n2), speye(n3)), format="csr")
return AvF
def getAverageMatrixE(h):
"""Construct the 3D averaging operator on cell edges."""
# Cell sizes in each direction
h1 = h[0]
h2 = h[1]
h3 = h[2]
# The number of cell centers in each direction
n1 = np.size(h1)
n2 = np.size(h2)
n3 = np.size(h3)
av1 = av(n1)
av2 = av(n2)
av3 = av(n3)
AvE = sparse.hstack(kron3(av3, av2, speye(n1)),
kron3(av3, speye(n2), av1),
kron3(speye(n3), av2, av1), format="csr")
return AvE
SimPEG/sputils.py
+12 -4
View File
@@ -3,16 +3,24 @@ from scipy import sparse
def ddx(n):
"""Define 1D derivatives"""
return sparse.spdiags((np.ones((n+1,1))*[-1,1]).T, [0,1], n, n+1)
return sparse.spdiags((np.ones((n+1,1))*[-1,1]).T, [0,1], n, n+1, format="csr")
def sdiag(h):
"""Sparse diagonal matrix"""
return sparse.spdiags(h, 0, np.size(h), np.size(h))
return sparse.spdiags(h, 0, np.size(h), np.size(h), format="csr")
def speye(n):
"""Sparse identity"""
return sparse.identity(n)
return sparse.identity(n, format="csr")
def kron3(A, B, C):
"""Two kron prods"""
return sparse.kron(sparse.kron(A, B), C)
return sparse.kron(sparse.kron(A, B), C, format="csr")
def spzeros(n1, n2):
"""spzeros"""
return sparse.coo_matrix((n1, n2)).tocsr()
def av(n):
"""Define 1D averaging operator"""
return sparse.spdiags((0.5*np.ones((n+1,1))*[1,1]).T, [0,1], n, n+1, format="csr")
SimPEG/tests/test_curl.py
+48
View File
@@ -0,0 +1,48 @@
import numpy as np
import sys
sys.path.append('../')
from TensorMesh import TensorMesh
from getDiffop import getCurlMatrix
err=0.
print '>> Test Curl operator'
for i in range(4):
icount=i+1
nc = 2**icount
# Define the mesh
h1 = np.ones((1,nc))/nc
h2 = np.ones((1,nc))/nc
h3 = np.ones((1,nc))/nc
h = [h1, h2, h3]
x0 = np.zeros((3, 1))
M = TensorMesh(h, x0)
#n = M.plotGrid()
# Generate DIV matrix
CURL = getCurlMatrix(h)
#Test function
fun = lambda x: np.cos(x) # i (cos(y)) + j (cos(z)) + k (cos(x))
sol = lambda x: np.sin(x) # i (sin(z)) + j (sin(x)) + k (sin(y))
Ex = fun(M.gridEx[:,1])
Ey = fun(M.gridEy[:,2])
Ez = fun(M.gridEz[:,0])
E = np.concatenate((Ex,Ey,Ez))
Fx = sol(M.gridFx[:,2])
Fy = sol(M.gridFy[:,0])
Fz = sol(M.gridFz[:,1])
curlE_anal = np.concatenate((Fx,Fy,Fz))
curlE = CURL*E
err = np.linalg.norm((curlE-curlE_anal), np.inf)
if icount == 1:
print 'h | inf norm | error ratio'
print '---------------------------------------'
print '%6.4f | %8.2e |'% (h1[0,0], err)
else:
print '%6.4f | %8.2e | %6.4f' % (h1[0,0], err, err_old/err)
err_old = err
SimPEG/tests/test_div.py
+15 -12
View File
@@ -5,16 +5,18 @@ sys.path.append('../')
from TensorMesh import TensorMesh
from getDIV import getDivMatrix, getarea, getvol
# Define the mesh
err=0.
print '>> Test face Divergence operator'
for i in range(4):
icount=i+1;
nc = 2*icount;
h1 = np.pi/nc*np.ones((1,nc))
h2 = np.pi/nc*np.ones((1,nc))
h3 = np.pi/nc*np.ones((1,nc))
icount=i+1
nc = 2**icount
# Define the mesh
h1 = np.ones((1,nc))/nc
h2 = np.ones((1,nc))/nc
h3 = np.ones((1,nc))/nc
h = [h1, h2, h3]
x0 = -np.pi/2*np.ones((3, 1))
x0 = np.zeros((3, 1))
M = TensorMesh(h, x0)
#n = M.plotGrid()
@@ -34,12 +36,13 @@ for i in range(4):
area = getarea(h)
vol = getvol(h)
err = np.linalg.norm((divF-divF_anal)*np.sqrt(vol), 2)
#err = np.linalg.norm((divF-divF_anal)*np.sqrt(vol), 2)
err = np.linalg.norm((divF-divF_anal), np.inf)
if icount == 1:
err1 = err
print 'h | 2 norm | error ratio'
print 'h | inf norm | error ratio'
print '---------------------------------------'
print '%6.4f | %8.2e |'% (h1[0,0], err)
else:
print '%6.4f | %8.2e | %6.4f' % (h1[0,0], err, err1/err)
print '%6.4f | %8.2e | %6.4f' % (h1[0,0], err, err_old/err)
err_old = err
SimPEG/tests/test_grad.py
+46
View File
@@ -0,0 +1,46 @@
import numpy as np
import sys
sys.path.append('../')
from TensorMesh import TensorMesh
from getDiffop import getGradMatrix
err=0.
print '>> Test nodal Gradient operator'
for i in range(4):
icount=i+1
nc = 2**icount
# Define the mesh
h1 = np.ones((1,nc))/nc
h2 = np.ones((1,nc))/nc
h3 = np.ones((1,nc))/nc
h = [h1, h2, h3]
x0 = np.zeros((3, 1))
M = TensorMesh(h, x0)
#n = M.plotGrid()
# Generate DIV matrix
GRAD = getGradMatrix(h)
#Test function
fun = lambda x, y, z: (np.cos(x)+np.cos(y)+np.cos(z))
sol = lambda x: -np.sin(x) # i (sin(x)) + j (sin(y)) + k (sin(z))
phi = fun(M.gridN[:,0], M.gridN[:,1], M.gridN[:,2])
gradE = GRAD*phi
Ex = sol(M.gridEx[:,0])
Ey = sol(M.gridEy[:,1])
Ez = sol(M.gridEz[:,2])
gradE_anal = np.concatenate((Ex,Ey,Ez))
err = np.linalg.norm((gradE-gradE_anal), np.inf)
if icount == 1:
print 'h | inf norm | error ratio'
print '---------------------------------------'
print '%6.4f | %8.2e |'% (h1[0,0], err)
else:
print '%6.4f | %8.2e | %6.4f' % (h1[0,0], err, err_old/err)
err_old = err