diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index 507127f7..d37bf038 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -1,4 +1,4 @@ -from SimPEG import Problem, Solver, Utils, np, sp +from SimPEG import Problem, Utils, np, sp, Solver as SimpegSolver from scipy.constants import mu_0 from SurveyFDEM import SurveyFDEM, DataFDEM, FieldsFDEM from simpegEM.Utils import Sources @@ -7,7 +7,7 @@ def omega(freq): """Change frequency to angular frequency, omega""" return 2.*np.pi*freq -class ProblemFDEM_e(Problem.BaseProblem): +class BaseProblemFDEM(Problem.BaseProblem): """ Frequency-Domain EM problem - E-formulation @@ -26,8 +26,8 @@ class ProblemFDEM_e(Problem.BaseProblem): surveyPair = SurveyFDEM dataPair = DataFDEM - solveOpts = {'factorize':False, 'backend':'scipy'} - + Solver = SimpegSolver + solverOpts = {'doDirect':True, 'options':{'factorize':False, 'backend':'scipy'}} #################################################### # Mass Matrices @@ -55,9 +55,14 @@ class ProblemFDEM_e(Problem.BaseProblem): #TODO: assuming constant mu self._MfMui = self.mesh.getFaceInnerProduct(1/mu_0) - #################################################### - # Internal Methods - #################################################### + +class ProblemFDEM_e(BaseProblemFDEM): + """ + Solving for e! + """ + def __init__(self, model, **kwargs): + BaseProblemFDEM.__init__(self, model, **kwargs) + def getA(self, freq): """ @@ -85,9 +90,11 @@ class ProblemFDEM_e(Problem.BaseProblem): SRCz = src(tx.loc, self.mesh.gridEz, 'z') rhs[i] = np.concatenate((SRCx, SRCy, SRCz)) - #TODO: this is completely wrong. b0 = self.mesh.edgeCurl*rhs but we are doing an e formulation... - j_s = np.concatenate(rhs).reshape((self.mesh.nE, len(Txs)), order='F') - return -1j*omega(freq)*self.Me*j_s + a = np.concatenate(rhs).reshape((self.mesh.nE, len(Txs)), order='F') + C = self.mesh.edgeCurl + j_s = C.T*self.MfMui*C*a + #TODO: self.Me* ?? + return -1j*omega(freq)*j_s def fields(self, m, useThisRhs=None): @@ -99,11 +106,13 @@ class ProblemFDEM_e(Problem.BaseProblem): for freq in self.survey.freqs: A = self.getA(freq) - b = self.getRHS(freq) - e = Solver(A, options=self.solveOpts).solve(b) + rhs = self.getRHS(freq) + solver = self.Solver(A, **self.solverOpts) + e = solver.solve(rhs) + + print np.linalg.norm(A*Utils.mkvc(e) - Utils.mkvc(rhs)) / np.linalg.norm(Utils.mkvc(rhs)) F[freq, 'e'] = e - #TODO: check if mass matrices needed: b = -1./(1j*omega(freq))*self.mesh.edgeCurl*e F[freq, 'b'] = b @@ -121,7 +130,7 @@ class ProblemFDEM_e(Problem.BaseProblem): for i, freq in enumerate(self.survey.freqs): e = u[freq, 'e'] A = self.getA(freq) - solver = Solver(A, options=self.solveOpts) + solver = self.Solver(A, **self.solverOpts) for txi, tx in enumerate(self.survey.getTransmitters(freq)): dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig, v=e[:,txi]) @@ -150,7 +159,7 @@ class ProblemFDEM_e(Problem.BaseProblem): for i, freq in enumerate(self.survey.freqs): e = u[freq, 'e'] AT = self.getA(freq).T - solver = Solver(AT, options=self.solveOpts) + solver = self.Solver(AT, **self.solverOpts) for txi, tx in enumerate(self.survey.getTransmitters(freq)): dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig, v=e[:,txi]) @@ -161,7 +170,119 @@ class ProblemFDEM_e(Problem.BaseProblem): return Jtv +class ProblemFDEM_b(BaseProblemFDEM): + """ + Solving for b! + """ + def __init__(self, model, **kwargs): + BaseProblemFDEM.__init__(self, model, **kwargs) + def getA(self, freq): + """ + :param float freq: Frequency + :rtype: scipy.sparse.csr_matrix + :return: A + """ + return self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui + 1j*omega(freq)*self.MfMui + + def getRHS(self, freq): + """ + :param float freq: Frequency + :rtype: numpy.ndarray (nE, nTx) + :return: RHS + """ + Txs = self.survey.getTransmitters(freq) + rhs = range(len(Txs)) + for i, tx in enumerate(Txs): + if tx.txType == 'VMD': + src = Sources.MagneticDipoleVectorPotential + else: + raise NotImplemented('%s txType is not implemented' % tx.txType) + SRCx = src(tx.loc, self.mesh.gridEx, 'x') + SRCy = src(tx.loc, self.mesh.gridEy, 'y') + SRCz = src(tx.loc, self.mesh.gridEz, 'z') + rhs[i] = np.concatenate((SRCx, SRCy, SRCz)) + + a = np.concatenate(rhs).reshape((self.mesh.nE, len(Txs)), order='F') + C = self.mesh.edgeCurl + b_0 = C*a + return -1j*omega(freq)*self.MfMui*b_0 + def fields(self, m, useThisRhs=None): + RHS = useThisRhs or self.getRHS + + self.makeMassMatrices(m) + + F = FieldsFDEM(self.mesh, self.survey) + + for freq in self.survey.freqs: + A = self.getA(freq) + rhs = self.getRHS(freq) + solver = self.Solver(A, **self.solverOpts) + #Note that we are solving for b_s + b = solver.solve(rhs) + + print np.linalg.norm(A*Utils.mkvc(b) - Utils.mkvc(rhs)) / np.linalg.norm(Utils.mkvc(rhs)) + + F[freq, 'b'] = b + e = self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui*b + F[freq, 'e'] = e + + return F + + + def Jvec(self, m, v, u=None): + if u is None: + u = self.fields(m) + + raise NotImplemented('') + + # Jv = self.dataPair(self.survey) + # sig = self.model.transform(m) + # dsig_dm = self.model.transformDeriv(m) + + # for i, freq in enumerate(self.survey.freqs): + # e = u[freq, 'e'] + # A = self.getA(freq) + # solver = self.Solver(A, **self.solverOpts) + + # for txi, tx in enumerate(self.survey.getTransmitters(freq)): + # dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig, v=e[:,txi]) + + # P = tx.projectFieldsDeriv(self.mesh, u) + # b = 1j*omega(freq) * ( dMe_dsig * ( dsig_dm * v ) ) + # Ainvb = solver.solve(b) + # Jv[tx] = -P*Ainvb + + # return Utils.mkvc(Jv) + + + def Jtvec(self, m, v, u=None): + if u is None: + u = self.fields(m) + + # Ensure v is a data object. + if not isinstance(v, self.dataPair): + v = self.dataPair(self.survey, v) + + raise NotImplemented('') + # Jtv = np.zeros(self.model.nP, dtype=complex) + + # sig = self.model.transform(m) + # dsig_dm = self.model.transformDeriv(m) + + # for i, freq in enumerate(self.survey.freqs): + # e = u[freq, 'e'] + # AT = self.getA(freq).T + # solver = self.Solver(AT, **self.solverOpts) + + # for txi, tx in enumerate(self.survey.getTransmitters(freq)): + # dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig, v=e[:,txi]) + + # P = tx.projectFieldsDeriv(self.mesh, u) + # w = solver.solve(P.T * v[tx]) + # Jtv += - 1j*omega(freq) * ( dsig_dm.T * ( dMe_dsig.T * w ) ) + + # return Jtv diff --git a/simpegEM/FDEM/__init__.py b/simpegEM/FDEM/__init__.py index f6583225..70124686 100644 --- a/simpegEM/FDEM/__init__.py +++ b/simpegEM/FDEM/__init__.py @@ -1,2 +1,2 @@ from SurveyFDEM import * -from FDEM import ProblemFDEM_e +from FDEM import ProblemFDEM_e, ProblemFDEM_b diff --git a/simpegEM/Tests/test_FDEM.py b/simpegEM/Tests/test_FDEM.py index 23c53528..3726427c 100644 --- a/simpegEM/Tests/test_FDEM.py +++ b/simpegEM/Tests/test_FDEM.py @@ -41,7 +41,7 @@ class FDEM_bDerivTests(unittest.TestCase): x0 = self.sigma def fun(x): return self.survey.dpred(x), lambda x: self.prb.Jvec(x0, x) - passed = Tests.checkDerivative(fun, x0, num=3, plotIt=False) + passed = Tests.checkDerivative(fun, x0, num=3, plotIt=False, eps=1e-18) self.assertTrue(passed) def test_Jtvec_adjointTest(self): diff --git a/simpegEM/Utils/Ana/FEM.py b/simpegEM/Utils/Ana/FEM.py new file mode 100644 index 00000000..df99ece2 --- /dev/null +++ b/simpegEM/Utils/Ana/FEM.py @@ -0,0 +1,17 @@ +import numpy as np +from scipy.constants import mu_0, pi +from scipy.special import erf + +def hzAnalyticDipoleF(r, freq, sigma): + """ + 4.56 in Ward and Hohmann + """ + r = np.abs(r) + k = np.sqrt(-1j*2.*np.pi*freq*mu_0*sigma) + + m = 1 + front = m / (2. * np.pi * (k**2) * (r**5) ) + back = 9 - ( 9 + 9j * k * r - 4 * (k**2) * (r**2) - 1j * (k**3) * (r**3)) * np.exp(-1j*k*r) + hz = front*back + hp =-1/(4*np.pi*r**3) + return hz-hp diff --git a/simpegEM/Utils/Ana/__init__.py b/simpegEM/Utils/Ana/__init__.py index 4553f927..e60286ca 100644 --- a/simpegEM/Utils/Ana/__init__.py +++ b/simpegEM/Utils/Ana/__init__.py @@ -1 +1,2 @@ -from TEM import hzAnalyticDipoleT \ No newline at end of file +from TEM import hzAnalyticDipoleT +from FEM import hzAnalyticDipoleF