From 7a15797fbeefbe5c51a91f605787e272cb334aae Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sun, 18 May 2014 16:39:10 -0700 Subject: [PATCH] Mesh and model updates. --- simpegEM/Base.py | 20 ++++--------------- simpegEM/FDEM/FDEM.py | 12 +++++------ simpegEM/TDEM/BaseTDEM.py | 18 +++++++++++++---- simpegEM/TDEM/TDEM_b.py | 4 ++-- simpegEM/Tests/test_FDEM_analytics.py | 7 ++++++- simpegEM/Tests/test_TDEM_b_DerivAdjoint.py | 9 +++++++-- .../Tests/test_TDEM_b_MultiTx_DerivAdjoint.py | 3 +-- simpegEM/Tests/test_TDEM_combos.py | 2 +- 8 files changed, 41 insertions(+), 34 deletions(-) diff --git a/simpegEM/Base.py b/simpegEM/Base.py index 2b6fd03f..373cee75 100644 --- a/simpegEM/Base.py +++ b/simpegEM/Base.py @@ -1,4 +1,4 @@ -from SimPEG import Survey, Problem, Utils, np, sp, Solver as SimpegSolver +from SimPEG import Survey, Problem, Utils, Models, np, sp, Solver as SimpegSolver from scipy.constants import mu_0 class BaseEMProblem(Problem.BaseProblem): @@ -38,7 +38,7 @@ class BaseEMProblem(Problem.BaseProblem): def MeSigma(self): #TODO: hardcoded to sigma as the model if getattr(self, '_MeSigma', None) is None: - sigma = self.curTModel + sigma = self.curModel.transform self._MeSigma = self.mesh.getEdgeInnerProduct(sigma) return self._MeSigma @@ -46,23 +46,11 @@ class BaseEMProblem(Problem.BaseProblem): def MeSigmaI(self): #TODO: hardcoded to sigma as the model if getattr(self, '_MeSigmaI', None) is None: - sigma = self.curTModel + sigma = self.curModel.transform self._MeSigmaI = self.mesh.getEdgeInnerProduct(sigma, invMat=True) return self._MeSigmaI - curModel = Utils.dependentProperty('_curModel', None, ['_MeSigma', '_MeSigmaI', '_curTModel', '_curTModelDeriv'], 'Sets the current model, and removes dependent mass matrices.') - - @property - def curTModel(self): - if getattr(self, '_curTModel', None) is None: - self._curTModel = self.mapping.transform(self.curModel) - return self._curTModel - - @property - def curTModelDeriv(self): - if getattr(self, '_curTModelDeriv', None) is None: - self._curTModelDeriv = self.mapping.transformDeriv(self.curModel) - return self._curTModelDeriv + deleteTheseOnModelUpdate = ['_MeSigma', '_MeSigmaI'] def fields(self, m): self.curModel = m diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index e6eb38ce..b201d403 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -142,8 +142,8 @@ class ProblemFDEM_e(BaseFDEMProblem): return C.T*mui*C + 1j*omega(freq)*sig def getADeriv(self, freq, u, v, adjoint=False): - sig = self.curTModel - dsig_dm = self.curTModelDeriv + sig = self.curModel.transform + dsig_dm = self.curModel.transformDeriv dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig, v=u) if adjoint: @@ -222,8 +222,8 @@ class ProblemFDEM_b(BaseFDEMProblem): mui = self.MfMui C = self.mesh.edgeCurl - sig = self.curTModel - dsig_dm = self.curTModelDeriv + sig = self.curModel.transform + dsig_dm = self.curModel.transformDeriv #TODO: This only works if diagonal (no tensors)... dMeSigmaI_dI = - self.MeSigmaI**2 @@ -275,8 +275,8 @@ class ProblemFDEM_b(BaseFDEMProblem): def calcFieldsDeriv(self, sol, freq, fieldType, v, adjoint=False): b = sol if fieldType == 'e': - sig = self.curTModel - dsig_dm = self.curTModelDeriv + sig = self.curModel.transform + dsig_dm = self.curModel.transformDeriv C = self.mesh.edgeCurl mui = self.MfMui diff --git a/simpegEM/TDEM/BaseTDEM.py b/simpegEM/TDEM/BaseTDEM.py index c2d4e71b..4cb63c41 100644 --- a/simpegEM/TDEM/BaseTDEM.py +++ b/simpegEM/TDEM/BaseTDEM.py @@ -18,14 +18,16 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem): _FieldsTDEM_pair = FieldsTDEM #: used for the forward calculation only def fields(self, m): + if self.verbose: print '%s\nCalculating fields(m)\n%s'%('*'*50,'*'*50) self.curModel = m # Create a fields storage object F = self._FieldsTDEM_pair(self.mesh, self.survey) for tx in self.survey.txList: # Set the initial conditions F[tx,:,0] = tx.getInitialFields(self.mesh) - return self.forward(m, self.getRHS, self.calcFields, F=F) - + F = self.forward(m, self.getRHS, self.calcFields, F=F) + if self.verbose: print '%s\nDone calculating fields(m)\n%s'%('*'*50,'*'*50) + return F def forward(self, m, RHS, CalcFields, F=None): self.curModel = m @@ -39,11 +41,13 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem): if Ainv is not None: Ainv.clean() A = self.getA(tInd) - if self.verbose: print 'Factoring... (dt = ' + str(dt) + ')' + if self.verbose: print 'Factoring... (dt = %e)'%dt Ainv = self.Solver(A, **self.solverOpts) if self.verbose: print 'Done' rhs = RHS(tInd, F) + if self.verbose: print ' Solving... (tInd = %d)'%tInd sol = Ainv * rhs + if self.verbose: print ' Done...' if sol.ndim == 1: sol.shape = (sol.size,1) F[:,:,tInd+1] = CalcFields(sol, tInd) @@ -62,11 +66,13 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem): if Ainv is not None: Ainv.clean() A = self.getA(tInd) - if self.verbose: print 'Factoring... (dt = ' + str(dt) + ')' + if self.verbose: print 'Factoring (Adjoint)... (dt = %e)'%dt Ainv = self.Solver(A, **self.solverOpts) if self.verbose: print 'Done' rhs = RHS(tInd, F) + if self.verbose: print ' Solving (Adjoint)... (tInd = %d)'%tInd sol = Ainv * rhs + if self.verbose: print ' Done...' if sol.ndim == 1: sol.shape = (sol.size,1) F[:,:,tInd+1] = CalcFields(sol, tInd) @@ -88,12 +94,14 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem): * Compute \\\(\\\\vec{w} = -\\\mathbf{Q} \\\\vec{y}\\\) """ + if self.verbose: print '%s\nCalculating J(v)\n%s'%('*'*50,'*'*50) self.curModel = m if u is None: u = self.fields(m) p = self.Gvec(m, v, u) y = self.solveAh(m, p) Jv = self.survey.projectFieldsDeriv(u, v=y) + if self.verbose: print '%s\nDone calculating J(v)\n%s'%('*'*50,'*'*50) return - mkvc(Jv) def Jtvec(self, m, v, u=None): @@ -111,6 +119,7 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem): * Compute \\\(\\\\vec{w} = -\\\mathbf{G}^\\\\top y\\\) """ + if self.verbose: print '%s\nCalculating J^T(v)\n%s'%('*'*50,'*'*50) self.curModel = m if u is None: u = self.fields(m) @@ -121,5 +130,6 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem): p = self.survey.projectFieldsDeriv(u, v=v, adjoint=True) y = self.solveAht(m, p) w = self.Gtvec(m, y, u) + if self.verbose: print '%s\nDone calculating J^T(v)\n%s'%('*'*50,'*'*50) return - mkvc(w) diff --git a/simpegEM/TDEM/TDEM_b.py b/simpegEM/TDEM/TDEM_b.py index 1a8fdbb2..bd553df3 100644 --- a/simpegEM/TDEM/TDEM_b.py +++ b/simpegEM/TDEM/TDEM_b.py @@ -97,7 +97,7 @@ class ProblemTDEM_b(BaseTDEMProblem): # fake initial 'e' fields p[:, 'e', 0] = 0.0 - c = self.mesh.getEdgeInnerProductDeriv(self.curTModel)*(self.curTModelDeriv*vec) + c = self.mesh.getEdgeInnerProductDeriv(self.curModel.transform)*(self.curModel.transformDeriv*vec) for i in range(1,self.nT+1): # TODO: G[1] may be dependent on the model # for a galvanic source (deriv of the dc problem) @@ -127,7 +127,7 @@ class ProblemTDEM_b(BaseTDEMProblem): if nTx > 1: vu = vu.sum(axis=1) tmp += vu - p = -mkvc(self.curTModelDeriv.T*(self.mesh.getEdgeInnerProductDeriv(self.curTModel).T*tmp)) + p = -mkvc(self.curModel.transformDeriv.T*(self.mesh.getEdgeInnerProductDeriv(self.curModel.transform).T*tmp)) return p def solveAh(self, m, p): diff --git a/simpegEM/Tests/test_FDEM_analytics.py b/simpegEM/Tests/test_FDEM_analytics.py index 16c75cea..2928c50c 100644 --- a/simpegEM/Tests/test_FDEM_analytics.py +++ b/simpegEM/Tests/test_FDEM_analytics.py @@ -28,7 +28,12 @@ class FDEM_analyticTests(unittest.TestCase): prb = EM.FDEM.ProblemFDEM_b(mesh, mapping=mapping) prb.pair(survey) - prb.Solver = SolverLU + + try: + from pymatsolver import MumpsSolver + prb.Solver = MumpsSolver + except ImportError, e: + prb.Solver = SolverLU sig = 1e-1 sigma = np.ones(mesh.nC)*sig diff --git a/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py b/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py index 9ce9175c..9f6050d5 100644 --- a/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py +++ b/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py @@ -18,8 +18,7 @@ class TDEM_bDerivTests(unittest.TestCase): active = mesh.vectorCCz<0. activeMap = Maps.ActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz) - mapping = Maps.ComboMap(mesh, - [Maps.ExpMap, Maps.Vertical1DMap, activeMap]) + mapping = Maps.ExpMap(mesh) * Maps.Vertical1DMap(mesh) * activeMap rxOffset = 40. rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20), 'bz') @@ -32,6 +31,12 @@ class TDEM_bDerivTests(unittest.TestCase): self.prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)] # self.prb.timeSteps = [(1e-05, 100)] + try: + from pymatsolver import MumpsSolver + self.prb.Solver = MumpsSolver + except ImportError, e: + self.prb.Solver = SolverLU + self.sigma = np.ones(mesh.nCz)*1e-8 self.sigma[mesh.vectorCCz<0] = 1e-1 self.sigma = np.log(self.sigma[active]) diff --git a/simpegEM/Tests/test_TDEM_b_MultiTx_DerivAdjoint.py b/simpegEM/Tests/test_TDEM_b_MultiTx_DerivAdjoint.py index e800e2b6..cd054d34 100644 --- a/simpegEM/Tests/test_TDEM_b_MultiTx_DerivAdjoint.py +++ b/simpegEM/Tests/test_TDEM_b_MultiTx_DerivAdjoint.py @@ -18,8 +18,7 @@ class TDEM_bDerivTests(unittest.TestCase): active = mesh.vectorCCz<0. activeMap = Maps.ActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz) - mapping = Maps.ComboMap(mesh, - [Maps.ExpMap, Maps.Vertical1DMap, activeMap]) + mapping = Maps.ExpMap(mesh) * Maps.Vertical1DMap(mesh) * activeMap rxOffset = 40. rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20), 'bz') diff --git a/simpegEM/Tests/test_TDEM_combos.py b/simpegEM/Tests/test_TDEM_combos.py index 02597cc8..286438ab 100644 --- a/simpegEM/Tests/test_TDEM_combos.py +++ b/simpegEM/Tests/test_TDEM_combos.py @@ -20,7 +20,7 @@ def getProb(meshType='CYL',rxTypes='bx,bz',nTx=1): active = mesh.vectorCCz<0. activeMap = Maps.ActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz) - mapping = Maps.ComboMap(mesh, [Maps.ExpMap, Maps.Vertical1DMap, activeMap]) + mapping = Maps.ExpMap(mesh) * Maps.Vertical1DMap(mesh) * activeMap rxOffset = 40.