From 9ebbe7613db9b72f4ce0053b99f03f1a20718744 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sat, 26 Apr 2014 23:15:31 -0700 Subject: [PATCH 1/8] 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() From 9f82ffff7bb43eeb21478ef44fc698404a37a5f0 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sat, 26 Apr 2014 23:22:02 -0700 Subject: [PATCH 2/8] column concat initial fields. --- simpegEM/TDEM/TDEM_b.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/simpegEM/TDEM/TDEM_b.py b/simpegEM/TDEM/TDEM_b.py index 2e8da6ca..e96525f2 100644 --- a/simpegEM/TDEM/TDEM_b.py +++ b/simpegEM/TDEM/TDEM_b.py @@ -37,7 +37,7 @@ class ProblemTDEM_b(BaseTDEMProblem): def getRHS(self, tInd, F): dt = self.timeSteps[tInd] - B_n = np.concatenate([F[tx,'b',tInd] for tx in self.survey.txList], axis=1) + B_n = np.c_[[F[tx,'b',tInd] for tx in self.survey.txList]].T RHS = (1.0/dt)*self.MfMui*B_n return RHS From 45884ba865c950162f5df72e3d2ab632ac6c8ccb Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sun, 27 Apr 2014 20:51:32 -0700 Subject: [PATCH 3/8] Updates to TDEM (Jtvec still not working.) --- simpegEM/TDEM/BaseTDEM.py | 4 +- simpegEM/TDEM/SurveyTDEM.py | 287 ++++++++------ simpegEM/TDEM/TDEM_b.py | 85 ++-- simpegEM/Tests/test_TDEM_b_DerivAdjoint.py | 396 ++++++++++--------- simpegEM/Tests/test_TDEM_forward_Analytic.py | 5 - 5 files changed, 404 insertions(+), 373 deletions(-) diff --git a/simpegEM/TDEM/BaseTDEM.py b/simpegEM/TDEM/BaseTDEM.py index 25a6d473..12911847 100644 --- a/simpegEM/TDEM/BaseTDEM.py +++ b/simpegEM/TDEM/BaseTDEM.py @@ -58,7 +58,7 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem): def adjoint(self, m, RHS, CalcFields, F=None): if F is None: - F = FieldsTDEM(self.mesh, self.survey.nTx, self.nT, store=self.storeTheseFields) + F = FieldsTDEM(self.mesh, self.survey) dtFact = None for tInd, dt in reversed(list(enumerate(self.timeSteps))): @@ -73,6 +73,6 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem): if sol.ndim == 1: sol.shape = (sol.size,1) newFields = CalcFields(sol, self.solType, tInd) - F.update(newFields, tInd) + F[:,:,tInd] = newFields return F diff --git a/simpegEM/TDEM/SurveyTDEM.py b/simpegEM/TDEM/SurveyTDEM.py index d24fc18b..3bc5fde6 100644 --- a/simpegEM/TDEM/SurveyTDEM.py +++ b/simpegEM/TDEM/SurveyTDEM.py @@ -37,15 +37,32 @@ class RxTDEM(Survey.BaseTimeRx): P = self.getP(mesh, timeMesh) if not adjoint: - return P * v + return P * Utils.mkvc(v[tx, self.projField, :]) elif adjoint: - return P.T * v + Ptv = P.T * v[tx, self] + return Ptv class FieldsTDEM(Survey.TimeFields): """Fancy Field Storage for a TDEM survey.""" knownFields = {'b': 'F', 'e': 'E'} + def tovec(self): + nTx, nF, nE = self.survey.nTx, self.mesh.nF, self.mesh.nE + u = np.empty(0 if nTx == 1 else (0, nTx)) + + for i in range(self.survey.prob.nT): + if 'b' in self: + b = self[:,'b',i+1] + else: + b = np.zeros(nF if nTx == 1 else (nF, nTx)) + + if 'e' in self: + e = self[:,'e',i+1] + else: + e = np.zeros(nE if nTx == 1 else (nE, nTx)) + u = np.r_[u, b, e] + return u class TxTDEM(Survey.BaseTx): rxPair = RxTDEM @@ -94,132 +111,148 @@ class SurveyTDEM(Survey.BaseSurvey): 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.') + def projectFieldsDeriv(self, u, v=None, adjoint=False): + assert v is not None, 'v to multiply must be provided.' - -class SurveyTDEM1D(BaseSurvey): - """ - docstring for SurveyTDEM1D - """ - - txLoc = None #: txLoc - txType = None #: txType - rxLoc = None #: rxLoc - rxType = None #: rxType - timeCh = None #: timeCh - nTx = 1 #: Number of transmitters - - @property - def nTimeCh(self): - """Number of time channels""" - return self.timeCh.size - - def __init__(self, **kwargs): - BaseSurvey.__init__(self, **kwargs) - Utils.setKwargs(self, **kwargs) - - def projectFields(self, u): - #TODO: this is hardcoded to 1Tx - return self.Qrx.dot(u.b[:,:,0].T).T - - def projectFieldsAdjoint(self, d): - # TODO: make the following self.nTimeCh - d = d.reshape((self.prob.nT, self.nTx), order='F') - #TODO: *Qtime.T need to multiply by a time projection. (outside for loop??) - ii = 0 - F = FieldsTDEM(self.prob.mesh, self.nTx, self.prob.nT, 'b') - for ii in range(self.prob.nT): - b = self.Qrx.T*d[ii,:] - F.set_b(b, ii) - F.set_e(np.zeros((self.prob.mesh.nE,self.nTx)), ii) - return F - - #################################################### - # Interpolation Matrices - #################################################### - - @property - def Qrx(self): - if self._Qrx is None: - if self.rxType == 'bz': - locType = 'Fz' - self._Qrx = self.prob.mesh.getInterpolationMat(self.rxLoc, locType=locType) - return self._Qrx - _Qrx = None - - -class FieldsTDEM_OLD(object): - """docstring for FieldsTDEM""" - - phi0 = None #: Initial electric potential - A0 = None #: Initial magnetic vector potential - e0 = None #: Initial electric field - b0 = None #: Initial magnetic flux density - j0 = None #: Initial current density - h0 = None #: Initial magnetic field - - phi = None #: Electric potential - A = None #: Magnetic vector potential - e = None #: Electric field - b = None #: Magnetic flux density - j = None #: Current density - h = None #: Magnetic field - - def __init__(self, mesh, nTx, nT, store='b'): - - self.nT = nT #: Number of times - self.nTx = nTx #: Number of transmitters - self.mesh = mesh - - def update(self, newFields, tInd): - self.set_b(newFields['b'], tInd) - self.set_e(newFields['e'], tInd) - - def fieldVec(self): - u = np.ndarray((0, self.nTx)) - for i in range(self.nT): - u = np.r_[u, self.get_b(i), self.get_e(i)] - if self.nTx == 1: - u = u.flatten() - return u - - #################################################### - # Get Methods - #################################################### - - def get_b(self, ind): - if ind == -1: - return self.b0 + if not adjoint: + data = Survey.Data(self) + for tx in self.txList: + for rx in tx.rxList: + data[tx, rx] = rx.projectFieldsDeriv(tx, self.mesh, self.prob.timeMesh, u, v) + return data else: - return self.b[ind,:,:] - - def get_e(self, ind): - if ind == -1: - return self.e0 - else: - return self.e[ind,:,:] - - #################################################### - # Set Methods - #################################################### - - def set_b(self, b, ind): - if self.b is None: - self.b = np.zeros((self.nT, np.sum(self.mesh.nF), self.nTx)) - self.b[:] = np.nan - if len(b.shape) == 1: - b = b[:, np.newaxis] - self.b[ind,:,:] = b - - def set_e(self, e, ind): - if self.e is None: - self.e = np.zeros((self.nT, np.sum(self.mesh.nE), self.nTx)) - self.e[:] = np.nan - if len(e.shape) == 1: - e = e[:, np.newaxis] - self.e[ind,:,:] = e + f = FieldsTDEM(self.mesh, self) + for tx in self.txList: + for rx in tx.rxList: + Ptv = rx.projectFieldsDeriv(tx, self.mesh, self.prob.timeMesh, u, v, adjoint=True) + Ptv = Ptv.reshape((-1, 1, self.prob.timeMesh.nN), order='F') + f[tx, rx.projField, :] = Ptv + return f - def __contains__(self, key): - return key in self.children + +# class SurveyTDEM1D(BaseSurvey): +# """ +# docstring for SurveyTDEM1D +# """ + +# txLoc = None #: txLoc +# txType = None #: txType +# rxLoc = None #: rxLoc +# rxType = None #: rxType +# timeCh = None #: timeCh +# nTx = 1 #: Number of transmitters + +# @property +# def nTimeCh(self): +# """Number of time channels""" +# return self.timeCh.size + +# def __init__(self, **kwargs): +# BaseSurvey.__init__(self, **kwargs) +# Utils.setKwargs(self, **kwargs) + +# def projectFields(self, u): +# #TODO: this is hardcoded to 1Tx +# return self.Qrx.dot(u.b[:,:,0].T).T + +# def projectFieldsAdjoint(self, d): +# # TODO: make the following self.nTimeCh +# d = d.reshape((self.prob.nT, self.nTx), order='F') +# #TODO: *Qtime.T need to multiply by a time projection. (outside for loop??) +# ii = 0 +# F = FieldsTDEM(self.prob.mesh, self.nTx, self.prob.nT, 'b') +# for ii in range(self.prob.nT): +# b = self.Qrx.T*d[ii,:] +# F.set_b(b, ii) +# F.set_e(np.zeros((self.prob.mesh.nE,self.nTx)), ii) +# return F + +# #################################################### +# # Interpolation Matrices +# #################################################### + +# @property +# def Qrx(self): +# if self._Qrx is None: +# if self.rxType == 'bz': +# locType = 'Fz' +# self._Qrx = self.prob.mesh.getInterpolationMat(self.rxLoc, locType=locType) +# return self._Qrx +# _Qrx = None + + +# class FieldsTDEM_OLD(object): +# """docstring for FieldsTDEM""" + +# phi0 = None #: Initial electric potential +# A0 = None #: Initial magnetic vector potential +# e0 = None #: Initial electric field +# b0 = None #: Initial magnetic flux density +# j0 = None #: Initial current density +# h0 = None #: Initial magnetic field + +# phi = None #: Electric potential +# A = None #: Magnetic vector potential +# e = None #: Electric field +# b = None #: Magnetic flux density +# j = None #: Current density +# h = None #: Magnetic field + +# def __init__(self, mesh, nTx, nT, store='b'): + +# self.nT = nT #: Number of times +# self.nTx = nTx #: Number of transmitters +# self.mesh = mesh + +# def update(self, newFields, tInd): +# self.set_b(newFields['b'], tInd) +# self.set_e(newFields['e'], tInd) + +# def fieldVec(self): +# u = np.ndarray((0, self.nTx)) +# for i in range(self.nT): +# u = np.r_[u, self.get_b(i), self.get_e(i)] +# if self.nTx == 1: +# u = u.flatten() +# return u + +# #################################################### +# # Get Methods +# #################################################### + +# def get_b(self, ind): +# if ind == -1: +# return self.b0 +# else: +# return self.b[ind,:,:] + +# def get_e(self, ind): +# if ind == -1: +# return self.e0 +# else: +# return self.e[ind,:,:] + +# #################################################### +# # Set Methods +# #################################################### + +# def set_b(self, b, ind): +# if self.b is None: +# self.b = np.zeros((self.nT, np.sum(self.mesh.nF), self.nTx)) +# self.b[:] = np.nan +# if len(b.shape) == 1: +# b = b[:, np.newaxis] +# self.b[ind,:,:] = b + +# def set_e(self, e, ind): +# if self.e is None: +# self.e = np.zeros((self.nT, np.sum(self.mesh.nE), self.nTx)) +# self.e[:] = np.nan +# if len(e.shape) == 1: +# e = e[:, np.newaxis] +# self.e[ind,:,:] = e + + +# def __contains__(self, key): +# return key in self.children diff --git a/simpegEM/TDEM/TDEM_b.py b/simpegEM/TDEM/TDEM_b.py index e96525f2..787c1442 100644 --- a/simpegEM/TDEM/TDEM_b.py +++ b/simpegEM/TDEM/TDEM_b.py @@ -51,12 +51,17 @@ class ProblemTDEM_b(BaseTDEMProblem): u = self.fields(m) p = self.Gvec(m, v, u) y = self.solveAh(m, p) - return self.survey.dpred(m, u=y) + Jv = self.survey.projectFieldsDeriv(u, v=y) + return mkvc(Jv) def Jtvec(self, m, v, u=None): if u is None: u = self.fields(m) - p = self.survey.projectFieldsAdjoint(v) + + if not isinstance(v, self.dataPair): + v = self.dataPair(self.survey, v) + + p = self.survey.projectFieldsDeriv(u, v=v, adjoint=True) y = self.solveAht(m, p) w = self.Gtvec(m, y, u) return w @@ -73,25 +78,25 @@ class ProblemTDEM_b(BaseTDEMProblem): """ if u is None: u = self.fields(m) - p = FieldsTDEM(self.mesh, 1, self.nT, 'b') + + p = FieldsTDEM(self.mesh, self.survey) + p[:, 'b', :] = 0.0 #np.zeros((self.mesh.nF, self.survey.nTx, self.prob.nT)) + p[:, 'e', 0] = 0.0 #np.zeros((self.mesh.nF, self.survey.nTx)) + # p = FieldsTDEM(self.mesh, 1, self.nT, 'b') curModel = self.mapping.transform(m) c = self.mesh.getEdgeInnerProductDeriv(curModel)*self.mapping.transformDeriv(m)*vec for i in range(self.nT): - ei = u.get_e(i) - pVal = np.empty_like(ei) - for j in range(ei.shape[1]): - pVal[:,j] = -ei[:,j]*c - - p.set_e(pVal,i) - p.set_b(np.zeros((self.mesh.nF,1)), i) + for tx in self.survey.txList: + p[tx, 'e', i+1] = -u[tx,'e',i+1]*c return p def Gtvec(self, m, v, u=None): if u is None: u = self.fields(m) - tmp = np.zeros((self.mesh.nE,self.survey.nTx)) - for i in range(self.nT): - tmp += v.get_e(i)*u.get_e(i) + nTx, nE = self.survey.nTx, self.mesh.nE + tmp = np.zeros(nE if nTx == 1 else (nE,nTx)) + for i in range(1,self.nT+1): + tmp += v[:,'e',i]*u[:,'e',i] curModel = self.mapping.transform(m) p = -mkvc(self.mapping.transformDeriv(m).T*self.mesh.getEdgeInnerProductDeriv(curModel).T*tmp) @@ -99,15 +104,17 @@ class ProblemTDEM_b(BaseTDEMProblem): def solveAh(self, m, p): def AhRHS(tInd, u): - rhs = self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*p.get_e(tInd) + p.get_b(tInd) + rhs = self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*p[:,'e',tInd+1] + p[:,'b',tInd+1] if tInd == 0: return rhs dt = self.timeSteps[tInd] - return rhs + 1.0/dt*self.MfMui*u.get_b(tInd-1) + return rhs + 1.0/dt*self.MfMui*u[:,'b',tInd] def AhCalcFields(sol, solType, tInd): b = sol - e = self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui*b - self.MeSigmaI*p.get_e(tInd) + if self.survey.nTx == 1: + b = mkvc(b) + e = self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui*b - self.MeSigmaI*p[:,'e',tInd+1] return {'b':b, 'e':e} self.curModel = m @@ -116,15 +123,17 @@ class ProblemTDEM_b(BaseTDEMProblem): def solveAht(self, m, p): def AhtRHS(tInd, u): - rhs = self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*p.get_e(tInd) + p.get_b(tInd) + rhs = self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*p[:,'e',tInd] + p[:,'b',tInd] if tInd == self.nT-1: return rhs dt = self.timeSteps[tInd+1] - return rhs + 1.0/dt*self.MfMui*u.get_b(tInd+1) + return rhs + 1.0/dt*self.MfMui*u[:,'b',tInd+1] def AhtCalcFields(sol, solType, tInd): b = sol - e = self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui*b - self.MeSigmaI*p.get_e(tInd) + if self.survey.nTx == 1: + b = mkvc(b) + e = self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui*b - self.MeSigmaI*p[:,'e',tInd] return {'b':b, 'e':e} self.curModel = m @@ -169,18 +178,14 @@ class ProblemTDEM_b(BaseTDEMProblem): """ 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) - f = FieldsTDEM(self.mesh, 1, self.nT, 'b') - f.set_b(b, 0) - f.set_e(e, 0) - for i in range(1,self.nT): - dt = self.timeSteps[i] - 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) + f = FieldsTDEM(self.mesh, self.survey) + for i in range(1,self.nT+1): + dt = self.timeSteps[i-1] + b = 1.0/dt*self.MfMui*vec[:,'b',i] + self.MfMui*self.mesh.edgeCurl*vec[:,'e',i] + if i > 1: + b = b - 1.0/dt*self.MfMui*vec[:,'b',i-1] + f[:,'b',i] = b + f[:,'e',i] = self.mesh.edgeCurl.T*self.MfMui*vec[:,'b',i] - self.MeSigma*vec[:,'e',i] return f def AhtVec(self, m, vec): @@ -217,17 +222,13 @@ class ProblemTDEM_b(BaseTDEMProblem): \\right] \\\\ """ 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) - 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) - N = self.nT - 1 - b = 1.0/self.timeSteps[N]*self.MfMui*vec.get_b(N) + self.MfMui*self.mesh.edgeCurl*vec.get_e(N) - e = self.mesh.edgeCurl.T*self.MfMui*vec.get_b(N) - self.MeSigma*vec.get_e(N) - f.set_b(b, N) - f.set_e(e, N) + f = FieldsTDEM(self.mesh, self.survey) + for i in range(1,self.nT+1): + b = 1.0/self.timeSteps[i-1]*self.MfMui*vec[:,'b',i] + self.MfMui*self.mesh.edgeCurl*vec[:,'e',i] + if i < self.nT: + b = b - 1.0/self.timeSteps[i]*self.MfMui*vec[:,'b',i+1] + f[:,'b', i] = b + f[:,'e', i] = self.mesh.edgeCurl.T*self.MfMui*vec[:,'b',i] - self.MeSigma*vec[:,'e',i] return f diff --git a/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py b/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py index 1eccc566..a7b8ea0f 100644 --- a/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py +++ b/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py @@ -21,14 +21,11 @@ class TDEM_bDerivTests(unittest.TestCase): mapping = Maps.ComboMap(mesh, [Maps.ExpMap, Maps.Vertical1DMap, activeMap]) + rxOffset = 40. + rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20), 'bz') + tx = EM.TDEM.TxTDEM(np.array([0., 0., 0.]), 'VMD_MVP', [rx]) - opts = {'txLoc':0., - 'txType': 'VMD_MVP', - 'rxLoc':np.r_[40., 0., 0.], - 'rxType':'bz', - 'timeCh':np.logspace(-4,-2,20), - } - self.dat = EM.TDEM.SurveyTDEM1D(**opts) + survey = EM.TDEM.SurveyTDEM([tx]) self.prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping) self.prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)] @@ -37,265 +34,270 @@ class TDEM_bDerivTests(unittest.TestCase): self.sigma[mesh.vectorCCz<0] = 1e-1 self.sigma = np.log(self.sigma[active]) - self.prb.pair(self.dat) + self.prb.pair(survey) self.mesh = mesh - def test_AhVec(self): - """ - Test that fields and AhVec produce consistent results - """ + # def test_AhVec(self): + # """ + # Test that fields and AhVec produce consistent results + # """ - prb = self.prb - sigma = self.sigma + # prb = self.prb + # sigma = self.sigma - u = prb.fields(sigma) - Ahu = prb.AhVec(sigma, u) + # u = prb.fields(sigma) + # Ahu = prb.AhVec(sigma, u) - V1 = Ahu.get_b(0) - V2 = 1./prb.timeSteps[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[:,'b',1] + # V2 = 1./prb.timeSteps[0]*prb.MfMui*u[:,'b',0] + # self.assertLess(np.linalg.norm(V1-V2)/np.linalg.norm(V2), 1.e-6) - V1 = Ahu.get_e(0) - self.assertTrue(np.linalg.norm(V1) < 1.e-6) + # V1 = Ahu[:,'e',1] + # self.assertLess(np.linalg.norm(V1), 1.e-6) - for i in range(1,u.nT): + # for i in range(2,prb.nT): - dt = prb.timeSteps[i] + # dt = prb.timeSteps[i] - V1 = Ahu.get_b(i) - V2 = 1/dt*prb.MfMui*u.get_b(i-1) - self.assertTrue(np.linalg.norm(V1)/np.linalg.norm(V2) < 1.e-6) + # V1 = Ahu[:,'b',i] + # V2 = 1.0/dt*prb.MfMui*u[:,'b', i-1] + # # print np.linalg.norm(V1), np.linalg.norm(V2) + # self.assertLess(np.linalg.norm(V1)/np.linalg.norm(V2), 1.e-6) - V1 = Ahu.get_e(i) - V2 = prb.MeSigma*u.get_e(i) - self.assertTrue(np.linalg.norm(V1)/np.linalg.norm(V2) < 1.e-6) + # V1 = Ahu[:,'e',i] + # V2 = prb.MeSigma*u[:,'e',i] + # # print np.linalg.norm(V1), np.linalg.norm(V2) + # self.assertLess(np.linalg.norm(V1)/np.linalg.norm(V2), 1.e-6) - def test_AhVecVSMat_OneTS(self): + # def test_AhVecVSMat_OneTS(self): - prb = self.prb - prb.timeSteps = [1e-05] - sigma = self.sigma - prb.curModel = sigma + # prb = self.prb + # prb.timeSteps = [1e-05] + # sigma = self.sigma + # prb.curModel = sigma - dt = prb.timeSteps[0] - a11 = 1/dt*prb.MfMui*sp.eye(prb.mesh.nF) - a12 = prb.MfMui*prb.mesh.edgeCurl - a21 = prb.mesh.edgeCurl.T*prb.MfMui - a22 = -prb.MeSigma - A = sp.bmat([[a11,a12],[a21,a22]]) + # dt = prb.timeSteps[0] + # a11 = 1/dt*prb.MfMui*sp.eye(prb.mesh.nF) + # a12 = prb.MfMui*prb.mesh.edgeCurl + # a21 = prb.mesh.edgeCurl.T*prb.MfMui + # a22 = -prb.MeSigma + # A = sp.bmat([[a11,a12],[a21,a22]]) - f = prb.fields(sigma) - u1 = A*f.fieldVec() - u2 = prb.AhVec(sigma,f).fieldVec() + # f = prb.fields(sigma) + # u1 = A*f.tovec() + # u2 = prb.AhVec(sigma,f).tovec() - self.assertTrue(np.linalg.norm(u1-u2)/np.linalg.norm(u1)<1e-12) + # self.assertTrue(np.linalg.norm(u1-u2)/np.linalg.norm(u1)<1e-12) - def test_solveAhVSMat_OneTS(self): - prb = self.prb + # def test_solveAhVSMat_OneTS(self): + # prb = self.prb - prb.timeSteps = [1e-05] + # prb.timeSteps = [1e-05] - sigma = self.sigma - prb.curModel = sigma + # sigma = self.sigma + # prb.curModel = sigma - dt = prb.timeSteps[0] - a11 = 1/dt*prb.MfMui*sp.eye(prb.mesh.nF) - a12 = prb.MfMui*prb.mesh.edgeCurl - a21 = prb.mesh.edgeCurl.T*prb.MfMui - a22 = -prb.MeSigma - A = sp.bmat([[a11,a12],[a21,a22]]) + # dt = prb.timeSteps[0] + # a11 = 1.0/dt*prb.MfMui*sp.eye(prb.mesh.nF) + # a12 = prb.MfMui*prb.mesh.edgeCurl + # a21 = prb.mesh.edgeCurl.T*prb.MfMui + # a22 = -prb.MeSigma + # A = sp.bmat([[a11,a12],[a21,a22]]) - f = prb.fields(sigma) - f.set_b(np.zeros((prb.mesh.nF,1)),0) - f.set_e(np.random.rand(prb.mesh.nE,1),0) + # f = prb.fields(sigma) + # f[:,:,0] = {'e':0,'b':0} + # f[:,'b',1] = 0 + # f[:,'e',1] = np.random.rand(prb.mesh.nE,1) - u1 = prb.solveAh(sigma,f).fieldVec().flatten() - u2 = sp.linalg.spsolve(A.tocsr(),f.fieldVec()) + # self.assertTrue(np.all(np.r_[f[:,'b',1],f[:,'e',1]] == f.tovec())) - self.assertTrue(np.linalg.norm(u1-u2)<1e-8) + # u1 = prb.solveAh(sigma,f).tovec().flatten() + # u2 = sp.linalg.spsolve(A.tocsr(),f.tovec()) - def test_solveAhVsAhVec(self): + # self.assertLess(np.linalg.norm(u1-u2),1e-8) - prb = self.prb - mesh = self.prb.mesh - sigma = self.sigma - self.prb.curModel = sigma + # def test_solveAhVsAhVec(self): - f = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nT, 'b') - for i in range(f.nT): - f.set_b(np.zeros((mesh.nF, 1)), i) - f.set_e(np.random.rand(mesh.nE, 1), i) + # prb = self.prb + # mesh = self.prb.mesh + # sigma = self.sigma + # self.prb.curModel = sigma - Ahf = prb.AhVec(sigma, f) - f_test = prb.solveAh(sigma, Ahf) + # f = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) + # f[:,'b',:] = 0.0 + # for i in range(prb.nT): + # f[:,'e', i] = np.random.rand(mesh.nE, 1) - u1 = f.fieldVec() - u2 = f_test.fieldVec() - self.assertTrue(np.linalg.norm(u1-u2)<1e-8) + # Ahf = prb.AhVec(sigma, f) + # f_test = prb.solveAh(sigma, Ahf) - def test_DerivG(self): - """ - Test the derivative of c with respect to sigma - """ + # u1 = f.tovec() + # u2 = f_test.tovec() + # self.assertTrue(np.linalg.norm(u1-u2)<1e-8) - # Random model and perturbation - sigma = np.random.rand(self.prb.mapping.nP) + # def test_DerivG(self): + # """ + # Test the derivative of c with respect to sigma + # """ - f = self.prb.fields(sigma) - dm = 1000*np.random.rand(self.prb.mapping.nP) - h = 0.01 + # # Random model and perturbation + # sigma = np.random.rand(self.prb.mapping.nP) - derChk = lambda m: [self.prb.AhVec(m, f).fieldVec(), lambda mx: self.prb.Gvec(sigma, mx, u=f).fieldVec()] - print '\ntest_DerivG' - passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20) - self.assertTrue(passed) + # f = self.prb.fields(sigma) + # dm = 1000*np.random.rand(self.prb.mapping.nP) + # h = 0.01 - def test_Deriv_dUdM(self): + # derChk = lambda m: [self.prb.AhVec(m, f).tovec(), lambda mx: self.prb.Gvec(sigma, mx, u=f).tovec()] + # print '\ntest_DerivG' + # passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20) + # self.assertTrue(passed) - prb = self.prb - prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)] - mesh = self.mesh - sigma = self.sigma + # def test_Deriv_dUdM(self): - dm = 10*np.random.rand(prb.mapping.nP) - f = prb.fields(sigma) + # prb = self.prb + # prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)] + # mesh = self.mesh + # sigma = self.sigma - derChk = lambda m: [self.prb.fields(m).fieldVec(), lambda mx: -prb.solveAh(sigma, prb.Gvec(sigma, mx, u=f)).fieldVec()] - print '\n' - print 'test_Deriv_dUdM' - passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20) - self.assertTrue(passed) + # dm = 10*np.random.rand(prb.mapping.nP) + # f = prb.fields(sigma) - def test_Deriv_J(self): + # derChk = lambda m: [self.prb.fields(m).tovec(), lambda mx: -prb.solveAh(sigma, prb.Gvec(sigma, mx, u=f)).tovec()] + # print '\n' + # print 'test_Deriv_dUdM' + # passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20) + # self.assertTrue(passed) - prb = self.prb - prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)] - mesh = self.mesh - sigma = self.sigma + # def test_Deriv_J(self): - # d_sig = 0.8*sigma #np.random.rand(mesh.nCz) - d_sig = 10*np.random.rand(prb.mapping.nP) + # prb = self.prb + # prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)] + # mesh = self.mesh + # sigma = self.sigma + + # # d_sig = 0.8*sigma #np.random.rand(mesh.nCz) + # d_sig = 10*np.random.rand(prb.mapping.nP) - derChk = lambda m: [prb.survey.dpred(m), lambda mx: -prb.Jvec(sigma, mx)] - print '\n' - print 'test_Deriv_J' - passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=4, eps=1e-20) - self.assertTrue(passed) + # derChk = lambda m: [prb.survey.dpred(m), lambda mx: -prb.Jvec(sigma, mx)] + # print '\n' + # print 'test_Deriv_J' + # passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=4, eps=1e-20) + # self.assertTrue(passed) - def test_projectAdjoint(self): - prb = self.prb - dat = self.dat - mesh = self.mesh + # def test_projectAdjoint(self): + # prb = self.prb + # survey = prb.survey + # mesh = self.mesh - # Generate random fields and data - f = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nT, 'b') - for i in range(f.nT): - f.set_b(np.random.rand(mesh.nF, 1), i) - f.set_e(np.random.rand(mesh.nE, 1), i) - d = np.random.rand(dat.prob.nT, dat.nTx) + # # Generate random fields and data + # f = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) + # for i in range(prb.nT): + # f[:,'b',i] = np.random.rand(mesh.nF, 1) + # f[:,'e',i] = np.random.rand(mesh.nE, 1) + # d_vec = np.random.rand(survey.nD, survey.nTx).flatten() + # d = Survey.Data(survey,v=d_vec) - # Check that d.T*Q*f = f.T*Q.T*d - V1 = d.T.dot(dat.projectFields(f)) - V2 = f.fieldVec().dot(dat.projectFieldsAdjoint(d).fieldVec()) + # # Check that d.T*Q*f = f.T*Q.T*d + # V1 = d_vec.dot(survey.projectFieldsDeriv(None, v=f).tovec()) + # V2 = f.tovec().dot(survey.projectFieldsDeriv(None, v=d, adjoint=True).tovec()) - self.assertLess((V1-V2)/np.abs(V1), 1e-6) + # self.assertLess((V1-V2)/np.abs(V1), 1e-6) - def test_adjointAhVsAht(self): - prb = self.prb - mesh = self.mesh - sigma = self.sigma + # def test_adjointAhVsAht(self): + # prb = self.prb + # mesh = self.mesh + # sigma = self.sigma - f1 = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nT, 'b') - for i in range(f1.nT): - f1.set_b(np.random.rand(mesh.nF, 1), i) - f1.set_e(np.random.rand(mesh.nE, 1), i) + # f1 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) + # for i in range(1,prb.nT+1): + # f1[:,'b',i] = np.random.rand(mesh.nF, 1) + # f1[:,'e',i] = np.random.rand(mesh.nE, 1) - f2 = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nT, 'b') - for i in range(f2.nT): - f2.set_b(np.random.rand(mesh.nF, 1), i) - f2.set_e(np.random.rand(mesh.nE, 1), i) + # f2 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) + # for i in range(1,prb.nT+1): + # f2[:,'b',i] = np.random.rand(mesh.nF, 1) + # f2[:,'e',i] = np.random.rand(mesh.nE, 1) - V1 = f2.fieldVec().dot(prb.AhVec(sigma, f1).fieldVec()) - V2 = f1.fieldVec().dot(prb.AhtVec(sigma, f2).fieldVec()) - self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6) + # V1 = f2.tovec().dot(prb.AhVec(sigma, f1).tovec()) + # V2 = f1.tovec().dot(prb.AhtVec(sigma, f2).tovec()) + # self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6) - def test_solveAhtVsAhtVec(self): - prb = self.prb - mesh = self.mesh - sigma = np.random.rand(prb.mapping.nP) + # def test_solveAhtVsAhtVec(self): + # prb = self.prb + # mesh = self.mesh + # sigma = np.random.rand(prb.mapping.nP) - f1 = EM.TDEM.FieldsTDEM(mesh, 1, prb.nT, 'b') - for i in range(prb.nT): - f1.set_b(np.random.rand(mesh.nF, 1), i) - f1.set_e(np.random.rand(mesh.nE, 1), i) + # f1 = EM.TDEM.FieldsTDEM(mesh, 1, prb.nT, 'b') + # for i in range(prb.nT): + # f1.set_b(np.random.rand(mesh.nF, 1), i) + # f1.set_e(np.random.rand(mesh.nE, 1), i) - f2 = prb.solveAht(sigma, f1) - f3 = prb.AhtVec(sigma, f2) + # f2 = prb.solveAht(sigma, f1) + # f3 = prb.AhtVec(sigma, f2) - if plotIt: - import matplotlib.pyplot as plt - plt.plot(f3.fieldVec()) - plt.plot(f1.fieldVec()) - plt.show() - 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) + # if plotIt: + # import matplotlib.pyplot as plt + # plt.plot(f3.tovec()) + # plt.plot(f1.tovec()) + # plt.show() + # V1 = np.linalg.norm(f3.tovec()-f1.tovec()) + # V2 = np.linalg.norm(f1.tovec()) + # print V1, V2 + # print 'I am gunna fail this one: boo. :(' + # self.assertLess(V1/V2, 1e-6) def test_adjointsolveAhVssolveAht(self): prb = self.prb mesh = self.mesh sigma = self.sigma - f1 = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nT, 'b') - for i in range(f1.nT): - f1.set_b(np.random.rand(mesh.nF, 1), i) - f1.set_e(np.random.rand(mesh.nE, 1), i) + f1 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) + for i in range(1,prb.nT+1): + f1[:,'b',i] = np.random.rand(mesh.nF, 1) + f1[:,'e',i] = np.random.rand(mesh.nE, 1) - f2 = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nT, 'b') - for i in range(f2.nT): - f2.set_b(np.random.rand(mesh.nF, 1), i) - f2.set_e(np.random.rand(mesh.nE, 1), i) + f2 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) + for i in range(1,prb.nT+1): + f2[:,'b',i] = np.random.rand(mesh.nF, 1) + f2[:,'e',i] = np.random.rand(mesh.nE, 1) - V1 = f2.fieldVec().dot(prb.solveAh(sigma, f1).fieldVec()) - V2 = f1.fieldVec().dot(prb.solveAht(sigma, f2).fieldVec()) + V1 = f2.tovec().dot(prb.solveAh(sigma, f1).tovec()) + V2 = f1.tovec().dot(prb.solveAht(sigma, f2).tovec()) self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6) - def test_adjointGvecVsGtvec(self): - mesh = self.mesh - prb = self.prb + # def test_adjointGvecVsGtvec(self): + # mesh = self.mesh + # prb = self.prb - m = np.random.rand(prb.mapping.nP) - sigma = np.random.rand(prb.mapping.nP) + # m = np.random.rand(prb.mapping.nP) + # sigma = np.random.rand(prb.mapping.nP) - u = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nT, 'b') - for i in range(u.nT): - u.set_b(np.random.rand(mesh.nF, 1), i) - u.set_e(np.random.rand(mesh.nE, 1), i) + # u = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) + # for i in range(prb.nT): + # u[:,'b',i] = np.random.rand(mesh.nF, 1) + # u[:,'e',i] = np.random.rand(mesh.nE, 1) - v = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nT, 'b') - for i in range(v.nT): - v.set_b(np.random.rand(mesh.nF, 1), i) - v.set_e(np.random.rand(mesh.nE, 1), i) + # v = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) + # for i in range(prb.nT): + # v[:,'b',i] = np.random.rand(mesh.nF, 1) + # v[:,'e',i] = np.random.rand(mesh.nE, 1) - V1 = m.dot(prb.Gtvec(sigma, v, u)) - V2 = v.fieldVec().dot(prb.Gvec(sigma, m, u).fieldVec()) - self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6) + # V1 = m.dot(prb.Gtvec(sigma, v, u)) + # V2 = v.tovec().dot(prb.Gvec(sigma, m, u).tovec()) + # self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6) - def test_adjointJvecVsJtvec(self): - mesh = self.mesh - prb = self.prb - sigma = self.sigma + # def test_adjointJvecVsJtvec(self): + # mesh = self.mesh + # prb = self.prb + # sigma = self.sigma - m = np.random.rand(prb.mapping.nP) - d = np.random.rand(prb.nT) + # m = np.random.rand(prb.mapping.nP) + # d = np.random.rand(prb.survey.nD) - V1 = d.dot(prb.Jvec(sigma, m)) - V2 = m.dot(prb.Jtvec(sigma, d)) - self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6) + # V1 = d.dot(prb.Jvec(sigma, m)) + # V2 = m.dot(prb.Jtvec(sigma, d)) + # self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6) diff --git a/simpegEM/Tests/test_TDEM_forward_Analytic.py b/simpegEM/Tests/test_TDEM_forward_Analytic.py index ed308906..84912b36 100644 --- a/simpegEM/Tests/test_TDEM_forward_Analytic.py +++ b/simpegEM/Tests/test_TDEM_forward_Analytic.py @@ -28,11 +28,6 @@ def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=[1e-5, survey = EM.TDEM.SurveyTDEM([tx]) prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping) prb.Solver = Utils.SolverUtils.DSolverWrap(sp.linalg.splu, factorize=True) - # try: - # from mumpsSCI import MumpsSolver - # prb.Solver = MumpsSolver - # except ImportError, e: - # pass prb.timeSteps = [(1e-06, 40), (5e-06, 40), (1e-05, 40), (5e-05, 40), (0.0001, 40), (0.0005, 40)] From d1f23081d648d72714099b9c683c855d25171c81 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sun, 27 Apr 2014 21:49:00 -0700 Subject: [PATCH 4/8] comment G vec and Gt vec --- simpegEM/TDEM/TDEM_b.py | 30 ++++++++++++++++++++++-------- 1 file changed, 22 insertions(+), 8 deletions(-) diff --git a/simpegEM/TDEM/TDEM_b.py b/simpegEM/TDEM/TDEM_b.py index 787c1442..2999f01e 100644 --- a/simpegEM/TDEM/TDEM_b.py +++ b/simpegEM/TDEM/TDEM_b.py @@ -74,29 +74,43 @@ class ProblemTDEM_b(BaseTDEMProblem): :rtype: simpegEM.TDEM.FieldsTDEM :return: f - Multiply G by a vector where + Multiply G by a vector """ if u is None: u = self.fields(m) + # Note: Fields has shape (nF/E, nTx, nT+1) + # However, p will only really fill (:,:,1:nT+1) + # meaning the 'initial fields' are zero (:,:,0) p = FieldsTDEM(self.mesh, self.survey) - p[:, 'b', :] = 0.0 #np.zeros((self.mesh.nF, self.survey.nTx, self.prob.nT)) - p[:, 'e', 0] = 0.0 #np.zeros((self.mesh.nF, self.survey.nTx)) - # p = FieldsTDEM(self.mesh, 1, self.nT, 'b') + p[:, 'b', :] = 0.0 # b at all times is zero. + p[:, 'e', 0] = 0.0 # fake initial fields curModel = self.mapping.transform(m) c = self.mesh.getEdgeInnerProductDeriv(curModel)*self.mapping.transformDeriv(m)*vec - for i in range(self.nT): + 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) for tx in self.survey.txList: - p[tx, 'e', i+1] = -u[tx,'e',i+1]*c + p[tx, 'e', i] = -u[tx,'e',i]*c # - diag(e) * MsigDeriv * v return p - def Gtvec(self, m, v, u=None): + def Gtvec(self, m, vec, u=None): + """ + :param numpy.array m: Conductivity model + :param numpy.array vec: vector (like a fields) + :param simpegEM.TDEM.FieldsTDEM u: Fields resulting from m + :rtype: np.ndarray (like a model) + :return: p + + Multiply G.T by a vector + """ if u is None: u = self.fields(m) nTx, nE = self.survey.nTx, self.mesh.nE tmp = np.zeros(nE if nTx == 1 else (nE,nTx)) + # Here we can do internal multiplications of Gt*v and then multiply by MsigDeriv.T in one go. for i in range(1,self.nT+1): - tmp += v[:,'e',i]*u[:,'e',i] + tmp += vec[:,'e',i]*u[:,'e',i] curModel = self.mapping.transform(m) p = -mkvc(self.mapping.transformDeriv(m).T*self.mesh.getEdgeInnerProductDeriv(curModel).T*tmp) From df7676f14c2dc8023d999ec8f6088e56f25e2979 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sun, 27 Apr 2014 22:17:55 -0700 Subject: [PATCH 5/8] working jvec --- docs/api_TDEM_derivation.rst | 24 +++------- simpegEM/TDEM/TDEM_b.py | 15 +++--- simpegEM/Tests/test_TDEM_b_DerivAdjoint.py | 54 +++++++++++----------- 3 files changed, 41 insertions(+), 52 deletions(-) diff --git a/docs/api_TDEM_derivation.rst b/docs/api_TDEM_derivation.rst index 364070e5..0d8bfcea 100644 --- a/docs/api_TDEM_derivation.rst +++ b/docs/api_TDEM_derivation.rst @@ -255,25 +255,8 @@ Multiplying **J** onto a vector can be broken into three steps \vec{p}_e^{(n)} = - \diag{\e^{(n)}} \Ace \diag{V} m \end{align} -First time step -.. math:: - - \begin{align} - \frac{1}{\delta t} \MfMui \vec{y}_{b}^{(1)} + \MfMui \dcurl \vec{y}_{e}^{(1)} = \vec{p}_b^{(1)} \\ - \dcurl^\top \MfMui \vec{y}_b^{(1)} - \MeSig \vec{y}_e^{(1)} = \vec{p}_e^{(1)} - \end{align} - - -.. math:: - - \begin{align} - \left( \MfMui \dcurl \MeSig^{-1} \dcurl^\top \MfMui + \frac{1}{\delta t} \MfMui \right) \vec{y}_{b}^{(1)} = \MfMui \dcurl \MeSig^{-1} \vec{p}_e^{(1)} + \vec{p}_b^{(1)} \\ - \vec{y}_e^{(1)} = \MeSig^{-1} \dcurl^\top \MfMui \vec{y}_b^{(1)} - \MeSig^{-1} \vec{p}_e^{(1)} - \end{align} - - -Remaining time steps: +For all time steps: .. math:: @@ -295,6 +278,11 @@ and \vec{y}_e^{(t+1)} = \MeSig^{-1} \dcurl^\top \MfMui \vec{y}_b^{(t+1)} - \MeSig^{-1} \vec{p}_e^{(t+1)} \end{align} +.. note:: + + For the first time step, \\\(t=0\\\), the term: \\\(\\frac{1}{\\delta t} \\MfMui \\vec{y}_b^{(0)}\\\) is zero. + + Implementing **J** transpose times a vector diff --git a/simpegEM/TDEM/TDEM_b.py b/simpegEM/TDEM/TDEM_b.py index 2999f01e..e60b9ec9 100644 --- a/simpegEM/TDEM/TDEM_b.py +++ b/simpegEM/TDEM/TDEM_b.py @@ -52,7 +52,7 @@ class ProblemTDEM_b(BaseTDEMProblem): p = self.Gvec(m, v, u) y = self.solveAh(m, p) Jv = self.survey.projectFieldsDeriv(u, v=y) - return mkvc(Jv) + return - mkvc(Jv) def Jtvec(self, m, v, u=None): if u is None: @@ -117,19 +117,20 @@ class ProblemTDEM_b(BaseTDEMProblem): return p def solveAh(self, m, p): - def AhRHS(tInd, u): + + def AhRHS(tInd, y): rhs = self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*p[:,'e',tInd+1] + p[:,'b',tInd+1] if tInd == 0: return rhs dt = self.timeSteps[tInd] - return rhs + 1.0/dt*self.MfMui*u[:,'b',tInd] + return rhs + 1.0/dt*self.MfMui*y[:,'b',tInd] def AhCalcFields(sol, solType, tInd): - b = sol + y_b = sol if self.survey.nTx == 1: - b = mkvc(b) - e = self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui*b - self.MeSigmaI*p[:,'e',tInd+1] - return {'b':b, 'e':e} + y_b = mkvc(y_b) + y_e = self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui*y_b - self.MeSigmaI*p[:,'e',tInd+1] + return {'b':y_b, 'e':y_e} self.curModel = m return self.forward(m, AhRHS, AhCalcFields) diff --git a/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py b/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py index a7b8ea0f..9f505ce3 100644 --- a/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py +++ b/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py @@ -168,22 +168,22 @@ class TDEM_bDerivTests(unittest.TestCase): # passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20) # self.assertTrue(passed) - # def test_Deriv_J(self): + def test_Deriv_J(self): - # prb = self.prb - # prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)] - # mesh = self.mesh - # sigma = self.sigma + prb = self.prb + prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)] + mesh = self.mesh + sigma = self.sigma - # # d_sig = 0.8*sigma #np.random.rand(mesh.nCz) - # d_sig = 10*np.random.rand(prb.mapping.nP) + # d_sig = 0.8*sigma #np.random.rand(mesh.nCz) + d_sig = 10*np.random.rand(prb.mapping.nP) - # derChk = lambda m: [prb.survey.dpred(m), lambda mx: -prb.Jvec(sigma, mx)] - # print '\n' - # print 'test_Deriv_J' - # passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=4, eps=1e-20) - # self.assertTrue(passed) + derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(sigma, mx)] + print '\n' + print 'test_Deriv_J' + passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=4, eps=1e-20) + self.assertTrue(passed) # def test_projectAdjoint(self): # prb = self.prb @@ -247,24 +247,24 @@ class TDEM_bDerivTests(unittest.TestCase): # print 'I am gunna fail this one: boo. :(' # self.assertLess(V1/V2, 1e-6) - def test_adjointsolveAhVssolveAht(self): - prb = self.prb - mesh = self.mesh - sigma = self.sigma + # def test_adjointsolveAhVssolveAht(self): + # prb = self.prb + # mesh = self.mesh + # sigma = self.sigma - f1 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) - for i in range(1,prb.nT+1): - f1[:,'b',i] = np.random.rand(mesh.nF, 1) - f1[:,'e',i] = np.random.rand(mesh.nE, 1) + # f1 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) + # for i in range(1,prb.nT+1): + # f1[:,'b',i] = np.random.rand(mesh.nF, 1) + # f1[:,'e',i] = np.random.rand(mesh.nE, 1) - f2 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) - for i in range(1,prb.nT+1): - f2[:,'b',i] = np.random.rand(mesh.nF, 1) - f2[:,'e',i] = np.random.rand(mesh.nE, 1) + # f2 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) + # for i in range(1,prb.nT+1): + # f2[:,'b',i] = np.random.rand(mesh.nF, 1) + # f2[:,'e',i] = np.random.rand(mesh.nE, 1) - V1 = f2.tovec().dot(prb.solveAh(sigma, f1).tovec()) - V2 = f1.tovec().dot(prb.solveAht(sigma, f2).tovec()) - self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6) + # V1 = f2.tovec().dot(prb.solveAh(sigma, f1).tovec()) + # V2 = f1.tovec().dot(prb.solveAht(sigma, f2).tovec()) + # self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6) # def test_adjointGvecVsGtvec(self): # mesh = self.mesh From 3514ad56d314406564a723d0e51ab475b694bd58 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sun, 27 Apr 2014 23:24:49 -0700 Subject: [PATCH 6/8] Adjoint solve not yet working. --- docs/api_TDEM_derivation.rst | 44 ++++------ simpegEM/TDEM/BaseTDEM.py | 8 +- simpegEM/TDEM/SurveyTDEM.py | 3 +- simpegEM/TDEM/TDEM_b.py | 23 +++-- simpegEM/Tests/test_TDEM_b_DerivAdjoint.py | 98 +++++++++++----------- 5 files changed, 85 insertions(+), 91 deletions(-) diff --git a/docs/api_TDEM_derivation.rst b/docs/api_TDEM_derivation.rst index 0d8bfcea..af3fc2fc 100644 --- a/docs/api_TDEM_derivation.rst +++ b/docs/api_TDEM_derivation.rst @@ -307,38 +307,21 @@ Multiplying \\(\\mathbf{J}^\\top\\) onto a vector can be broken into three steps \end{array} \right] -For the last time-step \\(t=N\\): +For the all time-steps (going backwards in time): + .. math:: - \begin{align} - \frac{1}{\delta t} \MfMui \vec{y}_{b}^{(N)} + \MfMui \dcurl \vec{y}_{e}^{(N)} = \vec{p}_b^{(N)} \\ - \dcurl^\top \MfMui \vec{y}_b^{(N)} - \MeSig \vec{y}_e^{(N)} = \vec{p}_e^{(N)} - \end{align} + A \vec{y}^{(t)} + B \vec{y}^{(t+1)} = \vec{p}^{(t)} .. math:: \begin{align} - \left( \MfMui \dcurl \MeSig^{-1} \dcurl^\top \MfMui + \frac{1}{\delta t} \MfMui \right) \vec{y}_{b}^{(N)} = \MfMui \dcurl \MeSig^{-1} \vec{p}_e^{(N)} + \vec{p}_b^{(N)} \\ - \vec{y}_e^{(N)} = \MeSig^{-1} \dcurl^\top \MfMui \vec{y}_b^{(N)} - \MeSig^{-1} \vec{p}_e^{(N)} - \end{align} - -For the rest of the time-steps (going backwards in time) - - -.. math:: - - A \vec{y}^{(t-1)} + B \vec{y}^{(t)} = \vec{p}^{(t-1)} - - -.. math:: - - \begin{align} - \frac{1}{\delta t} \MfMui\vec{y}_{b}^{(t-1)} + \MfMui\dcurl \vec{y}_{e}^{(t-1)} - - \frac{1}{\delta t} \MfMui \vec{y}_{b}^{(t)} - = \vec{p}_b^{(t-1)} \\ - \dcurl^\top \MfMui \vec{y}_b^{(t-1)} - \MeSig \vec{y}_e^{(t-1)} = \vec{p}_e^{(t-1)} + \frac{1}{\delta t} \MfMui\vec{y}_{b}^{(t)} + \MfMui\dcurl \vec{y}_{e}^{(t)} + - \frac{1}{\delta t} \MfMui \vec{y}_{b}^{(t+1)} + = \vec{p}_b^{(t)} \\ + \dcurl^\top \MfMui \vec{y}_b^{(t)} - \MeSig \vec{y}_e^{(t)} = \vec{p}_e^{(t)} \end{align} and @@ -346,8 +329,13 @@ and .. math:: \begin{align} - \left( \MfMui \dcurl \MeSig^{-1} \dcurl^\top \MfMui + \frac{1}{\delta t} \MfMui \right) \vec{y}_{b}^{(t-1)} = - \frac{1}{\delta t} \MfMui \vec{y}_b^{(t)} - + \MfMui \dcurl \MeSig^{-1} \vec{p}_e^{(t-1)} + \vec{p}_b^{(t-1)} \\ - \vec{y}_e^{(t-1)} = \MeSig^{-1} \dcurl^\top \MfMui \vec{y}_b^{(t-1)} - \MeSig^{-1} \vec{p}_e^{(t-1)} + \left( \MfMui \dcurl \MeSig^{-1} \dcurl^\top \MfMui + \frac{1}{\delta t} \MfMui \right) \vec{y}_{b}^{(t)} = + \frac{1}{\delta t} \MfMui \vec{y}_b^{(t+1)} + + \MfMui \dcurl \MeSig^{-1} \vec{p}_e^{(t)} + \vec{p}_b^{(t)} \\ + \vec{y}_e^{(t)} = \MeSig^{-1} \dcurl^\top \MfMui \vec{y}_b^{(t)} - \MeSig^{-1} \vec{p}_e^{(t)} \end{align} + + +.. note:: + + For the last time step, \\\(t=N\\\), the term: \\\(\\frac{1}{\\delta t} \\MfMui \\vec{y}_b^{(N+1)}\\\) is zero. diff --git a/simpegEM/TDEM/BaseTDEM.py b/simpegEM/TDEM/BaseTDEM.py index 12911847..cacf6eaa 100644 --- a/simpegEM/TDEM/BaseTDEM.py +++ b/simpegEM/TDEM/BaseTDEM.py @@ -57,8 +57,7 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem): return F def adjoint(self, m, RHS, CalcFields, F=None): - if F is None: - F = FieldsTDEM(self.mesh, self.survey) + F = F or FieldsTDEM(self.mesh, self.survey) dtFact = None for tInd, dt in reversed(list(enumerate(self.timeSteps))): @@ -66,13 +65,12 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem): dtFact = dt A = self.getA(tInd) # print 'Factoring... (dt = ' + str(dt) + ')' - Asolve = Solver(A, options=self.solverOpts) + 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[:,:,tInd] = newFields + F[:,:,tInd+1] = CalcFields(sol, self.solType, tInd) return F diff --git a/simpegEM/TDEM/SurveyTDEM.py b/simpegEM/TDEM/SurveyTDEM.py index 3bc5fde6..722846f5 100644 --- a/simpegEM/TDEM/SurveyTDEM.py +++ b/simpegEM/TDEM/SurveyTDEM.py @@ -39,8 +39,7 @@ class RxTDEM(Survey.BaseTimeRx): if not adjoint: return P * Utils.mkvc(v[tx, self.projField, :]) elif adjoint: - Ptv = P.T * v[tx, self] - return Ptv + return P.T * v[tx, self] class FieldsTDEM(Survey.TimeFields): diff --git a/simpegEM/TDEM/TDEM_b.py b/simpegEM/TDEM/TDEM_b.py index e60b9ec9..95304962 100644 --- a/simpegEM/TDEM/TDEM_b.py +++ b/simpegEM/TDEM/TDEM_b.py @@ -64,7 +64,7 @@ class ProblemTDEM_b(BaseTDEMProblem): p = self.survey.projectFieldsDeriv(u, v=v, adjoint=True) y = self.solveAht(m, p) w = self.Gtvec(m, y, u) - return w + return - w def Gvec(self, m, vec, u=None): """ @@ -138,18 +138,27 @@ class ProblemTDEM_b(BaseTDEMProblem): def solveAht(self, m, p): def AhtRHS(tInd, u): - rhs = self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*p[:,'e',tInd] + p[:,'b',tInd] + nTx, nF = self.survey.nTx, self.mesh.nF + rhs = np.zeros(nF if nTx == 1 else (nF, nTx)) + + if 'e' in p: + rhs += self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*p[:,'e',tInd] + if 'b' in p: + rhs += p[:,'b',tInd] + if tInd == self.nT-1: return rhs - dt = self.timeSteps[tInd+1] + dt = self.timeSteps[tInd] return rhs + 1.0/dt*self.MfMui*u[:,'b',tInd+1] def AhtCalcFields(sol, solType, tInd): - b = sol + y_b = sol if self.survey.nTx == 1: - b = mkvc(b) - e = self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui*b - self.MeSigmaI*p[:,'e',tInd] - return {'b':b, 'e':e} + y_b = mkvc(y_b) + y_e = self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui*y_b + if 'e' in p: + y_e += - self.MeSigmaI*p[:,'e',tInd] + return {'b':y_b, 'e':y_e} self.curModel = m return self.adjoint(m, AhtRHS, AhtCalcFields) diff --git a/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py b/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py index 9f505ce3..7b32f484 100644 --- a/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py +++ b/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py @@ -168,22 +168,22 @@ class TDEM_bDerivTests(unittest.TestCase): # passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20) # self.assertTrue(passed) - def test_Deriv_J(self): + # def test_Deriv_J(self): - prb = self.prb - prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)] - mesh = self.mesh - sigma = self.sigma + # prb = self.prb + # prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)] + # mesh = self.mesh + # sigma = self.sigma - # d_sig = 0.8*sigma #np.random.rand(mesh.nCz) - d_sig = 10*np.random.rand(prb.mapping.nP) + # # d_sig = 0.8*sigma #np.random.rand(mesh.nCz) + # d_sig = 10*np.random.rand(prb.mapping.nP) - derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(sigma, mx)] - print '\n' - print 'test_Deriv_J' - passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=4, eps=1e-20) - self.assertTrue(passed) + # derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(sigma, mx)] + # print '\n' + # print 'test_Deriv_J' + # passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=4, eps=1e-20) + # self.assertTrue(passed) # def test_projectAdjoint(self): # prb = self.prb @@ -223,48 +223,48 @@ class TDEM_bDerivTests(unittest.TestCase): # V2 = f1.tovec().dot(prb.AhtVec(sigma, f2).tovec()) # self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6) - # def test_solveAhtVsAhtVec(self): - # prb = self.prb - # mesh = self.mesh - # sigma = np.random.rand(prb.mapping.nP) + def test_solveAhtVsAhtVec(self): + prb = self.prb + mesh = self.mesh + sigma = np.random.rand(prb.mapping.nP) - # f1 = EM.TDEM.FieldsTDEM(mesh, 1, prb.nT, 'b') - # for i in range(prb.nT): - # f1.set_b(np.random.rand(mesh.nF, 1), i) - # f1.set_e(np.random.rand(mesh.nE, 1), i) + f1 = EM.TDEM.FieldsTDEM(mesh,prb.survey) + for i in range(prb.nT): + f1[:,'b',i] = np.random.rand(mesh.nF, 1) + f1[:,'e',i] = np.random.rand(mesh.nE, 1) - # f2 = prb.solveAht(sigma, f1) - # f3 = prb.AhtVec(sigma, f2) + f2 = prb.solveAht(sigma, f1) + f3 = prb.AhtVec(sigma, f2) - # if plotIt: - # import matplotlib.pyplot as plt - # plt.plot(f3.tovec()) - # plt.plot(f1.tovec()) - # plt.show() - # V1 = np.linalg.norm(f3.tovec()-f1.tovec()) - # V2 = np.linalg.norm(f1.tovec()) - # print V1, V2 - # print 'I am gunna fail this one: boo. :(' - # self.assertLess(V1/V2, 1e-6) + if plotIt: + import matplotlib.pyplot as plt + plt.plot(f3.tovec()) + plt.plot(f1.tovec()) + plt.show() + V1 = np.linalg.norm(f3.tovec()-f1.tovec()) + V2 = np.linalg.norm(f1.tovec()) + print V1, V2 + print 'I am gunna fail this one: boo. :(' + self.assertLess(V1/V2, 1e-6) - # def test_adjointsolveAhVssolveAht(self): - # prb = self.prb - # mesh = self.mesh - # sigma = self.sigma + def test_adjointsolveAhVssolveAht(self): + prb = self.prb + mesh = self.mesh + sigma = self.sigma - # f1 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) - # for i in range(1,prb.nT+1): - # f1[:,'b',i] = np.random.rand(mesh.nF, 1) - # f1[:,'e',i] = np.random.rand(mesh.nE, 1) + f1 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) + for i in range(1,prb.nT+1): + f1[:,'b',i] = np.random.rand(mesh.nF, 1) + f1[:,'e',i] = np.random.rand(mesh.nE, 1) - # f2 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) - # for i in range(1,prb.nT+1): - # f2[:,'b',i] = np.random.rand(mesh.nF, 1) - # f2[:,'e',i] = np.random.rand(mesh.nE, 1) + f2 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) + for i in range(1,prb.nT+1): + f2[:,'b',i] = np.random.rand(mesh.nF, 1) + f2[:,'e',i] = np.random.rand(mesh.nE, 1) - # V1 = f2.tovec().dot(prb.solveAh(sigma, f1).tovec()) - # V2 = f1.tovec().dot(prb.solveAht(sigma, f2).tovec()) - # self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6) + V1 = f2.tovec().dot(prb.solveAh(sigma, f1).tovec()) + V2 = f1.tovec().dot(prb.solveAht(sigma, f2).tovec()) + self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6) # def test_adjointGvecVsGtvec(self): # mesh = self.mesh @@ -274,12 +274,12 @@ class TDEM_bDerivTests(unittest.TestCase): # sigma = np.random.rand(prb.mapping.nP) # u = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) - # for i in range(prb.nT): + # for i in range(1,prb.nT+1): # u[:,'b',i] = np.random.rand(mesh.nF, 1) # u[:,'e',i] = np.random.rand(mesh.nE, 1) # v = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) - # for i in range(prb.nT): + # for i in range(1,prb.nT+1): # v[:,'b',i] = np.random.rand(mesh.nF, 1) # v[:,'e',i] = np.random.rand(mesh.nE, 1) From 4c15267c0ce9e58319fd25591b86851ab4bc14b3 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Mon, 28 Apr 2014 11:59:25 -0700 Subject: [PATCH 7/8] Adjoint test working. --- simpegEM/TDEM/TDEM_b.py | 18 +- simpegEM/Tests/test_TDEM_b_DerivAdjoint.py | 456 +++++++++++---------- 2 files changed, 243 insertions(+), 231 deletions(-) diff --git a/simpegEM/TDEM/TDEM_b.py b/simpegEM/TDEM/TDEM_b.py index 95304962..01c6a50b 100644 --- a/simpegEM/TDEM/TDEM_b.py +++ b/simpegEM/TDEM/TDEM_b.py @@ -137,19 +137,27 @@ class ProblemTDEM_b(BaseTDEMProblem): def solveAht(self, m, p): - def AhtRHS(tInd, u): + # fLoc 0 1 2 3 + # |-----|-----|-----| + # tInd 0 1 2 / / + # / __/ + # 2 (tInd=2 uses fields 3 and would use 4 but it doesn't exist) + # / __/ + # 1 (tInd=1 uses fields 2 and 3) + + def AhtRHS(tInd, y): nTx, nF = self.survey.nTx, self.mesh.nF rhs = np.zeros(nF if nTx == 1 else (nF, nTx)) if 'e' in p: - rhs += self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*p[:,'e',tInd] + rhs += self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*p[:,'e',tInd+1] if 'b' in p: - rhs += p[:,'b',tInd] + rhs += p[:,'b',tInd+1] if tInd == self.nT-1: return rhs - dt = self.timeSteps[tInd] - return rhs + 1.0/dt*self.MfMui*u[:,'b',tInd+1] + dt = self.timeSteps[tInd+1] + return rhs + 1.0/dt*self.MfMui*y[:,'b',tInd+2] def AhtCalcFields(sol, solType, tInd): y_b = sol diff --git a/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py b/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py index 7b32f484..841cef25 100644 --- a/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py +++ b/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py @@ -28,7 +28,9 @@ class TDEM_bDerivTests(unittest.TestCase): survey = EM.TDEM.SurveyTDEM([tx]) self.prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping) + # self.prb.timeSteps = [1e-5] self.prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)] + # self.prb.timeSteps = [(1e-05, 100)] self.sigma = np.ones(mesh.nCz)*1e-8 self.sigma[mesh.vectorCCz<0] = 1e-1 @@ -37,217 +39,174 @@ class TDEM_bDerivTests(unittest.TestCase): self.prb.pair(survey) self.mesh = mesh - # def test_AhVec(self): - # """ - # Test that fields and AhVec produce consistent results - # """ + def test_AhVec(self): + """ + Test that fields and AhVec produce consistent results + """ - # prb = self.prb - # sigma = self.sigma - - # u = prb.fields(sigma) - # Ahu = prb.AhVec(sigma, u) - - # V1 = Ahu[:,'b',1] - # V2 = 1./prb.timeSteps[0]*prb.MfMui*u[:,'b',0] - # self.assertLess(np.linalg.norm(V1-V2)/np.linalg.norm(V2), 1.e-6) - - # V1 = Ahu[:,'e',1] - # self.assertLess(np.linalg.norm(V1), 1.e-6) - - # for i in range(2,prb.nT): - - # dt = prb.timeSteps[i] - - # V1 = Ahu[:,'b',i] - # V2 = 1.0/dt*prb.MfMui*u[:,'b', i-1] - # # print np.linalg.norm(V1), np.linalg.norm(V2) - # self.assertLess(np.linalg.norm(V1)/np.linalg.norm(V2), 1.e-6) - - # V1 = Ahu[:,'e',i] - # V2 = prb.MeSigma*u[:,'e',i] - # # print np.linalg.norm(V1), np.linalg.norm(V2) - # self.assertLess(np.linalg.norm(V1)/np.linalg.norm(V2), 1.e-6) - - # def test_AhVecVSMat_OneTS(self): - - # prb = self.prb - # prb.timeSteps = [1e-05] - # sigma = self.sigma - # prb.curModel = sigma - - # dt = prb.timeSteps[0] - # a11 = 1/dt*prb.MfMui*sp.eye(prb.mesh.nF) - # a12 = prb.MfMui*prb.mesh.edgeCurl - # a21 = prb.mesh.edgeCurl.T*prb.MfMui - # a22 = -prb.MeSigma - # A = sp.bmat([[a11,a12],[a21,a22]]) - - # f = prb.fields(sigma) - # u1 = A*f.tovec() - # u2 = prb.AhVec(sigma,f).tovec() - - # self.assertTrue(np.linalg.norm(u1-u2)/np.linalg.norm(u1)<1e-12) - - # def test_solveAhVSMat_OneTS(self): - # prb = self.prb - - # prb.timeSteps = [1e-05] - - # sigma = self.sigma - # prb.curModel = sigma - - # dt = prb.timeSteps[0] - # a11 = 1.0/dt*prb.MfMui*sp.eye(prb.mesh.nF) - # a12 = prb.MfMui*prb.mesh.edgeCurl - # a21 = prb.mesh.edgeCurl.T*prb.MfMui - # a22 = -prb.MeSigma - # A = sp.bmat([[a11,a12],[a21,a22]]) - - # f = prb.fields(sigma) - # f[:,:,0] = {'e':0,'b':0} - # f[:,'b',1] = 0 - # f[:,'e',1] = np.random.rand(prb.mesh.nE,1) - - # self.assertTrue(np.all(np.r_[f[:,'b',1],f[:,'e',1]] == f.tovec())) - - # u1 = prb.solveAh(sigma,f).tovec().flatten() - # u2 = sp.linalg.spsolve(A.tocsr(),f.tovec()) - - # self.assertLess(np.linalg.norm(u1-u2),1e-8) - - # def test_solveAhVsAhVec(self): - - # prb = self.prb - # mesh = self.prb.mesh - # sigma = self.sigma - # self.prb.curModel = sigma - - # f = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) - # f[:,'b',:] = 0.0 - # for i in range(prb.nT): - # f[:,'e', i] = np.random.rand(mesh.nE, 1) - - # Ahf = prb.AhVec(sigma, f) - # f_test = prb.solveAh(sigma, Ahf) - - # u1 = f.tovec() - # u2 = f_test.tovec() - # self.assertTrue(np.linalg.norm(u1-u2)<1e-8) - - # def test_DerivG(self): - # """ - # Test the derivative of c with respect to sigma - # """ - - # # Random model and perturbation - # sigma = np.random.rand(self.prb.mapping.nP) - - # f = self.prb.fields(sigma) - # dm = 1000*np.random.rand(self.prb.mapping.nP) - # h = 0.01 - - # derChk = lambda m: [self.prb.AhVec(m, f).tovec(), lambda mx: self.prb.Gvec(sigma, mx, u=f).tovec()] - # print '\ntest_DerivG' - # passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20) - # self.assertTrue(passed) - - # def test_Deriv_dUdM(self): - - # prb = self.prb - # prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)] - # mesh = self.mesh - # sigma = self.sigma - - # dm = 10*np.random.rand(prb.mapping.nP) - # f = prb.fields(sigma) - - # derChk = lambda m: [self.prb.fields(m).tovec(), lambda mx: -prb.solveAh(sigma, prb.Gvec(sigma, mx, u=f)).tovec()] - # print '\n' - # print 'test_Deriv_dUdM' - # passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20) - # self.assertTrue(passed) - - # def test_Deriv_J(self): - - # prb = self.prb - # prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)] - # mesh = self.mesh - # sigma = self.sigma - - # # d_sig = 0.8*sigma #np.random.rand(mesh.nCz) - # d_sig = 10*np.random.rand(prb.mapping.nP) - - - # derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(sigma, mx)] - # print '\n' - # print 'test_Deriv_J' - # passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=4, eps=1e-20) - # self.assertTrue(passed) - - # def test_projectAdjoint(self): - # prb = self.prb - # survey = prb.survey - # mesh = self.mesh - - # # Generate random fields and data - # f = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) - # for i in range(prb.nT): - # f[:,'b',i] = np.random.rand(mesh.nF, 1) - # f[:,'e',i] = np.random.rand(mesh.nE, 1) - # d_vec = np.random.rand(survey.nD, survey.nTx).flatten() - # d = Survey.Data(survey,v=d_vec) - - # # Check that d.T*Q*f = f.T*Q.T*d - # V1 = d_vec.dot(survey.projectFieldsDeriv(None, v=f).tovec()) - # V2 = f.tovec().dot(survey.projectFieldsDeriv(None, v=d, adjoint=True).tovec()) - - # self.assertLess((V1-V2)/np.abs(V1), 1e-6) - - # def test_adjointAhVsAht(self): - # prb = self.prb - # mesh = self.mesh - # sigma = self.sigma - - # f1 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) - # for i in range(1,prb.nT+1): - # f1[:,'b',i] = np.random.rand(mesh.nF, 1) - # f1[:,'e',i] = np.random.rand(mesh.nE, 1) - - # f2 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) - # for i in range(1,prb.nT+1): - # f2[:,'b',i] = np.random.rand(mesh.nF, 1) - # f2[:,'e',i] = np.random.rand(mesh.nE, 1) - - # V1 = f2.tovec().dot(prb.AhVec(sigma, f1).tovec()) - # V2 = f1.tovec().dot(prb.AhtVec(sigma, f2).tovec()) - # self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6) - - def test_solveAhtVsAhtVec(self): prb = self.prb - mesh = self.mesh - sigma = np.random.rand(prb.mapping.nP) + sigma = self.sigma - f1 = EM.TDEM.FieldsTDEM(mesh,prb.survey) + u = prb.fields(sigma) + Ahu = prb.AhVec(sigma, u) + + V1 = Ahu[:,'b',1] + V2 = 1./prb.timeSteps[0]*prb.MfMui*u[:,'b',0] + self.assertLess(np.linalg.norm(V1-V2)/np.linalg.norm(V2), 1.e-6) + + V1 = Ahu[:,'e',1] + self.assertLess(np.linalg.norm(V1), 1.e-6) + + for i in range(2,prb.nT): + + dt = prb.timeSteps[i] + + V1 = Ahu[:,'b',i] + V2 = 1.0/dt*prb.MfMui*u[:,'b', i-1] + # print np.linalg.norm(V1), np.linalg.norm(V2) + self.assertLess(np.linalg.norm(V1)/np.linalg.norm(V2), 1.e-6) + + V1 = Ahu[:,'e',i] + V2 = prb.MeSigma*u[:,'e',i] + # print np.linalg.norm(V1), np.linalg.norm(V2) + self.assertLess(np.linalg.norm(V1)/np.linalg.norm(V2), 1.e-6) + + def test_AhVecVSMat_OneTS(self): + + prb = self.prb + prb.timeSteps = [1e-05] + sigma = self.sigma + prb.curModel = sigma + + dt = prb.timeSteps[0] + a11 = 1/dt*prb.MfMui*sp.eye(prb.mesh.nF) + a12 = prb.MfMui*prb.mesh.edgeCurl + a21 = prb.mesh.edgeCurl.T*prb.MfMui + a22 = -prb.MeSigma + A = sp.bmat([[a11,a12],[a21,a22]]) + + f = prb.fields(sigma) + u1 = A*f.tovec() + u2 = prb.AhVec(sigma,f).tovec() + + self.assertTrue(np.linalg.norm(u1-u2)/np.linalg.norm(u1)<1e-12) + + def test_solveAhVSMat_OneTS(self): + prb = self.prb + + prb.timeSteps = [1e-05] + + sigma = self.sigma + prb.curModel = sigma + + dt = prb.timeSteps[0] + a11 = 1.0/dt*prb.MfMui*sp.eye(prb.mesh.nF) + a12 = prb.MfMui*prb.mesh.edgeCurl + a21 = prb.mesh.edgeCurl.T*prb.MfMui + a22 = -prb.MeSigma + A = sp.bmat([[a11,a12],[a21,a22]]) + + f = prb.fields(sigma) + f[:,:,0] = {'e':0,'b':0} + f[:,'b',1] = 0 + f[:,'e',1] = np.random.rand(prb.mesh.nE,1) + + self.assertTrue(np.all(np.r_[f[:,'b',1],f[:,'e',1]] == f.tovec())) + + u1 = prb.solveAh(sigma,f).tovec().flatten() + u2 = sp.linalg.spsolve(A.tocsr(),f.tovec()) + + self.assertLess(np.linalg.norm(u1-u2),1e-8) + + def test_solveAhVsAhVec(self): + + prb = self.prb + mesh = self.prb.mesh + sigma = self.sigma + self.prb.curModel = sigma + + f = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) + f[:,'b',:] = 0.0 for i in range(prb.nT): - f1[:,'b',i] = np.random.rand(mesh.nF, 1) - f1[:,'e',i] = np.random.rand(mesh.nE, 1) + f[:,'e', i] = np.random.rand(mesh.nE, 1) - f2 = prb.solveAht(sigma, f1) - f3 = prb.AhtVec(sigma, f2) + Ahf = prb.AhVec(sigma, f) + f_test = prb.solveAh(sigma, Ahf) - if plotIt: - import matplotlib.pyplot as plt - plt.plot(f3.tovec()) - plt.plot(f1.tovec()) - plt.show() - V1 = np.linalg.norm(f3.tovec()-f1.tovec()) - V2 = np.linalg.norm(f1.tovec()) - print V1, V2 - print 'I am gunna fail this one: boo. :(' - self.assertLess(V1/V2, 1e-6) + u1 = f.tovec() + u2 = f_test.tovec() + self.assertTrue(np.linalg.norm(u1-u2)<1e-8) - def test_adjointsolveAhVssolveAht(self): + def test_DerivG(self): + """ + Test the derivative of c with respect to sigma + """ + + # Random model and perturbation + sigma = np.random.rand(self.prb.mapping.nP) + + f = self.prb.fields(sigma) + dm = 1000*np.random.rand(self.prb.mapping.nP) + h = 0.01 + + derChk = lambda m: [self.prb.AhVec(m, f).tovec(), lambda mx: self.prb.Gvec(sigma, mx, u=f).tovec()] + print '\ntest_DerivG' + passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20) + self.assertTrue(passed) + + def test_Deriv_dUdM(self): + + prb = self.prb + prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)] + mesh = self.mesh + sigma = self.sigma + + dm = 10*np.random.rand(prb.mapping.nP) + f = prb.fields(sigma) + + derChk = lambda m: [self.prb.fields(m).tovec(), lambda mx: -prb.solveAh(sigma, prb.Gvec(sigma, mx, u=f)).tovec()] + print '\n' + print 'test_Deriv_dUdM' + passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20) + self.assertTrue(passed) + + def test_Deriv_J(self): + + prb = self.prb + prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)] + mesh = self.mesh + sigma = self.sigma + + # d_sig = 0.8*sigma #np.random.rand(mesh.nCz) + d_sig = 10*np.random.rand(prb.mapping.nP) + + + derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(sigma, mx)] + print '\n' + print 'test_Deriv_J' + passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=4, eps=1e-20) + self.assertTrue(passed) + + def test_projectAdjoint(self): + prb = self.prb + survey = prb.survey + mesh = self.mesh + + # Generate random fields and data + f = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) + for i in range(prb.nT): + f[:,'b',i] = np.random.rand(mesh.nF, 1) + f[:,'e',i] = np.random.rand(mesh.nE, 1) + d_vec = np.random.rand(survey.nD, survey.nTx).flatten() + d = Survey.Data(survey,v=d_vec) + + # Check that d.T*Q*f = f.T*Q.T*d + V1 = d_vec.dot(survey.projectFieldsDeriv(None, v=f).tovec()) + V2 = f.tovec().dot(survey.projectFieldsDeriv(None, v=d, adjoint=True).tovec()) + + self.assertLess((V1-V2)/np.abs(V1), 1e-6) + + def test_adjointAhVsAht(self): prb = self.prb mesh = self.mesh sigma = self.sigma @@ -262,43 +221,88 @@ class TDEM_bDerivTests(unittest.TestCase): f2[:,'b',i] = np.random.rand(mesh.nF, 1) f2[:,'e',i] = np.random.rand(mesh.nE, 1) - V1 = f2.tovec().dot(prb.solveAh(sigma, f1).tovec()) - V2 = f1.tovec().dot(prb.solveAht(sigma, f2).tovec()) + V1 = f2.tovec().dot(prb.AhVec(sigma, f1).tovec()) + V2 = f1.tovec().dot(prb.AhtVec(sigma, f2).tovec()) self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6) - # def test_adjointGvecVsGtvec(self): - # mesh = self.mesh + # def test_solveAhtVsAhtVec(self): # prb = self.prb - - # m = np.random.rand(prb.mapping.nP) + # mesh = self.mesh # sigma = np.random.rand(prb.mapping.nP) - # u = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) + # f1 = EM.TDEM.FieldsTDEM(mesh,prb.survey) # for i in range(1,prb.nT+1): - # u[:,'b',i] = np.random.rand(mesh.nF, 1) - # u[:,'e',i] = np.random.rand(mesh.nE, 1) + # f1[:,'b',i] = np.random.rand(mesh.nF, 1) + # f1[:,'e',i] = np.random.rand(mesh.nE, 1) - # v = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) - # for i in range(1,prb.nT+1): - # v[:,'b',i] = np.random.rand(mesh.nF, 1) - # v[:,'e',i] = np.random.rand(mesh.nE, 1) + # f2 = prb.solveAht(sigma, f1) + # f3 = prb.AhtVec(sigma, f2) - # V1 = m.dot(prb.Gtvec(sigma, v, u)) - # V2 = v.tovec().dot(prb.Gvec(sigma, m, u).tovec()) - # self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6) + # if True: + # import matplotlib.pyplot as plt + # plt.plot(f3.tovec(),'b') + # plt.plot(f1.tovec(),'r') + # plt.show() + # V1 = np.linalg.norm(f3.tovec()-f1.tovec()) + # V2 = np.linalg.norm(f1.tovec()) + # print 'AhtVsAhtVec', V1, V2, f1.tovec() + # print 'I am gunna fail this one: boo. :(' + # self.assertLess(V1/V2, 1e-6) - # def test_adjointJvecVsJtvec(self): - # mesh = self.mesh + # def test_adjointsolveAhVssolveAht(self): # prb = self.prb + # mesh = self.mesh # sigma = self.sigma - # m = np.random.rand(prb.mapping.nP) - # d = np.random.rand(prb.survey.nD) + # f1 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) + # for i in range(1,prb.nT+1): + # f1[:,'b',i] = np.random.rand(mesh.nF, 1) + # f1[:,'e',i] = np.random.rand(mesh.nE, 1) - # V1 = d.dot(prb.Jvec(sigma, m)) - # V2 = m.dot(prb.Jtvec(sigma, d)) + # f2 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) + # for i in range(1,prb.nT+1): + # f2[:,'b',i] = np.random.rand(mesh.nF, 1) + # f2[:,'e',i] = np.random.rand(mesh.nE, 1) + + # V1 = f2.tovec().dot(prb.solveAh(sigma, f1).tovec()) + # V2 = f1.tovec().dot(prb.solveAht(sigma, f2).tovec()) + # print V1, V2 # self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6) + def test_adjointGvecVsGtvec(self): + mesh = self.mesh + prb = self.prb + + m = np.random.rand(prb.mapping.nP) + sigma = np.random.rand(prb.mapping.nP) + + u = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) + for i in range(1,prb.nT+1): + u[:,'b',i] = np.random.rand(mesh.nF, 1) + u[:,'e',i] = np.random.rand(mesh.nE, 1) + + v = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) + for i in range(1,prb.nT+1): + v[:,'b',i] = np.random.rand(mesh.nF, 1) + v[:,'e',i] = np.random.rand(mesh.nE, 1) + + V1 = m.dot(prb.Gtvec(sigma, v, u)) + V2 = v.tovec().dot(prb.Gvec(sigma, m, u).tovec()) + self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6) + + def test_adjointJvecVsJtvec(self): + mesh = self.mesh + prb = self.prb + sigma = self.sigma + + m = np.random.rand(prb.mapping.nP) + d = np.random.rand(prb.survey.nD) + + V1 = d.dot(prb.Jvec(sigma, m)) + V2 = m.dot(prb.Jtvec(sigma, d)) + print 'AdjointTest', V1, V2 + self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6) + if __name__ == '__main__': From 8e5c93accb73493242fb3bc13c25fd89c37421ae Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Mon, 28 Apr 2014 12:34:44 -0700 Subject: [PATCH 8/8] Updates to derivs for multiple Txs --- simpegEM/TDEM/SurveyTDEM.py | 4 +- simpegEM/TDEM/TDEM_b.py | 15 +- simpegEM/Tests/test_TDEM_b_DerivAdjoint.py | 2 +- .../Tests/test_TDEM_b_MultiTx_DerivAdjoint.py | 150 ++++++++++++++++++ 4 files changed, 165 insertions(+), 6 deletions(-) create mode 100644 simpegEM/Tests/test_TDEM_b_MultiTx_DerivAdjoint.py diff --git a/simpegEM/TDEM/SurveyTDEM.py b/simpegEM/TDEM/SurveyTDEM.py index 722846f5..0e04f10a 100644 --- a/simpegEM/TDEM/SurveyTDEM.py +++ b/simpegEM/TDEM/SurveyTDEM.py @@ -60,8 +60,8 @@ class FieldsTDEM(Survey.TimeFields): e = self[:,'e',i+1] else: e = np.zeros(nE if nTx == 1 else (nE, nTx)) - u = np.r_[u, b, e] - return u + u = np.concatenate((u, b, e)) + return Utils.mkvc(u) class TxTDEM(Survey.BaseTx): rxPair = RxTDEM diff --git a/simpegEM/TDEM/TDEM_b.py b/simpegEM/TDEM/TDEM_b.py index 01c6a50b..9622af14 100644 --- a/simpegEM/TDEM/TDEM_b.py +++ b/simpegEM/TDEM/TDEM_b.py @@ -64,7 +64,7 @@ class ProblemTDEM_b(BaseTDEMProblem): p = self.survey.projectFieldsDeriv(u, v=v, adjoint=True) y = self.solveAht(m, p) w = self.Gtvec(m, y, u) - return - w + return - mkvc(w) def Gvec(self, m, vec, u=None): """ @@ -107,10 +107,13 @@ class ProblemTDEM_b(BaseTDEMProblem): if u is None: u = self.fields(m) nTx, nE = self.survey.nTx, self.mesh.nE - tmp = np.zeros(nE if nTx == 1 else (nE,nTx)) + tmp = np.zeros(nE) # Here we can do internal multiplications of Gt*v and then multiply by MsigDeriv.T in one go. for i in range(1,self.nT+1): - tmp += vec[:,'e',i]*u[:,'e',i] + vu = vec[:,'e',i]*u[:,'e',i] + if nTx > 1: + vu = vu.sum(axis=1) + tmp += vu curModel = self.mapping.transform(m) p = -mkvc(self.mapping.transformDeriv(m).T*self.mesh.getEdgeInnerProductDeriv(curModel).T*tmp) @@ -137,6 +140,12 @@ class ProblemTDEM_b(BaseTDEMProblem): def solveAht(self, m, p): + # Mini Example: + # + # nT = 3, len(times) == 4, fields stored in F[:,:,1:4] + # + # 0 is held for initial conditions (this shifts the storage by +1) + # ^ # fLoc 0 1 2 3 # |-----|-----|-----| # tInd 0 1 2 / / diff --git a/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py b/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py index 841cef25..36176944 100644 --- a/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py +++ b/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py @@ -197,7 +197,7 @@ class TDEM_bDerivTests(unittest.TestCase): for i in range(prb.nT): f[:,'b',i] = np.random.rand(mesh.nF, 1) f[:,'e',i] = np.random.rand(mesh.nE, 1) - d_vec = np.random.rand(survey.nD, survey.nTx).flatten() + d_vec = np.random.rand(survey.nD) d = Survey.Data(survey,v=d_vec) # Check that d.T*Q*f = f.T*Q.T*d diff --git a/simpegEM/Tests/test_TDEM_b_MultiTx_DerivAdjoint.py b/simpegEM/Tests/test_TDEM_b_MultiTx_DerivAdjoint.py new file mode 100644 index 00000000..477952e1 --- /dev/null +++ b/simpegEM/Tests/test_TDEM_b_MultiTx_DerivAdjoint.py @@ -0,0 +1,150 @@ +import unittest +from SimPEG import * +import simpegEM as EM + +plotIt = False + +class TDEM_bDerivTests(unittest.TestCase): + + def setUp(self): + + cs = 5. + ncx = 20 + ncy = 6 + npad = 20 + hx = [(cs,ncx), (cs,npad,1.3)] + hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)] + mesh = Mesh.CylMesh([hx,1,hy], '00C') + + 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]) + + rxOffset = 40. + rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20), 'bz') + tx = EM.TDEM.TxTDEM(np.array([0., 0., 0.]), 'VMD_MVP', [rx]) + rx2 = EM.TDEM.RxTDEM(np.array([[rxOffset-10, 0., 0.]]), np.logspace(-5,-4, 25), 'bz') + tx2 = EM.TDEM.TxTDEM(np.array([0., 0., 0.]), 'VMD_MVP', [rx2]) + + survey = EM.TDEM.SurveyTDEM([tx,tx2]) + + self.prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping) + # self.prb.timeSteps = [1e-5] + self.prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)] + # self.prb.timeSteps = [(1e-05, 100)] + + self.sigma = np.ones(mesh.nCz)*1e-8 + self.sigma[mesh.vectorCCz<0] = 1e-1 + self.sigma = np.log(self.sigma[active]) + + self.prb.pair(survey) + self.mesh = mesh + + def test_DerivG(self): + """ + Test the derivative of c with respect to sigma + """ + + # Random model and perturbation + sigma = np.random.rand(self.prb.mapping.nP) + + f = self.prb.fields(sigma) + dm = 1000*np.random.rand(self.prb.mapping.nP) + h = 0.01 + + derChk = lambda m: [self.prb.AhVec(m, f).tovec(), lambda mx: self.prb.Gvec(sigma, mx, u=f).tovec()] + print '\ntest_DerivG' + passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20) + self.assertTrue(passed) + + def test_Deriv_dUdM(self): + + prb = self.prb + prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)] + mesh = self.mesh + sigma = self.sigma + + dm = 10*np.random.rand(prb.mapping.nP) + f = prb.fields(sigma) + + derChk = lambda m: [self.prb.fields(m).tovec(), lambda mx: -prb.solveAh(sigma, prb.Gvec(sigma, mx, u=f)).tovec()] + print '\n' + print 'test_Deriv_dUdM' + passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20) + self.assertTrue(passed) + + def test_Deriv_J(self): + + prb = self.prb + prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)] + mesh = self.mesh + sigma = self.sigma + + # d_sig = 0.8*sigma #np.random.rand(mesh.nCz) + d_sig = 10*np.random.rand(prb.mapping.nP) + + + derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(sigma, mx)] + print '\n' + print 'test_Deriv_J' + passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=4, eps=1e-20) + self.assertTrue(passed) + + def test_projectAdjoint(self): + prb = self.prb + survey = prb.survey + mesh = self.mesh + + # Generate random fields and data + f = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) + for i in range(prb.nT): + f[:,'b',i] = np.random.rand(mesh.nF, 1) + f[:,'e',i] = np.random.rand(mesh.nE, 1) + d_vec = np.random.rand(survey.nD) + d = Survey.Data(survey,v=d_vec) + + # Check that d.T*Q*f = f.T*Q.T*d + V1 = d_vec.dot(survey.projectFieldsDeriv(None, v=f).tovec()) + V2 = f.tovec().dot(survey.projectFieldsDeriv(None, v=d, adjoint=True).tovec()) + + self.assertLess((V1-V2)/np.abs(V1), 1e-6) + + def test_adjointGvecVsGtvec(self): + mesh = self.mesh + prb = self.prb + + m = np.random.rand(prb.mapping.nP) + sigma = np.random.rand(prb.mapping.nP) + + u = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) + for i in range(1,prb.nT+1): + u[:,'b',i] = np.random.rand(mesh.nF, 2) + u[:,'e',i] = np.random.rand(mesh.nE, 2) + + v = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) + for i in range(1,prb.nT+1): + v[:,'b',i] = np.random.rand(mesh.nF, 2) + v[:,'e',i] = np.random.rand(mesh.nE, 2) + + V1 = m.dot(prb.Gtvec(sigma, v, u)) + V2 = v.tovec().dot(prb.Gvec(sigma, m, u).tovec()) + self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6) + + def test_adjointJvecVsJtvec(self): + mesh = self.mesh + prb = self.prb + sigma = self.sigma + + m = np.random.rand(prb.mapping.nP) + d = np.random.rand(prb.survey.nD) + + V1 = d.dot(prb.Jvec(sigma, m)) + V2 = m.dot(prb.Jtvec(sigma, d)) + print 'AdjointTest', V1, V2 + self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6) + + + +if __name__ == '__main__': + unittest.main()