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

Update tests with new calls. Note that there seem to be problems with rotateLOM on some operator tests.

Browse Source
This commit is contained in:
Rowan Cockett
2013-08-06 15:09:01 -07:00
parent a510755926
commit 262aacce85
2 changed files with 79 additions and 95 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/tests/test_massMatrices.py
+22 -46
View File
@@ -64,33 +64,21 @@ class TestInnerProducts(OrderTest):
analytic = 69881./21600 # Found using matlab symbolic toolbox.
if self.location == 'edges':
if self.M._meshType == 'TENSOR':
Ex = call(ex, self.M.gridEx)
Ey = call(ey, self.M.gridEy)
Ez = call(ez, self.M.gridEz)
E = np.matrix(np.r_[Ex, Ey, Ez]).T
elif self.M._meshType == 'LOM':
cart = lambda g: np.c_[call(ex, g), call(ey, g), call(ez, g)]
Ec = np.vstack((cart(self.M.gridEx),
cart(self.M.gridEy),
cart(self.M.gridEz)))
E = np.matrix(self.M.projectEdgeVector(Ec))
cart = lambda g: np.c_[call(ex, g), call(ey, g), call(ez, g)]
Ec = np.vstack((cart(self.M.gridEx),
cart(self.M.gridEy),
cart(self.M.gridEz)))
E = self.M.projectEdgeVector(Ec)
A = self.M.getEdgeInnerProduct(sigma)
numeric = E.T*A*E
numeric = E.T.dot(A.dot(E))
elif self.location == 'faces':
if self.M._meshType == 'TENSOR':
Fx = call(ex, self.M.gridFx)
Fy = call(ey, self.M.gridFy)
Fz = call(ez, self.M.gridFz)
F = np.matrix(np.r_[Fx, Fy, Fz]).T
elif self.M._meshType == 'LOM':
cart = lambda g: np.c_[call(ex, g), call(ey, g), call(ez, g)]
Fc = np.vstack((cart(self.M.gridFx),
cart(self.M.gridFy),
cart(self.M.gridFz)))
F = np.matrix(self.M.projectFaceVector(Fc))
cart = lambda g: np.c_[call(ex, g), call(ey, g), call(ez, g)]
Fc = np.vstack((cart(self.M.gridFx),
cart(self.M.gridFy),
cart(self.M.gridFz)))
F = self.M.projectFaceVector(Fc)
A = self.M.getFaceInnerProduct(sigma)
numeric = F.T*A*F
numeric = F.T.dot(A.dot(F))
err = np.abs(numeric - analytic)
return err
@@ -164,31 +152,19 @@ class TestInnerProducts2D(OrderTest):
analytic = 781427./360 # Found using matlab symbolic toolbox. z=5
if self.location == 'edges':
if self.M._meshType == 'TENSOR':
Ex = call(ex, self.M.gridEx)
Ey = call(ey, self.M.gridEy)
E = np.matrix(np.r_[Ex, Ey]).T
elif self.M._meshType == 'LOM':
cart = lambda g: np.c_[call(ex, g), call(ey, g)]
Ec = np.vstack((cart(self.M.gridEx),
cart(self.M.gridEy)))
E = np.matrix(self.M.projectEdgeVector(Ec))
cart = lambda g: np.c_[call(ex, g), call(ey, g)]
Ec = np.vstack((cart(self.M.gridEx),
cart(self.M.gridEy)))
E = self.M.projectEdgeVector(Ec)
A = self.M.getEdgeInnerProduct(sigma)
numeric = E.T*A*E
numeric = E.T.dot(A.dot(E))
elif self.location == 'faces':
if self.M._meshType == 'TENSOR':
Fx = call(ex, self.M.gridFx)
Fy = call(ey, self.M.gridFy)
F = np.matrix(np.r_[Fx, Fy]).T
elif self.M._meshType == 'LOM':
cart = lambda g: np.c_[call(ex, g), call(ey, g)]
Fc = np.vstack((cart(self.M.gridFx),
cart(self.M.gridFy)))
F = np.matrix(self.M.projectFaceVector(Fc))
cart = lambda g: np.c_[call(ex, g), call(ey, g)]
Fc = np.vstack((cart(self.M.gridFx),
cart(self.M.gridFy)))
F = self.M.projectFaceVector(Fc)
A = self.M.getFaceInnerProduct(sigma)
numeric = F.T*A*F
numeric = F.T.dot(A.dot(F))
err = np.abs(numeric - analytic)
return err
SimPEG/tests/test_operators.py
+57 -49
View File
@@ -4,7 +4,15 @@ import sys
sys.path.append('../')
from OrderTest import OrderTest
MESHTYPES = ['uniformTensorMesh', 'uniformLOM'] # , 'rotateLOM'
MESHTYPES = ['uniformTensorMesh', 'uniformLOM', 'rotateLOM']
call2 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1])
call3 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2])
cart_row2 = lambda g, xfun, yfun: np.c_[call2(xfun, g), call2(yfun, g)]
cart_row3 = lambda g, xfun, yfun, zfun: np.c_[call3(xfun, g), call3(yfun, g), call3(zfun, g)]
cartF2 = lambda M, fx, fy: np.vstack((cart_row2(M.gridFx, fx, fy), cart_row2(M.gridFy, fx, fy)))
cartE2 = lambda M, ex, ey: np.vstack((cart_row2(M.gridEx, ex, ey), cart_row2(M.gridEy, ex, ey)))
cartF3 = lambda M, fx, fy, fz: np.vstack((cart_row3(M.gridFx, fx, fy, fz), cart_row3(M.gridFy, fx, fy, fz), cart_row3(M.gridFz, fx, fy, fz)))
cartE3 = lambda M, ex, ey, ez: np.vstack((cart_row3(M.gridEx, ex, ey, ez), cart_row3(M.gridEy, ex, ey, ez), cart_row3(M.gridEz, ex, ey, ez)))
class TestCurl(OrderTest):
@@ -12,24 +20,26 @@ class TestCurl(OrderTest):
meshTypes = MESHTYPES
def getError(self):
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))
# fun: i (cos(y)) + j (cos(z)) + k (cos(x))
# sol: i (sin(z)) + j (sin(x)) + k (sin(y))
Ex = fun(self.M.gridEx[:, 1])
Ey = fun(self.M.gridEy[:, 2])
Ez = fun(self.M.gridEz[:, 0])
E = np.concatenate((Ex, Ey, Ez))
funX = lambda x, y, z: np.cos(y)
funY = lambda x, y, z: np.cos(z)
funZ = lambda x, y, z: np.cos(x)
Fx = sol(self.M.gridFx[:, 2])
Fy = sol(self.M.gridFy[:, 0])
Fz = sol(self.M.gridFz[:, 1])
curlE_anal = np.concatenate((Fx, Fy, Fz))
solX = lambda x, y, z: np.sin(z)
solY = lambda x, y, z: np.sin(x)
solZ = lambda x, y, z: np.sin(y)
Ec = cartE3(self.M, funX, funY, funZ)
E = self.M.projectEdgeVector(Ec)
Fc = cartF3(self.M, solX, solY, solZ)
curlE_anal = self.M.projectFaceVector(Fc)
# Generate DIV matrix
CURL = self.M.edgeCurl
curlE = CURL*E
err = np.linalg.norm((curlE-curlE_anal), np.inf)
curlE = self.M.edgeCurl.dot(E)
err = np.linalg.norm((curlE - curlE_anal), np.inf)
return err
def test_order(self):
@@ -41,18 +51,17 @@ class TestFaceDiv(OrderTest):
meshTypes = MESHTYPES
def getError(self):
DIV = self.M.faceDiv
#Test function
fun = lambda x: np.sin(x)
Fx = fun(self.M.gridFx[:, 0])
Fy = fun(self.M.gridFy[:, 1])
Fz = fun(self.M.gridFz[:, 2])
F = np.concatenate((Fx, Fy, Fz))
divF = DIV*F
fx = lambda x, y, z: np.sin(x)
fy = lambda x, y, z: np.sin(y)
fz = lambda x, y, z: np.sin(z)
sol = lambda x, y, z: (np.cos(x)+np.cos(y)+np.cos(z))
divF_anal = sol(self.M.gridCC[:, 0], self.M.gridCC[:, 1], self.M.gridCC[:, 2])
Fc = cartF3(self.M, fx, fy, fz)
F = self.M.projectFaceVector(Fc)
divF = self.M.faceDiv.dot(F)
divF_anal = call3(sol, self.M.gridCC)
err = np.linalg.norm((divF-divF_anal), np.inf)
@@ -68,17 +77,16 @@ class TestFaceDiv2D(OrderTest):
meshDimension = 2
def getError(self):
DIV = self.M.faceDiv
#Test function
fun = lambda x: np.sin(x)
Fx = fun(self.M.gridFx[:, 0])
Fy = fun(self.M.gridFy[:, 1])
F = np.concatenate((Fx, Fy))
divF = DIV*F
fx = lambda x, y: np.sin(x)
fy = lambda x, y: np.sin(y)
sol = lambda x, y: (np.cos(x)+np.cos(y))
divF_anal = sol(self.M.gridCC[:, 0], self.M.gridCC[:, 1])
Fc = cartF2(self.M, fx, fy)
F = self.M.projectFaceVector(Fc)
divF = self.M.faceDiv.dot(F)
divF_anal = call2(sol, self.M.gridCC)
err = np.linalg.norm((divF-divF_anal), np.inf)
@@ -93,19 +101,19 @@ class TestNodalGrad(OrderTest):
meshTypes = MESHTYPES
def getError(self):
GRAD = self.M.nodalGrad
#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))
# i (sin(x)) + j (sin(y)) + k (sin(z))
solX = lambda x, y, z: -np.sin(x)
solY = lambda x, y, z: -np.sin(y)
solZ = lambda x, y, z: -np.sin(z)
phi = fun(self.M.gridN[:, 0], self.M.gridN[:, 1], self.M.gridN[:, 2])
gradE = GRAD*phi
phi = call3(fun, self.M.gridN)
gradE = self.M.nodalGrad.dot(phi)
Ex = sol(self.M.gridEx[:, 0])
Ey = sol(self.M.gridEy[:, 1])
Ez = sol(self.M.gridEz[:, 2])
Ec = cartE3(self.M, solX, solY, solZ)
gradE_anal = self.M.projectEdgeVector(Ec)
gradE_anal = np.concatenate((Ex, Ey, Ez))
err = np.linalg.norm((gradE-gradE_anal), np.inf)
return err
@@ -120,18 +128,18 @@ class TestNodalGrad2D(OrderTest):
meshDimension = 2
def getError(self):
GRAD = self.M.nodalGrad
#Test function
fun = lambda x, y: (np.cos(x)+np.cos(y))
sol = lambda x: -np.sin(x) # i (sin(x)) + j (sin(y)) + k (sin(z))
# i (sin(x)) + j (sin(y)) + k (sin(z))
solX = lambda x, y: -np.sin(x)
solY = lambda x, y: -np.sin(y)
phi = fun(self.M.gridN[:, 0], self.M.gridN[:, 1])
gradE = GRAD*phi
phi = call2(fun, self.M.gridN)
gradE = self.M.nodalGrad.dot(phi)
Ex = sol(self.M.gridEx[:, 0])
Ey = sol(self.M.gridEy[:, 1])
Ec = cartE2(self.M, solX, solY)
gradE_anal = self.M.projectEdgeVector(Ec)
gradE_anal = np.concatenate((Ex, Ey))
err = np.linalg.norm((gradE-gradE_anal), np.inf)
return err