diff --git a/simpegEM/Base.py b/simpegEM/Base.py new file mode 100644 index 00000000..da7c97fd --- /dev/null +++ b/simpegEM/Base.py @@ -0,0 +1,68 @@ +from SimPEG import Survey, Problem, Utils, np, sp, Solver as SimpegSolver +from scipy.constants import mu_0 + +class BaseEMProblem(Problem.BaseProblem): + + def __init__(self, mesh, **kwargs): + Problem.BaseProblem.__init__(self, mesh, **kwargs) + + solType = None + storeTheseFields = ['e', 'b'] + + surveyPair = Survey.BaseSurvey + dataPair = Survey.Data + + Solver = SimpegSolver + solverOpts = {} + + #################################################### + # Mass Matrices + #################################################### + + @property + def MfMui(self): + #TODO: assuming constant mu + if getattr(self, '_MfMui', None) is None: + self._MfMui = self.mesh.getFaceInnerProduct(1/mu_0) + return self._MfMui + + @property + def Me(self): + if getattr(self, '_Me', None) is None: + self._Me = self.mesh.getEdgeInnerProduct() + return self._Me + + @property + def MeSigma(self): + #TODO: hardcoded to sigma as the model + if getattr(self, '_MeSigma', None) is None: + sigma = self.curTModel + self._MeSigma = self.mesh.getEdgeInnerProduct(sigma) + return self._MeSigma + + @property + def MeSigmaI(self): + #TODO: hardcoded to sigma as the model + if getattr(self, '_MeSigmaI', None) is None: + sigma = self.curTModel + 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 + + def fields(self, m): + self.curModel = m + F = self.forward(m, self.getRHS, self.calcFields) + return F diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index f88af098..3fc68dc2 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -2,12 +2,13 @@ from SimPEG import Survey, Problem, Utils, np, sp, Solver as SimpegSolver from scipy.constants import mu_0 from SurveyFDEM import SurveyFDEM, FieldsFDEM from simpegEM.Utils import Sources +from simpegEM.Base import BaseEMProblem def omega(freq): """Change frequency to angular frequency, omega""" return 2.*np.pi*freq -class BaseProblemFDEM(Problem.BaseProblem): +class BaseProblemFDEM(BaseEMProblem): """ We start by looking at Maxwell's equations in the electric field \\(\\vec{E}\\) and the magnetic flux density \\(\\vec{B}\\): @@ -17,69 +18,8 @@ class BaseProblemFDEM(Problem.BaseProblem): \\nabla \\times \\mu^{-1} \\vec{B} - \\sigma \\vec{E} = \\vec{J_s} """ - def __init__(self, model, **kwargs): - Problem.BaseProblem.__init__(self, model, **kwargs) - - solType = None - storeTheseFields = ['e', 'b'] surveyPair = SurveyFDEM - dataPair = Survey.Data - - Solver = SimpegSolver - solverOpts = {} - - #################################################### - # Mass Matrices - #################################################### - - @property - def MfMui(self): - #TODO: assuming constant mu - if getattr(self, '_MfMui', None) is None: - self._MfMui = self.mesh.getFaceInnerProduct(1/mu_0) - return self._MfMui - - @property - def Me(self): - if getattr(self, '_Me', None) is None: - self._Me = self.mesh.getEdgeInnerProduct() - return self._Me - - @property - def MeSigma(self): - #TODO: hardcoded to sigma as the model - if getattr(self, '_MeSigma', None) is None: - sigma = self.curTModel - self._MeSigma = self.mesh.getEdgeInnerProduct(sigma) - return self._MeSigma - - @property - def MeSigmaI(self): - #TODO: hardcoded to sigma as the model - if getattr(self, '_MeSigmaI', None) is None: - sigma = self.curTModel - 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 - - def fields(self, m): - self.curModel = m - F = self.forward(m, self.getRHS, self.calcFields) - return F def forward(self, m, RHS, CalcFields): diff --git a/simpegEM/FDEM/RHSem.py b/simpegEM/FDEM/RHSem.py index 2d498e1b..fb1f8f0e 100644 --- a/simpegEM/FDEM/RHSem.py +++ b/simpegEM/FDEM/RHSem.py @@ -50,10 +50,6 @@ def path2edgeModel(mesh, pts): edgeModel = np.r_[edm_x, edm_y, edm_z] return edgeModel -def rho(x1, y1, x, y): - r = np.sqrt((x-x1)**2+(y-y1)**2) - return r - def MMRhalf(loc1, loc2, x, y): """ Anaytic function for MMR response (B^{1D}) - loc1=(x1,y1): x, y location for (+) charge @@ -67,6 +63,9 @@ def MMRhalf(loc1, loc2, x, y): y2=loc2[1] mu0 = 4*np.pi*1e-7 I = 1 + + rho = lambda x1, y1, x, y: np.sqrt((x-x1)**2+(y-y1)**2) + By =mu0*I/(4*np.pi)*np.array((x-x1)/rho(x1,y1,x,y)**2-(x-x2)/rho(x2,y2,x,y)**2) Bx =mu0*I/(4*np.pi)*np.array(-(y-y1)/rho(x1,y1,x,y)**2+(y-y2)/rho(x2,y2,x,y)**2) diff --git a/simpegEM/TDEM/BaseTDEM.py b/simpegEM/TDEM/BaseTDEM.py index d3b648da..d8c3d2ee 100644 --- a/simpegEM/TDEM/BaseTDEM.py +++ b/simpegEM/TDEM/BaseTDEM.py @@ -5,6 +5,7 @@ from SurveyTDEM import FieldsTDEM from scipy.constants import mu_0 from SimPEG.Utils import sdiag, mkvc from SimPEG import Utils, Mesh +from simpegEM.Base import BaseEMProblem import numpy as np @@ -46,44 +47,12 @@ class MixinInitialFieldCalc(object): return F -class ProblemBaseTDEM(MixinInitialFieldCalc, BaseTimeProblem): +class ProblemBaseTDEM(MixinInitialFieldCalc, BaseTimeProblem, BaseEMProblem): """docstring for ProblemTDEM1D""" def __init__(self, mesh, mapping=None, **kwargs): BaseTimeProblem.__init__(self, mesh, mapping=mapping, **kwargs) - #################################################### - # Physical Properties - #################################################### - - @property - def sigma(self): - return self._sigma - @sigma.setter - def sigma(self, value): - self._sigma = value - _sigma = None - - #################################################### - # Mass Matrices - #################################################### - - @property - def MfMui(self): return self._MfMui - - @property - def MeSigma(self): return self._MeSigma - - @property - def MeSigmaI(self): return self._MeSigmaI - - def makeMassMatrices(self, m): - sig = self.mapping.transform(m) - self._MeSigma = self.mesh.getEdgeInnerProduct(sig) - self._MeSigmaI = Utils.sdInv(self.MeSigma) - self._MfMui = self.mesh.getFaceInnerProduct(1.0/mu_0) - - def calcFields(self, sol, solType, tInd): if solType == 'b': @@ -100,7 +69,7 @@ class ProblemBaseTDEM(MixinInitialFieldCalc, BaseTimeProblem): solveOpts = {} def fields(self, m): - self.makeMassMatrices(m) + self.curModel = m F = self.getInitialFields() return self.forward(m, self.getRHS, self.calcFields, F=F) @@ -111,7 +80,7 @@ class ProblemBaseTDEM(MixinInitialFieldCalc, BaseTimeProblem): dtFact = None for tInd, dt in enumerate(self.timeSteps): - if dt!=dtFact: + if dt != dtFact: dtFact = dt A = self.getA(tInd) # print 'Factoring... (dt = ' + str(dt) + ')' @@ -131,7 +100,7 @@ class ProblemBaseTDEM(MixinInitialFieldCalc, BaseTimeProblem): dtFact = None for tInd, dt in reversed(list(enumerate(self.timeSteps))): - if dt!=dtFact: + if dt != dtFact: dtFact = dt A = self.getA(tInd) # print 'Factoring... (dt = ' + str(dt) + ')' diff --git a/simpegEM/TDEM/TDEM_b.py b/simpegEM/TDEM/TDEM_b.py index e41c5644..8bd60161 100644 --- a/simpegEM/TDEM/TDEM_b.py +++ b/simpegEM/TDEM/TDEM_b.py @@ -109,7 +109,7 @@ class ProblemTDEM_b(ProblemBaseTDEM): e = self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui*b - self.MeSigmaI*p.get_e(tInd) return {'b':b, 'e':e} - self.makeMassMatrices(m) + self.curModel = m return self.forward(m, AhRHS, AhCalcFields) def solveAht(self, m, p): @@ -126,16 +126,16 @@ class ProblemTDEM_b(ProblemBaseTDEM): e = self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui*b - self.MeSigmaI*p.get_e(tInd) return {'b':b, 'e':e} - self.makeMassMatrices(m) + self.curModel = m return self.adjoint(m, AhtRHS, AhtCalcFields) #################################################### # Functions for tests #################################################### - def AhVec(self, sigma, vec): + def AhVec(self, m, vec): """ - :param numpy.array sigma: Conductivity model + :param numpy.array m: Conductivity model :param simpegEM.TDEM.FieldsTDEM vec: Fields object :rtype: simpegEM.TDEM.FieldsTDEM :return: f @@ -167,7 +167,7 @@ class ProblemTDEM_b(ProblemBaseTDEM): \\right] \\\\ """ - self.makeMassMatrices(sigma) + self.curModel = m dt = self.timeSteps[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) @@ -182,9 +182,9 @@ class ProblemTDEM_b(ProblemBaseTDEM): f.set_e(e, i) return f - def AhtVec(self, sigma, vec): + def AhtVec(self, m, vec): """ - :param numpy.array sigma: Conductivity model + :param numpy.array m: Conductivity model :param simpegEM.TDEM.FieldsTDEM vec: Fields object :rtype: simpegEM.TDEM.FieldsTDEM :return: f @@ -215,7 +215,7 @@ class ProblemTDEM_b(ProblemBaseTDEM): \end{array} \\right] \\\\ """ - self.makeMassMatrices(sigma) + self.curModel = m f = FieldsTDEM(self.mesh, 1, self.nT, 'b') for i in range(self.nT-1): b = 1.0/self.timeSteps[i]*self.MfMui*vec.get_b(i) + self.MfMui*self.mesh.edgeCurl*vec.get_e(i) - 1.0/self.timeSteps[i+1]*self.MfMui*vec.get_b(i+1) @@ -257,9 +257,9 @@ if __name__ == '__main__': # prb.setTimes([1e-5, 5e-5, 2.5e-4], [10, 10, 10]) prb.timeSteps = [(1e-5, 10)] prb.pair(survey) - sigma = np.random.rand(mesh.nCz) + m = np.random.rand(mesh.nCz) - print survey.dpred(sigma) + print survey.dpred(m) diff --git a/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py b/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py index a1665324..1eccc566 100644 --- a/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py +++ b/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py @@ -17,7 +17,7 @@ class TDEM_bDerivTests(unittest.TestCase): mesh = Mesh.CylMesh([hx,1,hy], '00C') active = mesh.vectorCCz<0. - activeMap = Maps.ActiveCells(mesh, active, -8, nC=mesh.nCz) + activeMap = Maps.ActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz) mapping = Maps.ComboMap(mesh, [Maps.ExpMap, Maps.Vertical1DMap, activeMap]) @@ -76,7 +76,7 @@ class TDEM_bDerivTests(unittest.TestCase): prb = self.prb prb.timeSteps = [1e-05] sigma = self.sigma - prb.makeMassMatrices(sigma) + prb.curModel = sigma dt = prb.timeSteps[0] a11 = 1/dt*prb.MfMui*sp.eye(prb.mesh.nF) @@ -97,7 +97,7 @@ class TDEM_bDerivTests(unittest.TestCase): prb.timeSteps = [1e-05] sigma = self.sigma - prb.makeMassMatrices(sigma) + prb.curModel = sigma dt = prb.timeSteps[0] a11 = 1/dt*prb.MfMui*sp.eye(prb.mesh.nF) @@ -120,7 +120,7 @@ class TDEM_bDerivTests(unittest.TestCase): prb = self.prb mesh = self.prb.mesh sigma = self.sigma - self.prb.makeMassMatrices(sigma) + self.prb.curModel = sigma f = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nT, 'b') for i in range(f.nT): @@ -242,6 +242,7 @@ class TDEM_bDerivTests(unittest.TestCase): V1 = np.linalg.norm(f3.fieldVec()-f1.fieldVec()) V2 = np.linalg.norm(f1.fieldVec()) print V1, V2 + print 'I am gunna fail this one: boo. :(' self.assertLess(V1/V2, 1e-6) def test_adjointsolveAhVssolveAht(self): diff --git a/simpegEM/__init__.py b/simpegEM/__init__.py index 13667f2f..5beb775b 100644 --- a/simpegEM/__init__.py +++ b/simpegEM/__init__.py @@ -2,3 +2,4 @@ import Utils import TDEM import FDEM +import Base