From 12b0755dfc753dbde379100b2f8735b0f228daec Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 21 Feb 2014 17:28:10 -0800 Subject: [PATCH] start of the FDEM --- simpegEM/FDEM/DataFDEM.py | 50 ++++++++++++++++++++++ simpegEM/FDEM/FDEM.py | 83 ++++++++++++++++++++++++++----------- simpegEM/FDEM/FieldsFDEM.py | 52 +++++++++++++++++++++++ simpegEM/FDEM/__init__.py | 3 ++ simpegEM/TDEM/TDEM_b.py | 15 ++++--- simpegEM/__init__.py | 3 +- 6 files changed, 175 insertions(+), 31 deletions(-) create mode 100644 simpegEM/FDEM/DataFDEM.py create mode 100644 simpegEM/FDEM/FieldsFDEM.py diff --git a/simpegEM/FDEM/DataFDEM.py b/simpegEM/FDEM/DataFDEM.py new file mode 100644 index 00000000..bd9877b2 --- /dev/null +++ b/simpegEM/FDEM/DataFDEM.py @@ -0,0 +1,50 @@ +from SimPEG import Utils, np +from SimPEG.Data import BaseData +from FieldsFDEM import FieldsFDEM + +class DataFDEM(BaseData): + """ + docstring for DataFDEM + """ + + txLoc = None #: txLoc + txType = None #: txType + nTx = 1 #: Number of transmitters + rxLoc = None #: rxLoc + rxType = None #: rxType + freq = None #: freq + + + @property + def omega(self): + return 2*np.pi*self.freq + + @property + def nFreq(self): + """Number of frequencies""" + return self.freq.size + + def __init__(self, **kwargs): + BaseData.__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: fix this + pass + + #################################################### + # 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 diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index d930cf4a..6580bac9 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -1,7 +1,9 @@ -from SimPEG import Problem +from SimPEG import Problem, Solver import numpy as np from scipy.constants import mu_0 from SimPEG.Utils import sdiag, mkvc +from FieldsFDEM import FieldsFDEM +from DataFDEM import DataFDEM class ProblemFDEM_e(Problem.BaseProblem): @@ -18,11 +20,18 @@ class ProblemFDEM_e(Problem.BaseProblem): Problem.BaseProblem.__init__(self, mesh, model, **kwargs) solType = 'b' + storeTheseFields = 'e' - #TODO: - # j_s - # getOmega - # getFieldsObject + dataPair = DataFDEM + + solveOpts = {'factorize':False, 'backend':'scipy'} + + j_s = None + + + def getFieldsObject(self): + F = FieldsFDEM(self.mesh, self.data.nTx, self.data.nFreq, store=self.storeTheseFields) + return F #################################################### # Mass Matrices @@ -41,8 +50,10 @@ class ProblemFDEM_e(Problem.BaseProblem): def MeSigmaI(self): return self._MeSigmaI def makeMassMatrices(self, m): + #TODO: hardcoded to sigma as the model + sigma = self.model.transform(m) self._Me = self.mesh.getEdgeInnerProduct() - self._MeSigma = self.mesh.getEdgeInnerProduct(m) + self._MeSigma = self.mesh.getEdgeInnerProduct(sigma) # TODO: this will not work if tensor conductivity self._MeSigmaI = sdiag(1/self.MeSigma.diagonal()) #TODO: assuming constant mu @@ -52,17 +63,17 @@ class ProblemFDEM_e(Problem.BaseProblem): # Internal Methods #################################################### - def getA(self, omegaInd): + def getA(self, freqInd): """ :param int tInd: Time index :rtype: scipy.sparse.csr_matrix :return: A """ - omega = self.getOmega(omegaInd) + omega = self.data.omega[freqInd] return self.mesh.edgeCurl.T*self.MfMui*self.mesh.edgeCurl + 1j*omega*self.MeSigma - def getRHS(self, omegaInd): - omega = self.getOmega(omegaInd) + def getRHS(self, freqInd): + omega = self.data.omega[freqInd] return -1j*omega*self.Me*self.j_s @@ -73,8 +84,18 @@ class ProblemFDEM_e(Problem.BaseProblem): F = self.getFieldsObject() + for freqInd in range(self.data.nFreq): + A = self.getA(freqInd) + b = self.getRHS(freqInd) + e = Solver(A, options=self.solveOpts).solve(b) - return + F.set_e(e, freqInd) + omega = self.data.omega[freqInd] + #TODO: check if mass matrices needed: + b = -1./(1j*omega)*self.mesh.edgeCurl*e + F.set_b(b, freqInd) + + return F def Jvec(self, m, v, u=None): @@ -99,28 +120,42 @@ if __name__ == '__main__': import matplotlib.pyplot as plt cs = 5. - ncx = 20 + ncx = 6 ncy = 6 - npad = 20 - hx = Utils.meshTensors(((0,cs), (ncx,cs), (npad,cs))) + ncz = 6 + npad = 3 + hx = Utils.meshTensors(((npad,cs), (ncx,cs), (npad,cs))) hy = Utils.meshTensors(((npad,cs), (ncy,cs), (npad,cs))) - mesh = Mesh.Cyl1DMesh([hx,hy], -hy.sum()/2) - model = Model.Vertical1DModel(mesh) + hz = Utils.meshTensors(((npad,cs), (ncz,cs), (npad,cs))) + mesh = Mesh.TensorMesh([hx,hy,hz]) + + XY = Utils.ndgrid(np.linspace(20,50,3), np.linspace(20,50,3)) + rxLoc = np.c_[XY, np.ones(XY.shape[0])*40] + + model = Model.LogModel(mesh) opts = {'txLoc':0., 'txType':'VMD_MVP', - 'rxLoc':np.r_[150., 0.], + 'rxLoc': rxLoc, 'rxType':'bz', - 'timeCh':np.logspace(-4,-2,20), + 'freq': np.logspace(0,3,4), } - dat = EM.TDEM.DataTDEM1D(**opts) + dat = EM.FDEM.DataFDEM(**opts) - prb = EM.TDEM.ProblemTDEM_b(mesh, model) - # prb.setTimes([1e-5, 5e-5, 2.5e-4], [150, 150, 150]) - # prb.setTimes([1e-5, 5e-5, 2.5e-4], [10, 10, 10]) - prb.setTimes([1e-5], [1]) + prb = EM.FDEM.ProblemFDEM_e(mesh, model) prb.pair(dat) - sigma = np.random.rand(mesh.nCz) + + sigma = np.log(np.ones(mesh.nC)*1e-3) + + j_sx = np.zeros(mesh.vnEx) + j_sx[6,6,6] = 1 + j_s = np.r_[Utils.mkvc(j_sx),np.zeros(mesh.nEy+mesh.nEz)] + + prb.j_s = j_s + f = prb.fields(sigma) + + colorbar(mesh.plotSlice((f.get_e(3)), 'E', ind=11, normal='Z', view='real')[0]) + plt.show() diff --git a/simpegEM/FDEM/FieldsFDEM.py b/simpegEM/FDEM/FieldsFDEM.py new file mode 100644 index 00000000..93add159 --- /dev/null +++ b/simpegEM/FDEM/FieldsFDEM.py @@ -0,0 +1,52 @@ +import numpy as np + + +class FieldsFDEM(object): + """docstring for FieldsFDEM""" + + 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, nFreq, store='e'): + + self.nFreq = nFreq #: Number of times + self.nTx = nTx #: Number of transmitters + self.mesh = mesh + + def update(self, newFields, fInd): + self.set_b(newFields['b'], fInd) + self.set_e(newFields['e'], fInd) + + #################################################### + # Get Methods + #################################################### + + def get_b(self, ind): + return self.b[ind,:,:] + + def get_e(self, ind): + return self.e[ind,:,:] + + #################################################### + # Set Methods + #################################################### + + def set_b(self, b, ind): + if self.b is None: + self.b = np.zeros((self.nFreq, np.sum(self.mesh.nF), self.nTx), dtype=complex) + 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.nFreq, np.sum(self.mesh.nE), self.nTx), dtype=complex) + self.e[:] = np.nan + if len(e.shape) == 1: + e = e[:, np.newaxis] + self.e[ind,:,:] = e diff --git a/simpegEM/FDEM/__init__.py b/simpegEM/FDEM/__init__.py index e69de29b..2accd044 100644 --- a/simpegEM/FDEM/__init__.py +++ b/simpegEM/FDEM/__init__.py @@ -0,0 +1,3 @@ +from FieldsFDEM import FieldsFDEM +from DataFDEM import DataFDEM +from FDEM import ProblemFDEM_e diff --git a/simpegEM/TDEM/TDEM_b.py b/simpegEM/TDEM/TDEM_b.py index cbdf6721..6064120b 100644 --- a/simpegEM/TDEM/TDEM_b.py +++ b/simpegEM/TDEM/TDEM_b.py @@ -2,13 +2,14 @@ from BaseTDEM import ProblemBaseTDEM from FieldsTDEM import FieldsTDEM from SimPEG.Utils import mkvc import numpy as np +from DataTDEM import DataTDEM1D class ProblemTDEM_b(ProblemBaseTDEM): """ Time-Domain EM problem - B-formulation TDEM_b treats the following discretization of Maxwell's equations - + .. math:: \dcurl \e^{(t+1)} + \\frac{\\b^{(t+1)} - \\b^{(t)}}{\delta t} = 0 \\\\ \dcurl^\\top \MfMui \\b^{(t+1)} - \MeSig \e^{(t+1)} = \Me \j_s^{(t+1)} @@ -20,6 +21,8 @@ class ProblemTDEM_b(ProblemBaseTDEM): solType = 'b' + dataPair = DataTDEM1D + #################################################### # Internal Methods #################################################### @@ -66,7 +69,7 @@ class ProblemTDEM_b(ProblemBaseTDEM): :rtype: simpegEM.TDEM.FieldsTDEM :return: f - Multiply G by a vector where + Multiply G by a vector where """ if u is None: u = self.fields(sigma) @@ -108,14 +111,14 @@ class ProblemTDEM_b(ProblemBaseTDEM): return self.forward(m, AhRHS, AhCalcFields) 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) if tInd == self.nTimes-1: return rhs dt = self.getDt(tInd+1) return rhs + 1./dt*self.MfMui*u.get_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) @@ -159,7 +162,7 @@ class ProblemTDEM_b(ProblemBaseTDEM): -\\frac{1}{\delta t} \MfMui & 0 \\\\ 0 & 0 \end{array} - \\right] \\\\ + \\right] \\\\ """ self.makeMassMatrices(sigma) @@ -208,7 +211,7 @@ class ProblemTDEM_b(ProblemBaseTDEM): -\\frac{1}{\delta t} \MfMui & 0 \\\\ 0 & 0 \end{array} - \\right] \\\\ + \\right] \\\\ """ self.makeMassMatrices(sigma) f = FieldsTDEM(self.mesh, 1, self.times.size, 'b') diff --git a/simpegEM/__init__.py b/simpegEM/__init__.py index 4dbea793..13667f2f 100644 --- a/simpegEM/__init__.py +++ b/simpegEM/__init__.py @@ -1,3 +1,4 @@ # from EM import * import Utils -import TDEM \ No newline at end of file +import TDEM +import FDEM