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()