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)