From 904921581124cf9e5949183138b980b6d770f788 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 3 Mar 2015 11:29:10 -0800 Subject: [PATCH] add apparent resistivity tests --- .../Tests/test_ApparentResistivityAnalytic.py | 48 ++++++++ simpegMT/Tests/test_FieldsObject.py | 103 ------------------ 2 files changed, 48 insertions(+), 103 deletions(-) create mode 100644 simpegMT/Tests/test_ApparentResistivityAnalytic.py delete mode 100644 simpegMT/Tests/test_FieldsObject.py diff --git a/simpegMT/Tests/test_ApparentResistivityAnalytic.py b/simpegMT/Tests/test_ApparentResistivityAnalytic.py new file mode 100644 index 00000000..395e8ae0 --- /dev/null +++ b/simpegMT/Tests/test_ApparentResistivityAnalytic.py @@ -0,0 +1,48 @@ +import unittest +from SimPEG import * +import simpegMT as MT + +TOL = 1e-6 + +def appResPhs(freq,z): + app_res = ((1./(8e-7*np.pi**2))/freq)*np.abs(z)**2 + app_phs = np.arctan2(-z.imag,z.real)*(180/np.pi) + return app_res, app_phs + +def appResNorm(sigmaHalf): + nFreq = 26 + + m1d = Mesh.TensorMesh([[(100,5,1.5),(100.,10),(100,5,1.5)]], x0=['C']) + sigma = np.zeros(m1d.nC) + sigmaHalf + sigma[m1d.gridCC[:]>200] = 1e-8 + + # Calculate the analytic fields + freqs = np.logspace(4,-4,nFreq) + Z = [] + for freq in freqs: + Ed, Eu, Hd, Hu = MT.Utils.getEHfields(m1d,sigma,freq,np.array([200])) + Z.append((Ed + Eu)/(Hd + Hu)) + + Zarr = np.concatenate(Z) + + app_r, app_p = appResPhs(freqs,Zarr) + + return np.linalg.norm(np.abs(app_r - np.ones(nFreq)/sigmaHalf)) / np.log10(sigmaHalf) + + +class TestAnalytics(unittest.TestCase): + + def setUp(self): + pass + def test_appRes2en1(self):self.assertLess(appResNorm(2e-1), TOL) + def test_appRes2en2(self):self.assertLess(appResNorm(2e-2), TOL) + def test_appRes2en3(self):self.assertLess(appResNorm(2e-3), TOL) + def test_appRes2en4(self):self.assertLess(appResNorm(2e-4), TOL) + def test_appRes2en5(self):self.assertLess(appResNorm(2e-5), TOL) + def test_appRes2en6(self):self.assertLess(appResNorm(2e-6), TOL) + + + + +if __name__ == '__main__': + unittest.main() diff --git a/simpegMT/Tests/test_FieldsObject.py b/simpegMT/Tests/test_FieldsObject.py deleted file mode 100644 index eac83bad..00000000 --- a/simpegMT/Tests/test_FieldsObject.py +++ /dev/null @@ -1,103 +0,0 @@ -import unittest -from SimPEG import * -import simpegMT as MT - -class FieldsTest(unittest.TestCase): - - def setUp(self): - mesh = Mesh.TensorMesh([np.ones(n)*5 for n in [10,11,12]],[0,0,-30]) - x = np.linspace(5,10,3) - XYZ = Utils.ndgrid(x,x,np.r_[0.]) - txLoc = np.r_[0,0,0.] - rxList0 = MT.FDEM.RxFDEM(XYZ, 'exi') - Tx0 = MT.FDEM.TxFDEM(txLoc, 'VMD', 3., [rxList0]) - rxList1 = MT.FDEM.RxFDEM(XYZ, 'bxi') - Tx1 = MT.FDEM.TxFDEM(txLoc, 'VMD', 3., [rxList1]) - rxList2 = MT.FDEM.RxFDEM(XYZ, 'bxi') - Tx2 = MT.FDEM.TxFDEM(txLoc, 'VMD', 2., [rxList2]) - rxList3 = MT.FDEM.RxFDEM(XYZ, 'bxi') - Tx3 = MT.FDEM.TxFDEM(txLoc, 'VMD', 2., [rxList3]) - Tx4 = MT.FDEM.TxFDEM(txLoc, 'VMD', 1., [rxList0, rxList1, rxList2, rxList3]) - txList = [Tx0,Tx1,Tx2,Tx3,Tx4] - survey = MT.FDEM.SurveyFDEM(txList) - self.F = MT.FDEM.FieldsFDEM(mesh, survey) - self.Tx0 = Tx0 - self.Tx1 = Tx1 - self.mesh = mesh - self.XYZ = XYZ - - def test_SetGet(self): - F = self.F - for freq in F.survey.freqs: - nFreq = F.survey.nTxByFreq[freq] - Txs = F.survey.getTransmitters(freq) - e = np.random.rand(F.mesh.nE, nFreq) - F[Txs, 'e'] = e - b = np.random.rand(F.mesh.nF, nFreq) - F[Txs, 'b'] = b - if nFreq == 1: - F[Txs, 'b'] = Utils.mkvc(b) - if e.shape[1] == 1: - e, b = Utils.mkvc(e), Utils.mkvc(b) - self.assertTrue(np.all(F[Txs, 'e'] == e)) - self.assertTrue(np.all(F[Txs, 'b'] == b)) - F[Txs] = {'b':b,'e':e} - self.assertTrue(np.all(F[Txs, 'e'] == e)) - self.assertTrue(np.all(F[Txs, 'b'] == b)) - - lastFreq = F[Txs] - self.assertTrue(type(lastFreq) is dict) - self.assertTrue(sorted([k for k in lastFreq]) == ['b','e']) - self.assertTrue(np.all(lastFreq['b'] == b)) - self.assertTrue(np.all(lastFreq['e'] == e)) - - Tx_f3 = F.survey.getTransmitters(3.) - self.assertTrue(F[Tx_f3,'b'].shape == (F.mesh.nF, 2)) - - b = np.random.rand(F.mesh.nF, 2) - Tx_f0 = F.survey.getTransmitters(self.Tx0.freq) - F[Tx_f0,'b'] = b - self.assertTrue(F[self.Tx0]['b'].shape == (F.mesh.nF,)) - self.assertTrue(F[self.Tx0,'b'].shape == (F.mesh.nF,)) - self.assertTrue(np.all(F[self.Tx0,'b'] == b[:,0])) - self.assertTrue(np.all(F[self.Tx1,'b'] == b[:,1])) - - def test_assertions(self): - freq = self.F.survey.freqs[0] - Txs = self.F.survey.getTransmitters(freq) - bWrongSize = np.random.rand(self.F.mesh.nE, self.F.survey.nTxByFreq[freq]) - def fun(): self.F[Txs, 'b'] = bWrongSize - self.assertRaises(ValueError, fun) - def fun(): self.F[-999.] - self.assertRaises(KeyError, fun) - def fun(): self.F['notRight'] - self.assertRaises(KeyError, fun) - def fun(): self.F[Txs,'notThere'] - self.assertRaises(KeyError, fun) - - def test_FieldProjections(self): - F = self.F - for freq in F.survey.freqs: - nFreq = F.survey.nTxByFreq[freq] - Txs = F.survey.getTransmitters(freq) - e = np.random.rand(F.mesh.nE, nFreq) - b = np.random.rand(F.mesh.nF, nFreq) - F[Txs] = {'b':b,'e':e} - - Txs = F.survey.getTransmitters(freq) - for ii, tx in enumerate(Txs): - for jj, rx in enumerate(tx.rxList): - dat = rx.projectFields(tx, self.mesh, F) - self.assertTrue(dat.dtype == float) - fieldType = rx.projField - u = {'b':b[:,ii], 'e': e[:,ii]}[fieldType] - real_or_imag = rx.projComp - u = getattr(u, real_or_imag) - gloc = rx.projGLoc - d = self.mesh.getInterpolationMat(self.XYZ, gloc)*u - self.assertTrue(np.all(dat == d)) - - - -if __name__ == '__main__': - unittest.main()