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

Created an OrderTest class

Browse Source 
This creates a testable framework so that we can easily test order of convergence on things like operators by just writing the guts of the code.

Merged Seogi's code into this framework.
This commit is contained in:
Rowan Cockett
2013-07-22 12:26:17 -07:00
parent a6f2a512ae
commit d0efdb509c
5 changed files with 138 additions and 135 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/tests/OrderTest.py
+50
View File
@@ -0,0 +1,50 @@
import sys
sys.path.append('../../')
from SimPEG import TensorMesh
import numpy as np
import unittest
class OrderTest(unittest.TestCase):
"""Order test sets up the basics for testing order of decrease for a function on a mesh."""
name = "Order Test"
expectedOrder = 2
meshSizes = [4, 8, 16, 32]
def setupMesh(self, nc):
# Define the mesh
h1 = np.ones(nc)/nc
h2 = np.ones(nc)/nc
h3 = np.ones(nc)/nc
h = [h1, h2, h3]
self.M = TensorMesh(h)
def getError(self):
"""Overwrite this function with the guts of the test."""
return 1.
def orderTest(self):
order = []
err_old = 0.
nc_old = 0.
for ii, nc in enumerate(self.meshSizes):
self.setupMesh(nc)
err = self.getError()
if ii == 0:
print ''
print 'Testing order of: ' + self.name
print ' h | inf norm | ratio | order'
print '------------------------------------------'
print '%4i | %8.2e |' % (nc, err)
else:
order.append(np.log(err/err_old)/np.log(float(nc_old)/float(nc)))
print '%4i | %8.2e | %6.4f | %6.4f' % (nc, err, err_old/err, order[-1])
err_old = err
nc_old = nc
self.assertTrue(len(np.where(np.array(order) > 0.9*self.expectedOrder)[0]) > np.floor(0.75*len(order)))
if __name__ == '__main__':
unittest.main()
SimPEG/tests/test_curl.py
-46
View File
@@ -1,46 +0,0 @@
import numpy as np
import sys
sys.path.append('../')
from TensorMesh import TensorMesh
err=0.
print '>> Test Curl operator'
for i in range(4):
icount=i+1
nc = 2**icount
# Define the mesh
h1 = np.ones(nc)/nc
h2 = np.ones(nc)/nc
h3 = np.ones(nc)/nc
h = [h1, h2, h3]
M = TensorMesh(h)
#n = M.plotGrid()
# Generate DIV matrix
CURL = M.edgeCurl
#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], err)
else:
print '%6.4f | %8.2e | %6.4f' % (h1[0], err, err_old/err)
err_old = err
SimPEG/tests/test_div.py
-45
View File
@@ -1,45 +0,0 @@
import numpy as np
import sys
sys.path.append('../')
from TensorMesh import TensorMesh
err=0.
print '>> Test face Divergence operator'
for i in range(4):
icount=i+1
nc = 2**icount
# Define the mesh
h1 = np.ones(nc)/nc
h2 = np.ones(nc)/nc
h3 = np.ones(nc)/nc
h = [h1, h2, h3]
x0 = np.zeros((3, 1))
M = TensorMesh(h, x0)
#n = M.plotGrid()
# Generate DIV matrix
DIV = M.faceDiv
#Test function
fun = lambda x: np.sin(x)
Fx = fun(M.gridFx[:, 0])
Fy = fun(M.gridFy[:, 1])
Fz = fun(M.gridFz[:, 2])
F = np.concatenate((Fx,Fy,Fz))
divF = DIV*F
sol = lambda x, y, z: (np.cos(x)+np.cos(y)+np.cos(z))
divF_anal = sol(M.gridCC[:,0], M.gridCC[:,1], M.gridCC[:,2])
#err = np.linalg.norm((divF-divF_anal)*np.sqrt(vol), 2)
err = np.linalg.norm((divF-divF_anal), np.inf)
if icount == 1:
print 'h | inf norm | error ratio'
print '---------------------------------------'
print '%6.4f | %8.2e |'% (h1[0], err)
else:
print '%6.4f | %8.2e | %6.4f' % (h1[0], err, err_old/err)
err_old = err
SimPEG/tests/test_grad.py
-44
View File
@@ -1,44 +0,0 @@
import numpy as np
import sys
sys.path.append('../')
from TensorMesh import TensorMesh
err=0.
print '>> Test nodal Gradient operator'
for i in range(4):
icount=i+1
nc = 2**icount
# Define the mesh
h1 = np.ones(nc)/nc
h2 = np.ones(nc)/nc
h3 = np.ones(nc)/nc
h = [h1, h2, h3]
M = TensorMesh(h)
#n = M.plotGrid()
# Generate DIV matrix
GRAD = 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))
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], err)
else:
print '%6.4f | %8.2e | %6.4f' % (h1[0], err, err_old/err)
err_old = err
SimPEG/tests/test_operatorOrders.py
+88
View File
@@ -0,0 +1,88 @@
import numpy as np
from OrderTest import OrderTest
import unittest
class TestCurl(OrderTest):
name = "Curl"
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))
Ex = fun(self.M.gridEx[:, 1])
Ey = fun(self.M.gridEy[:, 2])
Ez = fun(self.M.gridEz[:, 0])
E = np.concatenate((Ex, Ey, Ez))
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))
# Generate DIV matrix
CURL = self.M.edgeCurl
curlE = CURL*E
err = np.linalg.norm((curlE-curlE_anal), np.inf)
return err
def test_order(self):
self.orderTest()
class TestFaceDiv(OrderTest):
name = "Face Divergence"
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
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])
err = np.linalg.norm((divF-divF_anal), np.inf)
return err
def test_order(self):
self.orderTest()
class TestNodalGrad(OrderTest):
name = "Nodal Gradient"
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))
phi = fun(self.M.gridN[:, 0], self.M.gridN[:, 1], self.M.gridN[:, 2])
gradE = GRAD*phi
Ex = sol(self.M.gridEx[:, 0])
Ey = sol(self.M.gridEy[:, 1])
Ez = sol(self.M.gridEz[:, 2])
gradE_anal = np.concatenate((Ex, Ey, Ez))
err = np.linalg.norm((gradE-gradE_anal), np.inf)
return err
def test_order(self):
self.orderTest()
if __name__ == '__main__':
unittest.main()