From 0a5a7aa5e3174cc14ae206face92da5448c75751 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Wed, 19 Mar 2014 23:58:44 -0700 Subject: [PATCH] Allows calculation of other fields, and derivatives. e.g. b from when solving e, and visa versa --- simpegEM/FDEM/FDEM.py | 94 +++++++++++++++++++++++++++---------- simpegEM/FDEM/SurveyFDEM.py | 21 +++++---- simpegEM/Tests/test_FDEM.py | 8 ++-- 3 files changed, 87 insertions(+), 36 deletions(-) diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index fd510abe..925a4315 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -89,7 +89,8 @@ class BaseProblemFDEM(Problem.BaseProblem): rhs = RHS(freq) solver = self.Solver(A, **self.solverOpts) sol = solver.solve(rhs) - F[freq] = CalcFields(sol, freq) + for fieldType in self.storeTheseFields: + F[freq, fieldType] = CalcFields(sol, freq, fieldType) return F @@ -106,10 +107,17 @@ class BaseProblemFDEM(Problem.BaseProblem): solver = self.Solver(A, **self.solverOpts) for tx in self.survey.getTransmitters(freq): - w = self.getADeriv(freq, u[tx, self.solType], v) + u_tx = u[tx, self.solType] + w = self.getADeriv(freq, u_tx, v) Ainvw = solver.solve(w) + fAinvw = self.calcFields(Ainvw, freq, tx.rxList.fieldType) P = tx.projectFieldsDeriv(self.mesh, u) - Jv[tx] = -P*Ainvw + + df_dm = self.calcFieldsDeriv(u_tx, freq, tx.rxList.fieldType, v) + if df_dm is None: + Jv[tx] = - P*fAinvw + else: + Jv[tx] = - P*fAinvw + P*df_dm return Utils.mkvc(Jv) @@ -131,8 +139,20 @@ class BaseProblemFDEM(Problem.BaseProblem): for tx in self.survey.getTransmitters(freq): P = tx.projectFieldsDeriv(self.mesh, u) - w = solver.solve( - P.T * v[tx]) - Jtv += self.getADeriv(freq, u[tx, self.solType], w, adjoint=True) + + PTv = P.T * v[tx] + fPTv = self.calcFields(PTv, freq, tx.rxList.fieldType, adjoint=True) + + w = solver.solve( fPTv ) + u_tx = u[tx, self.solType] + Jtv_tx = self.getADeriv(freq, u_tx, w, adjoint=True) + + df_dm = self.calcFieldsDeriv(u_tx, freq, tx.rxList.fieldType, PTv, adjoint=True) + + if df_dm is None: + Jtv += - Jtv_tx + else: + Jtv += - Jtv_tx + df_dm return Jtv @@ -192,19 +212,25 @@ class ProblemFDEM_e(BaseProblemFDEM): j_s = C.T*mui*C*a return -1j*omega(freq)*j_s - def calcFields(self, sol, freq): + def calcFields(self, sol, freq, fieldType, adjoint=False): e = sol + if fieldType == 'e': + return e + elif fieldType == 'b': + if not adjoint: + b = -1./(1j*omega(freq))*self.mesh.edgeCurl*e + else: + b = -1./(1j*omega(freq))*self.mesh.edgeCurl.T*e + return b + raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) - fDict = {} - - if 'e' in self.storeTheseFields: - fDict['e'] = e - - if 'b' in self.storeTheseFields: - b = -1./(1j*omega(freq))*self.mesh.edgeCurl*e - fDict['b'] = b - - return fDict + def calcFieldsDeriv(self, sol, freq, fieldType, v, adjoint=False): + e = sol + if fieldType == 'e': + return None + elif fieldType == 'b': + return None + raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) class ProblemFDEM_b(BaseProblemFDEM): @@ -270,17 +296,37 @@ class ProblemFDEM_b(BaseProblemFDEM): b_0 = C*a return -1j*omega(freq)*mui*b_0 - def calcFields(self, sol, freq): + def calcFields(self, sol, freq, fieldType, adjoint=False): b = sol + if fieldType == 'e': + if not adjoint: + e = self.MeSigmaI * ( self.mesh.edgeCurl.T * ( self.MfMui * b ) ) + else: + e = self.MfMui.T * ( self.mesh.edgeCurl * ( self.MeSigmaI.T * b ) ) + return e + elif fieldType == 'b': + return b + raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) - fDict = {} + def calcFieldsDeriv(self, sol, freq, fieldType, v, adjoint=False): + b = sol + if fieldType == 'e': + sig = self.curTModel + dsig_dm = self.curTModelDeriv - if 'b' in self.storeTheseFields: - fDict['b'] = b + C = self.mesh.edgeCurl + mui = self.MfMui - if 'e' in self.storeTheseFields: - e = self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui*b - fDict['e'] = e + #TODO: This only works if diagonal (no tensors)... + dMeSigmaI_dI = - self.MeSigmaI**2 - return fDict + vec = C.T * ( mui * b ) + dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig, v=vec) + if not adjoint: + return dMeSigmaI_dI * ( dMe_dsig * ( dsig_dm * v ) ) + else: + return dsig_dm.T * ( dMe_dsig.T * ( dMeSigmaI_dI.T * v ) ) + elif fieldType == 'b': + return None + raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) diff --git a/simpegEM/FDEM/SurveyFDEM.py b/simpegEM/FDEM/SurveyFDEM.py index 4bb4d1c7..253f2aaf 100644 --- a/simpegEM/FDEM/SurveyFDEM.py +++ b/simpegEM/FDEM/SurveyFDEM.py @@ -3,13 +3,13 @@ from SimPEG import Survey, Utils, np, sp class RxListFDEM(Survey.BaseRxList): knownRxTypes = { - 'Ex':'Ex', - 'Ey':'Ey', - 'Ez':'Ez', + 'ex':'Ex', + 'ey':'Ey', + 'ez':'Ez', - 'Bx':'Fx', - 'By':'Fy', - 'Bz':'Fz', + 'bx':'Fx', + 'by':'Fy', + 'bz':'Fz', } def __init__(self, locs, rxType): @@ -17,6 +17,11 @@ class RxListFDEM(Survey.BaseRxList): self._Ps = {} + @property + def fieldType(self): + #TODO: This assumes that it has the structure ex, by ... + return self.rxType[0] + def getP(self, mesh): if mesh not in self._Ps: self._Ps[mesh] = mesh.getInterpolationMat(self.locs, self.knownRxTypes[self.rxType]) @@ -43,9 +48,9 @@ class TxFDEM(Survey.BaseTx): def projectFields(self, mesh, u): - if self.rxList.rxType in ['Ex', 'Ey', 'Ez']: + if self.rxList.rxType in ['ex', 'ey', 'ez']: u_part = u[self, 'e'] - elif self.rxList.rxType in ['Bx', 'By', 'Bz']: + elif self.rxList.rxType in ['bx', 'by', 'bz']: u_part = u[self, 'b'] else: raise NotImplemented('Unknown receiver type.') diff --git a/simpegEM/Tests/test_FDEM.py b/simpegEM/Tests/test_FDEM.py index 489c9f2f..9212d79e 100644 --- a/simpegEM/Tests/test_FDEM.py +++ b/simpegEM/Tests/test_FDEM.py @@ -18,9 +18,9 @@ def getProblem(fdemType): x = np.linspace(5,10,3) XYZ = Utils.ndgrid(x,x,np.r_[0]) if fdemType == 'e': - rxList = EM.FDEM.RxListFDEM(XYZ, 'Ex') + rxList = EM.FDEM.RxListFDEM(XYZ, 'bx') elif fdemType == 'b': - rxList = EM.FDEM.RxListFDEM(XYZ, 'Bx') + rxList = EM.FDEM.RxListFDEM(XYZ, 'ex') else: raise NotImplementedError() @@ -29,9 +29,9 @@ def getProblem(fdemType): x = np.linspace(5,10,3) XYZ = Utils.ndgrid(x,x,np.r_[0]) if fdemType == 'e': - rxList = EM.FDEM.RxListFDEM(XYZ, 'Ey') + rxList = EM.FDEM.RxListFDEM(XYZ, 'ey') elif fdemType == 'b': - rxList = EM.FDEM.RxListFDEM(XYZ, 'By') + rxList = EM.FDEM.RxListFDEM(XYZ, 'ey') else: raise NotImplementedError()