diff --git a/SimPEG/EM/Analytics/DC.py b/SimPEG/EM/Analytics/DC.py new file mode 100644 index 00000000..3949b2e5 --- /dev/null +++ b/SimPEG/EM/Analytics/DC.py @@ -0,0 +1,119 @@ +import numpy as np +from scipy.constants import mu_0, pi +from scipy import special + +def DCAnalyticHalf(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 + +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, halfspace=False): +# 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 halfspace: + scale = 2 + else: + scale = 1 + + if flag == "sec": + return scale*(out-prim) + elif flag == "total": + return scale*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 5b7a8851..9df2aef7 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 DCAnalyticHalf, DCAnalyticSphere diff --git a/SimPEG/EM/Base.py b/SimPEG/EM/Base.py index a16cdb91..496f5227 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 @@ -127,7 +134,6 @@ class BaseEMProblem(Problem.BaseProblem): """ return self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma)(u) * self.curModel.sigmaDeriv - @property def MeSigmaI(self): """ @@ -150,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): """ @@ -183,7 +188,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.getFaceInnerProductDeriv(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): @@ -192,7 +203,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 +213,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/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/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/BoundaryUtils.py b/SimPEG/EM/Static/DC/BoundaryUtils.py new file mode 100644 index 00000000..3967eb46 --- /dev/null +++ b/SimPEG/EM/Static/DC/BoundaryUtils.py @@ -0,0 +1,160 @@ +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] + 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) + 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") + + fxm,fxp,fym,fyp = mesh.faceBoundaryInd + 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.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) + 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 = 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[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.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) + 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/FieldsDC.py b/SimPEG/EM/Static/DC/FieldsDC.py new file mode 100644 index 00000000..9999c56a --- /dev/null +++ b/SimPEG/EM/Static/DC/FieldsDC.py @@ -0,0 +1,111 @@ +import SimPEG +from SimPEG.Utils import Identity, Zero +import numpy as np + +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=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): + 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 = float) + + 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 = float) + + +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(self, src, v, adjoint = False): + return Identity()*v + + def _phiDeriv_m(self, src, v, adjoint = False): + return Zero() + + def _j(self, phiSolution, srcList): + raise NotImplementedError + + 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(self, src, v, adjoint = False): + return Identity()*v + + def _phiDeriv_m(self, 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/FieldsDC_2D.py b/SimPEG/EM/Static/DC/FieldsDC_2D.py new file mode 100644 index 00000000..da1fbf97 --- /dev/null +++ b/SimPEG/EM/Static/DC/FieldsDC_2D.py @@ -0,0 +1,146 @@ +import SimPEG +from SimPEG.Utils import Identity, Zero +import numpy as np + +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 + +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 new file mode 100644 index 00000000..0dd4f259 --- /dev/null +++ b/SimPEG/EM/Static/DC/ProblemDC.py @@ -0,0 +1,307 @@ +from SimPEG import Problem, Utils +from SimPEG.EM.Base import BaseEMProblem +from SurveyDC import Survey +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): + + 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) + 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): + + 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): + """ + 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 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) + self.setBC() + + def getA(self): + """ + + Make the A matrix for the cell centered DC resistivity problem + + A = D MfRhoI D^\\top V + + """ + + D = self.Div + G = self.Grad + # TODO: this won't work for full anisotropy + MfRhoI = self.MfRhoI + A = D * MfRhoI * G + + # I think we should deprecate this for DC problem. + # if self._makeASymmetric is True: + # return V.T * A + return A + + def getADeriv(self, u, v, adjoint= False): + + D = self.Div + G = self.Grad + MfRhoIDeriv = self.MfRhoIDeriv + + if adjoint: + # if self._makeASymmetric is True: + # v = V * 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( G * u ) * v) + + def getRHS(self): + """ + RHS for the DC problem + + q + """ + + RHS = self.getSourceTerm() + + # I think we should deprecate this for DC problem. + # 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() + + 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): + + _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. + + 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() + 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() + + + + diff --git a/SimPEG/EM/Static/DC/ProblemDC_2D.py b/SimPEG/EM/Static/DC/ProblemDC_2D.py new file mode 100644 index 00000000..2b5ed305 --- /dev/null +++ b/SimPEG/EM/Static/DC/ProblemDC_2D.py @@ -0,0 +1,349 @@ +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, Fields_ky_N +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 + kys = 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.kys[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 + Jv0 = self.dataPair(self.survey) + + # 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 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() + 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) + + # 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: + 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] + 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): + """ + 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_N + + 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 + vol = self.mesh.vol + MfRhoIDeriv = self.MfRhoIDeriv + 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): + """ + 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 + +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/RxDC.py b/SimPEG/EM/Static/DC/RxDC.py new file mode 100644 index 00000000..d2ec6098 --- /dev/null +++ b/SimPEG/EM/Static/DC/RxDC.py @@ -0,0 +1,129 @@ +import SimPEG +import numpy as np +from SimPEG.Utils import Zero, closestPoints + +class BaseRx(SimPEG.Survey.BaseRx): + locs = 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, locs, rxType, **kwargs): + SimPEG.Survey.BaseRx.__init__(self, 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) + comp + return f._GLoc(self.rxType) + + def eval(self, src, mesh, f): + P = self.getP(mesh, self.projGLoc(f)) + return P*f[src, self.projField] + + def evalDeriv(self, 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 + +# 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' + 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 + + +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, kys, src, mesh, f): + P = self.getP(mesh, self.projGLoc(f)) + Pf = P*f[src, self.projField,:] + return self.IntTrapezoidal(kys, 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, kys, Pf, y=0.): + phi = np.zeros(Pf.shape[0]) + nky = kys.size + dky = np.diff(kys) + dky = np.r_[dky[0], dky] + phi0 = 1./np.pi*Pf[:,0] + for iky in range(nky): + 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 new file mode 100644 index 00000000..02dae23e --- /dev/null +++ b/SimPEG/EM/Static/DC/SrcDC.py @@ -0,0 +1,53 @@ +import SimPEG +# from SimPEG.EM.Base import BaseEMSurvey +from SimPEG.Utils import Zero, closestPoints, mkvc +import numpy as np + +class BaseSrc(SimPEG.Survey.BaseSrc): + + current = 1.0 + loc = None + + def __init__(self, rxList, **kwargs): + SimPEG.Survey.BaseSrc.__init__(self, rxList, **kwargs) + + def eval(self, prob): + raise NotImplementedError + + def evalDeriv(self, prob): + return 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, gridLoc='CC') + q = np.zeros(prob.mesh.nC) + q[inds] = self.current * np.r_[1., -1.] + 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 = self.current * mkvc(qa+qb) + return q + +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': + q = prob.mesh.getInterpolationMat(self.loc, locType='N').todense() + q = self.current * mkvc(q) + return q + diff --git a/SimPEG/EM/Static/DC/SurveyDC.py b/SimPEG/EM/Static/DC/SurveyDC.py new file mode 100644 index 00000000..fb3d49a7 --- /dev/null +++ b/SimPEG/EM/Static/DC/SurveyDC.py @@ -0,0 +1,36 @@ +import SimPEG +from SimPEG.EM.Base import BaseEMSurvey +from SimPEG import sp, Survey +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) + +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) + kys = self.prob.kys + for src in self.srcList: + for rx in src.rxList: + data[src, rx] = rx.eval(kys, src, self.mesh, f) + return data diff --git a/SimPEG/EM/Static/DC/Utils.py b/SimPEG/EM/Static/DC/Utils.py new file mode 100644 index 00000000..626ae368 --- /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).reshape([1,-1]),getLoc(i,2).reshape([1,-1])) + 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 new file mode 100644 index 00000000..1080e391 --- /dev/null +++ b/SimPEG/EM/Static/DC/__init__.py @@ -0,0 +1,8 @@ +from ProblemDC import Problem3D_CC, Problem3D_N +from ProblemDC_2D import Problem2D_CC, Problem2D_N +from SurveyDC import Survey, Survey_ky +import SrcDC as Src #Pole +import RxDC as Rx +from FieldsDC import Fields_CC +from BoundaryUtils import getxBCyBC_CC +import Utils 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 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 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) diff --git a/SimPEG/Mesh/View.py b/SimPEG/Mesh/View.py index 8eb22098..569e0392 100644 --- a/SimPEG/Mesh/View.py +++ b/SimPEG/Mesh/View.py @@ -218,7 +218,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/__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_2D_jvecjtvecadj.py b/tests/em/static/test_DC_2D_jvecjtvecadj.py new file mode 100644 index 00000000..0740adc7 --- /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)*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): + + 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_N(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) + 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() diff --git a/tests/em/static/test_DC_analytic.py b/tests/em/static/test_DC_analytic.py new file mode 100644 index 00000000..53d494e2 --- /dev/null +++ b/tests/em/static/test_DC_analytic.py @@ -0,0 +1,71 @@ +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_Problem3D_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 + 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_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 passed" + else: + passed = False + print ">> DC analytic test for Problem3D_CC is failed" + self.assertTrue(passed) + +if __name__ == '__main__': + unittest.main() + diff --git a/tests/em/static/test_DC_jvecjtvecadj.py b/tests/em/static/test_DC_jvecjtvecadj.py new file mode 100644 index 00000000..227d4614 --- /dev/null +++ b/tests/em/static/test_DC_jvecjtvecadj.py @@ -0,0 +1,127 @@ +import unittest +from SimPEG import * +import SimPEG.EM.Static.DC as DC + + +class DCProblemTestsCC(unittest.TestCase): + + def setUp(self): + + aSpacing=2.5 + nElecs=5 + + 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, 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() diff --git a/tests/mesh/test_Mixed_boundaryPoisson.py b/tests/mesh/test_Mixed_boundaryPoisson.py new file mode 100644 index 00000000..3aa1dbcd --- /dev/null +++ b/tests/mesh/test_Mixed_boundaryPoisson.py @@ -0,0 +1,411 @@ +import numpy as np +import scipy.sparse as sp +import unittest +import matplotlib.pyplot as plt +from SimPEG import * + +MESHTYPES = ['uniformTensorMesh'] + +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] + 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) + 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") + + fxm,fxp,fym,fyp = mesh.faceBoundaryInd + 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.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) + 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 = 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[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.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) + 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_CC(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 + # 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 + + 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_CC(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[:,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) + 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_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) + 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_CC(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()