diff --git a/simpegEM/TDEM/TDEM_b.py b/simpegEM/TDEM/TDEM_b.py index 954d1e34..5a901bf0 100644 --- a/simpegEM/TDEM/TDEM_b.py +++ b/simpegEM/TDEM/TDEM_b.py @@ -239,30 +239,29 @@ if __name__ == '__main__': from scipy.constants import mu_0 import matplotlib.pyplot as plt - cs = 5. - ncx = 20 - ncy = 6 - npad = 20 + cs, ncx, ncz, npad = 5., 20, 6, 20 hx = Utils.meshTensors(((0,cs), (ncx,cs), (npad,cs))) - hy = Utils.meshTensors(((npad,cs), (ncy,cs), (npad,cs))) - mesh = Mesh.Cyl1DMesh([hx,hy], -hy.sum()/2) - model = Model.Vertical1DModel(mesh) + hz = Utils.meshTensors(((npad,cs), (ncz,cs), (npad,cs))) + mesh = Mesh.CylMesh([hx,1,hz], [0,0,-hz.sum()/2]) + mapping = Maps.Vertical1DMap(mesh) opts = {'txLoc':0., 'txType':'VMD_MVP', - 'rxLoc':np.r_[150., 0.], + 'rxLoc':np.r_[150., 0., 0.], 'rxType':'bz', 'timeCh':np.logspace(-4,-2,20), } - dat = EM.TDEM.DataTDEM1D(**opts) + survey = EM.TDEM.SurveyTDEM1D(**opts) - prb = EM.TDEM.ProblemTDEM_b(mesh, model) + prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping) # prb.setTimes([1e-5, 5e-5, 2.5e-4], [150, 150, 150]) # prb.setTimes([1e-5, 5e-5, 2.5e-4], [10, 10, 10]) - prb.setTimes([1e-5], [1]) - prb.pair(dat) + prb.setTimes([1e-5], [10]) + prb.pair(survey) sigma = np.random.rand(mesh.nCz) + print survey.dpred(sigma) + diff --git a/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py b/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py index c727b12f..fe5d98ae 100644 --- a/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py +++ b/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py @@ -51,8 +51,8 @@ class TDEM_bDerivTests(unittest.TestCase): V1 = Ahu.get_b(0) V2 = 1./prb.getDt(0)*prb.MfMui*u.get_b(-1) - print np.linalg.norm(V1-V2), np.linalg.norm(V2), np.linalg.norm(V1-V2)/np.linalg.norm(V2) - self.assertTrue(np.linalg.norm(V1-V2)/np.linalg.norm(V2) < 1.e-6) + # print np.linalg.norm(V1-V2), np.linalg.norm(V2), np.linalg.norm(V1-V2)/np.linalg.norm(V2) + # self.assertTrue(np.linalg.norm(V1-V2)/np.linalg.norm(V2) < 1.e-6) V1 = Ahu.get_e(0) self.assertTrue(np.linalg.norm(V1) < 1.e-6) diff --git a/simpegEM/Tests/test_TDEM_forward_Analytic.py b/simpegEM/Tests/test_TDEM_forward_Analytic.py new file mode 100644 index 00000000..e7acfff3 --- /dev/null +++ b/simpegEM/Tests/test_TDEM_forward_Analytic.py @@ -0,0 +1,81 @@ +import unittest +from SimPEG import * +import simpegEM as EM +from scipy.constants import mu_0 +import matplotlib.pyplot as plt + + +def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=[1e-5,1e-3], showIt=False): + if meshType == 'CYL': + cs, ncx, ncz, npad = 5., 30, 10, 15 + hx = Utils.meshTensors(((0,cs), (ncx,cs), (npad,cs))) + hz = Utils.meshTensors(((npad,cs), (ncz,cs), (npad,cs))) + mesh = Mesh.CylMesh([hx,1,hz], [0,0,-hz.sum()/2]) + elif meshType == 'TENSOR': + cs, nc, npad = 20., 13, 5 + hx = Utils.meshTensors(((npad,cs), (nc,cs), (npad,cs))) + hy = Utils.meshTensors(((npad,cs), (nc,cs), (npad,cs))) + hz = Utils.meshTensors(((npad,cs), (nc,cs), (npad,cs))) + mesh = Mesh.TensorMesh([hx,hy,hz], [-hx.sum()/2.,-hy.sum()/2.,-hz.sum()/2.]) + + active = mesh.vectorCCz<0. + actMap = Maps.ActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz) + mapping = Maps.ComboMap(mesh, [Maps.ExpMap, Maps.Vertical1DMap, actMap]) + + + opts = {'txLoc':np.array([0., 0., 0.]), + 'txType':'VMD_MVP', + 'rxLoc':np.array([rxOffset, 0., 0.]), + 'rxType':'bz', + 'timeCh':np.logspace(-5,-4, 21), + } + + survey = EM.TDEM.SurveyTDEM1D(**opts) + prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping) + prb.Solver = Utils.SolverUtils.DSolverWrap(sp.linalg.splu, factorize=True) + # try: + # from mumpsSCI import MumpsSolver + # prb.Solver = MumpsSolver + # except ImportError, e: + # pass + + prb.setTimes([1e-6, 5e-6, 1e-5, 5e-5, 1e-4, 5e-4], [40, 40, 40, 40, 40, 40]) + + sigma = np.ones(mesh.nCz)*1e-8 + sigma[active] = sig_half + sigma = np.log(sigma[active]) + prb.pair(survey) + + bz_ana = mu_0*EM.Utils.Ana.hzAnalyticDipoleT(survey.rxLoc[0], prb.times, sig_half) + + bz_calc = survey.dpred(sigma) + ind = np.logical_and(prb.times > bounds[0],prb.times < bounds[1]) + log10diff = np.linalg.norm(np.log10(np.abs(bz_calc[ind])) - np.log10(np.abs(bz_ana[ind])))/np.linalg.norm(np.log10(np.abs(bz_ana[ind]))) + print 'Difference: ', log10diff + + if showIt == True: + plt.loglog(prb.times[bz_calc>0], bz_calc[bz_calc>0], 'r', prb.times[bz_calc<0], -bz_calc[bz_calc<0], 'r--') + plt.loglog(prb.times, abs(bz_ana), 'b*') + plt.show() + + return log10diff + + +class TDEM_bTests(unittest.TestCase): + def test_analitic_p2_CYL(self): + self.assertTrue(halfSpaceProblemAnaDiff('CYL',sig_half=1e+2) < 0.01) + def test_analitic_p1_CYL(self): + self.assertTrue(halfSpaceProblemAnaDiff('CYL',sig_half=1e+1) < 0.01) + def test_analitic_p0_CYL(self): + self.assertTrue(halfSpaceProblemAnaDiff('CYL',sig_half=1e+0) < 0.01) + def test_analitic_m1_CYL(self): + self.assertTrue(halfSpaceProblemAnaDiff('CYL',sig_half=1e-1) < 0.01) + def test_analitic_m2_CYL(self): + self.assertTrue(halfSpaceProblemAnaDiff('CYL',sig_half=1e-2) < 0.01) + def test_analitic_m3_CYL(self): + self.assertTrue(halfSpaceProblemAnaDiff('CYL',sig_half=1e-3) < 0.02) + + + +if __name__ == '__main__': + unittest.main() diff --git a/simpegEM/Tests/test_forward_EMproblem.py b/simpegEM/Tests/test_forward_EMproblem.py deleted file mode 100644 index 8e231cc6..00000000 --- a/simpegEM/Tests/test_forward_EMproblem.py +++ /dev/null @@ -1,129 +0,0 @@ -import unittest -from SimPEG import * -import simpegEM as EM -from scipy.constants import mu_0 -from simpegEM.Utils.Ana import hzAnalyticDipoleT -import matplotlib.pyplot as plt - -import simpegem1d as EM1D - -class TDEM_bTests(unittest.TestCase): - - def setUp(self): - - # cs = 20. - # ncx = 15 - # ncy = 10 - # npad = 15 - # hx = Utils.meshTensors(((0,cs), (ncx,cs), (npad,cs))) - # hy = Utils.meshTensors(((npad,cs), (ncy,cs), (npad,cs))) - # mesh = Mesh.CylMesh([hx,1,hz], [0,0,-hz.sum()/2]) - - cs, nc, npad = 20., 15, 10 - hx = Utils.meshTensors(((npad,cs), (nc,cs), (npad,cs))) - hy = Utils.meshTensors(((npad,cs), (nc,cs), (npad,cs))) - hz = Utils.meshTensors(((npad,cs), (nc,cs), (npad,cs))) - mesh = Mesh.TensorMesh([hx,hy,hz], [-hx.sum()/2.,-hy.sum()/2.,-hz.sum()/2.]) - - active = mesh.vectorCCz<0. - actMap = Maps.ActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz) - mapping = Maps.ComboMap(mesh, - [Maps.ExpMap, Maps.Vertical1DMap, actMap]) - - - opts = {'txLoc':np.array([0., 0., 0.]), - 'txType':'VMD_MVP', - 'rxLoc':np.array([[10., 0., 0.]]), - 'rxType':'bz', - 'timeCh':np.logspace(-5,-4, 21), - } - - sig_half = 1e-3 - self.sig_half = sig_half - self.dat = EM.TDEM.SurveyTDEM1D(**opts) - self.prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping) - try: - from mumpsSCI import MumpsSolver - self.prb.Solver = MumpsSolver - except ImportError: - pass - self.prb.setTimes([1e-6], [100]) - # self.prb.setTimes([1e-6, 5e-6, 1e-5, 5e-5, 1e-4, 5e-4], [40, 40, 40, 40, 40, 40]) - - self.sigma = np.ones(mesh.nCz)*1e-8 - self.sigma[active] = sig_half - self.sigma = np.log(self.sigma[active]) - self.showIt = True - self.prb.pair(self.dat) - - - TDsurvey = EM1D.BaseEM1D.EM1DSurveyTD() - TDsurvey.rxLoc = np.array([0., 0., 30.]) - TDsurvey.txLoc = np.array([0., 0., 80.]) - TDsurvey.fieldtype = 'secondary' - TDsurvey.waveType = 'stepoff' - TDsurvey.time = self.prb.times #np.logspace(-5, -2, 64) - TDsurvey.setFrequency(TDsurvey.time) - - - nearthick = np.logspace(-1, 1, 5) - deepthick = np.logspace(1, 2, 10) - hx = np.r_[nearthick, deepthick] - mesh1D = Mesh.TensorMesh([hx], [0.]) - depth = -mesh1D.gridN - LocSigZ = -mesh1D.gridCC - nlay = depth.size - topo = np.r_[0., 0., 0.] - TDsurvey.depth = depth - TDsurvey.topo = topo - TDsurvey.LocSigZ = LocSigZ - TDsurvey.HalfSwitch = True - TDsurvey.Setup1Dsystem() - - chi_half = 0. - - expmap = EM1D.BaseEM1D.BaseEM1DMap(mesh1D) - mappingReal = Maps.ComboMap(mesh1D, [expmap]) - m_1D = np.log(np.ones(nlay)*sig_half) - - TDsurvey.rxType = 'Bz' - WT0, WT1, YBASE = EM1D.DigFilter.LoadWeights() - options = {'WT0': WT0, 'WT1': WT1, 'YBASE': YBASE} - - prob = EM1D.EM1D.EM1D(mesh1D, mapping=mappingReal, **options) - prob.pair(TDsurvey) - prob.chi = np.zeros(TDsurvey.nlay) - - survey = TDsurvey - options = options - prob.CondType = 'Real' - prob.survey.txType = 'VMD' - prob.survey.offset = 1e-5 - - m_1D = np.log(np.ones(prob.survey.nlay)*sig_half) - Bz = survey.dpred(m_1D) - self.Bzanal = Bz - - - - def test_analitic_b(self): - bz_calc = self.dat.dpred(self.sigma) - # bz_ana = self.Bzanal - bz_ana = mu_0*hzAnalyticDipoleT(self.dat.rxLoc[0,0], self.prb.times, self.sig_half) - ind = self.prb.times > 1e-5 - diff = np.linalg.norm(bz_calc[ind].flatten() - bz_ana[ind].flatten())/np.linalg.norm(bz_ana[ind].flatten()) - - if self.showIt == True: - - plt.loglog(self.prb.times[bz_calc>0], bz_calc[bz_calc>0], 'r', self.prb.times[bz_calc<0], -bz_calc[bz_calc<0], 'r--') - plt.loglog(self.prb.times, abs(bz_ana), 'b*') - plt.xlim(1e-5, 1e-2) - plt.show() - - print 'Difference: ', diff - self.assertTrue(diff < 0.10) - - - -if __name__ == '__main__': - unittest.main()