From 113c97ca87e5ecac2cd61bc96f902fde35536085 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Tue, 15 Apr 2014 10:40:40 -0700 Subject: [PATCH 1/3] 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) From eab87c892b2f919cb17056881e666764a5552db7 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Tue, 15 Apr 2014 12:07:58 -0700 Subject: [PATCH 2/3] Updates to for mapping changes in SimPEG --- simpegEM/FDEM/FDEM.py | 6 ++--- simpegEM/TDEM/BaseTDEM.py | 6 ++--- simpegEM/TDEM/TDEM_b.py | 23 +++++++++------- simpegEM/Tests/test_FDEM.py | 14 +++++++--- simpegEM/Tests/test_FDEM_analytics.py | 4 +-- simpegEM/Tests/test_TDEM_b_DerivAdjoint.py | 31 +++++++++++----------- simpegEM/Tests/test_forward_EMproblem.py | 24 ++++++++++------- 7 files changed, 61 insertions(+), 47 deletions(-) diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index 521e74c0..ea1cb0bc 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -66,13 +66,13 @@ class BaseProblemFDEM(Problem.BaseProblem): @property def curTModel(self): if getattr(self, '_curTModel', None) is None: - self._curTModel = self.model.transform(self.curModel) + self._curTModel = self.mapping.transform(self.curModel) return self._curTModel @property def curTModelDeriv(self): if getattr(self, '_curTModelDeriv', None) is None: - self._curTModelDeriv = self.model.transformDeriv(self.curModel) + self._curTModelDeriv = self.mapping.transformDeriv(self.curModel) return self._curTModelDeriv def fields(self, m): @@ -132,7 +132,7 @@ class BaseProblemFDEM(Problem.BaseProblem): if not isinstance(v, self.dataPair): v = self.dataPair(self.survey, v) - Jtv = np.zeros(self.model.nP) + Jtv = np.zeros(self.mapping.nP) for freq in self.survey.freqs: AT = self.getA(freq).T diff --git a/simpegEM/TDEM/BaseTDEM.py b/simpegEM/TDEM/BaseTDEM.py index db1fbf92..05a084a1 100644 --- a/simpegEM/TDEM/BaseTDEM.py +++ b/simpegEM/TDEM/BaseTDEM.py @@ -92,8 +92,8 @@ class MixinTimeStuff(object): class ProblemBaseTDEM(MixinTimeStuff, MixinInitialFieldCalc, BaseProblem): """docstring for ProblemTDEM1D""" - def __init__(self, model, **kwargs): - BaseProblem.__init__(self, model, **kwargs) + def __init__(self, mesh, mapping=None, **kwargs): + BaseProblem.__init__(self, mesh, mapping=mapping, **kwargs) #################################################### @@ -122,7 +122,7 @@ class ProblemBaseTDEM(MixinTimeStuff, MixinInitialFieldCalc, BaseProblem): def MeSigmaI(self): return self._MeSigmaI def makeMassMatrices(self, m): - sig = self.model.transform(m) + sig = self.mapping.transform(m) self._MeSigma = self.mesh.getEdgeInnerProduct(sig) self._MeSigmaI = sdiag(1/self.MeSigma.diagonal()) self._MfMui = self.mesh.getFaceInnerProduct(1/mu_0) diff --git a/simpegEM/TDEM/TDEM_b.py b/simpegEM/TDEM/TDEM_b.py index 7ba039d6..954d1e34 100644 --- a/simpegEM/TDEM/TDEM_b.py +++ b/simpegEM/TDEM/TDEM_b.py @@ -16,8 +16,8 @@ class ProblemTDEM_b(ProblemBaseTDEM): with \\\(\\b\\\) defined on cell faces and \\\(\e\\\) defined on edges. """ - def __init__(self, model, **kwargs): - ProblemBaseTDEM.__init__(self, model, **kwargs) + def __init__(self, mesh, mapping=None, **kwargs): + ProblemBaseTDEM.__init__(self, mesh, mapping=mapping, **kwargs) solType = 'b' @@ -35,11 +35,11 @@ class ProblemTDEM_b(ProblemBaseTDEM): """ dt = self.getDt(tInd) - return self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui + (1/dt)*self.MfMui + return self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui + (1.0/dt)*self.MfMui def getRHS(self, tInd, F): dt = self.getDt(tInd) - return (1/dt)*self.MfMui*F.get_b(tInd-1) + return (1.0/dt)*self.MfMui*F.get_b(tInd-1) #################################################### @@ -74,7 +74,8 @@ class ProblemTDEM_b(ProblemBaseTDEM): if u is None: u = self.fields(m) p = FieldsTDEM(self.mesh, 1, self.times.size, 'b') - c = self.mesh.getEdgeMassDeriv()*self.model.transformDeriv(m)*vec + curModel = self.mapping.transform(m) + c = self.mesh.getEdgeInnerProductDeriv(curModel)*self.mapping.transformDeriv(m)*vec for i in range(self.times.size): ei = u.get_e(i) pVal = np.empty_like(ei) @@ -91,7 +92,9 @@ class ProblemTDEM_b(ProblemBaseTDEM): tmp = np.zeros((self.mesh.nE,self.survey.nTx)) for i in range(self.nTimes): tmp += v.get_e(i)*u.get_e(i) - p = -mkvc(self.model.transformDeriv(m).T*self.mesh.getEdgeMassDeriv().T*tmp) + + curModel = self.mapping.transform(m) + p = -mkvc(self.mapping.transformDeriv(m).T*self.mesh.getEdgeInnerProductDeriv(curModel).T*tmp) return p def solveAh(self, m, p): @@ -100,7 +103,7 @@ class ProblemTDEM_b(ProblemBaseTDEM): if tInd == 0: return rhs dt = self.getDt(tInd) - return rhs + 1./dt*self.MfMui*u.get_b(tInd-1) + return rhs + 1.0/dt*self.MfMui*u.get_b(tInd-1) def AhCalcFields(sol, solType, tInd): b = sol @@ -117,7 +120,7 @@ class ProblemTDEM_b(ProblemBaseTDEM): if tInd == self.nTimes-1: return rhs dt = self.getDt(tInd+1) - return rhs + 1./dt*self.MfMui*u.get_b(tInd+1) + return rhs + 1.0/dt*self.MfMui*u.get_b(tInd+1) def AhtCalcFields(sol, solType, tInd): b = sol @@ -167,14 +170,14 @@ class ProblemTDEM_b(ProblemBaseTDEM): self.makeMassMatrices(sigma) dt = self.getDt(0) - b = 1/dt*self.MfMui*vec.get_b(0) + self.MfMui*self.mesh.edgeCurl*vec.get_e(0) + b = 1.0/dt*self.MfMui*vec.get_b(0) + self.MfMui*self.mesh.edgeCurl*vec.get_e(0) e = self.mesh.edgeCurl.T*self.MfMui*vec.get_b(0) - self.MeSigma*vec.get_e(0) f = FieldsTDEM(self.mesh, 1, self.times.size, 'b') f.set_b(b, 0) f.set_e(e, 0) for i in range(1,self.nTimes): dt = self.getDt(i) - b = 1/dt*self.MfMui*vec.get_b(i) + self.MfMui*self.mesh.edgeCurl*vec.get_e(i) - 1/dt*self.MfMui*vec.get_b(i-1) + b = 1.0/dt*self.MfMui*vec.get_b(i) + self.MfMui*self.mesh.edgeCurl*vec.get_e(i) - 1.0/dt*self.MfMui*vec.get_b(i-1) e = self.mesh.edgeCurl.T*self.MfMui*vec.get_b(i) - self.MeSigma*vec.get_e(i) f.set_b(b, i) f.set_e(e, i) diff --git a/simpegEM/Tests/test_FDEM.py b/simpegEM/Tests/test_FDEM.py index aa548845..3f53684a 100644 --- a/simpegEM/Tests/test_FDEM.py +++ b/simpegEM/Tests/test_FDEM.py @@ -15,7 +15,7 @@ def getProblem(fdemType, comp): hz = Utils.meshTensors(((npad,cs), (ncz,cs), (npad,cs))) mesh = Mesh.TensorMesh([hx,hy,hz],[-hx.sum()/2., -hy.sum()/2., -hz.sum()/2.]) - model = Model.LogModel(mesh) + mapping = Maps.ExpMap(mesh) x = np.linspace(-30,30,6) XYZ = Utils.ndgrid(x,x,np.r_[0]) @@ -25,13 +25,19 @@ def getProblem(fdemType, comp): survey = EM.FDEM.SurveyFDEM([Tx0]) if fdemType == 'e': - prb = EM.FDEM.ProblemFDEM_e(model) + prb = EM.FDEM.ProblemFDEM_e(mesh, mapping=mapping) elif fdemType == 'b': - prb = EM.FDEM.ProblemFDEM_b(model) + prb = EM.FDEM.ProblemFDEM_b(mesh, mapping=mapping) else: raise NotImplementedError() prb.pair(survey) + try: + from mumpsSCI import MumpsSolver + prb.Solver = MumpsSolver + except ImportError, e: + pass + return prb def adjointTest(fdemType, comp): @@ -42,7 +48,7 @@ def adjointTest(fdemType, comp): survey = prb.survey v = np.random.rand(survey.nD) - w = np.random.rand(prb.model.nP) + w = np.random.rand(prb.mapping.nP) u = prb.fields(m) vJw = v.dot(prb.Jvec(m, w, u=u)) diff --git a/simpegEM/Tests/test_FDEM_analytics.py b/simpegEM/Tests/test_FDEM_analytics.py index 29db5c78..d0fd3482 100644 --- a/simpegEM/Tests/test_FDEM_analytics.py +++ b/simpegEM/Tests/test_FDEM_analytics.py @@ -17,7 +17,7 @@ class FDEM_analyticTests(unittest.TestCase): hz = Utils.meshTensors(((npad,cs), (ncz,cs), (npad,cs))) mesh = Mesh.TensorMesh([hx,hy,hz], x0=[-hx.sum()/2.,-hy.sum()/2.,-hz.sum()/2.,]) - model = Model.LogModel(mesh) + mapping = Maps.ExpMap(mesh) x = np.linspace(-10,10,5) XYZ = Utils.ndgrid(x,np.r_[0],np.r_[0]) @@ -26,7 +26,7 @@ class FDEM_analyticTests(unittest.TestCase): survey = EM.FDEM.SurveyFDEM([Tx0]) - prb = EM.FDEM.ProblemFDEM_b(model) + prb = EM.FDEM.ProblemFDEM_b(mesh, mapping=mapping) prb.pair(survey) prb.Solver = Utils.SolverUtils.DSolverWrap(sp.linalg.splu, checkAccuracy=False) diff --git a/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py b/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py index 0dcd7122..c727b12f 100644 --- a/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py +++ b/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py @@ -12,23 +12,23 @@ 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.CylMesh([hx,hy], -hy.sum()/2) + mesh = Mesh.CylMesh([hx,1,hy], [0,0,-hy.sum()/2]) active = mesh.vectorCCz<0. - model = Model.ActiveModel(mesh, active, -8, nC=mesh.nCz) - model = Model.ComboModel(mesh, - [Model.LogModel, Model.Vertical1DModel, model]) + activeMap = Maps.ActiveCells(mesh, active, -8, nC=mesh.nCz) + mapping = Maps.ComboMap(mesh, + [Maps.ExpMap, Maps.Vertical1DMap, activeMap]) opts = {'txLoc':0., 'txType': 'VMD_MVP', - 'rxLoc':np.r_[150., 0.], + 'rxLoc':np.r_[40., 0., 0.], 'rxType':'bz', 'timeCh':np.logspace(-4,-2,20), } self.dat = EM.TDEM.SurveyTDEM1D(**opts) - self.prb = EM.TDEM.ProblemTDEM_b(model) + self.prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping) self.prb.setTimes([1e-5, 5e-5, 2.5e-4], [10, 10, 10]) self.sigma = np.ones(mesh.nCz)*1e-8 @@ -50,7 +50,8 @@ class TDEM_bDerivTests(unittest.TestCase): Ahu = prb.AhVec(sigma, u) V1 = Ahu.get_b(0) - V2 = 1/prb.getDt(0)*prb.MfMui*u.get_b(-1) + 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) V1 = Ahu.get_e(0) @@ -137,10 +138,10 @@ class TDEM_bDerivTests(unittest.TestCase): """ # Random model and perturbation - sigma = np.random.rand(self.prb.model.nP) + sigma = np.random.rand(self.prb.mapping.nP) f = self.prb.fields(sigma) - dm = 1000*np.random.rand(self.prb.model.nP) + dm = 1000*np.random.rand(self.prb.mapping.nP) h = 0.01 derChk = lambda m: [self.prb.AhVec(m, f).fieldVec(), lambda mx: self.prb.Gvec(sigma, mx, u=f).fieldVec()] @@ -155,7 +156,7 @@ class TDEM_bDerivTests(unittest.TestCase): mesh = self.mesh sigma = self.sigma - dm = 10*np.random.rand(prb.model.nP) + dm = 10*np.random.rand(prb.mapping.nP) f = prb.fields(sigma) derChk = lambda m: [self.prb.fields(m).fieldVec(), lambda mx: -prb.solveAh(sigma, prb.Gvec(sigma, mx, u=f)).fieldVec()] @@ -172,7 +173,7 @@ class TDEM_bDerivTests(unittest.TestCase): sigma = self.sigma # d_sig = 0.8*sigma #np.random.rand(mesh.nCz) - d_sig = 10*np.random.rand(prb.model.nP) + d_sig = 10*np.random.rand(prb.mapping.nP) derChk = lambda m: [prb.survey.dpred(m), lambda mx: -prb.Jvec(sigma, mx)] @@ -221,7 +222,7 @@ class TDEM_bDerivTests(unittest.TestCase): def test_solveAhtVsAhtVec(self): prb = self.prb mesh = self.mesh - sigma = np.random.rand(prb.model.nP) + sigma = np.random.rand(prb.mapping.nP) f1 = EM.TDEM.FieldsTDEM(mesh, 1, prb.nTimes, 'b') for i in range(f1.nTimes): @@ -258,8 +259,8 @@ class TDEM_bDerivTests(unittest.TestCase): mesh = self.mesh prb = self.prb - m = np.random.rand(prb.model.nP) - sigma = np.random.rand(prb.model.nP) + m = np.random.rand(prb.mapping.nP) + sigma = np.random.rand(prb.mapping.nP) u = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nTimes, 'b') for i in range(u.nTimes): @@ -280,7 +281,7 @@ class TDEM_bDerivTests(unittest.TestCase): prb = self.prb sigma = self.sigma - m = np.random.rand(prb.model.nP) + m = np.random.rand(prb.mapping.nP) d = np.random.rand(prb.nTimes) V1 = d.dot(prb.Jvec(sigma, m)) diff --git a/simpegEM/Tests/test_forward_EMproblem.py b/simpegEM/Tests/test_forward_EMproblem.py index f3b39e10..fa1aecf4 100644 --- a/simpegEM/Tests/test_forward_EMproblem.py +++ b/simpegEM/Tests/test_forward_EMproblem.py @@ -26,9 +26,10 @@ class TDEM_bTests(unittest.TestCase): 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]) + actMap = Maps.ActiveCells(mesh, active, -8, nC=mesh.nCz) + mapping = Maps.ComboMap(mesh, + [Maps.ExpMap, Maps.Vertical1DMap, actMap]) + opts = {'txLoc':np.array([0., 0., 0.]), 'txType':'VMD_MVP', @@ -40,11 +41,14 @@ class TDEM_bTests(unittest.TestCase): sig_half = 1e-3 self.sig_half = sig_half self.dat = EM.TDEM.SurveyTDEM1D(**opts) - self.prb = EM.TDEM.ProblemTDEM_b(model) - - from mumpsSCI import MumpsSolver - self.prb.Solver = MumpsSolver + 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 @@ -78,15 +82,15 @@ class TDEM_bTests(unittest.TestCase): chi_half = 0. - Logmodel = EM1D.BaseEM1D.BaseEM1DModel(mesh1D) - modelReal = Model.ComboModel(mesh1D, [Logmodel]) + 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(modelReal, **options) + prob = EM1D.EM1D.EM1D(mesh1D, mapping=mappingReal, **options) prob.pair(TDsurvey) prob.chi = np.zeros(TDsurvey.nlay) From 2d90799415e71da60a810bb0228754dba9c29650 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Tue, 15 Apr 2014 15:39:59 -0700 Subject: [PATCH 3/3] travis updates. --- README.md | 2 +- docs/index.rst | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 0b3aa035..1f8dc810 100644 --- a/README.md +++ b/README.md @@ -15,7 +15,7 @@ Tests: [https://travis-ci.org/simpeg/simpegem](https://travis-ci.org/simpeg/simpegem) Build Status: -[![Build Status](https://travis-ci.org/simpeg/simpegem.png)](https://travis-ci.org/simpeg/simpegem) +[![Build Status](https://travis-ci.org/simpeg/simpegem.svg?branch=master)](https://travis-ci.org/simpeg/simpegem) Bugs & Issues: [https://github.com/simpeg/simpegem/issues](https://github.com/simpeg/simpegem/issues) diff --git a/docs/index.rst b/docs/index.rst index 2c29fdfe..9666e5e8 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -52,7 +52,7 @@ Testing simpegEM ================ * Master Branch - .. image:: https://travis-ci.org/simpeg/simpegem.png?branch=master + .. image:: https://travis-ci.org/simpeg/simpegem.svg?branch=master :target: https://travis-ci.org/simpeg/simpegem :alt: Master Branch :align: center