From 9ebbe7613db9b72f4ce0053b99f03f1a20718744 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sat, 26 Apr 2014 23:15:31 -0700 Subject: [PATCH] Forward problem working. Not yet tested with multi Tx. --- simpegEM/FDEM/FDEM.py | 10 +- simpegEM/TDEM/BaseTDEM.py | 62 +++--------- simpegEM/TDEM/SurveyTDEM.py | 101 ++++++++++++++++++- simpegEM/TDEM/TDEM_b.py | 15 +-- simpegEM/TDEM/__init__.py | 4 +- simpegEM/Tests/test_TDEM_forward_Analytic.py | 20 ++-- 6 files changed, 134 insertions(+), 78 deletions(-) diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index 3fc68dc2..16571fc3 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -8,7 +8,7 @@ def omega(freq): """Change frequency to angular frequency, omega""" return 2.*np.pi*freq -class BaseProblemFDEM(BaseEMProblem): +class BaseFDEMProblem(BaseEMProblem): """ We start by looking at Maxwell's equations in the electric field \\(\\vec{E}\\) and the magnetic flux density \\(\\vec{B}\\): @@ -106,7 +106,7 @@ class BaseProblemFDEM(BaseEMProblem): return Jtv -class ProblemFDEM_e(BaseProblemFDEM): +class ProblemFDEM_e(BaseFDEMProblem): """ By eliminating the magnetic flux density using @@ -127,7 +127,7 @@ class ProblemFDEM_e(BaseProblemFDEM): solType = 'e' def __init__(self, model, **kwargs): - BaseProblemFDEM.__init__(self, model, **kwargs) + BaseFDEMProblem.__init__(self, model, **kwargs) def getA(self, freq): """ @@ -197,14 +197,14 @@ class ProblemFDEM_e(BaseProblemFDEM): raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) -class ProblemFDEM_b(BaseProblemFDEM): +class ProblemFDEM_b(BaseFDEMProblem): """ Solving for b! """ solType = 'b' def __init__(self, model, **kwargs): - BaseProblemFDEM.__init__(self, model, **kwargs) + BaseFDEMProblem.__init__(self, model, **kwargs) def getA(self, freq): """ diff --git a/simpegEM/TDEM/BaseTDEM.py b/simpegEM/TDEM/BaseTDEM.py index d8c3d2ee..25a6d473 100644 --- a/simpegEM/TDEM/BaseTDEM.py +++ b/simpegEM/TDEM/BaseTDEM.py @@ -1,7 +1,7 @@ from SimPEG import Solver from SimPEG.Problem import BaseTimeProblem from simpegEM.Utils import Sources -from SurveyTDEM import FieldsTDEM +from SurveyTDEM import FieldsTDEM, SurveyTDEM from scipy.constants import mu_0 from SimPEG.Utils import sdiag, mkvc from SimPEG import Utils, Mesh @@ -9,49 +9,12 @@ from simpegEM.Base import BaseEMProblem import numpy as np -class MixinInitialFieldCalc(object): - """docstring for MixinInitialFieldCalc""" - - storeTheseFields = 'b' - - def getInitialFields(self): - if self.survey.txType == 'VMD_MVP': - # Vertical magnetic dipole, magnetic vector potential - F = self._getInitialFields_VMD_MVP() - else: - exStr = 'Invalid txType: ' + str(self.survey.txType) - raise Exception(exStr) - return F - - def _getInitialFields_VMD_MVP(self): - 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.nT, store=self.storeTheseFields) - - # Set initial B - F.b0 = self.mesh.edgeCurl*MVP - - return F - - -class ProblemBaseTDEM(MixinInitialFieldCalc, BaseTimeProblem, BaseEMProblem): +class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem): """docstring for ProblemTDEM1D""" def __init__(self, mesh, mapping=None, **kwargs): BaseTimeProblem.__init__(self, mesh, mapping=mapping, **kwargs) + surveyPair = SurveyTDEM def calcFields(self, sol, solType, tInd): @@ -65,18 +28,18 @@ class ProblemBaseTDEM(MixinInitialFieldCalc, BaseTimeProblem, BaseEMProblem): return {'b':b, 'e':e} - Solver = Solver - solveOpts = {} - def fields(self, m): self.curModel = m - F = self.getInitialFields() + # Create a fields storage object + F = FieldsTDEM(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) def forward(self, m, RHS, CalcFields, F=None): - if F is None: - F = FieldsTDEM(self.mesh, self.survey.nTx, self.nT, store=self.storeTheseFields) + F = F or FieldsTDEM(self.mesh, self.survey) dtFact = None for tInd, dt in enumerate(self.timeSteps): @@ -84,14 +47,13 @@ class ProblemBaseTDEM(MixinInitialFieldCalc, BaseTimeProblem, BaseEMProblem): dtFact = dt A = self.getA(tInd) # print 'Factoring... (dt = ' + str(dt) + ')' - Asolve = self.Solver(A, **self.solveOpts) + Asolve = self.Solver(A, **self.solverOpts) # print 'Done' rhs = RHS(tInd, F) sol = Asolve.solve(rhs) if sol.ndim == 1: sol.shape = (sol.size,1) - newFields = CalcFields(sol, self.solType, tInd) - F.update(newFields, tInd) + F[:,:,tInd+1] = CalcFields(sol, self.solType, tInd) return F def adjoint(self, m, RHS, CalcFields, F=None): @@ -104,7 +66,7 @@ class ProblemBaseTDEM(MixinInitialFieldCalc, BaseTimeProblem, BaseEMProblem): dtFact = dt A = self.getA(tInd) # print 'Factoring... (dt = ' + str(dt) + ')' - Asolve = Solver(A, options=self.solveOpts) + Asolve = Solver(A, options=self.solverOpts) # print 'Done' rhs = RHS(tInd, F) sol = Asolve.solve(rhs) diff --git a/simpegEM/TDEM/SurveyTDEM.py b/simpegEM/TDEM/SurveyTDEM.py index 81c2421e..d24fc18b 100644 --- a/simpegEM/TDEM/SurveyTDEM.py +++ b/simpegEM/TDEM/SurveyTDEM.py @@ -1,5 +1,102 @@ -from SimPEG import Utils, np +from SimPEG import Utils, Survey, np from SimPEG.Survey import BaseSurvey +from simpegEM.Utils import Sources + + +class RxTDEM(Survey.BaseTimeRx): + + knownRxTypes = { + 'ex':['e', 'Ex'], + 'ey':['e', 'Ey'], + 'ez':['e', 'Ez'], + + 'bx':['b', 'Fx'], + 'by':['b', 'Fy'], + 'bz':['b', 'Fz'], + } + + def __init__(self, locs, times, rxType): + Survey.BaseTimeRx.__init__(self, locs, times, rxType) + + @property + def projField(self): + """Field Type projection (e.g. e b ...)""" + return self.knownRxTypes[self.rxType][0] + + @property + def projGLoc(self): + """Grid Location projection (e.g. Ex Fy ...)""" + return self.knownRxTypes[self.rxType][1] + + def projectFields(self, tx, mesh, timeMesh, u): + P = self.getP(mesh, timeMesh) + u_part = Utils.mkvc(u[tx, self.projField, :]) + return P*u_part + + def projectFieldsDeriv(self, tx, mesh, timeMesh, u, v, adjoint=False): + P = self.getP(mesh, timeMesh) + + if not adjoint: + return P * v + elif adjoint: + return P.T * v + + +class FieldsTDEM(Survey.TimeFields): + """Fancy Field Storage for a TDEM survey.""" + knownFields = {'b': 'F', 'e': 'E'} + + +class TxTDEM(Survey.BaseTx): + rxPair = RxTDEM + knownTxTypes = ['VMD_MVP'] + + def getInitialFields(self, mesh): + F0 = getattr(self, '_getInitialFields_' + self.txType)(mesh) + return F0 + + def _getInitialFields_VMD_MVP(self, mesh): + """Vertical magnetic dipole, magnetic vector potential""" + if mesh._meshType is 'CYL': + if mesh.isSymmetric: + MVP = Sources.MagneticDipoleVectorPotential(self.loc, mesh.gridEy, 'y') + else: + raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') + elif mesh._meshType is 'TENSOR': + MVPx = Sources.MagneticDipoleVectorPotential(self.loc, mesh.gridEx, 'x') + MVPy = Sources.MagneticDipoleVectorPotential(self.loc, mesh.gridEy, 'y') + MVPz = Sources.MagneticDipoleVectorPotential(self.loc, mesh.gridEz, 'z') + MVP = np.concatenate((MVPx, MVPy, MVPz)) + else: + raise Exception('Unknown mesh for VMD') + + return {"b": mesh.edgeCurl*MVP} + + def getJs(self, time): + return None + +class SurveyTDEM(Survey.BaseSurvey): + """ + docstring for SurveyTDEM + """ + + txPair = TxTDEM + + def __init__(self, txList, **kwargs): + # Sort these by frequency + self.txList = txList + Survey.BaseSurvey.__init__(self, **kwargs) + + def projectFields(self, u): + data = Survey.Data(self) + for tx in self.txList: + for rx in tx.rxList: + data[tx, rx] = rx.projectFields(tx, self.mesh, self.prob.timeMesh, u) + return data + + def projectFieldsDeriv(self, u): + raise Exception('Use Transmitters to project fields deriv.') + class SurveyTDEM1D(BaseSurvey): """ @@ -52,7 +149,7 @@ class SurveyTDEM1D(BaseSurvey): _Qrx = None -class FieldsTDEM(object): +class FieldsTDEM_OLD(object): """docstring for FieldsTDEM""" phi0 = None #: Initial electric potential diff --git a/simpegEM/TDEM/TDEM_b.py b/simpegEM/TDEM/TDEM_b.py index 8bd60161..2e8da6ca 100644 --- a/simpegEM/TDEM/TDEM_b.py +++ b/simpegEM/TDEM/TDEM_b.py @@ -1,9 +1,9 @@ -from BaseTDEM import ProblemBaseTDEM +from BaseTDEM import BaseTDEMProblem from SimPEG.Utils import mkvc import numpy as np -from SurveyTDEM import SurveyTDEM1D, FieldsTDEM +from SurveyTDEM import SurveyTDEM, FieldsTDEM -class ProblemTDEM_b(ProblemBaseTDEM): +class ProblemTDEM_b(BaseTDEMProblem): """ Time-Domain EM problem - B-formulation @@ -16,11 +16,11 @@ class ProblemTDEM_b(ProblemBaseTDEM): with \\\(\\b\\\) defined on cell faces and \\\(\e\\\) defined on edges. """ def __init__(self, mesh, mapping=None, **kwargs): - ProblemBaseTDEM.__init__(self, mesh, mapping=mapping, **kwargs) + BaseTDEMProblem.__init__(self, mesh, mapping=mapping, **kwargs) solType = 'b' - surveyPair = SurveyTDEM1D + surveyPair = SurveyTDEM #################################################### # Internal Methods @@ -32,13 +32,14 @@ class ProblemTDEM_b(ProblemBaseTDEM): :rtype: scipy.sparse.csr_matrix :return: A """ - dt = self.timeSteps[tInd] 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.timeSteps[tInd] - return (1.0/dt)*self.MfMui*F.get_b(tInd-1) + B_n = np.concatenate([F[tx,'b',tInd] for tx in self.survey.txList], axis=1) + RHS = (1.0/dt)*self.MfMui*B_n + return RHS #################################################### diff --git a/simpegEM/TDEM/__init__.py b/simpegEM/TDEM/__init__.py index 8a64e946..fed1a2dc 100644 --- a/simpegEM/TDEM/__init__.py +++ b/simpegEM/TDEM/__init__.py @@ -1,3 +1,3 @@ -from BaseTDEM import ProblemBaseTDEM -from SurveyTDEM import SurveyTDEM1D, FieldsTDEM +from SurveyTDEM import SurveyTDEM, FieldsTDEM, RxTDEM, TxTDEM +from BaseTDEM import BaseTDEMProblem from TDEM_b import ProblemTDEM_b diff --git a/simpegEM/Tests/test_TDEM_forward_Analytic.py b/simpegEM/Tests/test_TDEM_forward_Analytic.py index b2733d86..ed308906 100644 --- a/simpegEM/Tests/test_TDEM_forward_Analytic.py +++ b/simpegEM/Tests/test_TDEM_forward_Analytic.py @@ -22,15 +22,10 @@ def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=[1e-5, actMap = Maps.ActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz) mapping = Maps.ComboMap(mesh, [Maps.ExpMap, Maps.Vertical1DMap, actMap]) + rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-5,-4, 21), 'bz') + tx = EM.TDEM.TxTDEM(np.array([0., 0., 0.]), 'VMD_MVP', [rx]) - opts = {'txLoc':np.array([0., 0., 0.]), - 'txType':'VMD_MVP', - 'rxLoc':np.array([rxOffset, 0., 0.]), - 'rxType':'bz', - 'timeCh':np.logspace(-5,-4, 21), - } - - survey = EM.TDEM.SurveyTDEM1D(**opts) + survey = EM.TDEM.SurveyTDEM([tx]) prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping) prb.Solver = Utils.SolverUtils.DSolverWrap(sp.linalg.splu, factorize=True) # try: @@ -46,16 +41,17 @@ def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=[1e-5, sigma = np.log(sigma[active]) prb.pair(survey) - bz_ana = mu_0*EM.Utils.Ana.hzAnalyticDipoleT(survey.rxLoc[0]+1e-3, prb.times[1:], sig_half) + bz_ana = mu_0*EM.Utils.Ana.hzAnalyticDipoleT(rx.locs[0][0]+1e-3, rx.times, sig_half) bz_calc = survey.dpred(sigma) - ind = np.logical_and(prb.times[1:] > bounds[0],prb.times[1:] < bounds[1]) + + ind = np.logical_and(rx.times > bounds[0],rx.times < bounds[1]) log10diff = np.linalg.norm(np.log10(np.abs(bz_calc[ind])) - np.log10(np.abs(bz_ana[ind])))/np.linalg.norm(np.log10(np.abs(bz_ana[ind]))) print 'Difference: ', log10diff if showIt == True: - plt.loglog(prb.times[1:][bz_calc>0], bz_calc[bz_calc>0], 'r', prb.times[1:][bz_calc<0], -bz_calc[bz_calc<0], 'r--') - plt.loglog(prb.times[1:], abs(bz_ana), 'b*') + plt.loglog(rx.times[bz_calc>0], bz_calc[bz_calc>0], 'r', rx.times[bz_calc<0], -bz_calc[bz_calc<0], 'r--') + plt.loglog(rx.times, abs(bz_ana), 'b*') plt.title('sig_half = %e'%sig_half) plt.show()