From 113c97ca87e5ecac2cd61bc96f902fde35536085 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Tue, 15 Apr 2014 10:40:40 -0700 Subject: [PATCH] cylMesh Conventions.. A start --- simpegEM/TDEM/BaseTDEM.py | 23 ++++-- simpegEM/TDEM/SurveyTDEM.py | 2 +- simpegEM/Tests/test_TDEM_b_DerivAdjoint.py | 4 +- simpegEM/Tests/test_forward_EMproblem.py | 96 ++++++++++++++++++---- 4 files changed, 98 insertions(+), 27 deletions(-) diff --git a/simpegEM/TDEM/BaseTDEM.py b/simpegEM/TDEM/BaseTDEM.py index 7c0570f1..db1fbf92 100644 --- a/simpegEM/TDEM/BaseTDEM.py +++ b/simpegEM/TDEM/BaseTDEM.py @@ -22,13 +22,19 @@ class MixinInitialFieldCalc(object): return F def _getInitialFields_VMD_MVP(self): - if self.mesh._meshType is 'CYL1D': - MVP = Sources.MagneticDipoleVectorPotential(np.r_[0,0,self.survey.txLoc], np.c_[np.zeros(self.mesh.nN), self.mesh.gridN], 'x') + if self.mesh._meshType is 'CYL': + if self.mesh.isSymmetric: + MVP = Sources.MagneticDipoleVectorPotential(self.survey.txLoc, self.mesh.gridEy, 'y') + # MVP = Sources.MagneticDipoleVectorPotential(self.survey.txLoc, np.c_[np.zeros(self.mesh.nN), self.mesh.gridN], 'x') + else: + raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') elif self.mesh._meshType is 'TENSOR': MVPx = Sources.MagneticDipoleVectorPotential(self.survey.txLoc, self.mesh.gridEx, 'x') MVPy = Sources.MagneticDipoleVectorPotential(self.survey.txLoc, self.mesh.gridEy, 'y') MVPz = Sources.MagneticDipoleVectorPotential(self.survey.txLoc, self.mesh.gridEz, 'z') MVP = np.concatenate((MVPx, MVPy, MVPz)) + else: + raise Exception('Unknown mesh for VMD') # Initialize field object F = FieldsTDEM(self.mesh, 1, self.times.size, store=self.storeTheseFields) @@ -60,7 +66,7 @@ class MixinTimeStuff(object): nsteps = property(**nsteps()) def times(): - doc = "Modelling times" + doc = "Modeling times" def fget(self): t = np.r_[1:self.nsteps[0]+1]*self.dt[0] for i in range(1,self.dt.size): @@ -116,10 +122,10 @@ class ProblemBaseTDEM(MixinTimeStuff, MixinInitialFieldCalc, BaseProblem): def MeSigmaI(self): return self._MeSigmaI def makeMassMatrices(self, m): - m = self.model.transform(m) - self._MeSigma = self.mesh.getMass(m, loc='e') + sig = self.model.transform(m) + self._MeSigma = self.mesh.getEdgeInnerProduct(sig) self._MeSigmaI = sdiag(1/self.MeSigma.diagonal()) - self._MfMui = self.mesh.getMass(1/mu_0, loc='f') + self._MfMui = self.mesh.getFaceInnerProduct(1/mu_0) def calcFields(self, sol, solType, tInd): @@ -134,7 +140,8 @@ class ProblemBaseTDEM(MixinTimeStuff, MixinInitialFieldCalc, BaseProblem): return {'b':b, 'e':e} - solveOpts = {'factorize':True,'backend':'scipy'} + Solver = Solver + solveOpts = {} def fields(self, m): self.makeMassMatrices(m) @@ -153,7 +160,7 @@ class ProblemBaseTDEM(MixinTimeStuff, MixinInitialFieldCalc, BaseProblem): dtFact = dt A = self.getA(tInd) # print 'Factoring... (dt = ' + str(dt) + ')' - Asolve = Solver(A, options=self.solveOpts) + Asolve = self.Solver(A, **self.solveOpts) # print 'Done' rhs = RHS(tInd, F) sol = Asolve.solve(rhs) diff --git a/simpegEM/TDEM/SurveyTDEM.py b/simpegEM/TDEM/SurveyTDEM.py index dfaf8a91..4b1a91b3 100644 --- a/simpegEM/TDEM/SurveyTDEM.py +++ b/simpegEM/TDEM/SurveyTDEM.py @@ -47,7 +47,7 @@ class SurveyTDEM1D(BaseSurvey): def Qrx(self): if self._Qrx is None: if self.rxType == 'bz': - locType = 'fz' + locType = 'Fz' self._Qrx = self.prob.mesh.getInterpolationMat(self.rxLoc, locType=locType) return self._Qrx _Qrx = None diff --git a/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py b/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py index f65a46d3..0dcd7122 100644 --- a/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py +++ b/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py @@ -12,7 +12,7 @@ class TDEM_bDerivTests(unittest.TestCase): npad = 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) + mesh = Mesh.CylMesh([hx,hy], -hy.sum()/2) active = mesh.vectorCCz<0. model = Model.ActiveModel(mesh, active, -8, nC=mesh.nCz) @@ -21,7 +21,7 @@ class TDEM_bDerivTests(unittest.TestCase): opts = {'txLoc':0., - 'txType':'VMD_MVP', + 'txType': 'VMD_MVP', 'rxLoc':np.r_[150., 0.], 'rxType':'bz', 'timeCh':np.logspace(-4,-2,20), diff --git a/simpegEM/Tests/test_forward_EMproblem.py b/simpegEM/Tests/test_forward_EMproblem.py index 63dd864f..f3b39e10 100644 --- a/simpegEM/Tests/test_forward_EMproblem.py +++ b/simpegEM/Tests/test_forward_EMproblem.py @@ -5,54 +5,118 @@ 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 = 10. - ncx = 15 - ncy = 10 - npad = 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) + # 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. model = Model.ActiveModel(mesh, active, -8, nC=mesh.nCz) model = Model.ComboModel(mesh, [Model.LogModel, Model.Vertical1DModel, model]) - opts = {'txLoc':0., + opts = {'txLoc':np.array([0., 0., 0.]), 'txType':'VMD_MVP', - 'rxLoc':np.r_[30., 0.], + 'rxLoc':np.array([[10., 0., 0.]]), 'rxType':'bz', - 'timeCh':np.logspace(-4,-2.5, 21), + '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(model) - self.prb.setTimes([1e-6, 5e-6, 1e-5, 5e-5, 1e-4, 5e-4], [40, 40, 40, 40, 40, 40]) + + from mumpsSCI import MumpsSolver + self.prb.Solver = MumpsSolver + self.prb.setTimes([1e-6], [100]) self.sigma = np.ones(mesh.nCz)*1e-8 - self.sigma[mesh.vectorCCz<0] = 1e-3 + self.sigma[active] = sig_half self.sigma = np.log(self.sigma[active]) - self.showIt = False + 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. + + Logmodel = EM1D.BaseEM1D.BaseEM1DModel(mesh1D) + modelReal = Model.ComboModel(mesh1D, [Logmodel]) + 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(modelReal, **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 = mu_0*hzAnalyticDipoleT(self.dat.rxLoc[0], self.prb.times, np.exp(self.sigma[0])) + # 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], 'b', self.prb.times[bz_calc<0], -bz_calc[bz_calc<0], 'b--') + 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 diff + print 'Difference: ', diff self.assertTrue(diff < 0.10)