From 35bac38c8bdaf49235d79adffef23b95d3d0485f Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Thu, 14 Apr 2016 22:41:47 -0700 Subject: [PATCH 01/25] working on mixed BC --- SimPEG/Mesh/DiffOperators.py | 60 ++++++++++++++++++++++++++++++++++++ 1 file changed, 60 insertions(+) diff --git a/SimPEG/Mesh/DiffOperators.py b/SimPEG/Mesh/DiffOperators.py index d0363001..b4a4b1ba 100644 --- a/SimPEG/Mesh/DiffOperators.py +++ b/SimPEG/Mesh/DiffOperators.py @@ -584,7 +584,67 @@ class DiffOperators(object): return Pbc, Pin, Pout + def getBCProjWF_simple(self, discretization='CC'): + """ + The weak form boundary condition projection matrices + when mixed boundary condition is used + + + """ + + if discretization is not 'CC': + raise NotImplementedError('Boundary conditions only implemented for CC discretization.') + + def projBC(n): + ij = ([0,n], [0,1]) + vals = [0,0] + vals[0] = 1 + vals[1] = 1 + return sp.csr_matrix((vals, ij), shape=(n+1,2)) + + def projDirichlet(n, bc): + bc = checkBC(bc) + ij = ([0,n], [0,1]) + vals = [0,0] + if(bc[0] == 'dirichlet'): + vals[0] = -1 + if(bc[1] == 'dirichlet'): + vals[1] = 1 + return sp.csr_matrix((vals, ij), shape=(n+1,2)) + + BC = [['dirichlet','dirichlet'],['dirichlet','dirichlet'],['dirichlet','dirichlet']] + n = self.vnC + indF = self.faceBoundaryInd + if(self.dim == 1): + Pbc = projDirichlet(n[0], BC[0]) + B = projBC(n[0]) + indF = indF[0] | indF[1] + Pbc = Pbc*sdiag(self.area[indF]) + + elif(self.dim == 2): + Pbc1 = sp.kron(speye(n[1]), projDirichlet(n[0], BC[0])) + Pbc2 = sp.kron(projDirichlet(n[1], BC[1]), speye(n[0])) + Pbc = sp.block_diag((Pbc1, Pbc2), format="csr") + B1 = sp.kron(speye(n[1]), projBC(n[0])) + B2 = sp.kron(projBC(n[1]), speye(n[0])) + B = sp.block_diag((B1, B2), format="csr") + indF = np.r_[(indF[0] | indF[1]), (indF[2] | indF[3])] + Pbc = Pbc*sdiag(self.area[indF]) + + elif(self.dim == 3): + Pbc1 = kron3(speye(n[2]), speye(n[1]), projDirichlet(n[0], BC[0])) + Pbc2 = kron3(speye(n[2]), projDirichlet(n[1], BC[1]), speye(n[0])) + Pbc3 = kron3(projDirichlet(n[2], BC[2]), speye(n[1]), speye(n[0])) + Pbc = sp.block_diag((Pbc1, Pbc2, Pbc3), format="csr") + B1 = kron3(speye(n[2]), speye(n[1]), projBC(n[0])) + B2 = kron3(speye(n[2]), projBC(n[1]), speye(n[0])) + B3 = kron3(projBC(n[2]), speye(n[1]), speye(n[0])) + B = sp.block_diag((B1, B2, B3), format="csr") + indF = np.r_[(indF[0] | indF[1]), (indF[2] | indF[3]), (indF[4] | indF[5])] + Pbc = Pbc*sdiag(self.area[indF]) + + return Pbc, B.T # --------------- Averaging --------------------- @property From 119bc801c7d880e1bd5223e6b31a3801cf866478 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 15 Apr 2016 14:51:19 -0700 Subject: [PATCH 02/25] 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 From 8739ba0f2072daabfe7e1a820beffa64ef00cb62 Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Fri, 15 Apr 2016 16:53:43 -0700 Subject: [PATCH 03/25] working dc fwd --- SimPEG/EM/FDEM/FDEM.py | 1 + SimPEG/EM/Static/DC/ProblemDC.py | 28 +++++++++++++++++++++++++++- SimPEG/EM/Static/DC/RxDC.py | 25 ++++++++++++++----------- SimPEG/EM/Static/DC/SrcDC.py | 1 + SimPEG/EM/Static/__init__.py | 1 + SimPEG/EM/TDEM/SurveyTDEM.py | 12 ++++++------ SimPEG/EM/__init__.py | 1 + 7 files changed, 51 insertions(+), 18 deletions(-) create mode 100644 SimPEG/EM/Static/__init__.py diff --git a/SimPEG/EM/FDEM/FDEM.py b/SimPEG/EM/FDEM/FDEM.py index caca7602..1ac494b9 100644 --- a/SimPEG/EM/FDEM/FDEM.py +++ b/SimPEG/EM/FDEM/FDEM.py @@ -167,6 +167,7 @@ class BaseFDEMProblem(BaseEMProblem): for i, src in enumerate(Srcs): smi, sei = src.eval(self) + #Why are you adding? s_m[:,i] = s_m[:,i] + smi s_e[:,i] = s_e[:,i] + sei diff --git a/SimPEG/EM/Static/DC/ProblemDC.py b/SimPEG/EM/Static/DC/ProblemDC.py index f8f80234..f712e036 100644 --- a/SimPEG/EM/Static/DC/ProblemDC.py +++ b/SimPEG/EM/Static/DC/ProblemDC.py @@ -2,6 +2,7 @@ from SimPEG import Problem from SimPEG.EM.Base import BaseEMProblem from SurveyDC import Survey from FieldsDC import Fields, Fields_CC +import numpy as np class BaseDCProblem(BaseEMProblem): @@ -29,9 +30,28 @@ class BaseDCProblem(BaseEMProblem): """ takes concept of source and turns it into a matrix """ - raise NotImplementedError + """ + Evaluates the sources for a given frequency and puts them in matrix form + :param float freq: Frequency + :rtype: (numpy.ndarray, numpy.ndarray) + :return: s_m, s_e (nE or nF, nSrc) + """ + Srcs = self.survey.srcList + + if self._formulation is 'EB': + n = self.mesh.nN + # return NotImplementedError + + elif self._formulation is 'HJ': + n = self.mesh.nC + + q = np.zeros((n, len(Srcs))) + + for i, src in enumerate(Srcs): + q[:,i] = src.eval(self) + return q class Problem3D_CC(BaseDCProblem): @@ -63,6 +83,8 @@ class Problem3D_CC(BaseDCProblem): return V.T * A return A + def getADeriv(): + raise NotImplementedError def getRHS(self): """ @@ -76,4 +98,8 @@ class Problem3D_CC(BaseDCProblem): return self.Vol.T * RHS return RHS + def getRHSDeriv(): + raise NotImplementedError + + diff --git a/SimPEG/EM/Static/DC/RxDC.py b/SimPEG/EM/Static/DC/RxDC.py index f9a2546d..9fbd05cb 100644 --- a/SimPEG/EM/Static/DC/RxDC.py +++ b/SimPEG/EM/Static/DC/RxDC.py @@ -3,7 +3,7 @@ import SimPEG from SimPEG.Utils import Zero, closestPoints class BaseRx(SimPEG.Survey.BaseRx): - loc = None + locs = None rxType = None knownRxTypes = { @@ -16,8 +16,9 @@ class BaseRx(SimPEG.Survey.BaseRx): 'jz':['j','z'], } - def __init__(self, **kwargs): - SimPEG.Survey.BaseRx.__init__(locs, rxType, **kwargs) + def __init__(self, locs, rxType, **kwargs): + SimPEG.Survey.BaseRx.__init__(self, locs, rxType, **kwargs) + @property def projField(self): @@ -28,11 +29,11 @@ class BaseRx(SimPEG.Survey.BaseRx): """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]) + return f._GLoc(self.rxType) + comp + return f._GLoc(self.rxType) def eval(self, src, mesh, f): - P = self.getP(self.prob.mesh) + P = self.getP(mesh, self.projGLoc(f)) return P*f[src, self.projField] # DC.Rx.Dipole(locs) @@ -40,24 +41,26 @@ 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) + locs = [locsM, locsN] + # We may not need this ... + BaseRx.__init__(self, locs, rxType) @property def nD(self): """Number of data in the receiver.""" return self.locs[0].shape[0] - def getP(self,mesh): + def getP(self, mesh, Gloc): 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) + P0 = mesh.getInterpolationMat(self.locs[0], Gloc) + P1 = mesh.getInterpolationMat(self.locs[1], Gloc) P = P0 - P1 if self.storeProjections: self._Ps[mesh] = P + return P diff --git a/SimPEG/EM/Static/DC/SrcDC.py b/SimPEG/EM/Static/DC/SrcDC.py index 61f490fa..cc855118 100644 --- a/SimPEG/EM/Static/DC/SrcDC.py +++ b/SimPEG/EM/Static/DC/SrcDC.py @@ -1,6 +1,7 @@ import SimPEG # from SimPEG.EM.Base import BaseEMSurvey from SimPEG.Utils import Zero, closestPoints +import numpy as np class BaseSrc(SimPEG.Survey.BaseSrc): diff --git a/SimPEG/EM/Static/__init__.py b/SimPEG/EM/Static/__init__.py new file mode 100644 index 00000000..6ebc9df2 --- /dev/null +++ b/SimPEG/EM/Static/__init__.py @@ -0,0 +1 @@ +import DC diff --git a/SimPEG/EM/TDEM/SurveyTDEM.py b/SimPEG/EM/TDEM/SurveyTDEM.py index 4d04f0ae..9c9ed4b9 100644 --- a/SimPEG/EM/TDEM/SurveyTDEM.py +++ b/SimPEG/EM/TDEM/SurveyTDEM.py @@ -87,7 +87,7 @@ class SrcTDEM_VMD_MVP(SrcTDEM): def getInitialFields(self, mesh): """Vertical magnetic dipole, magnetic vector potential""" if self.waveformType == "STEPOFF": - print ">> Step waveform: Non-zero initial condition" + print ">> Step waveform: Non-zero initial condition" if mesh._meshType is 'CYL': if mesh.isSymmetric: MVP = MagneticDipoleVectorPotential(self.loc, mesh, 'Ey') @@ -96,8 +96,8 @@ class SrcTDEM_VMD_MVP(SrcTDEM): elif mesh._meshType is 'TENSOR': MVP = MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez']) else: - raise Exception('Unknown mesh for VMD') - return {"b": mesh.edgeCurl*MVP} + raise Exception('Unknown mesh for VMD') + return {"b": mesh.edgeCurl*MVP} elif self.waveformType == "GENERAL": print ">> General waveform: Zero initial condition" return {"b": np.zeros(mesh.nF)} @@ -113,7 +113,7 @@ class SrcTDEM_VMD_MVP(SrcTDEM): elif mesh._meshType is 'TENSOR': MVP = MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez']) else: - raise Exception('Unknown mesh for VMD') + raise Exception('Unknown mesh for VMD') return mesh.edgeCurl.T*MfMui*mesh.edgeCurl*MVP @@ -122,7 +122,7 @@ class SrcTDEM_CircularLoop_MVP(SrcTDEM): self.loc = loc self.radius = radius self.waveformType = waveformType - SrcTDEM.__init__(self,rxList) + SrcTDEM.__init__(self,rxList) def getInitialFields(self, mesh): """Circular Loop, magnetic vector potential""" @@ -153,7 +153,7 @@ class SrcTDEM_CircularLoop_MVP(SrcTDEM): elif mesh._meshType is 'TENSOR': MVP = MagneticLoopVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'], self.radius) else: - raise Exception('Unknown mesh for CircularLoop') + raise Exception('Unknown mesh for CircularLoop') return mesh.edgeCurl.T*MfMui*mesh.edgeCurl*MVP diff --git a/SimPEG/EM/__init__.py b/SimPEG/EM/__init__.py index 565f63a8..e42afbc5 100644 --- a/SimPEG/EM/__init__.py +++ b/SimPEG/EM/__init__.py @@ -1,5 +1,6 @@ import TDEM import FDEM +import Static import Base import Analytics import Utils From 28005dde45b4469fb6c9202c22b3b7b85dba37d6 Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Fri, 15 Apr 2016 17:17:18 -0700 Subject: [PATCH 04/25] change minor bug for meshIO --- SimPEG/Mesh/MeshIO.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SimPEG/Mesh/MeshIO.py b/SimPEG/Mesh/MeshIO.py index 7501a66f..16d1b457 100644 --- a/SimPEG/Mesh/MeshIO.py +++ b/SimPEG/Mesh/MeshIO.py @@ -21,7 +21,7 @@ class TensorMeshIO(object): if '*' in seg: st = seg sp = seg.split('*') - re = np.array(sp[0],dtype=int)*(' ' + sp[1]) + re = int(sp[0])*(' ' + sp[1]) line = line.replace(st,re.strip()) return np.array(line.split(),dtype=float) From cd5339322eb545c0ef0759e56e3de756a86506ea Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Thu, 21 Apr 2016 11:00:20 -0700 Subject: [PATCH 05/25] sketch of DC --- SimPEG/EM/Static/DC/ProblemDC.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/SimPEG/EM/Static/DC/ProblemDC.py b/SimPEG/EM/Static/DC/ProblemDC.py index f712e036..cc9c1ebb 100644 --- a/SimPEG/EM/Static/DC/ProblemDC.py +++ b/SimPEG/EM/Static/DC/ProblemDC.py @@ -83,8 +83,8 @@ class Problem3D_CC(BaseDCProblem): return V.T * A return A - def getADeriv(): - raise NotImplementedError + def getADeriv(self, u, v, adjoint= False): + def getRHS(self): """ From 2c09be9fc12a8046feadcc240f25b74dd2d439ae Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Thu, 21 Apr 2016 14:44:37 -0700 Subject: [PATCH 06/25] Working Mixed boundary conditions and testing ... --- tests/mesh/test_Mixed_boundaryPoisson.py | 395 +++++++++++++++++++++++ 1 file changed, 395 insertions(+) create mode 100644 tests/mesh/test_Mixed_boundaryPoisson.py diff --git a/tests/mesh/test_Mixed_boundaryPoisson.py b/tests/mesh/test_Mixed_boundaryPoisson.py new file mode 100644 index 00000000..7a15f36d --- /dev/null +++ b/tests/mesh/test_Mixed_boundaryPoisson.py @@ -0,0 +1,395 @@ +import numpy as np +import scipy.sparse as sp +import unittest +import matplotlib.pyplot as plt +from SimPEG import * + +MESHTYPES = ['uniformTensorMesh'] + +def getxBCyBC(mesh, alpha, beta, gamma): +# def getxBCyBC(mesh, alpha, beta, gamma): + """ + + """ + if mesh.dim == 1: #1D + if (len(alpha) != 2 or len(beta) != 2 or len(gamma) != 2): + raise Exception("Lenght of list, alpha should be 2") + fCCxm,fCCxp = mesh.cellBoundaryInd + nBC = fCCxm.sum()+fCCxp.sum() + h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp] + + alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0] + alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1] + + h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp] + + a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm) + b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm) + a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp) + b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp) + + xBC_xm = 0.5*a_xm + xBC_xp = 0.5*a_xp/b_xp + yBC_xm = 0.5*(1.-b_xm) + yBC_xp = 0.5*(1.-1./b_xp) + + xBC = np.r_[xBC_xm, xBC_xp] + yBC = np.r_[yBC_xm, yBC_xp] + + elif mesh.dim == 2: #2D + if (len(alpha) != 4 or len(beta) != 4 or len(gamma) != 4): + raise Exception("Lenght of list, alpha should be 4") + + fCCxm,fCCxp,fCCym,fCCyp = mesh.cellBoundaryInd + fxm,fxp,fym,fyp = mesh.faceBoundaryInd + nBC = fCCxm.sum()+fCCxp.sum()+fCCxm.sum()+fCCxp.sum() + h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp] + h_ym, h_yp = mesh.gridCC[fCCym], mesh.gridCC[fCCyp] + + alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0] + alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1] + alpha_ym, beta_ym, gamma_ym = alpha[2], beta[2], gamma[2] + alpha_yp, beta_yp, gamma_yp = alpha[3], beta[3], gamma[3] + + h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0] + h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1] + + a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm) + b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm) + a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp) + b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp) + + a_ym = gamma_ym/(0.5*alpha_ym-beta_ym/h_ym) + b_ym = (0.5*alpha_ym+beta_ym/h_ym)/(0.5*alpha_ym-beta_ym/h_ym) + a_yp = gamma_yp/(0.5*alpha_yp-beta_yp/h_yp) + b_yp = (0.5*alpha_yp+beta_yp/h_yp)/(0.5*alpha_yp-beta_yp/h_yp) + + xBC_xm = 0.5*a_xm + xBC_xp = 0.5*a_xp/b_xp + yBC_xm = 0.5*(1.-b_xm) + yBC_xp = 0.5*(1.-1./b_xp) + xBC_ym = 0.5*a_ym + xBC_yp = 0.5*a_yp/b_yp + yBC_ym = 0.5*(1.-b_ym) + yBC_yp = 0.5*(1.-1./b_yp) + + sortindsfx = np.argsort(np.r_[np.arange(mesh.nFx)[fxm], np.arange(mesh.nFx)[fxp]]) + sortindsfy = np.argsort(np.r_[np.arange(mesh.nFy)[fym], np.arange(mesh.nFy)[fyp]]) + + xBC_x = np.r_[xBC_xm, xBC_xp][sortindsfx] + xBC_y = np.r_[xBC_ym, xBC_yp][sortindsfy] + yBC_x = np.r_[yBC_xm, yBC_xp][sortindsfx] + yBC_y = np.r_[yBC_ym, yBC_yp][sortindsfy] + + xBC = np.r_[xBC_x, xBC_y] + yBC = np.r_[yBC_x, yBC_y] + + elif mesh.dim == 3: #3D + if (len(alpha) != 6 or len(beta) != 6 or len(gamma) != 6): + raise Exception("Lenght of list, alpha should be 6") + fCCxm,fCCxp,fCCym,fCCyp,fCCzm,fCCzp = mesh.cellBoundaryInd + fxm,fxp,fym,fyp,fzm,fzp = mesh.faceBoundaryInd + nBC = fCCxm.sum()+fCCxp.sum()+fCCxm.sum()+fCCxp.sum() + h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp] + h_ym, h_yp = mesh.gridCC[fCCym], mesh.gridCC[fCCyp] + h_zm, h_zp = mesh.gridCC[fCCzm], mesh.gridCC[fCCzp] + + alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0] + alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1] + alpha_ym, beta_ym, gamma_ym = alpha[2], beta[2], gamma[2] + alpha_yp, beta_yp, gamma_yp = alpha[3], beta[3], gamma[3] + alpha_zm, beta_zm, gamma_zm = alpha[2], beta[2], gamma[2] + alpha_zp, beta_zp, gamma_zp = alpha[3], beta[3], gamma[3] + + h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0] + h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1] + h_zm, h_zp = mesh.gridCC[fCCzm,2], mesh.gridCC[fCCzp,2] + + a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm) + b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm) + a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp) + b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp) + + a_ym = gamma_ym/(0.5*alpha_ym-beta_ym/h_ym) + b_ym = (0.5*alpha_ym+beta_ym/h_ym)/(0.5*alpha_ym-beta_ym/h_ym) + a_yp = gamma_yp/(0.5*alpha_yp-beta_yp/h_yp) + b_yp = (0.5*alpha_yp+beta_yp/h_yp)/(0.5*alpha_yp-beta_yp/h_yp) + + a_zm = gamma_zm/(0.5*alpha_zm-beta_zm/h_zm) + b_zm = (0.5*alpha_zm+beta_zm/h_zm)/(0.5*alpha_zm-beta_zm/h_zm) + a_zp = gamma_zp/(0.5*alpha_zp-beta_zp/h_zp) + b_zp = (0.5*alpha_zp+beta_zp/h_zp)/(0.5*alpha_zp-beta_zp/h_zp) + + xBC_xm = 0.5*a_xm + xBC_xp = 0.5*a_xp/b_xp + yBC_xm = 0.5*(1.-b_xm) + yBC_xp = 0.5*(1.-1./b_xp) + xBC_ym = 0.5*a_ym + xBC_yp = 0.5*a_yp/b_yp + yBC_ym = 0.5*(1.-b_ym) + yBC_yp = 0.5*(1.-1./b_yp) + xBC_zm = 0.5*a_zm + xBC_zp = 0.5*a_zp/b_zp + yBC_zm = 0.5*(1.-b_zm) + yBC_zp = 0.5*(1.-1./b_zp) + + sortindsfx = np.argsort(np.r_[np.arange(mesh.nFx)[fxm], np.arange(mesh.nFx)[fxp]]) + sortindsfy = np.argsort(np.r_[np.arange(mesh.nFy)[fym], np.arange(mesh.nFy)[fyp]]) + sortindsfz = np.argsort(np.r_[np.arange(mesh.nFz)[fzm], np.arange(mesh.nFz)[fzp]]) + + xBC_x = np.r_[xBC_xm, xBC_xp][sortindsfx] + xBC_y = np.r_[xBC_ym, xBC_yp][sortindsfy] + xBC_z = np.r_[xBC_zm, xBC_zp][sortindsfz] + + yBC_x = np.r_[yBC_xm, yBC_xp][sortindsfx] + yBC_y = np.r_[yBC_ym, yBC_yp][sortindsfy] + yBC_z = np.r_[yBC_zm, yBC_zp][sortindsfz] + + xBC = np.r_[xBC_x, xBC_y, xBC_z] + yBC = np.r_[yBC_x, yBC_y, yBC_z] + + return xBC, yBC + +class Test1D_InhomogeneousMixed(Tests.OrderTest): + name = "1D - Mixed" + meshTypes = MESHTYPES + meshDimension = 1 + expectedOrders = 2 + meshSizes = [4, 8, 16, 32] + + def getError(self): + #Test function + phi_fun = lambda x: np.cos(np.pi*x) + j_fun = lambda x: np.pi*np.sin(np.pi*x) + phi_deriv = lambda x: -j_fun(x) + q_fun = lambda x: (np.pi**2)*np.cos(np.pi*x) + + xc_ana = phi_fun(self.M.gridCC) + q_ana = q_fun(self.M.gridCC) + j_ana = j_fun(self.M.gridFx) + + # Get boundary locations + vecN = self.M.vectorNx + vecC = self.M.vectorCCx + + # Setup Mixed B.C (alpha, beta, gamma) + alpha_xm, alpha_xp = 1., 1. + beta_xm, beta_xp = 1., 1. + alpha = np.r_[alpha_xm, alpha_xp] + beta = np.r_[beta_xm, beta_xp] + vecN = self.M.vectorNx + vecC = self.M.vectorCCx + phi_bc = phi_fun(vecN[[0,-1]]) + phi_deriv_bc = phi_deriv(vecN[[0,-1]]) + gamma = alpha*phi_bc + beta*phi_deriv_bc + x_BC, y_BC = getxBCyBC(self.M, alpha, beta, gamma) + + + sigma = np.ones(self.M.nC) + Mfrho = self.M.getFaceInnerProduct(1./sigma) + MfrhoI = self.M.getFaceInnerProduct(1./sigma, invMat=True) + V = Utils.sdiag(self.M.vol) + Div = V*self.M.faceDiv + P_BC, B = self.M.getBCProjWF_simple() + q = q_fun(self.M.gridCC) + M = B*self.M.aveCC2F + G = Div.T - P_BC*Utils.sdiag(y_BC)*M + rhs = V*q + Div*MfrhoI*P_BC*x_BC + A = Div*MfrhoI*G + + if self.myTest == 'xc': + #TODO: fix the null space + Ainv = Solver(A) + xc = Ainv*rhs + err = np.linalg.norm((xc-xc_ana), np.inf) + else: + NotImplementedError + return err + + + def test_order(self): + print "==== Testing Mixed boudary conduction for CC-problem ====" + self.name = "1D" + self.myTest = 'xc' + self.orderTest() + +class Test2D_InhomogeneousMixed(Tests.OrderTest): + name = "2D - Mixed" + meshTypes = MESHTYPES + meshDimension = 2 + expectedOrders = 2 + meshSizes = [4, 8, 16, 32] + + def getError(self): + #Test function + phi_fun = lambda x: np.cos(np.pi*x[:,0])*np.cos(np.pi*x[:,1]) + j_funX = lambda x: +np.pi*np.sin(np.pi*x[:,0])*np.cos(np.pi*x[:,1]) + j_funY = lambda x: +np.pi*np.cos(np.pi*x[:,0])*np.sin(np.pi*x[:,1]) + phideriv_funX = lambda x: -j_funX(x) + phideriv_funY = lambda x: -j_funY(x) + q_fun = lambda x: +2*(np.pi**2)*phi_fun(x) + + xc_ana = phi_fun(self.M.gridCC) + q_ana = q_fun(self.M.gridCC) + jX_ana = j_funX(self.M.gridFx) + jY_ana = j_funY(self.M.gridFy) + j_ana = np.r_[jX_ana,jY_ana] + + # Get boundary locations + fxm,fxp,fym,fyp = self.M.faceBoundaryInd + gBFxm = self.M.gridFx[fxm,:] + gBFxp = self.M.gridFx[fxp,:] + gBFym = self.M.gridFy[fym,:] + gBFyp = self.M.gridFy[fyp,:] + + # Setup Mixed B.C (alpha, beta, gamma) + alpha_xm, alpha_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0]) + beta_xm, beta_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0]) + alpha_ym, alpha_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1]) + beta_ym, beta_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1]) + + phi_bc_xm, phi_bc_xp = phi_fun(gBFxm), phi_fun(gBFxp) + phi_bc_ym, phi_bc_yp = phi_fun(gBFym), phi_fun(gBFyp) + + phiderivX_bc_xm, phiderivX_bc_xp = phideriv_funX(gBFxm), phideriv_funX(gBFxp) + phiderivY_bc_ym, phiderivY_bc_yp = phideriv_funY(gBFym), phideriv_funY(gBFyp) + + gamma_fun = lambda alpha, beta, phi, phi_deriv: alpha*phi + beta*phi_deriv + gamma_xm = gamma_fun(alpha_xm, beta_xm, phi_bc_xm, phiderivX_bc_xm) + gamma_xp = gamma_fun(alpha_xp, beta_xp, phi_bc_xp, phiderivX_bc_xp) + gamma_ym = gamma_fun(alpha_ym, beta_ym, phi_bc_ym, phiderivY_bc_ym) + gamma_yp = gamma_fun(alpha_yp, beta_yp, phi_bc_yp, phiderivY_bc_yp) + + alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp] + beta = [beta_xm, beta_xp, beta_ym, beta_yp] + gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp] + + x_BC, y_BC = getxBCyBC(self.M, alpha, beta, gamma) + + + sigma = np.ones(self.M.nC) + Mfrho = self.M.getFaceInnerProduct(1./sigma) + MfrhoI = self.M.getFaceInnerProduct(1./sigma, invMat=True) + V = Utils.sdiag(self.M.vol) + Div = V*self.M.faceDiv + P_BC, B = self.M.getBCProjWF_simple() + q = q_fun(self.M.gridCC) + M = B*self.M.aveCC2F + G = Div.T - P_BC*Utils.sdiag(y_BC)*M + rhs = V*q + Div*MfrhoI*P_BC*x_BC + A = Div*MfrhoI*G + + if self.myTest == 'xc': + Ainv = Solver(A) + xc = Ainv*rhs + err = np.linalg.norm((xc-xc_ana), np.inf) + else: + NotImplementedError + return err + + + def test_order(self): + print "==== Testing Mixed boudary conduction for CC-problem ====" + self.name = "2D" + self.myTest = 'xc' + self.orderTest() + +class Test3D_InhomogeneousMixed(Tests.OrderTest): + name = "3D - Mixed" + meshTypes = MESHTYPES + meshDimension = 3 + expectedOrders = 2 + meshSizes = [4, 8, 16] + + def getError(self): + #Test function + phi_fun = lambda x: np.cos(np.pi*x[:,0])*np.cos(np.pi*x[:,1])*np.cos(np.pi*x[:,2]) + j_funX = lambda x: +np.pi*np.sin(np.pi*x[:,0])*np.cos(np.pi*x[:,1])*np.cos(np.pi*x[:,2]) + j_funY = lambda x: +np.pi*np.cos(np.pi*x[:,0])*np.sin(np.pi*x[:,1])*np.cos(np.pi*x[:,2]) + j_funZ = lambda x: +np.pi*np.cos(np.pi*x[:,0])*np.cos(np.pi*x[:,1])*np.sin(np.pi*x[:,2]) + + phideriv_funX = lambda x: -j_funX(x) + phideriv_funY = lambda x: -j_funY(x) + phideriv_funZ = lambda x: -j_funZ(x) + + q_fun = lambda x: 3*(np.pi**2)*phi_fun(x) + + xc_ana = phi_fun(self.M.gridCC) + q_ana = q_fun(self.M.gridCC) + jX_ana = j_funX(self.M.gridFx) + jY_ana = j_funY(self.M.gridFy) + j_ana = np.r_[jX_ana,jY_ana,jY_ana] + + # Get boundary locations + fxm,fxp,fym,fyp,fzm,fzp = self.M.faceBoundaryInd + gBFxm = self.M.gridFx[fxm,:] + gBFxp = self.M.gridFx[fxp,:] + gBFym = self.M.gridFy[fym,:] + gBFyp = self.M.gridFy[fyp,:] + gBFzm = self.M.gridFz[fzm,:] + gBFzp = self.M.gridFz[fzp,:] + + # Setup Mixed B.C (alpha, beta, gamma) + alpha_xm, alpha_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0]) + beta_xm, beta_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0]) + alpha_ym, alpha_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1]) + beta_ym, beta_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1]) + alpha_zm, alpha_zp = np.ones_like(gBFzm[:,1]), np.ones_like(gBFzp[:,1]) + beta_zm, beta_zp = np.ones_like(gBFzm[:,1]), np.ones_like(gBFzp[:,1]) + + + phi_bc_xm, phi_bc_xp = phi_fun(gBFxm), phi_fun(gBFxp) + phi_bc_ym, phi_bc_yp = phi_fun(gBFym), phi_fun(gBFyp) + phi_bc_zm, phi_bc_zp = phi_fun(gBFzm), phi_fun(gBFzp) + + phiderivX_bc_xm, phiderivX_bc_xp = phideriv_funX(gBFxm), phideriv_funX(gBFxp) + phiderivY_bc_ym, phiderivY_bc_yp = phideriv_funY(gBFym), phideriv_funY(gBFyp) + phiderivY_bc_zm, phiderivY_bc_zp = phideriv_funY(gBFzm), phideriv_funY(gBFzp) + + gamma_fun = lambda alpha, beta, phi, phi_deriv: alpha*phi + beta*phi_deriv + gamma_xm = gamma_fun(alpha_xm, beta_xm, phi_bc_xm, phiderivX_bc_xm) + gamma_xp = gamma_fun(alpha_xp, beta_xp, phi_bc_xp, phiderivX_bc_xp) + gamma_ym = gamma_fun(alpha_ym, beta_ym, phi_bc_ym, phiderivY_bc_ym) + gamma_yp = gamma_fun(alpha_yp, beta_yp, phi_bc_yp, phiderivY_bc_yp) + gamma_zm = gamma_fun(alpha_zm, beta_zm, phi_bc_zm, phiderivY_bc_zm) + gamma_zp = gamma_fun(alpha_zp, beta_zp, phi_bc_zp, phiderivY_bc_zp) + + alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp, alpha_zm, alpha_zp] + beta = [beta_xm, beta_xp, beta_ym, beta_yp, beta_zm, beta_zp] + gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp, gamma_zm, gamma_zp] + + x_BC, y_BC = getxBCyBC(self.M, alpha, beta, gamma) + + + sigma = np.ones(self.M.nC) + Mfrho = self.M.getFaceInnerProduct(1./sigma) + MfrhoI = self.M.getFaceInnerProduct(1./sigma, invMat=True) + V = Utils.sdiag(self.M.vol) + Div = V*self.M.faceDiv + P_BC, B = self.M.getBCProjWF_simple() + q = q_fun(self.M.gridCC) + M = B*self.M.aveCC2F + G = Div.T - P_BC*Utils.sdiag(y_BC)*M + rhs = V*q + Div*MfrhoI*P_BC*x_BC + A = Div*MfrhoI*G + + if self.myTest == 'xc': + #TODO: fix the null space + Ainv = Solver(A) + xc = Ainv*rhs + err = np.linalg.norm((xc-xc_ana), np.inf) + else: + NotImplementedError + return err + + + def test_order(self): + print "==== Testing Mixed boudary conduction for CC-problem ====" + self.name = "3D" + self.myTest = 'xc' + self.orderTest() + + + +if __name__ == '__main__': + unittest.main() From 5ec6e79a39ab85549f71f03394411953e41df8a4 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 22 Apr 2016 17:47:36 -0700 Subject: [PATCH 07/25] MfRhoIDeriv --- SimPEG/EM/Base.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/SimPEG/EM/Base.py b/SimPEG/EM/Base.py index dd349b43..0df18cdb 100644 --- a/SimPEG/EM/Base.py +++ b/SimPEG/EM/Base.py @@ -190,7 +190,13 @@ class BaseEMProblem(Problem.BaseProblem): """ Derivative of :code:`MfRhoI` with respect to the model. """ - return self.mesh.getFaceInnerProductDeriv(self.curModel.rho, invMat=True)(u) * self.curModel.rhoDeriv + + dMfRhoI_dI = -self.MfRhoI**2 + dMf_drho = self.mesh.getEdgeInnerProductDeriv(self.curModel.rho)(u) + drho_dm = self.curModel.rhoDeriv + return dMfRhoI_dI * ( dMf_drho * ( drho_dm)) + + # return self.mesh.getFaceInnerProductDeriv(self.curModel.rho, invMat=True)(u) * self.curModel.rhoDeriv class BaseEMSurvey(Survey.BaseSurvey): From 8775364d8f4b56d92e801c2df99cca4614a3a74a Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 22 Apr 2016 17:48:26 -0700 Subject: [PATCH 08/25] start of the sketch of Jvec (not to be trusted yet!) --- SimPEG/EM/Static/DC/ProblemDC.py | 42 +++++++++++++++-- SimPEG/EM/Static/DC/RxDC.py | 2 +- SimPEG/EM/Static/DC/Utils.py | 38 ++++++++++++++++ SimPEG/EM/Static/DC/__init__.py | 1 + tests/em/static/__init__.py | 12 +++++ tests/em/static/test_DC.py | 77 ++++++++++++++++++++++++++++++++ tests/em/static/test_DC_deriv.py | 0 7 files changed, 168 insertions(+), 4 deletions(-) create mode 100644 SimPEG/EM/Static/DC/Utils.py create mode 100644 tests/em/static/__init__.py create mode 100644 tests/em/static/test_DC.py create mode 100644 tests/em/static/test_DC_deriv.py diff --git a/SimPEG/EM/Static/DC/ProblemDC.py b/SimPEG/EM/Static/DC/ProblemDC.py index cc9c1ebb..602ce3d8 100644 --- a/SimPEG/EM/Static/DC/ProblemDC.py +++ b/SimPEG/EM/Static/DC/ProblemDC.py @@ -3,6 +3,7 @@ from SimPEG.EM.Base import BaseEMProblem from SurveyDC import Survey from FieldsDC import Fields, Fields_CC import numpy as np +from SimPEG.Utils import Zero class BaseDCProblem(BaseEMProblem): @@ -21,7 +22,30 @@ class BaseDCProblem(BaseEMProblem): return f def Jvec(self, m, v, f=None): - raise NotImplementedError + + if f is None: + f = self.fields(m) + + self.curModel = m + + Jv = self.dataPair(self.survey) #same size as the data + + A = self.getA() + Ainv = self.Solver(A, **self.solverOpts) + + for src in self.survey.srcList: + u_src = f[src, self._solutionType] # solution vector + dA_dm_v = self.getADeriv(u_src, v) + dRHS_dm_v = self.getRHSDeriv(src, v) + print type(dA_dm_v + dRHS_dm_v), (dA_dm_v + dRHS_dm_v).shape + du_dm_v = Ainv * ( - dA_dm_v + dRHS_dm_v ) + + for rx in src.rxList: + df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None) + df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False) + Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v) + Ainv.clean() + return Utils.mkvc(Jv) def Jtvec(self, m, v, f=None): raise NotImplementedError @@ -85,6 +109,18 @@ class Problem3D_CC(BaseDCProblem): def getADeriv(self, u, v, adjoint= False): + D = self.mesh.faceDiv + MfRhoIDeriv = self.MfRhoIDeriv + V = self.Vol + + if adjoint: + if self._makeASymmetric is True: + v = V * v + return V.T * ( D * ( MfRhoIDeriv(D * v) ) ) + + if self._makeASymmetric is True: + return V.T * ( D * ( MfRhoIDeriv( * D.T * ( V * u ) ) * v ) ) + return D * ( MfRhoIDeriv( D.T * ( V * v ) ) ) def getRHS(self): """ @@ -98,8 +134,8 @@ class Problem3D_CC(BaseDCProblem): return self.Vol.T * RHS return RHS - def getRHSDeriv(): - raise NotImplementedError + def getRHSDeriv(self, src, v, adjoint=False): + return Zero() diff --git a/SimPEG/EM/Static/DC/RxDC.py b/SimPEG/EM/Static/DC/RxDC.py index 9fbd05cb..779ffde0 100644 --- a/SimPEG/EM/Static/DC/RxDC.py +++ b/SimPEG/EM/Static/DC/RxDC.py @@ -48,7 +48,7 @@ class Dipole(BaseRx): @property def nD(self): """Number of data in the receiver.""" - return self.locs[0].shape[0] + return int(self.locs[0].size / 2) def getP(self, mesh, Gloc): if mesh in self._Ps: diff --git a/SimPEG/EM/Static/DC/Utils.py b/SimPEG/EM/Static/DC/Utils.py new file mode 100644 index 00000000..647590cc --- /dev/null +++ b/SimPEG/EM/Static/DC/Utils.py @@ -0,0 +1,38 @@ +import numpy as np + +def WennerSrcList(nElecs, aSpacing, in2D=False, plotIt=False): + + import SimPEG.EM.Static.DC as DC + + elocs = np.arange(0,aSpacing*nElecs,aSpacing) + elocs -= (nElecs*aSpacing - aSpacing)/2 + space = 1 + WENNER = np.zeros((0,),dtype=int) + for ii in range(nElecs): + for jj in range(nElecs): + test = np.r_[jj,jj+space,jj+space*2,jj+space*3] + if np.any(test >= nElecs): + break + WENNER = np.r_[WENNER, test] + space += 1 + WENNER = WENNER.reshape((-1,4)) + + + if plotIt: + for i, s in enumerate('rbkg'): + plt.plot(elocs[WENNER[:,i]],s+'.') + plt.show() + + # Create sources and receivers + i = 0 + if in2D: + getLoc = lambda ii, abmn: np.r_[elocs[WENNER[ii,abmn]],0] + else: + getLoc = lambda ii, abmn: np.r_[elocs[WENNER[ii,abmn]],0, 0] + srcList = [] + for i in range(WENNER.shape[0]): + rx = DC.Rx.Dipole(getLoc(i,1),getLoc(i,2)) + src = DC.Src.Dipole([rx], getLoc(i,0),getLoc(i,3)) + srcList += [src] + + return srcList diff --git a/SimPEG/EM/Static/DC/__init__.py b/SimPEG/EM/Static/DC/__init__.py index a3e1eba5..dcc40764 100644 --- a/SimPEG/EM/Static/DC/__init__.py +++ b/SimPEG/EM/Static/DC/__init__.py @@ -3,3 +3,4 @@ from SurveyDC import Survey import SrcDC as Src #Pole import RxDC as Rx from FieldsDC import Fields_CC +import Utils diff --git a/tests/em/static/__init__.py b/tests/em/static/__init__.py new file mode 100644 index 00000000..420388ef --- /dev/null +++ b/tests/em/static/__init__.py @@ -0,0 +1,12 @@ +import os +import glob +import unittest + +if __name__ == '__main__': + test_file_strings = glob.glob('test_*.py') + module_strings = [str[0:len(str)-3] for str in test_file_strings] + suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str + in module_strings] + testSuite = unittest.TestSuite(suites) + + unittest.TextTestRunner(verbosity=2).run(testSuite) diff --git a/tests/em/static/test_DC.py b/tests/em/static/test_DC.py new file mode 100644 index 00000000..a25dcf67 --- /dev/null +++ b/tests/em/static/test_DC.py @@ -0,0 +1,77 @@ +import unittest +from SimPEG import * +import SimPEG.EM.Static.DC as DC + + +class DCProblemTests(unittest.TestCase): + + def setUp(self): + + aSpacing=2.5 + nElecs=10 + + surveySize = nElecs*aSpacing - aSpacing + cs = surveySize/nElecs/4 + + mesh = Mesh.TensorMesh([ + [(cs,10, -1.3),(cs,surveySize/cs),(cs,10, 1.3)], + [(cs,3, -1.3),(cs,3,1.3)], + # [(cs,5, -1.3),(cs,10)] + ],'CN') + + srcList = DC.Utils.WennerSrcList(nElecs, aSpacing, in2D=True) + survey = DC.Survey(srcList) + problem = DC.Problem3D_CC(mesh, mapping=[('rho', Maps.IdentityMap(mesh))]) + problem.pair(survey) + + mSynth = np.ones(mesh.nC) + survey.makeSyntheticData(mSynth) + + # Now set up the problem to do some minimization + dmis = DataMisfit.l2_DataMisfit(survey) + reg = Regularization.Tikhonov(mesh) + opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6) + invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4) + inv = Inversion.BaseInversion(invProb) + + self.inv = inv + self.reg = reg + self.p = problem + self.mesh = mesh + self.m0 = mSynth + self.survey = survey + self.dmis = dmis + + def test_misfit(self): + derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)] + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False) + self.assertTrue(passed) + + def test_adjoint(self): + # Adjoint Test + u = np.random.rand(self.mesh.nC*self.survey.nSrc) + v = np.random.rand(self.mesh.nC) + w = np.random.rand(self.survey.dobs.shape[0]) + wtJv = w.dot(self.p.Jvec(self.m0, v)) + vtJtw = v.dot(self.p.Jtvec(self.m0, w)) + passed = np.abs(wtJv - vtJtw) < 1e-10 + print 'Adjoint Test', np.abs(wtJv - vtJtw), passed + self.assertTrue(passed) + + def test_dataObj(self): + derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)] + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False) + self.assertTrue(passed) + + + def test_massMatrices(self): + Gu = np.random.rand(self.mesh.nF) + def derChk(m): + self.p.curModel = m + return [self.p.Msig * Gu, self.p.dMdsig(Gu)] + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False) + self.assertTrue(passed) + + +if __name__ == '__main__': + unittest.main() diff --git a/tests/em/static/test_DC_deriv.py b/tests/em/static/test_DC_deriv.py new file mode 100644 index 00000000..e69de29b From 64b94861a01b2fcfe3da7a2c76df6ac97ff0b41c Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Fri, 22 Apr 2016 23:09:31 -0700 Subject: [PATCH 09/25] working on DC problem CC and N --- SimPEG/EM/Analytics/DC.py | 34 +++++++++++ SimPEG/EM/Analytics/__init__.py | 1 + SimPEG/EM/Static/DC/FieldsDC.py | 37 ++++++++++++ SimPEG/EM/Static/DC/ProblemDC.py | 100 +++++++++++++++++++++++++------ SimPEG/Mesh/MeshIO.py | 2 +- 5 files changed, 154 insertions(+), 20 deletions(-) create mode 100644 SimPEG/EM/Analytics/DC.py diff --git a/SimPEG/EM/Analytics/DC.py b/SimPEG/EM/Analytics/DC.py new file mode 100644 index 00000000..69d17090 --- /dev/null +++ b/SimPEG/EM/Analytics/DC.py @@ -0,0 +1,34 @@ +import numpy as np +from scipy.constants import mu_0, pi + +def DCAnalytic(txloc, rxlocs, sigma, flag="wholespace"): + """ + Analytic solution for electric potential from a postive pole + + Input variables: + + txloc = a xyz location of A (+) electrode (np.r_[xa, ya, za]) + + rxlocs = [M, N] + M: xyz locations of M (+) electrode (np.c_[xmlocs, ymlocs, zmlocs]) + N: xyz locations of N (-) electrode (np.c_[xnlocs, ynlocs, znlocs]) + + sigma = conductivity (either float or complex) + flag = "wholsespace" or "halfspace" + + """ + M = rxlocs[0] + N = rxlocs[1] + + rM = np.sqrt( (M[:,0]-txloc[0])**2 + (M[:,1]-txloc[1])**2 + (M[:,2]-txloc[1])**2 ) + rN = np.sqrt( (N[:,0]-txloc[0])**2 + (N[:,1]-txloc[1])**2 + (N[:,2]-txloc[1])**2 ) + + phiM = 1./(4*np.pi*rM*sigma) + phiN = 1./(4*np.pi*rN*sigma) + phi = phiM - phiN + + if flag == "halfspace": + phi *= 2 + + return phi + diff --git a/SimPEG/EM/Analytics/__init__.py b/SimPEG/EM/Analytics/__init__.py index 5b7a8851..d251f205 100644 --- a/SimPEG/EM/Analytics/__init__.py +++ b/SimPEG/EM/Analytics/__init__.py @@ -1,3 +1,4 @@ from TDEM import hzAnalyticDipoleT from FDEM import hzAnalyticDipoleF from FDEMcasing import * +from DC import DCAnalytic diff --git a/SimPEG/EM/Static/DC/FieldsDC.py b/SimPEG/EM/Static/DC/FieldsDC.py index 9ca221a0..91b243c1 100644 --- a/SimPEG/EM/Static/DC/FieldsDC.py +++ b/SimPEG/EM/Static/DC/FieldsDC.py @@ -69,4 +69,41 @@ class Fields_CC(Fields): def _e(self, phiSolution, srcList): raise NotImplementedError +class Fields_N(Fields): + knownFields = {'phiSolution':'N'} + aliasFields = { + 'phi': ['phiSolution','N','_phi'], + 'j' : ['phiSolution','E','_j'], + 'e' : ['phiSolution','E','_e'], + } + # primary - secondary + # N 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 'N' + elif fieldType == 'e' or fieldType == 'j': + return 'E' + 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 index f712e036..c41c8f62 100644 --- a/SimPEG/EM/Static/DC/ProblemDC.py +++ b/SimPEG/EM/Static/DC/ProblemDC.py @@ -1,7 +1,8 @@ from SimPEG import Problem from SimPEG.EM.Base import BaseEMProblem from SurveyDC import Survey -from FieldsDC import Fields, Fields_CC +from FieldsDC import Fields, Fields_CC, Fields_N +from SimPEG.Utils import sdiag import numpy as np class BaseDCProblem(BaseEMProblem): @@ -53,16 +54,15 @@ class BaseDCProblem(BaseEMProblem): q[:,i] = src.eval(self) return q -class Problem3D_CC(BaseDCProblem): +class Problem3D_N(BaseDCProblem): _solutionType = 'phiSolution' - _formulation = 'HJ' # CC potentials means J is on faces - fieldsPair = Fields_CC + _formulation = 'EB' # N potentials means B is on faces + fieldsPair = Fields_N def __init__(self, mesh, **kwargs): BaseDCProblem.__init__(self, mesh, **kwargs) - def getA(self): """ @@ -73,18 +73,21 @@ class Problem3D_CC(BaseDCProblem): """ # 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 + MeSigma = self.MeSigma + Grad = self.mesh.nodalGrad + A = Grad.T * MeSigma * Grad + # if self._makeASymmetric is True: + # return V.T * A return A - def getADeriv(): - raise NotImplementedError + def getADeriv(self, u, v, adoint=False): + """ + + Product of the derivative of our system matrix with respect to the model and a vector + + """ + return Div*self.MfRhoIDeriv(Div.T*u) + def getRHS(self): """ @@ -94,12 +97,71 @@ class Problem3D_CC(BaseDCProblem): """ RHS = self.getSourceTerm() - if self._makeASymmetric is True: - return self.Vol.T * RHS + # if self._makeASymmetric is True: + # return self.Vol.T * RHS return RHS - def getRHSDeriv(): - raise NotImplementedError + def getRHSDeriv(self, src, v, adjoint=False): + """ + Derivative of the right hand side with respect to the model + """ + qDeriv = src.evalDeriv(self, adjoint=adjoint) + return qDeriv + +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 + # V = self.Vol + # Div = V*self.mesh.faceDiv + MfRhoI = self.MfRhoI + A = self.Div * MfRhoI * self.Div.T + # if self._makeASymmetric is True: + # return V.T * A + return A + + def getADeriv(self, u, v, adoint=False): + """ + + Product of the derivative of our system matrix with respect to the model and a vector + + """ + return Div*self.MfRhoIDeriv(Div.T*u) + + + def getRHS(self): + """ + RHS for the DC problem + + q + """ + + RHS = self.getSourceTerm() + # if self._makeASymmetric is True: + # return self.Vol.T * RHS + return RHS + + def getRHSDeriv(self, src, v, adjoint=False): + """ + Derivative of the right hand side with respect to the model + """ + qDeriv = src.evalDeriv(self, adjoint=adjoint) + return qDeriv diff --git a/SimPEG/Mesh/MeshIO.py b/SimPEG/Mesh/MeshIO.py index 7501a66f..16d1b457 100644 --- a/SimPEG/Mesh/MeshIO.py +++ b/SimPEG/Mesh/MeshIO.py @@ -21,7 +21,7 @@ class TensorMeshIO(object): if '*' in seg: st = seg sp = seg.split('*') - re = np.array(sp[0],dtype=int)*(' ' + sp[1]) + re = int(sp[0])*(' ' + sp[1]) line = line.replace(st,re.strip()) return np.array(line.split(),dtype=float) From 73001abfc584edb1526ca1fbcee985a434cd528a Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Sat, 23 Apr 2016 01:24:31 -0700 Subject: [PATCH 10/25] fix bug --- SimPEG/EM/Base.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SimPEG/EM/Base.py b/SimPEG/EM/Base.py index ac28b3a1..491dcf71 100644 --- a/SimPEG/EM/Base.py +++ b/SimPEG/EM/Base.py @@ -192,7 +192,7 @@ class BaseEMProblem(Problem.BaseProblem): """ dMfRhoI_dI = -self.MfRhoI**2 - dMf_drho = self.mesh.getFaceInnerProduct(self.curModel.rho)(u) + dMf_drho = self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u) drho_dm = self.curModel.rhoDeriv return dMfRhoI_dI * ( dMf_drho * ( drho_dm)) From 8cac166fba57d361887f85fdc566cb557468f634 Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Sat, 23 Apr 2016 10:40:47 -0700 Subject: [PATCH 11/25] Working Jvec: Getting closer to understand how modular EM code is working ... --- SimPEG/EM/Static/DC/FieldsDC.py | 3 ++- SimPEG/EM/Static/DC/ProblemDC.py | 3 +-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/SimPEG/EM/Static/DC/FieldsDC.py b/SimPEG/EM/Static/DC/FieldsDC.py index 46a779f5..07b81806 100644 --- a/SimPEG/EM/Static/DC/FieldsDC.py +++ b/SimPEG/EM/Static/DC/FieldsDC.py @@ -12,6 +12,7 @@ class Fields(SimPEG.Problem.Fields): if adjoint: return self._phiDeriv_u(src, v, adjoint=adjoint), self._phiDeriv_m(src, v, adjoint=adjoint) + return np.array(self._phiDeriv_u(src, du_dm_v, adjoint) + self._phiDeriv_m(src, v, adjoint), dtype = float) def _eDeriv(self, src, du_dm_v, v, adjoint=False): @@ -59,7 +60,7 @@ class Fields_CC(Fields): return phiSolution def _phiDeriv_u(self, src, v, adjoint = False): - return Identity() + return Identity()*v def _phiDeriv_m(self, src, v, adjoint = False): return Zero() diff --git a/SimPEG/EM/Static/DC/ProblemDC.py b/SimPEG/EM/Static/DC/ProblemDC.py index 25ac7fcf..ff3fff0f 100644 --- a/SimPEG/EM/Static/DC/ProblemDC.py +++ b/SimPEG/EM/Static/DC/ProblemDC.py @@ -1,4 +1,4 @@ -from SimPEG import Problem +from SimPEG import Problem, Utils from SimPEG.EM.Base import BaseEMProblem from SurveyDC import Survey from FieldsDC import Fields, Fields_CC, Fields_N @@ -38,7 +38,6 @@ class BaseDCProblem(BaseEMProblem): u_src = f[src, self._solutionType] # solution vector dA_dm_v = self.getADeriv(u_src, v) dRHS_dm_v = self.getRHSDeriv(src, v) - print type(dA_dm_v + dRHS_dm_v), (dA_dm_v + dRHS_dm_v).shape du_dm_v = Ainv * ( - dA_dm_v + dRHS_dm_v ) for rx in src.rxList: From 0e16645b67d85ac2f0a9c5e653a07898a8dd3ac6 Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Sat, 23 Apr 2016 11:06:24 -0700 Subject: [PATCH 12/25] working Jtvec --- SimPEG/EM/Static/DC/ProblemDC.py | 43 ++++++++++++++++++++++++++------ SimPEG/EM/Static/DC/RxDC.py | 8 +++++- tests/em/static/test_DC.py | 38 +++++++++++----------------- 3 files changed, 56 insertions(+), 33 deletions(-) diff --git a/SimPEG/EM/Static/DC/ProblemDC.py b/SimPEG/EM/Static/DC/ProblemDC.py index ff3fff0f..a31731a1 100644 --- a/SimPEG/EM/Static/DC/ProblemDC.py +++ b/SimPEG/EM/Static/DC/ProblemDC.py @@ -44,22 +44,47 @@ class BaseDCProblem(BaseEMProblem): df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None) df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False) Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v) - Ainv.clean() return Utils.mkvc(Jv) def Jtvec(self, m, v, f=None): - raise NotImplementedError + if f is None: + f = self.fields(m) + + self.curModel = m + + # Ensure v is a data object. + if not isinstance(v, self.dataPair): + v = self.dataPair(self.survey, v) + + Jtv = np.zeros(m.size) + AT = self.getA().T + ATinv = self.Solver(AT, **self.solverOpts) + + for src in self.survey.srcList: + u_src = f[src, self._solutionType] + for rx in src.rxList: + PTv = rx.evalDeriv(src, self.mesh, f, v[src, rx], adjoint=True) # wrt f, need possibility wrt m + df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None) + df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True) + + ATinvdf_duT = ATinv * df_duT + + dA_dmT = self.getADeriv(u_src, ATinvdf_duT, adjoint=True) + dRHS_dmT = self.getRHSDeriv(src, ATinvdf_duT, adjoint=True) + du_dmT = -dA_dmT + dRHS_dmT + Jtv += df_dmT + du_dmT + + return Utils.mkvc(Jtv) def getSourceTerm(self): """ takes concept of source and turns it into a matrix """ """ - Evaluates the sources for a given frequency and puts them in matrix form + Evaluates the sources, and puts them in matrix form - :param float freq: Frequency :rtype: (numpy.ndarray, numpy.ndarray) - :return: s_m, s_e (nE or nF, nSrc) + :return: q (nC or nN, nSrc) """ Srcs = self.survey.srcList @@ -128,8 +153,10 @@ class Problem3D_N(BaseDCProblem): """ Derivative of the right hand side with respect to the model """ - qDeriv = src.evalDeriv(self, adjoint=adjoint) - return qDeriv + # TODO: add qDeriv for RHS depending on m + # qDeriv = src.evalDeriv(self, adjoint=adjoint) + # return qDeriv + return Zero() class Problem3D_CC(BaseDCProblem): @@ -169,7 +196,7 @@ class Problem3D_CC(BaseDCProblem): if adjoint: # if self._makeASymmetric is True: # v = V * v - return D * MfRhoIDeriv(D * v) + return( MfRhoIDeriv( D.T * u ).T) * ( D.T * v) # I think we should deprecate this for DC problem. # if self._makeASymmetric is True: diff --git a/SimPEG/EM/Static/DC/RxDC.py b/SimPEG/EM/Static/DC/RxDC.py index 7692c880..7d50392e 100644 --- a/SimPEG/EM/Static/DC/RxDC.py +++ b/SimPEG/EM/Static/DC/RxDC.py @@ -38,7 +38,13 @@ class BaseRx(SimPEG.Survey.BaseRx): def evalDeriv(self, src, mesh, f, v, adjoint=False): P = self.getP(mesh, self.projGLoc(f)) - return P*v + if not adjoint: + return P*v + elif adjoint: + return P.T*v + + + # DC.Rx.Dipole(locs) class Dipole(BaseRx): diff --git a/tests/em/static/test_DC.py b/tests/em/static/test_DC.py index 7de2ea45..6c777acf 100644 --- a/tests/em/static/test_DC.py +++ b/tests/em/static/test_DC.py @@ -47,31 +47,21 @@ class DCProblemTests(unittest.TestCase): passed = Tests.checkDerivative(derChk, self.m0, plotIt=False) self.assertTrue(passed) - # def test_adjoint(self): - # # Adjoint Test - # u = np.random.rand(self.mesh.nC*self.survey.nSrc) - # v = np.random.rand(self.mesh.nC) - # w = np.random.rand(self.survey.dobs.shape[0]) - # wtJv = w.dot(self.p.Jvec(self.m0, v)) - # vtJtw = v.dot(self.p.Jtvec(self.m0, w)) - # passed = np.abs(wtJv - vtJtw) < 1e-10 - # print 'Adjoint Test', np.abs(wtJv - vtJtw), passed - # self.assertTrue(passed) - - # def test_dataObj(self): - # derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)] - # passed = Tests.checkDerivative(derChk, self.m0, plotIt=False) - # self.assertTrue(passed) - - - # def test_massMatrices(self): - # Gu = np.random.rand(self.mesh.nF) - # def derChk(m): - # self.p.curModel = m - # return [self.p.Msig * Gu, self.p.dMdsig(Gu)] - # passed = Tests.checkDerivative(derChk, self.m0, plotIt=False) - # self.assertTrue(passed) + def test_adjoint(self): + # Adjoint Test + u = np.random.rand(self.mesh.nC*self.survey.nSrc) + v = np.random.rand(self.mesh.nC) + w = np.random.rand(self.survey.dobs.shape[0]) + wtJv = w.dot(self.p.Jvec(self.m0, v)) + vtJtw = v.dot(self.p.Jtvec(self.m0, w)) + passed = np.abs(wtJv - vtJtw) < 1e-10 + print 'Adjoint Test', np.abs(wtJv - vtJtw), passed + self.assertTrue(passed) + def test_dataObj(self): + derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)] + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False) + self.assertTrue(passed) if __name__ == '__main__': unittest.main() From 0bb001973cfd77d869f38c72e04902e8f97f693a Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Sun, 24 Apr 2016 13:01:03 -0700 Subject: [PATCH 13/25] workking nodal discretizations --- SimPEG/EM/Static/DC/BoundaryUtils.py | 158 +++++++++++++++++++++++++++ SimPEG/EM/Static/DC/ProblemDC.py | 25 +++-- SimPEG/EM/Static/DC/SrcDC.py | 12 +- SimPEG/EM/Static/DC/__init__.py | 2 +- tests/em/static/test_DC.py | 61 ++++++++++- 5 files changed, 242 insertions(+), 16 deletions(-) create mode 100644 SimPEG/EM/Static/DC/BoundaryUtils.py diff --git a/SimPEG/EM/Static/DC/BoundaryUtils.py b/SimPEG/EM/Static/DC/BoundaryUtils.py new file mode 100644 index 00000000..5ec3282c --- /dev/null +++ b/SimPEG/EM/Static/DC/BoundaryUtils.py @@ -0,0 +1,158 @@ +import numpy as np + +def getxBCyBC_CC(mesh, alpha, beta, gamma): +# def getxBCyBC(mesh, alpha, beta, gamma): + """ + This is a subfunction generating mixed-boundary condition: + + .. math:: + + \nabla \cdot \vec{j} = -\nabla \cdot \vec{j}_s = q + + \rho \vec{j} = -\nabla \phi \phi + + \alpha \phi + \beta \frac{\partial \phi}{\partial r} = \gamma \ at \ r = \partial \Omega + + xBC = f_1(\alpha, \beta, \gamma) + yBC = f(\alpha, \beta, \gamma) + + Computes xBC and yBC for cell-centered discretizations + """ + if mesh.dim == 1: #1D + if (len(alpha) != 2 or len(beta) != 2 or len(gamma) != 2): + raise Exception("Lenght of list, alpha should be 2") + fCCxm,fCCxp = mesh.cellBoundaryInd + nBC = fCCxm.sum()+fCCxp.sum() + h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp] + + alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0] + alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1] + + h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp] + + a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm) + b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm) + a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp) + b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp) + + xBC_xm = 0.5*a_xm + xBC_xp = 0.5*a_xp/b_xp + yBC_xm = 0.5*(1.-b_xm) + yBC_xp = 0.5*(1.-1./b_xp) + + xBC = np.r_[xBC_xm, xBC_xp] + yBC = np.r_[yBC_xm, yBC_xp] + + elif mesh.dim == 2: #2D + if (len(alpha) != 4 or len(beta) != 4 or len(gamma) != 4): + raise Exception("Lenght of list, alpha should be 4") + + fCCxm,fCCxp,fCCym,fCCyp = mesh.cellBoundaryInd + fxm,fxp,fym,fyp = mesh.faceBoundaryInd + nBC = fCCxm.sum()+fCCxp.sum()+fCCxm.sum()+fCCxp.sum() + h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp] + h_ym, h_yp = mesh.gridCC[fCCym], mesh.gridCC[fCCyp] + + alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0] + alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1] + alpha_ym, beta_ym, gamma_ym = alpha[2], beta[2], gamma[2] + alpha_yp, beta_yp, gamma_yp = alpha[3], beta[3], gamma[3] + + h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0] + h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1] + + a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm) + b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm) + a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp) + b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp) + + a_ym = gamma_ym/(0.5*alpha_ym-beta_ym/h_ym) + b_ym = (0.5*alpha_ym+beta_ym/h_ym)/(0.5*alpha_ym-beta_ym/h_ym) + a_yp = gamma_yp/(0.5*alpha_yp-beta_yp/h_yp) + b_yp = (0.5*alpha_yp+beta_yp/h_yp)/(0.5*alpha_yp-beta_yp/h_yp) + + xBC_xm = 0.5*a_xm + xBC_xp = 0.5*a_xp/b_xp + yBC_xm = 0.5*(1.-b_xm) + yBC_xp = 0.5*(1.-1./b_xp) + xBC_ym = 0.5*a_ym + xBC_yp = 0.5*a_yp/b_yp + yBC_ym = 0.5*(1.-b_ym) + yBC_yp = 0.5*(1.-1./b_yp) + + sortindsfx = np.argsort(np.r_[np.arange(mesh.nFx)[fxm], np.arange(mesh.nFx)[fxp]]) + sortindsfy = np.argsort(np.r_[np.arange(mesh.nFy)[fym], np.arange(mesh.nFy)[fyp]]) + + xBC_x = np.r_[xBC_xm, xBC_xp][sortindsfx] + xBC_y = np.r_[xBC_ym, xBC_yp][sortindsfy] + yBC_x = np.r_[yBC_xm, yBC_xp][sortindsfx] + yBC_y = np.r_[yBC_ym, yBC_yp][sortindsfy] + + xBC = np.r_[xBC_x, xBC_y] + yBC = np.r_[yBC_x, yBC_y] + + elif mesh.dim == 3: #3D + if (len(alpha) != 6 or len(beta) != 6 or len(gamma) != 6): + raise Exception("Lenght of list, alpha should be 6") + fCCxm,fCCxp,fCCym,fCCyp,fCCzm,fCCzp = mesh.cellBoundaryInd + fxm,fxp,fym,fyp,fzm,fzp = mesh.faceBoundaryInd + nBC = fCCxm.sum()+fCCxp.sum()+fCCxm.sum()+fCCxp.sum() + h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp] + h_ym, h_yp = mesh.gridCC[fCCym], mesh.gridCC[fCCyp] + h_zm, h_zp = mesh.gridCC[fCCzm], mesh.gridCC[fCCzp] + + alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0] + alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1] + alpha_ym, beta_ym, gamma_ym = alpha[2], beta[2], gamma[2] + alpha_yp, beta_yp, gamma_yp = alpha[3], beta[3], gamma[3] + alpha_zm, beta_zm, gamma_zm = alpha[2], beta[2], gamma[2] + alpha_zp, beta_zp, gamma_zp = alpha[3], beta[3], gamma[3] + + h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0] + h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1] + h_zm, h_zp = mesh.gridCC[fCCzm,2], mesh.gridCC[fCCzp,2] + + a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm) + b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm) + a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp) + b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp) + + a_ym = gamma_ym/(0.5*alpha_ym-beta_ym/h_ym) + b_ym = (0.5*alpha_ym+beta_ym/h_ym)/(0.5*alpha_ym-beta_ym/h_ym) + a_yp = gamma_yp/(0.5*alpha_yp-beta_yp/h_yp) + b_yp = (0.5*alpha_yp+beta_yp/h_yp)/(0.5*alpha_yp-beta_yp/h_yp) + + a_zm = gamma_zm/(0.5*alpha_zm-beta_zm/h_zm) + b_zm = (0.5*alpha_zm+beta_zm/h_zm)/(0.5*alpha_zm-beta_zm/h_zm) + a_zp = gamma_zp/(0.5*alpha_zp-beta_zp/h_zp) + b_zp = (0.5*alpha_zp+beta_zp/h_zp)/(0.5*alpha_zp-beta_zp/h_zp) + + xBC_xm = 0.5*a_xm + xBC_xp = 0.5*a_xp/b_xp + yBC_xm = 0.5*(1.-b_xm) + yBC_xp = 0.5*(1.-1./b_xp) + xBC_ym = 0.5*a_ym + xBC_yp = 0.5*a_yp/b_yp + yBC_ym = 0.5*(1.-b_ym) + yBC_yp = 0.5*(1.-1./b_yp) + xBC_zm = 0.5*a_zm + xBC_zp = 0.5*a_zp/b_zp + yBC_zm = 0.5*(1.-b_zm) + yBC_zp = 0.5*(1.-1./b_zp) + + sortindsfx = np.argsort(np.r_[np.arange(mesh.nFx)[fxm], np.arange(mesh.nFx)[fxp]]) + sortindsfy = np.argsort(np.r_[np.arange(mesh.nFy)[fym], np.arange(mesh.nFy)[fyp]]) + sortindsfz = np.argsort(np.r_[np.arange(mesh.nFz)[fzm], np.arange(mesh.nFz)[fzp]]) + + xBC_x = np.r_[xBC_xm, xBC_xp][sortindsfx] + xBC_y = np.r_[xBC_ym, xBC_yp][sortindsfy] + xBC_z = np.r_[xBC_zm, xBC_zp][sortindsfz] + + yBC_x = np.r_[yBC_xm, yBC_xp][sortindsfx] + yBC_y = np.r_[yBC_ym, yBC_yp][sortindsfy] + yBC_z = np.r_[yBC_zm, yBC_zp][sortindsfz] + + xBC = np.r_[xBC_x, xBC_y, xBC_z] + yBC = np.r_[yBC_x, yBC_y, yBC_z] + + return xBC, yBC diff --git a/SimPEG/EM/Static/DC/ProblemDC.py b/SimPEG/EM/Static/DC/ProblemDC.py index a31731a1..40d4183f 100644 --- a/SimPEG/EM/Static/DC/ProblemDC.py +++ b/SimPEG/EM/Static/DC/ProblemDC.py @@ -10,9 +10,14 @@ class BaseDCProblem(BaseEMProblem): surveyPair = Survey fieldsPair = Fields + Ainv = None def fields(self, m): self.curModel = m + + if not self.Ainv == None: + self.Ainv.clean() + f = self.fieldsPair(self.mesh, self.survey) A = self.getA() self.Ainv = self.Solver(A, **self.solverOpts) @@ -32,13 +37,12 @@ class BaseDCProblem(BaseEMProblem): Jv = self.dataPair(self.survey) #same size as the data A = self.getA() - Ainv = self.Solver(A, **self.solverOpts) for src in self.survey.srcList: u_src = f[src, self._solutionType] # solution vector dA_dm_v = self.getADeriv(u_src, v) dRHS_dm_v = self.getRHSDeriv(src, v) - du_dm_v = Ainv * ( - dA_dm_v + dRHS_dm_v ) + du_dm_v = self.Ainv * ( - dA_dm_v + dRHS_dm_v ) for rx in src.rxList: df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None) @@ -57,8 +61,8 @@ class BaseDCProblem(BaseEMProblem): v = self.dataPair(self.survey, v) Jtv = np.zeros(m.size) - AT = self.getA().T - ATinv = self.Solver(AT, **self.solverOpts) + AT = self.getA() + for src in self.survey.srcList: u_src = f[src, self._solutionType] @@ -67,7 +71,7 @@ class BaseDCProblem(BaseEMProblem): df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None) df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True) - ATinvdf_duT = ATinv * df_duT + ATinvdf_duT = self.Ainv * df_duT dA_dmT = self.getADeriv(u_src, ATinvdf_duT, adjoint=True) dRHS_dmT = self.getRHSDeriv(src, ATinvdf_duT, adjoint=True) @@ -128,13 +132,18 @@ class Problem3D_N(BaseDCProblem): # return V.T * A return A - def getADeriv(self, u, v, adoint=False): + def getADeriv(self, u, v, adjoint=False): """ Product of the derivative of our system matrix with respect to the model and a vector """ - return Div*self.MfRhoIDeriv(Div.T*u) + MeSigma = self.MeSigma + Grad = self.mesh.nodalGrad + if not adjoint: + return Grad.T*(self.MeSigmaDeriv(Grad*u)*v) + elif adjoint: + return self.MeSigmaDeriv(Grad*u).T * (Grad*v) def getRHS(self): @@ -196,7 +205,7 @@ class Problem3D_CC(BaseDCProblem): if adjoint: # if self._makeASymmetric is True: # v = V * v - return( MfRhoIDeriv( D.T * u ).T) * ( D.T * v) + return(MfRhoIDeriv( D.T * u ).T) * ( D.T * v) # I think we should deprecate this for DC problem. # if self._makeASymmetric is True: diff --git a/SimPEG/EM/Static/DC/SrcDC.py b/SimPEG/EM/Static/DC/SrcDC.py index cc855118..3ed6067d 100644 --- a/SimPEG/EM/Static/DC/SrcDC.py +++ b/SimPEG/EM/Static/DC/SrcDC.py @@ -31,9 +31,9 @@ class Dipole(BaseSrc): 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 + inds = closestPoints(prob.mesh, self.loc) + q = np.zeros(prob.mesh.nN) + q[inds] = self.current * np.r_[1., -1.] return q # def bc_contribution @@ -52,9 +52,9 @@ class Pole(BaseSrc): 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 + inds = closestPoints(prob.mesh, self.loc) + q = np.zeros(prob.mesh.nN) + q[inds] = self.current * np.r_[1.] return q # def bc_contribution diff --git a/SimPEG/EM/Static/DC/__init__.py b/SimPEG/EM/Static/DC/__init__.py index dcc40764..57da57e8 100644 --- a/SimPEG/EM/Static/DC/__init__.py +++ b/SimPEG/EM/Static/DC/__init__.py @@ -1,4 +1,4 @@ -from ProblemDC import Problem3D_CC +from ProblemDC import Problem3D_CC, Problem3D_N from SurveyDC import Survey import SrcDC as Src #Pole import RxDC as Rx diff --git a/tests/em/static/test_DC.py b/tests/em/static/test_DC.py index 6c777acf..3472a27e 100644 --- a/tests/em/static/test_DC.py +++ b/tests/em/static/test_DC.py @@ -3,7 +3,7 @@ from SimPEG import * import SimPEG.EM.Static.DC as DC -class DCProblemTests(unittest.TestCase): +class DCProblemTestsCC(unittest.TestCase): def setUp(self): @@ -63,5 +63,64 @@ class DCProblemTests(unittest.TestCase): passed = Tests.checkDerivative(derChk, self.m0, plotIt=False) self.assertTrue(passed) +class DCProblemTestsN(unittest.TestCase): + + def setUp(self): + + aSpacing=2.5 + nElecs=10 + + surveySize = nElecs*aSpacing - aSpacing + cs = surveySize/nElecs/4 + + mesh = Mesh.TensorMesh([ + [(cs,10, -1.3),(cs,surveySize/cs),(cs,10, 1.3)], + [(cs,3, -1.3),(cs,3,1.3)], + # [(cs,5, -1.3),(cs,10)] + ],'CN') + + srcList = DC.Utils.WennerSrcList(nElecs, aSpacing, in2D=True) + survey = DC.Survey(srcList) + problem = DC.Problem3D_N(mesh, mapping=[('rho', Maps.IdentityMap(mesh))]) + problem.pair(survey) + + mSynth = np.ones(mesh.nC) + survey.makeSyntheticData(mSynth) + + # Now set up the problem to do some minimization + dmis = DataMisfit.l2_DataMisfit(survey) + reg = Regularization.Tikhonov(mesh) + opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6) + invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4) + inv = Inversion.BaseInversion(invProb) + + self.inv = inv + self.reg = reg + self.p = problem + self.mesh = mesh + self.m0 = mSynth + self.survey = survey + self.dmis = dmis + + def test_misfit(self): + derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)] + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False) + self.assertTrue(passed) + + def test_adjoint(self): + # Adjoint Test + u = np.random.rand(self.mesh.nC*self.survey.nSrc) + v = np.random.rand(self.mesh.nC) + w = np.random.rand(self.survey.dobs.shape[0]) + wtJv = w.dot(self.p.Jvec(self.m0, v)) + vtJtw = v.dot(self.p.Jtvec(self.m0, w)) + passed = np.abs(wtJv - vtJtw) < 1e-10 + print 'Adjoint Test', np.abs(wtJv - vtJtw), passed + self.assertTrue(passed) + + def test_dataObj(self): + derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)] + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False) + self.assertTrue(passed) if __name__ == '__main__': unittest.main() From a48224ed8b6eca2233f2501d55c9ff08b29e5059 Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Sun, 24 Apr 2016 13:18:36 -0700 Subject: [PATCH 14/25] no message --- tests/mesh/test_Mixed_boundaryPoisson.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/mesh/test_Mixed_boundaryPoisson.py b/tests/mesh/test_Mixed_boundaryPoisson.py index 7a15f36d..d5dcb6fa 100644 --- a/tests/mesh/test_Mixed_boundaryPoisson.py +++ b/tests/mesh/test_Mixed_boundaryPoisson.py @@ -194,6 +194,7 @@ class Test1D_InhomogeneousMixed(Tests.OrderTest): q = q_fun(self.M.gridCC) M = B*self.M.aveCC2F G = Div.T - P_BC*Utils.sdiag(y_BC)*M + # Mrhoj = D.T V phi + P_BC*Utils.sdiag(y_BC)*M phi - P_BC*x_BC rhs = V*q + Div*MfrhoI*P_BC*x_BC A = Div*MfrhoI*G From eeee594f09a8f13d592848f5d9d8dec712c523ef Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Sun, 24 Apr 2016 13:29:38 -0700 Subject: [PATCH 15/25] Problem3D_N is tested! --- SimPEG/EM/Static/DC/FieldsDC.py | 2 +- tests/em/static/test_DC.py | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/SimPEG/EM/Static/DC/FieldsDC.py b/SimPEG/EM/Static/DC/FieldsDC.py index 07b81806..9999c56a 100644 --- a/SimPEG/EM/Static/DC/FieldsDC.py +++ b/SimPEG/EM/Static/DC/FieldsDC.py @@ -99,7 +99,7 @@ class Fields_N(Fields): return phiSolution def _phiDeriv_u(self, src, v, adjoint = False): - return Identity() + return Identity()*v def _phiDeriv_m(self, src, v, adjoint = False): return Zero() diff --git a/tests/em/static/test_DC.py b/tests/em/static/test_DC.py index 3472a27e..1d83d32c 100644 --- a/tests/em/static/test_DC.py +++ b/tests/em/static/test_DC.py @@ -104,7 +104,7 @@ class DCProblemTestsN(unittest.TestCase): def test_misfit(self): derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)] - passed = Tests.checkDerivative(derChk, self.m0, plotIt=False) + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) self.assertTrue(passed) def test_adjoint(self): @@ -114,13 +114,13 @@ class DCProblemTestsN(unittest.TestCase): w = np.random.rand(self.survey.dobs.shape[0]) wtJv = w.dot(self.p.Jvec(self.m0, v)) vtJtw = v.dot(self.p.Jtvec(self.m0, w)) - passed = np.abs(wtJv - vtJtw) < 1e-10 + passed = np.abs(wtJv - vtJtw) < 1e-8 print 'Adjoint Test', np.abs(wtJv - vtJtw), passed self.assertTrue(passed) def test_dataObj(self): derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)] - passed = Tests.checkDerivative(derChk, self.m0, plotIt=False) + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) self.assertTrue(passed) if __name__ == '__main__': unittest.main() From fcc2b8b22a9d09a837c91536f5f51bf1ecfcf603 Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Sun, 24 Apr 2016 13:32:44 -0700 Subject: [PATCH 16/25] Handling null space of A --- SimPEG/EM/Static/DC/ProblemDC.py | 4 ++++ tests/em/static/test_DC.py | 4 ++-- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/SimPEG/EM/Static/DC/ProblemDC.py b/SimPEG/EM/Static/DC/ProblemDC.py index 40d4183f..04a07c57 100644 --- a/SimPEG/EM/Static/DC/ProblemDC.py +++ b/SimPEG/EM/Static/DC/ProblemDC.py @@ -128,6 +128,10 @@ class Problem3D_N(BaseDCProblem): MeSigma = self.MeSigma Grad = self.mesh.nodalGrad A = Grad.T * MeSigma * Grad + + # Handling ... singularity + A[0,0] = A[0,0] + 1. + # if self._makeASymmetric is True: # return V.T * A return A diff --git a/tests/em/static/test_DC.py b/tests/em/static/test_DC.py index 1d83d32c..99b384aa 100644 --- a/tests/em/static/test_DC.py +++ b/tests/em/static/test_DC.py @@ -104,7 +104,7 @@ class DCProblemTestsN(unittest.TestCase): def test_misfit(self): derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)] - passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False) self.assertTrue(passed) def test_adjoint(self): @@ -120,7 +120,7 @@ class DCProblemTestsN(unittest.TestCase): def test_dataObj(self): derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)] - passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False) self.assertTrue(passed) if __name__ == '__main__': unittest.main() From 1936a046833d3d5b8a073d2e6a87c2206a1c9c97 Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Sun, 24 Apr 2016 13:56:37 -0700 Subject: [PATCH 17/25] Working on implementing mixed BC to DC problem --- SimPEG/EM/Static/DC/ProblemDC.py | 139 ++++++++++++++++--------------- 1 file changed, 74 insertions(+), 65 deletions(-) diff --git a/SimPEG/EM/Static/DC/ProblemDC.py b/SimPEG/EM/Static/DC/ProblemDC.py index 04a07c57..7ec32919 100644 --- a/SimPEG/EM/Static/DC/ProblemDC.py +++ b/SimPEG/EM/Static/DC/ProblemDC.py @@ -5,6 +5,7 @@ from FieldsDC import Fields, Fields_CC, Fields_N from SimPEG.Utils import sdiag import numpy as np from SimPEG.Utils import Zero +from BoundaryUtils import getxBCyBC_CC class BaseDCProblem(BaseEMProblem): @@ -106,71 +107,6 @@ class BaseDCProblem(BaseEMProblem): q[:,i] = src.eval(self) return q -class Problem3D_N(BaseDCProblem): - - _solutionType = 'phiSolution' - _formulation = 'EB' # N potentials means B is on faces - fieldsPair = Fields_N - - 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 - MeSigma = self.MeSigma - Grad = self.mesh.nodalGrad - A = Grad.T * MeSigma * Grad - - # Handling ... singularity - A[0,0] = A[0,0] + 1. - - # if self._makeASymmetric is True: - # return V.T * A - return A - - def getADeriv(self, u, v, adjoint=False): - """ - - Product of the derivative of our system matrix with respect to the model and a vector - - """ - MeSigma = self.MeSigma - Grad = self.mesh.nodalGrad - if not adjoint: - return Grad.T*(self.MeSigmaDeriv(Grad*u)*v) - elif adjoint: - return self.MeSigmaDeriv(Grad*u).T * (Grad*v) - - - def getRHS(self): - """ - RHS for the DC problem - - q - """ - - RHS = self.getSourceTerm() - # if self._makeASymmetric is True: - # return self.Vol.T * RHS - return RHS - - def getRHSDeriv(self, src, v, adjoint=False): - """ - Derivative of the right hand side with respect to the model - """ - # TODO: add qDeriv for RHS depending on m - # qDeriv = src.evalDeriv(self, adjoint=adjoint) - # return qDeriv - return Zero() - class Problem3D_CC(BaseDCProblem): _solutionType = 'phiSolution' @@ -180,6 +116,12 @@ class Problem3D_CC(BaseDCProblem): def __init__(self, mesh, **kwargs): BaseDCProblem.__init__(self, mesh, **kwargs) + def setBC(self): + self.Div = V * self.mesh.faceDiv + P_BC, B = self.mesh.getBCProjWF_simple() + M = B*self.mesh.aveCC2F + Grad = Div.T - P_BC*Utils.sdiag(y_BC)*M + def getA(self): """ @@ -241,4 +183,71 @@ class Problem3D_CC(BaseDCProblem): return Zero() +class Problem3D_N(BaseDCProblem): + + _solutionType = 'phiSolution' + _formulation = 'EB' # N potentials means B is on faces + fieldsPair = Fields_N + + 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 + MeSigma = self.MeSigma + Grad = self.mesh.nodalGrad + A = Grad.T * MeSigma * Grad + + # Handling Null space of A + A[0,0] = A[0,0] + 1. + + # if self._makeASymmetric is True: + # return V.T * A + return A + + def getADeriv(self, u, v, adjoint=False): + """ + + Product of the derivative of our system matrix with respect to the model and a vector + + """ + MeSigma = self.MeSigma + Grad = self.mesh.nodalGrad + if not adjoint: + return Grad.T*(self.MeSigmaDeriv(Grad*u)*v) + elif adjoint: + return self.MeSigmaDeriv(Grad*u).T * (Grad*v) + + + def getRHS(self): + """ + RHS for the DC problem + + q + """ + + RHS = self.getSourceTerm() + # if self._makeASymmetric is True: + # return self.Vol.T * RHS + return RHS + + def getRHSDeriv(self, src, v, adjoint=False): + """ + Derivative of the right hand side with respect to the model + """ + # TODO: add qDeriv for RHS depending on m + # qDeriv = src.evalDeriv(self, adjoint=adjoint) + # return qDeriv + return Zero() + + + From dcd4fbf97311762376c3fe2698cb980580d0f086 Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Sun, 24 Apr 2016 15:23:14 -0700 Subject: [PATCH 18/25] Implemented mixed B.C. to CC problem. Fix bugs in get fuction getxBCyBC_CC --- SimPEG/EM/Static/DC/BoundaryUtils.py | 36 +++++----- SimPEG/EM/Static/DC/ProblemDC.py | 84 ++++++++++++++++++++---- SimPEG/EM/Static/DC/SrcDC.py | 2 +- tests/em/static/test_DC.py | 9 +-- tests/mesh/test_Mixed_boundaryPoisson.py | 63 +++++++++++------- 5 files changed, 135 insertions(+), 59 deletions(-) diff --git a/SimPEG/EM/Static/DC/BoundaryUtils.py b/SimPEG/EM/Static/DC/BoundaryUtils.py index 5ec3282c..3967eb46 100644 --- a/SimPEG/EM/Static/DC/BoundaryUtils.py +++ b/SimPEG/EM/Static/DC/BoundaryUtils.py @@ -28,7 +28,8 @@ def getxBCyBC_CC(mesh, alpha, beta, gamma): alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0] alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1] - h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp] + # h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp] + h_xm, h_xp = mesh.hx[0], mesh.hx[-1] a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm) b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm) @@ -47,19 +48,19 @@ def getxBCyBC_CC(mesh, alpha, beta, gamma): if (len(alpha) != 4 or len(beta) != 4 or len(gamma) != 4): raise Exception("Lenght of list, alpha should be 4") - fCCxm,fCCxp,fCCym,fCCyp = mesh.cellBoundaryInd fxm,fxp,fym,fyp = mesh.faceBoundaryInd - nBC = fCCxm.sum()+fCCxp.sum()+fCCxm.sum()+fCCxp.sum() - h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp] - h_ym, h_yp = mesh.gridCC[fCCym], mesh.gridCC[fCCyp] + nBC = fxm.sum()+fxp.sum()+fxm.sum()+fxp.sum() alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0] alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1] alpha_ym, beta_ym, gamma_ym = alpha[2], beta[2], gamma[2] alpha_yp, beta_yp, gamma_yp = alpha[3], beta[3], gamma[3] - h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0] - h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1] + # h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0] + # h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1] + + h_xm, h_xp = mesh.hx[0]*np.ones_like(alpha_xm), mesh.hx[-1]*np.ones_like(alpha_xp) + h_ym, h_yp = mesh.hy[0]*np.ones_like(alpha_ym), mesh.hy[-1]*np.ones_like(alpha_yp) a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm) b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm) @@ -94,23 +95,24 @@ def getxBCyBC_CC(mesh, alpha, beta, gamma): elif mesh.dim == 3: #3D if (len(alpha) != 6 or len(beta) != 6 or len(gamma) != 6): raise Exception("Lenght of list, alpha should be 6") - fCCxm,fCCxp,fCCym,fCCyp,fCCzm,fCCzp = mesh.cellBoundaryInd + # fCCxm,fCCxp,fCCym,fCCyp,fCCzm,fCCzp = mesh.cellBoundaryInd fxm,fxp,fym,fyp,fzm,fzp = mesh.faceBoundaryInd - nBC = fCCxm.sum()+fCCxp.sum()+fCCxm.sum()+fCCxp.sum() - h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp] - h_ym, h_yp = mesh.gridCC[fCCym], mesh.gridCC[fCCyp] - h_zm, h_zp = mesh.gridCC[fCCzm], mesh.gridCC[fCCzp] + nBC = fxm.sum()+fxp.sum()+fxm.sum()+fxp.sum() alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0] alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1] alpha_ym, beta_ym, gamma_ym = alpha[2], beta[2], gamma[2] alpha_yp, beta_yp, gamma_yp = alpha[3], beta[3], gamma[3] - alpha_zm, beta_zm, gamma_zm = alpha[2], beta[2], gamma[2] - alpha_zp, beta_zp, gamma_zp = alpha[3], beta[3], gamma[3] + alpha_zm, beta_zm, gamma_zm = alpha[4], beta[4], gamma[4] + alpha_zp, beta_zp, gamma_zp = alpha[5], beta[5], gamma[5] - h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0] - h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1] - h_zm, h_zp = mesh.gridCC[fCCzm,2], mesh.gridCC[fCCzp,2] + # h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0] + # h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1] + # h_zm, h_zp = mesh.gridCC[fCCzm,2], mesh.gridCC[fCCzp,2] + + h_xm, h_xp = mesh.hx[0]*np.ones_like(alpha_xm), mesh.hx[-1]*np.ones_like(alpha_xp) + h_ym, h_yp = mesh.hy[0]*np.ones_like(alpha_ym), mesh.hy[-1]*np.ones_like(alpha_yp) + h_zm, h_zp = mesh.hz[0]*np.ones_like(alpha_zm), mesh.hz[-1]*np.ones_like(alpha_zp) a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm) b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm) diff --git a/SimPEG/EM/Static/DC/ProblemDC.py b/SimPEG/EM/Static/DC/ProblemDC.py index 7ec32919..9f05d876 100644 --- a/SimPEG/EM/Static/DC/ProblemDC.py +++ b/SimPEG/EM/Static/DC/ProblemDC.py @@ -115,12 +115,7 @@ class Problem3D_CC(BaseDCProblem): def __init__(self, mesh, **kwargs): BaseDCProblem.__init__(self, mesh, **kwargs) - - def setBC(self): - self.Div = V * self.mesh.faceDiv - P_BC, B = self.mesh.getBCProjWF_simple() - M = B*self.mesh.aveCC2F - Grad = Div.T - P_BC*Utils.sdiag(y_BC)*M + self.setBC() def getA(self): """ @@ -131,11 +126,11 @@ class Problem3D_CC(BaseDCProblem): """ - V = self.Vol - D = V * self.mesh.faceDiv + D = self.Div + G = self.Grad # TODO: this won't work for full anisotropy MfRhoI = self.MfRhoI - A = D * MfRhoI * D.T + A = D * MfRhoI * G # I think we should deprecate this for DC problem. # if self._makeASymmetric is True: @@ -144,19 +139,19 @@ class Problem3D_CC(BaseDCProblem): def getADeriv(self, u, v, adjoint= False): - V = self.Vol - D = V * self.mesh.faceDiv + D = self.Div + G = self.Grad MfRhoIDeriv = self.MfRhoIDeriv if adjoint: # if self._makeASymmetric is True: # v = V * v - return(MfRhoIDeriv( D.T * u ).T) * ( D.T * v) + return(MfRhoIDeriv( G * u ).T) * ( D.T * v) # I think we should deprecate this for DC problem. # if self._makeASymmetric is True: # return V.T * ( D * ( MfRhoIDeriv( D.T * ( V * u ) ) * v ) ) - return D * (MfRhoIDeriv( D.T * u ) * v) + return D * (MfRhoIDeriv( G * u ) * v) def getRHS(self): """ @@ -182,6 +177,69 @@ class Problem3D_CC(BaseDCProblem): # return qDeriv return Zero() + def setBC(self): + if self.mesh.dim==3: + fxm,fxp,fym,fyp,fzm,fzp = self.mesh.faceBoundaryInd + gBFxm = self.mesh.gridFx[fxm,:] + gBFxp = self.mesh.gridFx[fxp,:] + gBFym = self.mesh.gridFy[fym,:] + gBFyp = self.mesh.gridFy[fyp,:] + gBFzm = self.mesh.gridFz[fzm,:] + gBFzp = self.mesh.gridFz[fzp,:] + + # Setup Mixed B.C (alpha, beta, gamma) + temp_xm, temp_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0]) + temp_ym, temp_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1]) + temp_zm, temp_zp = np.ones_like(gBFzm[:,2]), np.ones_like(gBFzp[:,2]) + + alpha_xm, alpha_xp = temp_xm*0., temp_xp*0. + alpha_ym, alpha_yp = temp_ym*0., temp_yp*0. + alpha_zm, alpha_zp = temp_zm*0., temp_zp*0. + + beta_xm, beta_xp = temp_xm, temp_xp + beta_ym, beta_yp = temp_ym, temp_yp + beta_zm, beta_zp = temp_zm, temp_zp + + gamma_xm, gamma_xp = temp_xm*0., temp_xp*0. + gamma_ym, gamma_yp = temp_ym*0., temp_yp*0. + gamma_zm, gamma_zp = temp_zm*0., temp_zp*0. + + alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp, alpha_zm, alpha_zp] + beta = [beta_xm, beta_xp, beta_ym, beta_yp, beta_zm, beta_zp] + gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp, gamma_zm, gamma_zp] + + elif self.mesh.dim==2: + + fxm,fxp,fym,fyp = self.mesh.faceBoundaryInd + gBFxm = self.mesh.gridFx[fxm,:] + gBFxp = self.mesh.gridFx[fxp,:] + gBFym = self.mesh.gridFy[fym,:] + gBFyp = self.mesh.gridFy[fyp,:] + + # Setup Mixed B.C (alpha, beta, gamma) + temp_xm, temp_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0]) + temp_ym, temp_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1]) + + alpha_xm, alpha_xp = temp_xm*0., temp_xp*0. + alpha_ym, alpha_yp = temp_ym*0., temp_yp*0. + + beta_xm, beta_xp = temp_xm, temp_xp + beta_ym, beta_yp = temp_ym, temp_yp + + gamma_xm, gamma_xp = temp_xm*0., temp_xp*0. + gamma_ym, gamma_yp = temp_ym*0., temp_yp*0. + + alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp] + beta = [beta_xm, beta_xp, beta_ym, beta_yp] + gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp] + + x_BC, y_BC = getxBCyBC_CC(self.mesh, alpha, beta, gamma) + V = self.Vol + self.Div = V * self.mesh.faceDiv + P_BC, B = self.mesh.getBCProjWF_simple() + M = B*self.mesh.aveCC2F + self.Grad = self.Div.T - P_BC*Utils.sdiag(y_BC)*M + class Problem3D_N(BaseDCProblem): diff --git a/SimPEG/EM/Static/DC/SrcDC.py b/SimPEG/EM/Static/DC/SrcDC.py index 3ed6067d..1e64835e 100644 --- a/SimPEG/EM/Static/DC/SrcDC.py +++ b/SimPEG/EM/Static/DC/SrcDC.py @@ -15,7 +15,7 @@ class BaseSrc(SimPEG.Survey.BaseSrc): raise NotImplementedError def evalDeriv(self, prob): - Zero() + return Zero() class Dipole(BaseSrc): diff --git a/tests/em/static/test_DC.py b/tests/em/static/test_DC.py index 99b384aa..227d4614 100644 --- a/tests/em/static/test_DC.py +++ b/tests/em/static/test_DC.py @@ -8,7 +8,7 @@ class DCProblemTestsCC(unittest.TestCase): def setUp(self): aSpacing=2.5 - nElecs=10 + nElecs=5 surveySize = nElecs*aSpacing - aSpacing cs = surveySize/nElecs/4 @@ -16,7 +16,7 @@ class DCProblemTestsCC(unittest.TestCase): mesh = Mesh.TensorMesh([ [(cs,10, -1.3),(cs,surveySize/cs),(cs,10, 1.3)], [(cs,3, -1.3),(cs,3,1.3)], - # [(cs,5, -1.3),(cs,10)] + # [(cs,5, -1.3),(cs,10)] ],'CN') srcList = DC.Utils.WennerSrcList(nElecs, aSpacing, in2D=True) @@ -44,7 +44,7 @@ class DCProblemTestsCC(unittest.TestCase): def test_misfit(self): derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)] - passed = Tests.checkDerivative(derChk, self.m0, plotIt=False) + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) self.assertTrue(passed) def test_adjoint(self): @@ -60,7 +60,7 @@ class DCProblemTestsCC(unittest.TestCase): def test_dataObj(self): derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)] - passed = Tests.checkDerivative(derChk, self.m0, plotIt=False) + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) self.assertTrue(passed) class DCProblemTestsN(unittest.TestCase): @@ -122,5 +122,6 @@ class DCProblemTestsN(unittest.TestCase): derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)] passed = Tests.checkDerivative(derChk, self.m0, plotIt=False) self.assertTrue(passed) + if __name__ == '__main__': unittest.main() diff --git a/tests/mesh/test_Mixed_boundaryPoisson.py b/tests/mesh/test_Mixed_boundaryPoisson.py index d5dcb6fa..3aa1dbcd 100644 --- a/tests/mesh/test_Mixed_boundaryPoisson.py +++ b/tests/mesh/test_Mixed_boundaryPoisson.py @@ -6,10 +6,23 @@ from SimPEG import * MESHTYPES = ['uniformTensorMesh'] -def getxBCyBC(mesh, alpha, beta, gamma): +def getxBCyBC_CC(mesh, alpha, beta, gamma): # def getxBCyBC(mesh, alpha, beta, gamma): """ + This is a subfunction generating mixed-boundary condition: + .. math:: + + \nabla \cdot \vec{j} = -\nabla \cdot \vec{j}_s = q + + \rho \vec{j} = -\nabla \phi \phi + + \alpha \phi + \beta \frac{\partial \phi}{\partial r} = \gamma \ at \ r = \partial \Omega + + xBC = f_1(\alpha, \beta, \gamma) + yBC = f(\alpha, \beta, \gamma) + + Computes xBC and yBC for cell-centered discretizations """ if mesh.dim == 1: #1D if (len(alpha) != 2 or len(beta) != 2 or len(gamma) != 2): @@ -21,7 +34,8 @@ def getxBCyBC(mesh, alpha, beta, gamma): alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0] alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1] - h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp] + # h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp] + h_xm, h_xp = mesh.hx[0], mesh.hx[-1] a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm) b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm) @@ -40,19 +54,19 @@ def getxBCyBC(mesh, alpha, beta, gamma): if (len(alpha) != 4 or len(beta) != 4 or len(gamma) != 4): raise Exception("Lenght of list, alpha should be 4") - fCCxm,fCCxp,fCCym,fCCyp = mesh.cellBoundaryInd fxm,fxp,fym,fyp = mesh.faceBoundaryInd - nBC = fCCxm.sum()+fCCxp.sum()+fCCxm.sum()+fCCxp.sum() - h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp] - h_ym, h_yp = mesh.gridCC[fCCym], mesh.gridCC[fCCyp] + nBC = fxm.sum()+fxp.sum()+fxm.sum()+fxp.sum() alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0] alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1] alpha_ym, beta_ym, gamma_ym = alpha[2], beta[2], gamma[2] alpha_yp, beta_yp, gamma_yp = alpha[3], beta[3], gamma[3] - h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0] - h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1] + # h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0] + # h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1] + + h_xm, h_xp = mesh.hx[0]*np.ones_like(alpha_xm), mesh.hx[-1]*np.ones_like(alpha_xp) + h_ym, h_yp = mesh.hy[0]*np.ones_like(alpha_ym), mesh.hy[-1]*np.ones_like(alpha_yp) a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm) b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm) @@ -87,23 +101,24 @@ def getxBCyBC(mesh, alpha, beta, gamma): elif mesh.dim == 3: #3D if (len(alpha) != 6 or len(beta) != 6 or len(gamma) != 6): raise Exception("Lenght of list, alpha should be 6") - fCCxm,fCCxp,fCCym,fCCyp,fCCzm,fCCzp = mesh.cellBoundaryInd + # fCCxm,fCCxp,fCCym,fCCyp,fCCzm,fCCzp = mesh.cellBoundaryInd fxm,fxp,fym,fyp,fzm,fzp = mesh.faceBoundaryInd - nBC = fCCxm.sum()+fCCxp.sum()+fCCxm.sum()+fCCxp.sum() - h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp] - h_ym, h_yp = mesh.gridCC[fCCym], mesh.gridCC[fCCyp] - h_zm, h_zp = mesh.gridCC[fCCzm], mesh.gridCC[fCCzp] + nBC = fxm.sum()+fxp.sum()+fxm.sum()+fxp.sum() alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0] alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1] alpha_ym, beta_ym, gamma_ym = alpha[2], beta[2], gamma[2] alpha_yp, beta_yp, gamma_yp = alpha[3], beta[3], gamma[3] - alpha_zm, beta_zm, gamma_zm = alpha[2], beta[2], gamma[2] - alpha_zp, beta_zp, gamma_zp = alpha[3], beta[3], gamma[3] + alpha_zm, beta_zm, gamma_zm = alpha[4], beta[4], gamma[4] + alpha_zp, beta_zp, gamma_zp = alpha[5], beta[5], gamma[5] - h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0] - h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1] - h_zm, h_zp = mesh.gridCC[fCCzm,2], mesh.gridCC[fCCzp,2] + # h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0] + # h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1] + # h_zm, h_zp = mesh.gridCC[fCCzm,2], mesh.gridCC[fCCzp,2] + + h_xm, h_xp = mesh.hx[0]*np.ones_like(alpha_xm), mesh.hx[-1]*np.ones_like(alpha_xp) + h_ym, h_yp = mesh.hy[0]*np.ones_like(alpha_ym), mesh.hy[-1]*np.ones_like(alpha_yp) + h_zm, h_zp = mesh.hz[0]*np.ones_like(alpha_zm), mesh.hz[-1]*np.ones_like(alpha_zp) a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm) b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm) @@ -182,7 +197,7 @@ class Test1D_InhomogeneousMixed(Tests.OrderTest): phi_bc = phi_fun(vecN[[0,-1]]) phi_deriv_bc = phi_deriv(vecN[[0,-1]]) gamma = alpha*phi_bc + beta*phi_deriv_bc - x_BC, y_BC = getxBCyBC(self.M, alpha, beta, gamma) + x_BC, y_BC = getxBCyBC_CC(self.M, alpha, beta, gamma) sigma = np.ones(self.M.nC) @@ -265,7 +280,7 @@ class Test2D_InhomogeneousMixed(Tests.OrderTest): beta = [beta_xm, beta_xp, beta_ym, beta_yp] gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp] - x_BC, y_BC = getxBCyBC(self.M, alpha, beta, gamma) + x_BC, y_BC = getxBCyBC_CC(self.M, alpha, beta, gamma) sigma = np.ones(self.M.nC) @@ -335,8 +350,8 @@ class Test3D_InhomogeneousMixed(Tests.OrderTest): beta_xm, beta_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0]) alpha_ym, alpha_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1]) beta_ym, beta_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1]) - alpha_zm, alpha_zp = np.ones_like(gBFzm[:,1]), np.ones_like(gBFzp[:,1]) - beta_zm, beta_zp = np.ones_like(gBFzm[:,1]), np.ones_like(gBFzp[:,1]) + alpha_zm, alpha_zp = np.ones_like(gBFzm[:,2]), np.ones_like(gBFzp[:,2]) + beta_zm, beta_zp = np.ones_like(gBFzm[:,2]), np.ones_like(gBFzp[:,2]) phi_bc_xm, phi_bc_xp = phi_fun(gBFxm), phi_fun(gBFxp) @@ -345,7 +360,7 @@ class Test3D_InhomogeneousMixed(Tests.OrderTest): phiderivX_bc_xm, phiderivX_bc_xp = phideriv_funX(gBFxm), phideriv_funX(gBFxp) phiderivY_bc_ym, phiderivY_bc_yp = phideriv_funY(gBFym), phideriv_funY(gBFyp) - phiderivY_bc_zm, phiderivY_bc_zp = phideriv_funY(gBFzm), phideriv_funY(gBFzp) + phiderivY_bc_zm, phiderivY_bc_zp = phideriv_funZ(gBFzm), phideriv_funZ(gBFzp) gamma_fun = lambda alpha, beta, phi, phi_deriv: alpha*phi + beta*phi_deriv gamma_xm = gamma_fun(alpha_xm, beta_xm, phi_bc_xm, phiderivX_bc_xm) @@ -359,7 +374,7 @@ class Test3D_InhomogeneousMixed(Tests.OrderTest): beta = [beta_xm, beta_xp, beta_ym, beta_yp, beta_zm, beta_zp] gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp, gamma_zm, gamma_zp] - x_BC, y_BC = getxBCyBC(self.M, alpha, beta, gamma) + x_BC, y_BC = getxBCyBC_CC(self.M, alpha, beta, gamma) sigma = np.ones(self.M.nC) From f944f9b76b81ad85ab7477f18184cf6f0c6d60b6 Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Mon, 25 Apr 2016 10:58:54 -0700 Subject: [PATCH 19/25] 1. Add distributed source for nodal discretization 2. Add Analytic tests 3. Fix simple bug in PlotSlice for nodal variable 4. Add more analytic function (sphere) --- SimPEG/EM/Analytics/DC.py | 81 ++++++++++++++++++- SimPEG/EM/Analytics/__init__.py | 2 +- SimPEG/EM/Static/DC/SrcDC.py | 10 +-- SimPEG/Mesh/View.py | 2 +- tests/em/static/test_DC_analytic.py | 68 ++++++++++++++++ tests/em/static/test_DC_deriv.py | 0 .../{test_DC.py => test_DC_jvecjtvecadj.py} | 0 7 files changed, 155 insertions(+), 8 deletions(-) create mode 100644 tests/em/static/test_DC_analytic.py delete mode 100644 tests/em/static/test_DC_deriv.py rename tests/em/static/{test_DC.py => test_DC_jvecjtvecadj.py} (100%) diff --git a/SimPEG/EM/Analytics/DC.py b/SimPEG/EM/Analytics/DC.py index 69d17090..4fb03fd3 100644 --- a/SimPEG/EM/Analytics/DC.py +++ b/SimPEG/EM/Analytics/DC.py @@ -1,7 +1,7 @@ import numpy as np from scipy.constants import mu_0, pi -def DCAnalytic(txloc, rxlocs, sigma, flag="wholespace"): +def DCAnalyticHalf(txloc, rxlocs, sigma, flag="wholespace"): """ Analytic solution for electric potential from a postive pole @@ -32,3 +32,82 @@ def DCAnalytic(txloc, rxlocs, sigma, flag="wholespace"): return phi +deg2rad = lambda deg: deg/180.*np.pi +rad2deg = lambda rad: rad*180./np.pi + +def DCAnalyticSphere(txloc, rxloc, xc, radius, sigma, sigma1, \ + flag = "sec", order=12): +# def DCSpherePointCurrent(txloc, rxloc, xc, radius, rho, rho1, \ +# flag = "sec", order=12): + """ + + Parameters: + + txloc (array) : current electrode location (x,y,z) + xc (float) : x center of depressed sphere + rxloc (array) : electrode locations + (Nx3 array, # of electrodes) + radius (float): radius of the sphere (m) + rho (float) : resistivity of the background (ohm-m) + rho1 (float) : resistivity of the sphere + flag (string) : "sec", "total", "prim" + (default="sec") + "sec": secondary potential only due to sphere + "prim": primary potential from the point source + "total": "sec"+"prim" + order (float) : maximum order of Legendre polynomial + (default=12) + + Written by Seogi Kang (skang@eos.ubc.ca) + Ph.D. Candidate of University of British Columbia, Canada + + """ + + Pleg = [] + # Compute Legendre Polynomial + for i in range(order): + Pleg.append(special.legendre(i, monic=0)) + + + rho = 1./sigma + rho1 = 1./sigma1 + + # Center of the sphere should be aligned in txloc in y-direction + yc = txloc[1] + xyz = np.c_[rxloc[:,0]-xc, rxloc[:,1]-yc, rxloc[:,2]] + r = np.sqrt( (xyz**2).sum(axis=1) ) + + x0 = abs(txloc[0]-xc) + + costheta = xyz[:,0]/r * (txloc[0]-xc)/x0 + phi = np.zeros_like(r) + R = (r**2+x0**2.-2.*r*x0*costheta)**0.5 + # primary potential in a whole space + prim = rho*1./(4*np.pi*R) + + if flag =="prim": + return prim + + sphind = r < radius + out = np.zeros_like(r) + for n in range(order): + An, Bn = AnBnfun(n, radius, x0, rho, rho1) + dumout = An*r[~sphind]**(-n-1.)*Pleg[n](costheta[~sphind]) + out[~sphind] += dumout + dumin = Bn*r[sphind]**(n)*Pleg[n](costheta[sphind]) + out[sphind] += dumin + + out[~sphind] += prim[~sphind] + + if flag == "sec": + return out-prim + elif flag == "total": + return out + +def AnBnfun(n, radius, x0, rho, rho1, I=1.): + const = I*rho/(4*np.pi) + bunmo = n*rho + (n+1)*rho1 + An = const * radius**(2*n+1) / x0 ** (n+1.) * n * \ + (rho1-rho) / bunmo + Bn = const * 1. / x0 ** (n+1.) * (2*n+1) * (rho1) / bunmo + return An, Bn diff --git a/SimPEG/EM/Analytics/__init__.py b/SimPEG/EM/Analytics/__init__.py index d251f205..9df2aef7 100644 --- a/SimPEG/EM/Analytics/__init__.py +++ b/SimPEG/EM/Analytics/__init__.py @@ -1,4 +1,4 @@ from TDEM import hzAnalyticDipoleT from FDEM import hzAnalyticDipoleF from FDEMcasing import * -from DC import DCAnalytic +from DC import DCAnalyticHalf, DCAnalyticSphere diff --git a/SimPEG/EM/Static/DC/SrcDC.py b/SimPEG/EM/Static/DC/SrcDC.py index 1e64835e..4eda7fe2 100644 --- a/SimPEG/EM/Static/DC/SrcDC.py +++ b/SimPEG/EM/Static/DC/SrcDC.py @@ -1,6 +1,6 @@ import SimPEG # from SimPEG.EM.Base import BaseEMSurvey -from SimPEG.Utils import Zero, closestPoints +from SimPEG.Utils import Zero, closestPoints, mkvc import numpy as np class BaseSrc(SimPEG.Survey.BaseSrc): @@ -27,13 +27,13 @@ class Dipole(BaseSrc): def eval(self, prob): if prob._formulation == 'HJ': - inds = closestPoints(prob.mesh, self.loc) + inds = closestPoints(prob.mesh, self.loc, gridLoc='CC') q = np.zeros(prob.mesh.nC) q[inds] = self.current * np.r_[1., -1.] elif prob._formulation == 'EB': - inds = closestPoints(prob.mesh, self.loc) - q = np.zeros(prob.mesh.nN) - q[inds] = self.current * np.r_[1., -1.] + qa = prob.mesh.getInterpolationMat(self.loc[0], locType='N').todense() + qb = -prob.mesh.getInterpolationMat(self.loc[1], locType='N').todense() + q = mkvc(qa+qb) return q # def bc_contribution diff --git a/SimPEG/Mesh/View.py b/SimPEG/Mesh/View.py index 089d7d9a..6f009dc9 100644 --- a/SimPEG/Mesh/View.py +++ b/SimPEG/Mesh/View.py @@ -206,7 +206,7 @@ class TensorView(object): return out viewOpts = ['real','imag','abs','vec'] normalOpts = ['X', 'Y', 'Z'] - vTypeOpts = ['CC', 'CCv','F','E','Fx','Fy','Fz','E','Ex','Ey','Ez'] + vTypeOpts = ['CC', 'CCv','N','F','E','Fx','Fy','Fz','E','Ex','Ey','Ez'] # Some user error checking assert vType in vTypeOpts, "vType must be in ['%s']" % "','".join(vTypeOpts) diff --git a/tests/em/static/test_DC_analytic.py b/tests/em/static/test_DC_analytic.py new file mode 100644 index 00000000..3755c6ba --- /dev/null +++ b/tests/em/static/test_DC_analytic.py @@ -0,0 +1,68 @@ +import unittest +from SimPEG import Mesh, Utils, EM, Maps, np +import SimPEG.EM.Static.DC as DC + +class DCProblemAnalyticTests(unittest.TestCase): + + def setUp(self): + + cs = 25. + hx = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)] + hy = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)] + hz = [(cs,7, -1.3),(cs,20)] + mesh = Mesh.TensorMesh([hx, hy, hz],x0="CCN") + sigma = np.ones(mesh.nC)*1e-2 + + x = mesh.vectorCCx[(mesh.vectorCCx>-155.)&(mesh.vectorCCx<155.)] + y = mesh.vectorCCx[(mesh.vectorCCy>-155.)&(mesh.vectorCCy<155.)] + Aloc = np.r_[-200., 0., 0.] + Bloc = np.r_[200., 0., 0.] + M = Utils.ndgrid(x-25.,y, np.r_[0.]) + N = Utils.ndgrid(x+25.,y, np.r_[0.]) + phiA = EM.Analytics.DCAnalyticHalf(Aloc, [M,N], 1e-2, flag="halfspace") + phiB = EM.Analytics.DCAnalyticHalf(Bloc, [M,N], 1e-2, flag="halfspace") + data_anal = phiA-phiB + + rx = DC.Rx.Dipole(M, N) + src = DC.Src.Dipole([rx], Aloc, Bloc) + survey = DC.Survey([src]) + + self.survey = survey + self.mesh = mesh + self.sigma = sigma + self.data_anal = data_anal + + try: + from pymatsolver import MumpsSolver + self.Solver = MumpsSolver + except ImportError, e: + self.Solver = SolverLU + + def test_N(self): + problem = DC.Problem3D_N(self.mesh) + problem.Solver = self.Solver + problem.pair(self.survey) + data = self.survey.dpred(self.sigma) + err= np.linalg.norm(data-self.data_anal)/np.linalg.norm(self.data_anal) + if err < 0.2: + passed = True + else: + passed = False + self.assertTrue(passed) + + def test_CC(self): + problem = DC.Problem3D_N(self.mesh) + problem.Solver = self.Solver + problem.pair(self.survey) + data = self.survey.dpred(self.sigma) + err= np.linalg.norm(data-self.data_anal)/np.linalg.norm(self.data_anal) + if err < 0.2: + passed = True + print ">> DC analytic test for Problem3D_CC is pased" + else: + passed = False + self.assertTrue(passed) + +if __name__ == '__main__': + unittest.main() + diff --git a/tests/em/static/test_DC_deriv.py b/tests/em/static/test_DC_deriv.py deleted file mode 100644 index e69de29b..00000000 diff --git a/tests/em/static/test_DC.py b/tests/em/static/test_DC_jvecjtvecadj.py similarity index 100% rename from tests/em/static/test_DC.py rename to tests/em/static/test_DC_jvecjtvecadj.py From 6a064c5f96c106e54dc80412039f22ebf5b6d5d7 Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Mon, 25 Apr 2016 11:00:32 -0700 Subject: [PATCH 20/25] Minor changes --- tests/em/static/test_DC_analytic.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/tests/em/static/test_DC_analytic.py b/tests/em/static/test_DC_analytic.py index 3755c6ba..53d494e2 100644 --- a/tests/em/static/test_DC_analytic.py +++ b/tests/em/static/test_DC_analytic.py @@ -38,7 +38,7 @@ class DCProblemAnalyticTests(unittest.TestCase): except ImportError, e: self.Solver = SolverLU - def test_N(self): + def test_Problem3D_N(self): problem = DC.Problem3D_N(self.mesh) problem.Solver = self.Solver problem.pair(self.survey) @@ -46,21 +46,24 @@ class DCProblemAnalyticTests(unittest.TestCase): err= np.linalg.norm(data-self.data_anal)/np.linalg.norm(self.data_anal) if err < 0.2: passed = True + print ">> DC analytic test for Problem3D_N is passed" else: passed = False + print ">> DC analytic test for Problem3D_N is failed" self.assertTrue(passed) - def test_CC(self): - problem = DC.Problem3D_N(self.mesh) + def test_Problem3D_CC(self): + problem = DC.Problem3D_CC(self.mesh) problem.Solver = self.Solver problem.pair(self.survey) data = self.survey.dpred(self.sigma) err= np.linalg.norm(data-self.data_anal)/np.linalg.norm(self.data_anal) if err < 0.2: passed = True - print ">> DC analytic test for Problem3D_CC is pased" + print ">> DC analytic test for Problem3D_CC is passed" else: passed = False + print ">> DC analytic test for Problem3D_CC is failed" self.assertTrue(passed) if __name__ == '__main__': From 92e2fd67de0752a97806139c21c0cfc3e346644d Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Wed, 27 Apr 2016 15:25:40 -0700 Subject: [PATCH 21/25] minor fixes. --- SimPEG/EM/Analytics/DC.py | 12 +++++++++--- SimPEG/EM/Static/DC/SrcDC.py | 4 ++-- SimPEG/EM/Static/DC/__init__.py | 1 + 3 files changed, 12 insertions(+), 5 deletions(-) diff --git a/SimPEG/EM/Analytics/DC.py b/SimPEG/EM/Analytics/DC.py index 4fb03fd3..3949b2e5 100644 --- a/SimPEG/EM/Analytics/DC.py +++ b/SimPEG/EM/Analytics/DC.py @@ -1,5 +1,6 @@ import numpy as np from scipy.constants import mu_0, pi +from scipy import special def DCAnalyticHalf(txloc, rxlocs, sigma, flag="wholespace"): """ @@ -36,7 +37,7 @@ deg2rad = lambda deg: deg/180.*np.pi rad2deg = lambda rad: rad*180./np.pi def DCAnalyticSphere(txloc, rxloc, xc, radius, sigma, sigma1, \ - flag = "sec", order=12): + flag = "sec", order=12, halfspace=False): # def DCSpherePointCurrent(txloc, rxloc, xc, radius, rho, rho1, \ # flag = "sec", order=12): """ @@ -99,10 +100,15 @@ def DCAnalyticSphere(txloc, rxloc, xc, radius, sigma, sigma1, \ out[~sphind] += prim[~sphind] + if halfspace: + scale = 2 + else: + scale = 1 + if flag == "sec": - return out-prim + return scale*(out-prim) elif flag == "total": - return out + return scale*out def AnBnfun(n, radius, x0, rho, rho1, I=1.): const = I*rho/(4*np.pi) diff --git a/SimPEG/EM/Static/DC/SrcDC.py b/SimPEG/EM/Static/DC/SrcDC.py index 4eda7fe2..7879ba0c 100644 --- a/SimPEG/EM/Static/DC/SrcDC.py +++ b/SimPEG/EM/Static/DC/SrcDC.py @@ -5,7 +5,7 @@ import numpy as np class BaseSrc(SimPEG.Survey.BaseSrc): - current = 1 + current = 1.0 loc = None def __init__(self, rxList, **kwargs): @@ -33,7 +33,7 @@ class Dipole(BaseSrc): elif prob._formulation == 'EB': qa = prob.mesh.getInterpolationMat(self.loc[0], locType='N').todense() qb = -prob.mesh.getInterpolationMat(self.loc[1], locType='N').todense() - q = mkvc(qa+qb) + q = self.current * mkvc(qa+qb) return q # def bc_contribution diff --git a/SimPEG/EM/Static/DC/__init__.py b/SimPEG/EM/Static/DC/__init__.py index 57da57e8..82cf76d5 100644 --- a/SimPEG/EM/Static/DC/__init__.py +++ b/SimPEG/EM/Static/DC/__init__.py @@ -3,4 +3,5 @@ from SurveyDC import Survey import SrcDC as Src #Pole import RxDC as Rx from FieldsDC import Fields_CC +from BoundaryUtils import getxBCyBC_CC import Utils From d14cd444ac59b7594ded7b3635fc5f114f1d7bb4 Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Wed, 27 Apr 2016 23:04:28 -0700 Subject: [PATCH 22/25] working 2.5D fwd problem. --- SimPEG/EM/Static/DC/FieldsDC_2D.py | 106 +++++++++++++ SimPEG/EM/Static/DC/ProblemDC_2D.py | 235 ++++++++++++++++++++++++++++ SimPEG/EM/Static/DC/RxDC.py | 57 ++++++- SimPEG/EM/Static/DC/SrcDC.py | 2 - SimPEG/EM/Static/DC/SurveyDC.py | 22 ++- SimPEG/EM/Static/DC/__init__.py | 3 +- 6 files changed, 415 insertions(+), 10 deletions(-) create mode 100644 SimPEG/EM/Static/DC/FieldsDC_2D.py create mode 100644 SimPEG/EM/Static/DC/ProblemDC_2D.py diff --git a/SimPEG/EM/Static/DC/FieldsDC_2D.py b/SimPEG/EM/Static/DC/FieldsDC_2D.py new file mode 100644 index 00000000..77e3199e --- /dev/null +++ b/SimPEG/EM/Static/DC/FieldsDC_2D.py @@ -0,0 +1,106 @@ +import SimPEG +import Utils, numpy as np, scipy.sparse as sp + +class Fields_ky(SimPEG.Problem.TimeFields): + + """ + + Fancy Field Storage for a 2.5D code. + + u[:,'phi', kyInd] = phi + print u[src0,'phi'] + + Only one field type is stored for + each problem, the rest are computed. The fields obejct acts like an array and is indexed by + .. code-block:: python + f = problem.fields(m) + e = f[srcList,'e'] + j = f[srcList,'j'] + + If accessing all sources for a given field, use the :code:`:` + .. code-block:: python + f = problem.fields(m) + phi = f[:,'phi'] + e = f[:,'e'] + b = f[:,'b'] + The array returned will be size (nE or nF, nSrcs :math:`\\times` nFrequencies) + """ + + knownFields = {} + dtype = float + + def _phiDeriv(self,kyInd, 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(kyInd, src, v, adjoint=adjoint), self._phiDeriv_m(kyInd, src, v, adjoint=adjoint) + + return np.array(self._phiDeriv_u(kyInd, src, du_dm_v, adjoint) + self._phiDeriv_m(kyInd, src, v, adjoint), dtype = float) + + def _eDeriv(self,kyInd, 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(kyInd, src, v, adjoint), self._eDeriv_m(kyInd, src, v, adjoint) + return np.array(self._eDeriv_u(kyInd, src, du_dm_v, adjoint) + self._eDeriv_m(kyInd, src, v, adjoint), dtype = float) + + def _jDeriv(self,kyInd, 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(kyInd, src, v, adjoint), self._jDeriv_m(kyInd, src, v, adjoint) + return np.array(self._jDeriv_u(kyInd, src, du_dm_v, adjoint) + self._jDeriv_m(kyInd, src, v, adjoint), dtype = float) + + + # def _eDeriv(self, tInd, src, dun_dm_v, v, adjoint=False): + # if adjoint is True: + # return self._eDeriv_u(tInd, src, v, adjoint), self._eDeriv_m(tInd, src, v, adjoint) + # return self._eDeriv_u(tInd, src, dun_dm_v) + self._eDeriv_m(tInd, src, v) + + # def _bDeriv(self, tInd, src, dun_dm_v, v, adjoint=False): + # if adjoint is True: + # return self._bDeriv_u(tInd, src, v, adjoint), self._bDeriv_m(tInd, src, v, adjoint) + # return self._bDeriv_u(tInd, src, dun_dm_v) + self._bDeriv_m(tInd, src, v) + + +class Fields_ky_CC(Fields_ky): + 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_ky.__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, src, kyInd): + return phiSolution + + def _phiDeriv_u(self, kyInd, src, v, adjoint = False): + return Identity()*v + + def _phiDeriv_m(self, kyInd, src, v, adjoint = False): + return Zero() + + def _j(self, phiSolution, srcList): + raise NotImplementedError + + def _e(self, phiSolution, srcList): + raise NotImplementedError diff --git a/SimPEG/EM/Static/DC/ProblemDC_2D.py b/SimPEG/EM/Static/DC/ProblemDC_2D.py new file mode 100644 index 00000000..ee04a560 --- /dev/null +++ b/SimPEG/EM/Static/DC/ProblemDC_2D.py @@ -0,0 +1,235 @@ +from SimPEG import Problem, Utils +from SimPEG.EM.Base import BaseEMProblem +from SurveyDC import Survey, Survey_ky +from FieldsDC_2D import Fields_ky, Fields_ky_CC +from SimPEG.Utils import sdiag +import numpy as np +from SimPEG.Utils import Zero +from BoundaryUtils import getxBCyBC_CC + +class BaseDCProblem_2D(BaseEMProblem): + + surveyPair = Survey_ky + fieldsPair = Fields_ky + nky = 15 + ky = np.logspace(-4, 1, nky) + Ainv = [None for i in range(nky)] + nT = nky # Only for using TimeFields + + def fields(self, m): + self.curModel = m + + if not self.Ainv[0] == None: + for i in range(self.nky): + self.Ainv[i].clean() + + f = self.fieldsPair(self.mesh, self.survey) + Srcs = self.survey.srcList + for iky in range(self.nky): + ky = self.ky[iky] + A = self.getA(ky) + self.Ainv[iky] = self.Solver(A, **self.solverOpts) + RHS = self.getRHS(ky) + u = self.Ainv[iky] * RHS + f[Srcs, self._solutionType, iky] = u + return f + + # def Jvec(self, m, v, f=None): + + # if f is None: + # f = self.fields(m) + + # self.curModel = m + + # Jv = self.dataPair(self.survey) #same size as the data + + # A = self.getA() + + # for src in self.survey.srcList: + # u_src = f[src, self._solutionType] # solution vector + # dA_dm_v = self.getADeriv(u_src, v) + # dRHS_dm_v = self.getRHSDeriv(src, v) + # du_dm_v = self.Ainv * ( - dA_dm_v + dRHS_dm_v ) + + # for rx in src.rxList: + # df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None) + # df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False) + # Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v) + # return Utils.mkvc(Jv) + + # def Jtvec(self, m, v, f=None): + # if f is None: + # f = self.fields(m) + + # self.curModel = m + + # # Ensure v is a data object. + # if not isinstance(v, self.dataPair): + # v = self.dataPair(self.survey, v) + + # Jtv = np.zeros(m.size) + # AT = self.getA() + + + # for src in self.survey.srcList: + # u_src = f[src, self._solutionType] + # for rx in src.rxList: + # PTv = rx.evalDeriv(src, self.mesh, f, v[src, rx], adjoint=True) # wrt f, need possibility wrt m + # df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None) + # df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True) + + # ATinvdf_duT = self.Ainv * df_duT + + # dA_dmT = self.getADeriv(u_src, ATinvdf_duT, adjoint=True) + # dRHS_dmT = self.getRHSDeriv(src, ATinvdf_duT, adjoint=True) + # du_dmT = -dA_dmT + dRHS_dmT + # Jtv += df_dmT + du_dmT + + # return Utils.mkvc(Jtv) + + def getSourceTerm(self, ky): + """ + takes concept of source and turns it into a matrix + """ + """ + Evaluates the sources, and puts them in matrix form + + :rtype: (numpy.ndarray, numpy.ndarray) + :return: q (nC or nN, nSrc) + """ + + Srcs = self.survey.srcList + + if self._formulation is 'EB': + n = self.mesh.nN + # return NotImplementedError + + elif self._formulation is 'HJ': + n = self.mesh.nC + + q = np.zeros((n, len(Srcs))) + + for i, src in enumerate(Srcs): + q[:,i] = src.eval(self) + return q + +class Problem2D_CC(BaseDCProblem_2D): + + _solutionType = 'phiSolution' + _formulation = 'HJ' # CC potentials means J is on faces + fieldsPair = Fields_ky_CC + + def __init__(self, mesh, **kwargs): + BaseDCProblem_2D.__init__(self, mesh, **kwargs) + self.setBC() + + def getA(self, ky): + """ + + Make the A matrix for the cell centered DC resistivity problem + + A = D MfRhoI D^\\top V + + """ + + D = self.Div + G = self.Grad + vol = self.mesh.vol + # TODO: this won't work for full anisotropy + MfRhoI = self.MfRhoI + # Get resistivity rho + rho = self.curModel.rho + A = D * MfRhoI * G + Utils.sdiag(ky**2*vol/rho) + return A + + def getADeriv(self, ky, u, v, adjoint= False): + + D = self.Div + G = self.Grad + MfRhoIDeriv = self.MfRhoIDeriv + + if adjoint: + return(MfRhoIDeriv( G * u ).T) * ( D.T * v) + Utils.sdiag(ky**2*mesh.vol)*v + return D * ((MfRhoIDeriv( G * u )) * v) + Utils.sdiag(ky**2*mesh.vol)*v + + def getRHS(self, ky): + """ + RHS for the DC problem + + q + """ + + RHS = self.getSourceTerm(ky) + return RHS + + def getRHSDeriv(self, ky, src, v, adjoint=False): + """ + Derivative of the right hand side with respect to the model + """ + # TODO: add qDeriv for RHS depending on m + # qDeriv = src.evalDeriv(self, ky, adjoint=adjoint) + # return qDeriv + return Zero() + + def setBC(self): + if self.mesh.dim==3: + fxm,fxp,fym,fyp,fzm,fzp = self.mesh.faceBoundaryInd + gBFxm = self.mesh.gridFx[fxm,:] + gBFxp = self.mesh.gridFx[fxp,:] + gBFym = self.mesh.gridFy[fym,:] + gBFyp = self.mesh.gridFy[fyp,:] + gBFzm = self.mesh.gridFz[fzm,:] + gBFzp = self.mesh.gridFz[fzp,:] + + # Setup Mixed B.C (alpha, beta, gamma) + temp_xm, temp_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0]) + temp_ym, temp_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1]) + temp_zm, temp_zp = np.ones_like(gBFzm[:,2]), np.ones_like(gBFzp[:,2]) + + alpha_xm, alpha_xp = temp_xm*0., temp_xp*0. + alpha_ym, alpha_yp = temp_ym*0., temp_yp*0. + alpha_zm, alpha_zp = temp_zm*0., temp_zp*0. + + beta_xm, beta_xp = temp_xm, temp_xp + beta_ym, beta_yp = temp_ym, temp_yp + beta_zm, beta_zp = temp_zm, temp_zp + + gamma_xm, gamma_xp = temp_xm*0., temp_xp*0. + gamma_ym, gamma_yp = temp_ym*0., temp_yp*0. + gamma_zm, gamma_zp = temp_zm*0., temp_zp*0. + + alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp, alpha_zm, alpha_zp] + beta = [beta_xm, beta_xp, beta_ym, beta_yp, beta_zm, beta_zp] + gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp, gamma_zm, gamma_zp] + + elif self.mesh.dim==2: + + fxm,fxp,fym,fyp = self.mesh.faceBoundaryInd + gBFxm = self.mesh.gridFx[fxm,:] + gBFxp = self.mesh.gridFx[fxp,:] + gBFym = self.mesh.gridFy[fym,:] + gBFyp = self.mesh.gridFy[fyp,:] + + # Setup Mixed B.C (alpha, beta, gamma) + temp_xm, temp_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0]) + temp_ym, temp_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1]) + + alpha_xm, alpha_xp = temp_xm*0., temp_xp*0. + alpha_ym, alpha_yp = temp_ym*0., temp_yp*0. + + beta_xm, beta_xp = temp_xm, temp_xp + beta_ym, beta_yp = temp_ym, temp_yp + + gamma_xm, gamma_xp = temp_xm*0., temp_xp*0. + gamma_ym, gamma_yp = temp_ym*0., temp_yp*0. + + alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp] + beta = [beta_xm, beta_xp, beta_ym, beta_yp] + gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp] + + x_BC, y_BC = getxBCyBC_CC(self.mesh, alpha, beta, gamma) + V = self.Vol + self.Div = V * self.mesh.faceDiv + P_BC, B = self.mesh.getBCProjWF_simple() + M = B*self.mesh.aveCC2F + self.Grad = self.Div.T - P_BC*Utils.sdiag(y_BC)*M diff --git a/SimPEG/EM/Static/DC/RxDC.py b/SimPEG/EM/Static/DC/RxDC.py index 7d50392e..f7c7d352 100644 --- a/SimPEG/EM/Static/DC/RxDC.py +++ b/SimPEG/EM/Static/DC/RxDC.py @@ -1,5 +1,5 @@ import SimPEG -# from SimPEG.EM.Base import BaseEMSurvey +import numpy as np from SimPEG.Utils import Zero, closestPoints class BaseRx(SimPEG.Survey.BaseRx): @@ -43,9 +43,6 @@ class BaseRx(SimPEG.Survey.BaseRx): elif adjoint: return P.T*v - - - # DC.Rx.Dipole(locs) class Dipole(BaseRx): @@ -77,6 +74,56 @@ class Dipole(BaseRx): return P -# class Pole(BaseRx): +class Dipole_ky(BaseRx): + + def __init__(self, locsM, locsN, rxType = 'phi', **kwargs): + assert locsM.shape == locsN.shape, 'locsM and locsN need to be the same size' + locs = [locsM, locsN] + # We may not need this ... + BaseRx.__init__(self, locs, rxType) + + @property + def nD(self): + """Number of data in the receiver.""" + return self.locs[0].shape[0] + + # Not sure why ... + # return int(self.locs[0].size / 2) + + def getP(self, mesh, Gloc): + if mesh in self._Ps: + return self._Ps[mesh] + + P0 = mesh.getInterpolationMat(self.locs[0], Gloc) + P1 = mesh.getInterpolationMat(self.locs[1], Gloc) + P = P0 - P1 + if self.storeProjections: + self._Ps[mesh] = P + return P + + def eval(self, ky, src, mesh, f): + P = self.getP(mesh, self.projGLoc(f)) + Pf = P*f[src, self.projField,:] + return self.IntTrapezoidal(ky, Pf, y=0.) + + def evalDeriv(self, ky, src, mesh, f, v, adjoint=False): + P = self.getP(mesh, self.projGLoc(f)) + if not adjoint: + return P*v + elif adjoint: + return P.T*v + + def IntTrapezoidal(self, ky, Pf, y=0.): + phi = np.zeros(Pf.shape[0]) + nky = ky.size + dky = np.diff(ky) + dky = np.r_[dky[0], dky] + phi0 = Pf[:,0] + for iky in range(nky): + phi1 = 2./np.pi*Pf[:,iky]/2. + phi += phi1*dky[iky]/2.*np.cos(ky[iky]*y) + phi += phi0*dky[iky]/2.*np.cos(ky[iky]*y) + phi0 = phi1.copy() + return phi diff --git a/SimPEG/EM/Static/DC/SrcDC.py b/SimPEG/EM/Static/DC/SrcDC.py index 7879ba0c..60a53dff 100644 --- a/SimPEG/EM/Static/DC/SrcDC.py +++ b/SimPEG/EM/Static/DC/SrcDC.py @@ -57,5 +57,3 @@ class Pole(BaseSrc): q[inds] = self.current * np.r_[1.] return q - # def bc_contribution - diff --git a/SimPEG/EM/Static/DC/SurveyDC.py b/SimPEG/EM/Static/DC/SurveyDC.py index 3b631bef..62e2922a 100644 --- a/SimPEG/EM/Static/DC/SurveyDC.py +++ b/SimPEG/EM/Static/DC/SurveyDC.py @@ -1,6 +1,6 @@ import SimPEG from SimPEG.EM.Base import BaseEMSurvey -from SimPEG import sp +from SimPEG import sp, Survey from SimPEG.Utils import Zero, Identity from RxDC import BaseRx from SrcDC import BaseSrc @@ -13,6 +13,24 @@ class Survey(BaseEMSurvey): self.srcList = srcList BaseEMSurvey.__init__(self, srcList, **kwargs) +class Survey_ky(BaseEMSurvey): + rxPair = BaseRx + srcPair = BaseSrc + def __init__(self, srcList, **kwargs): + self.srcList = srcList + BaseEMSurvey.__init__(self, srcList, **kwargs) - + def eval(self, f): + """ + Project fields to receiver locations + :param Fields u: fields object + :rtype: numpy.ndarray + :return: data + """ + data = SimPEG.Survey.Data(self) + ky = self.prob.ky + for src in self.srcList: + for rx in src.rxList: + data[src, rx] = rx.eval(ky, src, self.mesh, f) + return data diff --git a/SimPEG/EM/Static/DC/__init__.py b/SimPEG/EM/Static/DC/__init__.py index 82cf76d5..8c790084 100644 --- a/SimPEG/EM/Static/DC/__init__.py +++ b/SimPEG/EM/Static/DC/__init__.py @@ -1,5 +1,6 @@ from ProblemDC import Problem3D_CC, Problem3D_N -from SurveyDC import Survey +from ProblemDC_2D import Problem2D_CC +from SurveyDC import Survey, Survey_ky import SrcDC as Src #Pole import RxDC as Rx from FieldsDC import Fields_CC From 0610289fdf049196b0a539593727920bd8b13060 Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Thu, 28 Apr 2016 11:18:37 -0700 Subject: [PATCH 23/25] Working Jvec for 2.5D DC code --- SimPEG/EM/Static/DC/FieldsDC_2D.py | 3 +- SimPEG/EM/Static/DC/ProblemDC_2D.py | 62 ++++++---- SimPEG/EM/Static/DC/RxDC.py | 18 +-- SimPEG/EM/Static/DC/SrcDC.py | 5 - SimPEG/EM/Static/DC/SurveyDC.py | 4 +- tests/em/static/test_DC_2D_jvecjtvecadj.py | 127 +++++++++++++++++++++ 6 files changed, 180 insertions(+), 39 deletions(-) create mode 100644 tests/em/static/test_DC_2D_jvecjtvecadj.py diff --git a/SimPEG/EM/Static/DC/FieldsDC_2D.py b/SimPEG/EM/Static/DC/FieldsDC_2D.py index 77e3199e..5b75031d 100644 --- a/SimPEG/EM/Static/DC/FieldsDC_2D.py +++ b/SimPEG/EM/Static/DC/FieldsDC_2D.py @@ -1,5 +1,6 @@ import SimPEG -import Utils, numpy as np, scipy.sparse as sp +from SimPEG.Utils import Identity, Zero +import numpy as np class Fields_ky(SimPEG.Problem.TimeFields): diff --git a/SimPEG/EM/Static/DC/ProblemDC_2D.py b/SimPEG/EM/Static/DC/ProblemDC_2D.py index ee04a560..8ebb9c67 100644 --- a/SimPEG/EM/Static/DC/ProblemDC_2D.py +++ b/SimPEG/EM/Static/DC/ProblemDC_2D.py @@ -12,7 +12,7 @@ class BaseDCProblem_2D(BaseEMProblem): surveyPair = Survey_ky fieldsPair = Fields_ky nky = 15 - ky = np.logspace(-4, 1, nky) + kys = np.logspace(-4, 1, nky) Ainv = [None for i in range(nky)] nT = nky # Only for using TimeFields @@ -26,7 +26,7 @@ class BaseDCProblem_2D(BaseEMProblem): f = self.fieldsPair(self.mesh, self.survey) Srcs = self.survey.srcList for iky in range(self.nky): - ky = self.ky[iky] + ky = self.kys[iky] A = self.getA(ky) self.Ainv[iky] = self.Solver(A, **self.solverOpts) RHS = self.getRHS(ky) @@ -34,28 +34,44 @@ class BaseDCProblem_2D(BaseEMProblem): f[Srcs, self._solutionType, iky] = u return f - # def Jvec(self, m, v, f=None): + def Jvec(self, m, v, f=None): - # if f is None: - # f = self.fields(m) + if f is None: + f = self.fields(m) - # self.curModel = m + self.curModel = m - # Jv = self.dataPair(self.survey) #same size as the data + Jv = self.dataPair(self.survey) #same size as the data + Jv0 = self.dataPair(self.survey) - # A = self.getA() + # Assume y=0. + # This needs some thoughts to implement in general when src is dipole + dky = np.diff(self.kys) + dky = np.r_[dky[0], dky] + y = 0. - # for src in self.survey.srcList: - # u_src = f[src, self._solutionType] # solution vector - # dA_dm_v = self.getADeriv(u_src, v) - # dRHS_dm_v = self.getRHSDeriv(src, v) - # du_dm_v = self.Ainv * ( - dA_dm_v + dRHS_dm_v ) - - # for rx in src.rxList: - # df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None) - # df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False) - # Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v) - # return Utils.mkvc(Jv) + for iky in range(self.nky): + ky = self.kys[iky] + A = self.getA(ky) + for src in self.survey.srcList: + u_src = f[src, self._solutionType, iky] # solution vector + dA_dm_v = self.getADeriv(ky, u_src, v) + dRHS_dm_v = self.getRHSDeriv(ky, src, v) + du_dm_v = self.Ainv[iky] * ( - dA_dm_v + dRHS_dm_v ) + for rx in src.rxList: + df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None) + df_dm_v = df_dmFun(iky, src, du_dm_v, v, adjoint=False) + # Trapezoidal intergration + Jv1_temp = 1./np.pi*rx.evalDeriv(ky, src, self.mesh, f, df_dm_v) + if iky==0: + #First assigment + Jv[src, rx] = Jv1_temp*dky[iky]*np.cos(ky*y) + else: + Jv[src, rx] += Jv1_temp*dky[iky] /2.*np.cos(ky*y) + Jv[src, rx] += Jv0[src, rx]*dky[iky]/2.*np.cos(ky*y) + Jv0[src, rx] = Jv1_temp.copy() + JV[iky,isrc,:] = Jv1_temp.copy() + return Utils.mkvc(Jv) # def Jtvec(self, m, v, f=None): # if f is None: @@ -146,11 +162,13 @@ class Problem2D_CC(BaseDCProblem_2D): D = self.Div G = self.Grad + vol = self.mesh.vol MfRhoIDeriv = self.MfRhoIDeriv - + rho = self.curModel.rho if adjoint: - return(MfRhoIDeriv( G * u ).T) * ( D.T * v) + Utils.sdiag(ky**2*mesh.vol)*v - return D * ((MfRhoIDeriv( G * u )) * v) + Utils.sdiag(ky**2*mesh.vol)*v + return(MfRhoIDeriv( G * u ).T) * ( D.T * v) + ky**2*Utils.sdiag(u.flatten()*vol*(-1./rho**2))*v + + return D * ((MfRhoIDeriv( G * u )) * v) + ky**2*Utils.sdiag(u.flatten()*vol*(-1./rho**2))*v def getRHS(self, ky): """ diff --git a/SimPEG/EM/Static/DC/RxDC.py b/SimPEG/EM/Static/DC/RxDC.py index f7c7d352..d2ec6098 100644 --- a/SimPEG/EM/Static/DC/RxDC.py +++ b/SimPEG/EM/Static/DC/RxDC.py @@ -102,10 +102,10 @@ class Dipole_ky(BaseRx): self._Ps[mesh] = P return P - def eval(self, ky, src, mesh, f): + def eval(self, kys, src, mesh, f): P = self.getP(mesh, self.projGLoc(f)) Pf = P*f[src, self.projField,:] - return self.IntTrapezoidal(ky, Pf, y=0.) + return self.IntTrapezoidal(kys, Pf, y=0.) def evalDeriv(self, ky, src, mesh, f, v, adjoint=False): P = self.getP(mesh, self.projGLoc(f)) @@ -114,16 +114,16 @@ class Dipole_ky(BaseRx): elif adjoint: return P.T*v - def IntTrapezoidal(self, ky, Pf, y=0.): + def IntTrapezoidal(self, kys, Pf, y=0.): phi = np.zeros(Pf.shape[0]) - nky = ky.size - dky = np.diff(ky) + nky = kys.size + dky = np.diff(kys) dky = np.r_[dky[0], dky] - phi0 = Pf[:,0] + phi0 = 1./np.pi*Pf[:,0] for iky in range(nky): - phi1 = 2./np.pi*Pf[:,iky]/2. - phi += phi1*dky[iky]/2.*np.cos(ky[iky]*y) - phi += phi0*dky[iky]/2.*np.cos(ky[iky]*y) + phi1 = 1./np.pi*Pf[:,iky] + phi += phi1*dky[iky]/2.*np.cos(kys[iky]*y) + phi += phi0*dky[iky]/2.*np.cos(kys[iky]*y) phi0 = phi1.copy() return phi diff --git a/SimPEG/EM/Static/DC/SrcDC.py b/SimPEG/EM/Static/DC/SrcDC.py index 60a53dff..5f4ac0e2 100644 --- a/SimPEG/EM/Static/DC/SrcDC.py +++ b/SimPEG/EM/Static/DC/SrcDC.py @@ -36,11 +36,6 @@ class Dipole(BaseSrc): q = self.current * mkvc(qa+qb) return q - # def bc_contribution - - -# How to treat boundary conditions here - class Pole(BaseSrc): def __init__(self, rxList, loc, **kwargs): diff --git a/SimPEG/EM/Static/DC/SurveyDC.py b/SimPEG/EM/Static/DC/SurveyDC.py index 62e2922a..fb3d49a7 100644 --- a/SimPEG/EM/Static/DC/SurveyDC.py +++ b/SimPEG/EM/Static/DC/SurveyDC.py @@ -29,8 +29,8 @@ class Survey_ky(BaseEMSurvey): :return: data """ data = SimPEG.Survey.Data(self) - ky = self.prob.ky + kys = self.prob.kys for src in self.srcList: for rx in src.rxList: - data[src, rx] = rx.eval(ky, src, self.mesh, f) + data[src, rx] = rx.eval(kys, src, self.mesh, f) return data diff --git a/tests/em/static/test_DC_2D_jvecjtvecadj.py b/tests/em/static/test_DC_2D_jvecjtvecadj.py new file mode 100644 index 00000000..ad7198e9 --- /dev/null +++ b/tests/em/static/test_DC_2D_jvecjtvecadj.py @@ -0,0 +1,127 @@ +import unittest +from SimPEG import * +import SimPEG.EM.Static.DC as DC + + +class DCProblem_2DTestsCC(unittest.TestCase): + + def setUp(self): + + cs = 12.5 + hx = [(cs,7, -1.3),(cs,61),(cs,7, 1.3)] + hy = [(cs,7, -1.3),(cs,20)] + mesh = Mesh.TensorMesh([hx, hy],x0="CN") + x = np.linspace(-135, 250., 20) + M = Utils.ndgrid(x-12.5, np.r_[0.]) + N = Utils.ndgrid(x+12.5, np.r_[0.]) + A0loc = np.r_[-150, 0.] + A1loc = np.r_[-130, 0.] + rxloc = [np.c_[M, np.zeros(20)], np.c_[N, np.zeros(20)]] + rx = DC.Rx.Dipole_ky(M, N) + src0 = DC.Src.Pole([rx], A0loc) + src1 = DC.Src.Pole([rx], A1loc) + survey = DC.Survey_ky([src0, src1]) + problem = DC.Problem2D_CC(mesh, mapping=[('rho', Maps.IdentityMap(mesh))]) + problem.pair(survey) + + mSynth = np.ones(mesh.nC) + survey.makeSyntheticData(mSynth) + + # Now set up the problem to do some minimization + dmis = DataMisfit.l2_DataMisfit(survey) + reg = Regularization.Tikhonov(mesh) + opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6) + invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4) + inv = Inversion.BaseInversion(invProb) + + self.inv = inv + self.reg = reg + self.p = problem + self.mesh = mesh + self.m0 = mSynth + self.survey = survey + self.dmis = dmis + + def test_misfit(self): + derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)] + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) + self.assertTrue(passed) + + # def test_adjoint(self): + # # Adjoint Test + # u = np.random.rand(self.mesh.nC*self.survey.nSrc) + # v = np.random.rand(self.mesh.nC) + # w = np.random.rand(self.survey.dobs.shape[0]) + # wtJv = w.dot(self.p.Jvec(self.m0, v)) + # vtJtw = v.dot(self.p.Jtvec(self.m0, w)) + # passed = np.abs(wtJv - vtJtw) < 1e-10 + # print 'Adjoint Test', np.abs(wtJv - vtJtw), passed + # self.assertTrue(passed) + + # def test_dataObj(self): + # derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)] + # passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) + # self.assertTrue(passed) + +# class DCProblemTestsN(unittest.TestCase): + +# def setUp(self): + +# aSpacing=2.5 +# nElecs=10 + +# surveySize = nElecs*aSpacing - aSpacing +# cs = surveySize/nElecs/4 + +# mesh = Mesh.TensorMesh([ +# [(cs,10, -1.3),(cs,surveySize/cs),(cs,10, 1.3)], +# [(cs,3, -1.3),(cs,3,1.3)], +# # [(cs,5, -1.3),(cs,10)] +# ],'CN') + +# srcList = DC.Utils.WennerSrcList(nElecs, aSpacing, in2D=True) +# survey = DC.Survey(srcList) +# problem = DC.Problem3D_N(mesh, mapping=[('rho', Maps.IdentityMap(mesh))]) +# problem.pair(survey) + +# mSynth = np.ones(mesh.nC) +# survey.makeSyntheticData(mSynth) + +# # Now set up the problem to do some minimization +# dmis = DataMisfit.l2_DataMisfit(survey) +# reg = Regularization.Tikhonov(mesh) +# opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6) +# invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4) +# inv = Inversion.BaseInversion(invProb) + +# self.inv = inv +# self.reg = reg +# self.p = problem +# self.mesh = mesh +# self.m0 = mSynth +# self.survey = survey +# self.dmis = dmis + +# def test_misfit(self): +# derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)] +# passed = Tests.checkDerivative(derChk, self.m0, plotIt=False) +# self.assertTrue(passed) + +# def test_adjoint(self): +# # Adjoint Test +# u = np.random.rand(self.mesh.nC*self.survey.nSrc) +# v = np.random.rand(self.mesh.nC) +# w = np.random.rand(self.survey.dobs.shape[0]) +# wtJv = w.dot(self.p.Jvec(self.m0, v)) +# vtJtw = v.dot(self.p.Jtvec(self.m0, w)) +# passed = np.abs(wtJv - vtJtw) < 1e-8 +# print 'Adjoint Test', np.abs(wtJv - vtJtw), passed +# self.assertTrue(passed) + +# def test_dataObj(self): +# derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)] +# passed = Tests.checkDerivative(derChk, self.m0, plotIt=False) +# self.assertTrue(passed) + +if __name__ == '__main__': + unittest.main() From ef602eaab1f03dfbfb5bba8b55170f4968f4dca4 Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Thu, 28 Apr 2016 18:13:18 -0700 Subject: [PATCH 24/25] working Jtvec --- SimPEG/EM/Static/DC/ProblemDC_2D.py | 64 ++++++++++++++-------- tests/em/static/test_DC_2D_jvecjtvecadj.py | 32 +++++------ 2 files changed, 57 insertions(+), 39 deletions(-) diff --git a/SimPEG/EM/Static/DC/ProblemDC_2D.py b/SimPEG/EM/Static/DC/ProblemDC_2D.py index 8ebb9c67..30c10d97 100644 --- a/SimPEG/EM/Static/DC/ProblemDC_2D.py +++ b/SimPEG/EM/Static/DC/ProblemDC_2D.py @@ -70,38 +70,57 @@ class BaseDCProblem_2D(BaseEMProblem): Jv[src, rx] += Jv1_temp*dky[iky] /2.*np.cos(ky*y) Jv[src, rx] += Jv0[src, rx]*dky[iky]/2.*np.cos(ky*y) Jv0[src, rx] = Jv1_temp.copy() - JV[iky,isrc,:] = Jv1_temp.copy() return Utils.mkvc(Jv) - # def Jtvec(self, m, v, f=None): - # if f is None: - # f = self.fields(m) + def Jtvec(self, m, v, f=None): + if f is None: + f = self.fields(m) - # self.curModel = m + self.curModel = m - # # Ensure v is a data object. - # if not isinstance(v, self.dataPair): - # v = self.dataPair(self.survey, v) + # Ensure v is a data object. + if not isinstance(v, self.dataPair): + v = self.dataPair(self.survey, v) - # Jtv = np.zeros(m.size) - # AT = self.getA() + Jtv = np.zeros(m.size) + + # Assume y=0. + # This needs some thoughts to implement in general when src is dipole + dky = np.diff(self.kys) + dky = np.r_[dky[0], dky] + y = 0. - # for src in self.survey.srcList: - # u_src = f[src, self._solutionType] - # for rx in src.rxList: - # PTv = rx.evalDeriv(src, self.mesh, f, v[src, rx], adjoint=True) # wrt f, need possibility wrt m - # df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None) - # df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True) + for src in self.survey.srcList: - # ATinvdf_duT = self.Ainv * df_duT + for rx in src.rxList: - # dA_dmT = self.getADeriv(u_src, ATinvdf_duT, adjoint=True) - # dRHS_dmT = self.getRHSDeriv(src, ATinvdf_duT, adjoint=True) - # du_dmT = -dA_dmT + dRHS_dmT - # Jtv += df_dmT + du_dmT + Jtv_temp1 = np.zeros(m.size) + Jtv_temp0 = np.zeros(m.size) - # return Utils.mkvc(Jtv) + for iky in range(self.nky): + u_src = f[src, self._solutionType, iky] + ky = self.kys[iky] + AT = self.getA(ky) + PTv = rx.evalDeriv(ky, src, self.mesh, f, v[src, rx], adjoint=True) # wrt f, need possibility wrt m + df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None) + df_duT, df_dmT = df_duTFun(iky, src, None, PTv, adjoint=True) + + ATinvdf_duT = self.Ainv[iky] * df_duT + + dA_dmT = self.getADeriv(ky, u_src, ATinvdf_duT, adjoint=True) + dRHS_dmT = self.getRHSDeriv(ky, src, ATinvdf_duT, adjoint=True) + du_dmT = -dA_dmT + dRHS_dmT + Jtv_temp1 = 1./np.pi*(df_dmT + du_dmT) + # Trapezoidal intergration + if iky==0: + #First assigment + Jtv += Jtv_temp1*dky[iky]*np.cos(ky*y) + else: + Jtv += Jtv_temp1*dky[iky]/2.*np.cos(ky*y) + Jtv += Jtv_temp0*dky[iky]/2.*np.cos(ky*y) + Jtv_temp0 = Jtv_temp1.copy() + return Utils.mkvc(Jtv) def getSourceTerm(self, ky): """ @@ -167,7 +186,6 @@ class Problem2D_CC(BaseDCProblem_2D): rho = self.curModel.rho if adjoint: return(MfRhoIDeriv( G * u ).T) * ( D.T * v) + ky**2*Utils.sdiag(u.flatten()*vol*(-1./rho**2))*v - return D * ((MfRhoIDeriv( G * u )) * v) + ky**2*Utils.sdiag(u.flatten()*vol*(-1./rho**2))*v def getRHS(self, ky): diff --git a/tests/em/static/test_DC_2D_jvecjtvecadj.py b/tests/em/static/test_DC_2D_jvecjtvecadj.py index ad7198e9..6beee640 100644 --- a/tests/em/static/test_DC_2D_jvecjtvecadj.py +++ b/tests/em/static/test_DC_2D_jvecjtvecadj.py @@ -24,14 +24,14 @@ class DCProblem_2DTestsCC(unittest.TestCase): problem = DC.Problem2D_CC(mesh, mapping=[('rho', Maps.IdentityMap(mesh))]) problem.pair(survey) - mSynth = np.ones(mesh.nC) + mSynth = np.ones(mesh.nC)*1. survey.makeSyntheticData(mSynth) # Now set up the problem to do some minimization dmis = DataMisfit.l2_DataMisfit(survey) reg = Regularization.Tikhonov(mesh) opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6) - invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4) + invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e0) inv = Inversion.BaseInversion(invProb) self.inv = inv @@ -47,21 +47,21 @@ class DCProblem_2DTestsCC(unittest.TestCase): passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) self.assertTrue(passed) - # def test_adjoint(self): - # # Adjoint Test - # u = np.random.rand(self.mesh.nC*self.survey.nSrc) - # v = np.random.rand(self.mesh.nC) - # w = np.random.rand(self.survey.dobs.shape[0]) - # wtJv = w.dot(self.p.Jvec(self.m0, v)) - # vtJtw = v.dot(self.p.Jtvec(self.m0, w)) - # passed = np.abs(wtJv - vtJtw) < 1e-10 - # print 'Adjoint Test', np.abs(wtJv - vtJtw), passed - # self.assertTrue(passed) + def test_adjoint(self): + # Adjoint Test + u = np.random.rand(self.mesh.nC*self.survey.nSrc) + v = np.random.rand(self.mesh.nC) + w = np.random.rand(self.survey.dobs.shape[0]) + wtJv = w.dot(self.p.Jvec(self.m0, v)) + vtJtw = v.dot(self.p.Jtvec(self.m0, w)) + passed = np.abs(wtJv - vtJtw) < 1e-10 + print 'Adjoint Test', np.abs(wtJv - vtJtw), passed + self.assertTrue(passed) - # def test_dataObj(self): - # derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)] - # passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) - # self.assertTrue(passed) + def test_dataObj(self): + derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)] + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) + self.assertTrue(passed) # class DCProblemTestsN(unittest.TestCase): From 38aef03f9df81a8660f08873eaf4f0269475a21a Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Fri, 29 Apr 2016 09:35:03 -0700 Subject: [PATCH 25/25] Working 2.5D fwd (nodal discretization) On going Jvec and Jtvec --- SimPEG/EM/Base.py | 2 - SimPEG/EM/Static/DC/FieldsDC_2D.py | 39 +++++ SimPEG/EM/Static/DC/ProblemDC.py | 4 - SimPEG/EM/Static/DC/ProblemDC_2D.py | 90 ++++++++++- SimPEG/EM/Static/DC/SrcDC.py | 5 +- SimPEG/EM/Static/DC/__init__.py | 2 +- tests/em/static/test_DC_2D_jvecjtvecadj.py | 166 ++++++++++----------- 7 files changed, 209 insertions(+), 99 deletions(-) diff --git a/SimPEG/EM/Base.py b/SimPEG/EM/Base.py index 491dcf71..496f5227 100644 --- a/SimPEG/EM/Base.py +++ b/SimPEG/EM/Base.py @@ -134,7 +134,6 @@ class BaseEMProblem(Problem.BaseProblem): """ return self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma)(u) * self.curModel.sigmaDeriv - @property def MeSigmaI(self): """ @@ -157,7 +156,6 @@ class BaseEMProblem(Problem.BaseProblem): return dMeSigmaI_dI * ( dMe_dsig * ( dsig_dm)) # return self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma, invMat=True)(u) - @property def MfRho(self): """ diff --git a/SimPEG/EM/Static/DC/FieldsDC_2D.py b/SimPEG/EM/Static/DC/FieldsDC_2D.py index 5b75031d..da1fbf97 100644 --- a/SimPEG/EM/Static/DC/FieldsDC_2D.py +++ b/SimPEG/EM/Static/DC/FieldsDC_2D.py @@ -105,3 +105,42 @@ class Fields_ky_CC(Fields_ky): def _e(self, phiSolution, srcList): raise NotImplementedError + +class Fields_ky_N(Fields_ky): + knownFields = {'phiSolution':'N'} + aliasFields = { + 'phi': ['phiSolution','N','_phi'], + 'j' : ['phiSolution','E','_j'], + 'e' : ['phiSolution','E','_e'], + } + # primary - secondary + # CC variables + + def __init__(self, mesh, survey, **kwargs): + Fields_ky.__init__(self, mesh, survey, **kwargs) + + def startup(self): + self.prob = self.survey.prob + + def _GLoc(self, fieldType): + if fieldType == 'phi': + return 'N' + elif fieldType == 'e' or fieldType == 'j': + return 'E' + else: + raise Exception('Field type must be phi, e, j') + + def _phi(self, phiSolution, src, kyInd): + return phiSolution + + def _phiDeriv_u(self, kyInd, src, v, adjoint = False): + return Identity()*v + + def _phiDeriv_m(self, kyInd, src, v, adjoint = False): + 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 index 9f05d876..0dd4f259 100644 --- a/SimPEG/EM/Static/DC/ProblemDC.py +++ b/SimPEG/EM/Static/DC/ProblemDC.py @@ -267,8 +267,6 @@ class Problem3D_N(BaseDCProblem): # Handling Null space of A A[0,0] = A[0,0] + 1. - # if self._makeASymmetric is True: - # return V.T * A return A def getADeriv(self, u, v, adjoint=False): @@ -293,8 +291,6 @@ class Problem3D_N(BaseDCProblem): """ RHS = self.getSourceTerm() - # if self._makeASymmetric is True: - # return self.Vol.T * RHS return RHS def getRHSDeriv(self, src, v, adjoint=False): diff --git a/SimPEG/EM/Static/DC/ProblemDC_2D.py b/SimPEG/EM/Static/DC/ProblemDC_2D.py index 30c10d97..2b5ed305 100644 --- a/SimPEG/EM/Static/DC/ProblemDC_2D.py +++ b/SimPEG/EM/Static/DC/ProblemDC_2D.py @@ -1,7 +1,7 @@ from SimPEG import Problem, Utils from SimPEG.EM.Base import BaseEMProblem from SurveyDC import Survey, Survey_ky -from FieldsDC_2D import Fields_ky, Fields_ky_CC +from FieldsDC_2D import Fields_ky, Fields_ky_CC, Fields_ky_N from SimPEG.Utils import sdiag import numpy as np from SimPEG.Utils import Zero @@ -90,14 +90,10 @@ class BaseDCProblem_2D(BaseEMProblem): dky = np.r_[dky[0], dky] y = 0. - for src in self.survey.srcList: - for rx in src.rxList: - Jtv_temp1 = np.zeros(m.size) Jtv_temp0 = np.zeros(m.size) - for iky in range(self.nky): u_src = f[src, self._solutionType, iky] ky = self.kys[iky] @@ -152,7 +148,7 @@ class Problem2D_CC(BaseDCProblem_2D): _solutionType = 'phiSolution' _formulation = 'HJ' # CC potentials means J is on faces - fieldsPair = Fields_ky_CC + fieldsPair = Fields_ky_N def __init__(self, mesh, **kwargs): BaseDCProblem_2D.__init__(self, mesh, **kwargs) @@ -269,3 +265,85 @@ class Problem2D_CC(BaseDCProblem_2D): P_BC, B = self.mesh.getBCProjWF_simple() M = B*self.mesh.aveCC2F self.Grad = self.Div.T - P_BC*Utils.sdiag(y_BC)*M + +class Problem2D_N(BaseDCProblem_2D): + + _solutionType = 'phiSolution' + _formulation = 'EB' # CC potentials means J is on faces + fieldsPair = Fields_ky_N + + def __init__(self, mesh, **kwargs): + BaseDCProblem_2D.__init__(self, mesh, **kwargs) + # self.setBC() + + @property + def MnSigma(self): + """ + Node inner product matrix for \\(\\sigma\\). Used in the E-B formulation + """ + # TODO: only works isotropic sigma + sigma = self.curModel.sigma + vol = self.mesh.vol + MnSigma = Utils.sdiag(self.mesh.aveN2CC.T*(Utils.sdiag(vol)*sigma)) + + return MnSigma + + def MnSigmaDeriv(self, u): + """ + Derivative of MnSigma with respect to the model + """ + sigma = self.curModel.sigma + sigmaderiv = self.curModel.sigmaDeriv + vol = self.mesh.vol + return Utils.sdiag(u)*self.mesh.aveN2CC.T*Utils.sdiag(vol) * self.curModel.sigmaDeriv + + def getA(self, ky): + """ + + 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 + MeSigma = self.MeSigma + MnSigma = self.MnSigma + Grad = self.mesh.nodalGrad + # Get conductivity sigma + sigma = self.curModel.sigma + A = Grad.T * MeSigma * Grad + ky**2*MnSigma + + # Handling Null space of A + A[0,0] = A[0,0] + 1. + return A + + def getADeriv(self, ky, u, v, adjoint= False): + + MeSigma = self.MeSigma + Grad = self.mesh.nodalGrad + sigma = self.curModel.sigma + vol = self.mesh.vol + + if adjoint: + return Grad.T*(self.MeSigmaDeriv(Grad*u)*v) + ky**2*self.MnSigmaDeriv(u)*v + return self.MeSigmaDeriv(Grad*u).T * (Grad*v) + ky**2*self.MnSigmaDeriv(u)*v + + def getRHS(self, ky): + """ + RHS for the DC problem + + q + """ + + RHS = self.getSourceTerm(ky) + return RHS + + def getRHSDeriv(self, ky, src, v, adjoint=False): + """ + Derivative of the right hand side with respect to the model + """ + # TODO: add qDeriv for RHS depending on m + # qDeriv = src.evalDeriv(self, ky, adjoint=adjoint) + # return qDeriv + return Zero() diff --git a/SimPEG/EM/Static/DC/SrcDC.py b/SimPEG/EM/Static/DC/SrcDC.py index 5f4ac0e2..02dae23e 100644 --- a/SimPEG/EM/Static/DC/SrcDC.py +++ b/SimPEG/EM/Static/DC/SrcDC.py @@ -47,8 +47,7 @@ class Pole(BaseSrc): q = np.zeros(prob.mesh.nC) q[inds] = self.current * np.r_[1.] elif prob._formulation == 'EB': - inds = closestPoints(prob.mesh, self.loc) - q = np.zeros(prob.mesh.nN) - q[inds] = self.current * np.r_[1.] + q = prob.mesh.getInterpolationMat(self.loc, locType='N').todense() + q = self.current * mkvc(q) return q diff --git a/SimPEG/EM/Static/DC/__init__.py b/SimPEG/EM/Static/DC/__init__.py index 8c790084..1080e391 100644 --- a/SimPEG/EM/Static/DC/__init__.py +++ b/SimPEG/EM/Static/DC/__init__.py @@ -1,5 +1,5 @@ from ProblemDC import Problem3D_CC, Problem3D_N -from ProblemDC_2D import Problem2D_CC +from ProblemDC_2D import Problem2D_CC, Problem2D_N from SurveyDC import Survey, Survey_ky import SrcDC as Src #Pole import RxDC as Rx diff --git a/tests/em/static/test_DC_2D_jvecjtvecadj.py b/tests/em/static/test_DC_2D_jvecjtvecadj.py index 6beee640..0740adc7 100644 --- a/tests/em/static/test_DC_2D_jvecjtvecadj.py +++ b/tests/em/static/test_DC_2D_jvecjtvecadj.py @@ -3,7 +3,67 @@ from SimPEG import * import SimPEG.EM.Static.DC as DC -class DCProblem_2DTestsCC(unittest.TestCase): +# class DCProblem_2DTestsCC(unittest.TestCase): + +# def setUp(self): + +# cs = 12.5 +# hx = [(cs,7, -1.3),(cs,61),(cs,7, 1.3)] +# hy = [(cs,7, -1.3),(cs,20)] +# mesh = Mesh.TensorMesh([hx, hy],x0="CN") +# x = np.linspace(-135, 250., 20) +# M = Utils.ndgrid(x-12.5, np.r_[0.]) +# N = Utils.ndgrid(x+12.5, np.r_[0.]) +# A0loc = np.r_[-150, 0.] +# A1loc = np.r_[-130, 0.] +# rxloc = [np.c_[M, np.zeros(20)], np.c_[N, np.zeros(20)]] +# rx = DC.Rx.Dipole_ky(M, N) +# src0 = DC.Src.Pole([rx], A0loc) +# src1 = DC.Src.Pole([rx], A1loc) +# survey = DC.Survey_ky([src0, src1]) +# problem = DC.Problem2D_CC(mesh, mapping=[('rho', Maps.IdentityMap(mesh))]) +# problem.pair(survey) + +# mSynth = np.ones(mesh.nC)*1. +# survey.makeSyntheticData(mSynth) + +# # Now set up the problem to do some minimization +# dmis = DataMisfit.l2_DataMisfit(survey) +# reg = Regularization.Tikhonov(mesh) +# opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6) +# invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e0) +# inv = Inversion.BaseInversion(invProb) + +# self.inv = inv +# self.reg = reg +# self.p = problem +# self.mesh = mesh +# self.m0 = mSynth +# self.survey = survey +# self.dmis = dmis + +# def test_misfit(self): +# derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)] +# passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) +# self.assertTrue(passed) + +# def test_adjoint(self): +# # Adjoint Test +# u = np.random.rand(self.mesh.nC*self.survey.nSrc) +# v = np.random.rand(self.mesh.nC) +# w = np.random.rand(self.survey.dobs.shape[0]) +# wtJv = w.dot(self.p.Jvec(self.m0, v)) +# vtJtw = v.dot(self.p.Jtvec(self.m0, w)) +# passed = np.abs(wtJv - vtJtw) < 1e-10 +# print 'Adjoint Test', np.abs(wtJv - vtJtw), passed +# self.assertTrue(passed) + +# def test_dataObj(self): +# derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)] +# passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) +# self.assertTrue(passed) + +class DCProblemTestsN(unittest.TestCase): def setUp(self): @@ -21,7 +81,7 @@ class DCProblem_2DTestsCC(unittest.TestCase): src0 = DC.Src.Pole([rx], A0loc) src1 = DC.Src.Pole([rx], A1loc) survey = DC.Survey_ky([src0, src1]) - problem = DC.Problem2D_CC(mesh, mapping=[('rho', Maps.IdentityMap(mesh))]) + problem = DC.Problem2D_N(mesh, mapping=[('rho', Maps.IdentityMap(mesh))]) problem.pair(survey) mSynth = np.ones(mesh.nC)*1. @@ -34,94 +94,34 @@ class DCProblem_2DTestsCC(unittest.TestCase): invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e0) inv = Inversion.BaseInversion(invProb) - self.inv = inv - self.reg = reg - self.p = problem - self.mesh = mesh - self.m0 = mSynth + self.inv = inv + self.reg = reg + self.p = problem + self.mesh = mesh + self.m0 = mSynth self.survey = survey - self.dmis = dmis + self.dmis = dmis def test_misfit(self): derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)] - passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False) self.assertTrue(passed) - def test_adjoint(self): - # Adjoint Test - u = np.random.rand(self.mesh.nC*self.survey.nSrc) - v = np.random.rand(self.mesh.nC) - w = np.random.rand(self.survey.dobs.shape[0]) - wtJv = w.dot(self.p.Jvec(self.m0, v)) - vtJtw = v.dot(self.p.Jtvec(self.m0, w)) - passed = np.abs(wtJv - vtJtw) < 1e-10 - print 'Adjoint Test', np.abs(wtJv - vtJtw), passed - self.assertTrue(passed) + # def test_adjoint(self): + # # Adjoint Test + # u = np.random.rand(self.mesh.nC*self.survey.nSrc) + # v = np.random.rand(self.mesh.nC) + # w = np.random.rand(self.survey.dobs.shape[0]) + # wtJv = w.dot(self.p.Jvec(self.m0, v)) + # vtJtw = v.dot(self.p.Jtvec(self.m0, w)) + # passed = np.abs(wtJv - vtJtw) < 1e-8 + # print 'Adjoint Test', np.abs(wtJv - vtJtw), passed + # self.assertTrue(passed) - def test_dataObj(self): - derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)] - passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) - self.assertTrue(passed) - -# class DCProblemTestsN(unittest.TestCase): - -# def setUp(self): - -# aSpacing=2.5 -# nElecs=10 - -# surveySize = nElecs*aSpacing - aSpacing -# cs = surveySize/nElecs/4 - -# mesh = Mesh.TensorMesh([ -# [(cs,10, -1.3),(cs,surveySize/cs),(cs,10, 1.3)], -# [(cs,3, -1.3),(cs,3,1.3)], -# # [(cs,5, -1.3),(cs,10)] -# ],'CN') - -# srcList = DC.Utils.WennerSrcList(nElecs, aSpacing, in2D=True) -# survey = DC.Survey(srcList) -# problem = DC.Problem3D_N(mesh, mapping=[('rho', Maps.IdentityMap(mesh))]) -# problem.pair(survey) - -# mSynth = np.ones(mesh.nC) -# survey.makeSyntheticData(mSynth) - -# # Now set up the problem to do some minimization -# dmis = DataMisfit.l2_DataMisfit(survey) -# reg = Regularization.Tikhonov(mesh) -# opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6) -# invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4) -# inv = Inversion.BaseInversion(invProb) - -# self.inv = inv -# self.reg = reg -# self.p = problem -# self.mesh = mesh -# self.m0 = mSynth -# self.survey = survey -# self.dmis = dmis - -# def test_misfit(self): -# derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)] -# passed = Tests.checkDerivative(derChk, self.m0, plotIt=False) -# self.assertTrue(passed) - -# def test_adjoint(self): -# # Adjoint Test -# u = np.random.rand(self.mesh.nC*self.survey.nSrc) -# v = np.random.rand(self.mesh.nC) -# w = np.random.rand(self.survey.dobs.shape[0]) -# wtJv = w.dot(self.p.Jvec(self.m0, v)) -# vtJtw = v.dot(self.p.Jtvec(self.m0, w)) -# passed = np.abs(wtJv - vtJtw) < 1e-8 -# print 'Adjoint Test', np.abs(wtJv - vtJtw), passed -# self.assertTrue(passed) - -# def test_dataObj(self): -# derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)] -# passed = Tests.checkDerivative(derChk, self.m0, plotIt=False) -# self.assertTrue(passed) + # def test_dataObj(self): + # derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)] + # passed = Tests.checkDerivative(derChk, self.m0, plotIt=False) + # self.assertTrue(passed) if __name__ == '__main__': unittest.main()