From 64b94861a01b2fcfe3da7a2c76df6ac97ff0b41c Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Fri, 22 Apr 2016 23:09:31 -0700 Subject: [PATCH] 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)