From d0efdb509c094f093b8a5059c71f5ded80a66434 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 22 Jul 2013 12:26:17 -0700 Subject: [PATCH] Created an OrderTest class 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. --- SimPEG/tests/OrderTest.py | 50 ++++++++++++++++ SimPEG/tests/test_curl.py | 46 --------------- SimPEG/tests/test_div.py | 45 --------------- SimPEG/tests/test_grad.py | 44 --------------- SimPEG/tests/test_operatorOrders.py | 88 +++++++++++++++++++++++++++++ 5 files changed, 138 insertions(+), 135 deletions(-) create mode 100644 SimPEG/tests/OrderTest.py delete mode 100644 SimPEG/tests/test_curl.py delete mode 100644 SimPEG/tests/test_div.py delete mode 100644 SimPEG/tests/test_grad.py create mode 100644 SimPEG/tests/test_operatorOrders.py diff --git a/SimPEG/tests/OrderTest.py b/SimPEG/tests/OrderTest.py new file mode 100644 index 00000000..68837f7a --- /dev/null +++ b/SimPEG/tests/OrderTest.py @@ -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() diff --git a/SimPEG/tests/test_curl.py b/SimPEG/tests/test_curl.py deleted file mode 100644 index e90f8cfa..00000000 --- a/SimPEG/tests/test_curl.py +++ /dev/null @@ -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 \ No newline at end of file diff --git a/SimPEG/tests/test_div.py b/SimPEG/tests/test_div.py deleted file mode 100644 index 53ca506b..00000000 --- a/SimPEG/tests/test_div.py +++ /dev/null @@ -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 diff --git a/SimPEG/tests/test_grad.py b/SimPEG/tests/test_grad.py deleted file mode 100644 index fbc36aeb..00000000 --- a/SimPEG/tests/test_grad.py +++ /dev/null @@ -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 - diff --git a/SimPEG/tests/test_operatorOrders.py b/SimPEG/tests/test_operatorOrders.py new file mode 100644 index 00000000..9bb9fa03 --- /dev/null +++ b/SimPEG/tests/test_operatorOrders.py @@ -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()