From 119bc801c7d880e1bd5223e6b31a3801cf866478 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 15 Apr 2016 14:51:19 -0700 Subject: [PATCH] start of DC refactor in EM/Static --- SimPEG/EM/Base.py | 13 ++++-- SimPEG/EM/FDEM/FieldsFDEM.py | 2 +- SimPEG/EM/FDEM/SrcFDEM.py | 1 - SimPEG/EM/FDEM/SurveyFDEM.py | 4 +- SimPEG/EM/Static/DC/FieldsDC.py | 72 +++++++++++++++++++++++++++++ SimPEG/EM/Static/DC/ProblemDC.py | 79 ++++++++++++++++++++++++++++++++ SimPEG/EM/Static/DC/RxDC.py | 66 ++++++++++++++++++++++++++ SimPEG/EM/Static/DC/SrcDC.py | 60 ++++++++++++++++++++++++ SimPEG/EM/Static/DC/SurveyDC.py | 18 ++++++++ SimPEG/EM/Static/DC/__init__.py | 5 ++ 10 files changed, 313 insertions(+), 7 deletions(-) create mode 100644 SimPEG/EM/Static/DC/FieldsDC.py create mode 100644 SimPEG/EM/Static/DC/ProblemDC.py create mode 100644 SimPEG/EM/Static/DC/RxDC.py create mode 100644 SimPEG/EM/Static/DC/SrcDC.py create mode 100644 SimPEG/EM/Static/DC/SurveyDC.py create mode 100644 SimPEG/EM/Static/DC/__init__.py diff --git a/SimPEG/EM/Base.py b/SimPEG/EM/Base.py index a16cdb91..dd349b43 100644 --- a/SimPEG/EM/Base.py +++ b/SimPEG/EM/Base.py @@ -1,6 +1,7 @@ from SimPEG import Survey, Problem, Utils, Models, Maps, PropMaps, np, sp, Solver as SimpegSolver from scipy.constants import mu_0 + class EMPropMap(Maps.PropMap): """ Property Map for EM Problems. The electrical conductivity (\\(\\sigma\\)) is the default inversion property, and the default value of the magnetic permeability is that of free space (\\(\\mu = 4\\pi\\times 10^{-7} \\) H/m) @@ -70,6 +71,12 @@ class BaseEMProblem(Problem.BaseProblem): self._Mf = self.mesh.getFaceInnerProduct() return self._Mf + @property + def Vol(self): + if getattr(self, '_Vol', None) is None: + self._Vol = Utils.sdiag(self.mesh.vol) + return self._Vol + # ----- Magnetic Permeability ----- # @property @@ -192,7 +199,7 @@ class BaseEMSurvey(Survey.BaseSurvey): self.srcList = srcList Survey.BaseSurvey.__init__(self, **kwargs) - def eval(self, u): + def eval(self, f): """ Project fields to receiver locations :param Fields u: fields object @@ -202,8 +209,8 @@ class BaseEMSurvey(Survey.BaseSurvey): data = Survey.Data(self) for src in self.srcList: for rx in src.rxList: - data[src, rx] = rx.eval(src, self.mesh, u) + data[src, rx] = rx.eval(src, self.mesh, f) return data - def evalDeriv(self, u): + def evalDeriv(self, f): raise Exception('Use Receivers to project fields deriv.') diff --git a/SimPEG/EM/FDEM/FieldsFDEM.py b/SimPEG/EM/FDEM/FieldsFDEM.py index e2193973..add2dbc9 100644 --- a/SimPEG/EM/FDEM/FieldsFDEM.py +++ b/SimPEG/EM/FDEM/FieldsFDEM.py @@ -181,7 +181,7 @@ class Fields_e(Fields): } def __init__(self, mesh, survey, **kwargs): - Fields.__init__(self,mesh,survey,**kwargs) + Fields.__init__(self, mesh, survey, **kwargs) def startup(self): self.prob = self.survey.prob diff --git a/SimPEG/EM/FDEM/SrcFDEM.py b/SimPEG/EM/FDEM/SrcFDEM.py index 87967dd5..f1dfef44 100644 --- a/SimPEG/EM/FDEM/SrcFDEM.py +++ b/SimPEG/EM/FDEM/SrcFDEM.py @@ -14,7 +14,6 @@ class BaseSrc(Survey.BaseSrc): def eval(self, prob): """ - Evaluate the source terms. - :math:`s_m` : magnetic source term - :math:`s_e` : electric source term diff --git a/SimPEG/EM/FDEM/SurveyFDEM.py b/SimPEG/EM/FDEM/SurveyFDEM.py index 1552a12c..9659ded6 100644 --- a/SimPEG/EM/FDEM/SurveyFDEM.py +++ b/SimPEG/EM/FDEM/SurveyFDEM.py @@ -63,9 +63,9 @@ class Rx(SimPEG.Survey.BaseRx): """Component projection (real/imag)""" return self.knownRxTypes[self.rxType][2] - def projGLoc(self, u): + def projGLoc(self, f): """Grid Location projection (e.g. Ex Fy ...)""" - return u._GLoc(self.rxType[0]) + self.knownRxTypes[self.rxType][1] + return f._GLoc(self.rxType[0]) + self.knownRxTypes[self.rxType][1] def eval(self, src, mesh, f): """ diff --git a/SimPEG/EM/Static/DC/FieldsDC.py b/SimPEG/EM/Static/DC/FieldsDC.py new file mode 100644 index 00000000..9ca221a0 --- /dev/null +++ b/SimPEG/EM/Static/DC/FieldsDC.py @@ -0,0 +1,72 @@ +import SimPEG +from SimPEG.Utils import Identity, Zero + +class Fields(SimPEG.Problem.Fields): + knownFields = {} + dtype = float + + def _phiDeriv(self, src, du_dm_v, v, adjoint=False): + if getattr(self, '_phiDeriv_u', None) is None or getattr(self, '_phiDeriv_m', None) is None: + raise NotImplementedError ('Getting phiDerivs from %s is not implemented' %self.knownFields.keys()[0]) + + if adjoint: + return self._phiDeriv_u(src, v, adjoint), self._phiDeriv_m(src, v, adjoint) + return np.array(self._phiDeriv_u(src, du_dm_v, adjoint) + self._phiDeriv_m(src, v, adjoint), dtype = complex) + + def _eDeriv(self, src, du_dm_v, v, adjoint=False): + if getattr(self, '_eDeriv_u', None) is None or getattr(self, '_eDeriv_m', None) is None: + raise NotImplementedError ('Getting eDerivs from %s is not implemented' %self.knownFields.keys()[0]) + + if adjoint: + return self._eDeriv_u(src, v, adjoint), self._eDeriv_m(src, v, adjoint) + return np.array(self._eDeriv_u(src, du_dm_v, adjoint) + self._eDeriv_m(src, v, adjoint), dtype = complex) + + def _jDeriv(self, src, du_dm_v, v, adjoint=False): + if getattr(self, '_jDeriv_u', None) is None or getattr(self, '_jDeriv_m', None) is None: + raise NotImplementedError ('Getting jDerivs from %s is not implemented' %self.knownFields.keys()[0]) + + if adjoint: + return self._jDeriv_u(src, v, adjoint), self._jDeriv_m(src, v, adjoint) + return np.array(self._jDeriv_u(src, du_dm_v, adjoint) + self._jDeriv_m(src, v, adjoint), dtype = complex) + + +class Fields_CC(Fields): + knownFields = {'phiSolution':'CC'} + aliasFields = { + 'phi': ['phiSolution','CC','_phi'], + 'j' : ['phiSolution','F','_j'], + 'e' : ['phiSolution','F','_e'], + } + # primary - secondary + # CC variables + + def __init__(self, mesh, survey, **kwargs): + Fields.__init__(self, mesh, survey, **kwargs) + + def startup(self): + self.prob = self.survey.prob + + def _GLoc(self, fieldType): + if fieldType == 'phi': + return 'CC' + elif fieldType == 'e' or fieldType == 'j': + return 'F' + else: + raise Exception('Field type must be phi, e, j') + + def _phi(self, phiSolution, srcList): + return phiSolution + + def _phiDeriv_u(): + return Identity() + + def _phiDeriv_m(): + return Zero() + + def _j(self, phiSolution, srcList): + raise NotImplementedError + + def _e(self, phiSolution, srcList): + raise NotImplementedError + + diff --git a/SimPEG/EM/Static/DC/ProblemDC.py b/SimPEG/EM/Static/DC/ProblemDC.py new file mode 100644 index 00000000..f8f80234 --- /dev/null +++ b/SimPEG/EM/Static/DC/ProblemDC.py @@ -0,0 +1,79 @@ +from SimPEG import Problem +from SimPEG.EM.Base import BaseEMProblem +from SurveyDC import Survey +from FieldsDC import Fields, Fields_CC + +class BaseDCProblem(BaseEMProblem): + + surveyPair = Survey + fieldsPair = Fields + + def fields(self, m): + self.curModel = m + f = self.fieldsPair(self.mesh, self.survey) + A = self.getA() + self.Ainv = self.Solver(A, **self.solverOpts) + RHS = self.getRHS() + u = self.Ainv * RHS + Srcs = self.survey.srcList + f[Srcs, self._solutionType] = u + return f + + def Jvec(self, m, v, f=None): + raise NotImplementedError + + def Jtvec(self, m, v, f=None): + raise NotImplementedError + + def getSourceTerm(self): + """ + takes concept of source and turns it into a matrix + """ + raise NotImplementedError + + + +class Problem3D_CC(BaseDCProblem): + + _solutionType = 'phiSolution' + _formulation = 'HJ' # CC potentials means J is on faces + fieldsPair = Fields_CC + + def __init__(self, mesh, **kwargs): + BaseDCProblem.__init__(self, mesh, **kwargs) + + + def getA(self): + """ + + Make the A matrix for the cell centered DC resistivity problem + + A = D MfRhoI D^\\top V + + """ + + # TODO: this won't work for full anisotropy + + D = self.mesh.faceDiv + MfRhoI = self.MfRhoI + V = self.Vol + A = D * ( MfRhoI * ( D.T * V ) ) + + if self._makeASymmetric is True: + return V.T * A + return A + + + def getRHS(self): + """ + RHS for the DC problem + + q + """ + + RHS = self.getSourceTerm() + if self._makeASymmetric is True: + return self.Vol.T * RHS + return RHS + + diff --git a/SimPEG/EM/Static/DC/RxDC.py b/SimPEG/EM/Static/DC/RxDC.py new file mode 100644 index 00000000..f9a2546d --- /dev/null +++ b/SimPEG/EM/Static/DC/RxDC.py @@ -0,0 +1,66 @@ +import SimPEG +# from SimPEG.EM.Base import BaseEMSurvey +from SimPEG.Utils import Zero, closestPoints + +class BaseRx(SimPEG.Survey.BaseRx): + loc = None + rxType = None + + knownRxTypes = { + 'phi':['phi',None], + 'ex':['e','x'], + 'ey':['e','y'], + 'ez':['e','z'], + 'jx':['j','x'], + 'jy':['j','y'], + 'jz':['j','z'], + } + + def __init__(self, **kwargs): + SimPEG.Survey.BaseRx.__init__(locs, rxType, **kwargs) + + @property + def projField(self): + """Field Type projection (e.g. e b ...)""" + return self.knownRxTypes[self.rxType][0] + + def projGLoc(self, f): + """Grid Location projection (e.g. Ex Fy ...)""" + comp = self.knownRxTypes[self.rxType][1] + if comp is not None: + return f._GLoc(self.rxType[0]) + comp + return f._GLoc(self.rxType[0]) + + def eval(self, src, mesh, f): + P = self.getP(self.prob.mesh) + return P*f[src, self.projField] + +# DC.Rx.Dipole(locs) +class Dipole(BaseRx): + + def __init__(self, locsM, locsN, rxType = 'phi', **kwargs): + assert locsM.shape == locsN.shape, 'locsM and locsN need to be the same size' + self.locs = [locsM, locsN] + BaseRx.__init__(self) + + @property + def nD(self): + """Number of data in the receiver.""" + return self.locs[0].shape[0] + + def getP(self,mesh): + if mesh in self._Ps: + return self._Ps[mesh] + + P0 = mesh.getInterpolationMat(self.locs[0], self.projGLoc) + P1 = mesh.getInterpolationMat(self.locs[1], self.projGLoc) + P = P0 - P1 + + if self.storeProjections: + self._Ps[mesh] = P + return P + + +# class Pole(BaseRx): + + diff --git a/SimPEG/EM/Static/DC/SrcDC.py b/SimPEG/EM/Static/DC/SrcDC.py new file mode 100644 index 00000000..61f490fa --- /dev/null +++ b/SimPEG/EM/Static/DC/SrcDC.py @@ -0,0 +1,60 @@ +import SimPEG +# from SimPEG.EM.Base import BaseEMSurvey +from SimPEG.Utils import Zero, closestPoints + +class BaseSrc(SimPEG.Survey.BaseSrc): + + current = 1 + loc = None + + def __init__(self, rxList, **kwargs): + SimPEG.Survey.BaseSrc.__init__(self, rxList, **kwargs) + + def eval(self, prob): + raise NotImplementedError + + def evalDeriv(self, prob): + Zero() + + +class Dipole(BaseSrc): + + def __init__(self, rxList, locA, locB, **kwargs): + assert locA.shape == locB.shape, 'Shape of locA and locB should be the same' + self.loc = [locA, locB] + BaseSrc.__init__(self, rxList, **kwargs) + + def eval(self, prob): + if prob._formulation == 'HJ': + inds = closestPoints(prob.mesh, self.loc) + q = np.zeros(prob.mesh.nC) + q[inds] = self.current * np.r_[1., -1.] + elif prob._formulation == 'EB': + # TODO: there is probably a faster way to do this + # Utils.cellNodes , Utils.cellFaces, Utils.cellEdges + raise NotImplementedError + return q + + # def bc_contribution + + +# How to treat boundary conditions here + +class Pole(BaseSrc): + + def __init__(self, rxList, loc, **kwargs): + BaseSrc.__init__(self, rxList, loc=loc, **kwargs) + + def eval(self, prob): + if prob._formulation == 'HJ': + inds = closestPoints(prob.mesh, self.loc) + q = np.zeros(prob.mesh.nC) + q[inds] = self.current * np.r_[1.] + elif prob._formulation == 'EB': + # TODO: there is probably a faster way to do this + # Utils.cellNodes , Utils.cellFaces, Utils.cellEdges + raise NotImplementedError + return q + + # def bc_contribution + diff --git a/SimPEG/EM/Static/DC/SurveyDC.py b/SimPEG/EM/Static/DC/SurveyDC.py new file mode 100644 index 00000000..3b631bef --- /dev/null +++ b/SimPEG/EM/Static/DC/SurveyDC.py @@ -0,0 +1,18 @@ +import SimPEG +from SimPEG.EM.Base import BaseEMSurvey +from SimPEG import sp +from SimPEG.Utils import Zero, Identity +from RxDC import BaseRx +from SrcDC import BaseSrc + +class Survey(BaseEMSurvey): + rxPair = BaseRx + srcPair = BaseSrc + + def __init__(self, srcList, **kwargs): + self.srcList = srcList + BaseEMSurvey.__init__(self, srcList, **kwargs) + + + + diff --git a/SimPEG/EM/Static/DC/__init__.py b/SimPEG/EM/Static/DC/__init__.py new file mode 100644 index 00000000..a3e1eba5 --- /dev/null +++ b/SimPEG/EM/Static/DC/__init__.py @@ -0,0 +1,5 @@ +from ProblemDC import Problem3D_CC +from SurveyDC import Survey +import SrcDC as Src #Pole +import RxDC as Rx +from FieldsDC import Fields_CC