wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-07-08 04:02:43 +08:00
Code Issues Packages Projects Releases Wiki Activity

start of getting any field from any formulation

Browse Source
This commit is contained in:
Lindsey
2015-12-08 18:11:01 -08:00
parent 7c0be4ec50
commit 884d27b541
6 changed files with 235 additions and 76 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/EM/FDEM/FieldsFDEM.py
+94 -2
View File
@@ -3,7 +3,7 @@ import scipy.sparse as sp
import SimPEG
from SimPEG import Utils
from SimPEG.EM.Utils import omega
from SimPEG.Utils import Zero, Identity
from SimPEG.Utils import Zero, Identity, sdiag
class Fields(SimPEG.Problem.Fields):
@@ -19,7 +19,9 @@ class Fields_e(Fields):
'eSecondary' : ['eSolution','E','_eSecondary'],
'b' : ['eSolution','F','_b'],
'bPrimary' : ['eSolution','F','_bPrimary'],
'bSecondary' : ['eSolution','F','_bSecondary']
'bSecondary' : ['eSolution','F','_bSecondary'],
'j' : ['eSolution','CC','_j'],
'h' : ['eSolution','CC','_h'],
}
def __init__(self,mesh,survey,**kwargs):
@@ -28,6 +30,21 @@ class Fields_e(Fields):
def startup(self):
self.prob = self.survey.prob
self._edgeCurl = self.survey.prob.mesh.edgeCurl
self._aveE2CCV = self.survey.prob.mesh.aveE2CCV
self._aveF2CCV = self.survey.prob.mesh.aveF2CCV
self._sigma = self.survey.prob.curModel.sigma
self._sigmaDeriv = self.survey.prob.curModel.sigmaDeriv
self._nC = self.survey.prob.mesh.nC
def _GLoc(self,fieldType):
if fieldType == 'e':
return 'E'
elif fieldType == 'b':
return 'F'
elif (fieldType == 'h') or (fieldType == 'j'):
return 'CC'
else:
raise Exception('Field type must be e, b, h, j')
def _ePrimary(self, eSolution, srcList):
ePrimary = np.zeros_like(eSolution)
@@ -87,6 +104,41 @@ class Fields_e(Fields):
# Assuming the primary does not depend on the model
return self._bSecondaryDeriv_m(src, v, adjoint)
def _j(self, eSolution, srcList):
sigma = self._sigma
aveE2CCV = self._aveE2CCV
n = int(aveE2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy
VI = sdiag(1./np.kron(np.ones(n), self.prob.mesh.vol))
Sigma = sdiag(np.kron(np.ones(n), sigma))
e = self._e(eSolution, srcList)
return Sigma * (aveE2CCV * e)
def _jDeriv_u(self, src, v, adjoint=False):
raise NotImplementedError
sigma = self._sigma
aveE2CCV = self._aveE2CCV
n = int(aveE2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy
Sigma = sdiag(sp.kron(np.ones(n), sigma))
if not adjoint:
return Sigma * (aveE2CCV * (v + self._eDeriv_u(src, v, adjoint)))
return aveE2CCV.T * Sigma.T * v
def _jDeriv_m(self, src, v, adjoint=False):
raise NotImplementedError
sigma = self._sigma
aveE2CCV = self._aveE2CCV
n = int(aveE2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy
Sigma = sdiag(sp.kron(np.ones(n), sigma))
if not adjoint:
dsigma_dm = self._sigmaDeriv(v)
dSigma_dm = sdiag(sp.kron(np.ones(n), dsigma_dm))
class Fields_b(Fields):
knownFields = {'bSolution':'F'}
@@ -110,6 +162,16 @@ class Fields_b(Fields):
self._MeSigmaIDeriv = self.survey.prob.MeSigmaIDeriv
self._Me = self.survey.prob.Me
def _GLoc(self,fieldType):
if fieldType == 'e':
return 'E'
elif fieldType == 'b':
return 'F'
elif (fieldType == 'h') or (fieldType == 'j'):
return'CC'
else:
raise Exception('Field type must be e, b, h, j')
def _bPrimary(self, bSolution, srcList):
bPrimary = np.zeros_like(bSolution)
for i, src in enumerate(srcList):
@@ -193,6 +255,7 @@ class Fields_j(Fields):
'h' : ['jSolution','E','_h'],
'hPrimary' : ['jSolution','E','_hPrimary'],
'hSecondary' : ['jSolution','E','_hSecondary'],
'e' : ['jSolution','C','_e'],
}
def __init__(self,mesh,survey,**kwargs):
@@ -205,6 +268,19 @@ class Fields_j(Fields):
self._MfRho = self.survey.prob.MfRho
self._MfRhoDeriv = self.survey.prob.MfRhoDeriv
self._Me = self.survey.prob.Me
self._rho = self.survey.prob.curModel.rho
self._aveF2CCV = self.survey.prob.mesh.aveF2CCV
self._nC = self.survey.prob.mesh.nC
def _GLoc(self,fieldType):
if fieldType == 'h':
return 'E'
elif fieldType == 'j':
return 'F'
elif (fieldType == 'e') or (fieldType == 'b'):
return 'CC'
else:
raise Exception('Field type must be e, b, h, j')
def _jPrimary(self, jSolution, srcList):
jPrimary = np.zeros_like(jSolution,dtype = complex)
@@ -281,6 +357,22 @@ class Fields_j(Fields):
# assuming the primary doesn't depend on the model
return self._hSecondaryDeriv_m(src, v, adjoint)
def _e(self, jSolution, srcList):
rho = self._rho
aveF2CCV = self._aveF2CCV
n = int(aveF2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy
Rho = sdiag(np.kron(np.ones(n), rho))
j = self._j(jSolution, srcList)
return Rho * (aveF2CCV * j)
def _eDeriv_u(self, src, v, adjoint=False):
raise NotImplementedError
def _eDeriv_m(self, src, v, adjoint=False):
raise NotImplementedError
class Fields_h(Fields):
knownFields = {'hSolution':'E'}
SimPEG/EM/FDEM/SurveyFDEM.py
+63 -30
View File
@@ -3,6 +3,7 @@ from SimPEG.EM.Utils import *
from scipy.constants import mu_0
from SimPEG.Utils import Zero, Identity
import SrcFDEM as Src
from SimPEG import sp
####################################################
@@ -12,33 +13,33 @@ import SrcFDEM as Src
class Rx(SimPEG.Survey.BaseRx):
knownRxTypes = {
'exr':['e', 'Ex', 'real'],
'eyr':['e', 'Ey', 'real'],
'ezr':['e', 'Ez', 'real'],
'exi':['e', 'Ex', 'imag'],
'eyi':['e', 'Ey', 'imag'],
'ezi':['e', 'Ez', 'imag'],
'exr':['e', 'x', 'real'],
'eyr':['e', 'y', 'real'],
'ezr':['e', 'z', 'real'],
'exi':['e', 'x', 'imag'],
'eyi':['e', 'y', 'imag'],
'ezi':['e', 'z', 'imag'],
'bxr':['b', 'Fx', 'real'],
'byr':['b', 'Fy', 'real'],
'bzr':['b', 'Fz', 'real'],
'bxi':['b', 'Fx', 'imag'],
'byi':['b', 'Fy', 'imag'],
'bzi':['b', 'Fz', 'imag'],
'bxr':['b', 'x', 'real'],
'byr':['b', 'y', 'real'],
'bzr':['b', 'z', 'real'],
'bxi':['b', 'x', 'imag'],
'byi':['b', 'y', 'imag'],
'bzi':['b', 'z', 'imag'],
'jxr':['j', 'Fx', 'real'],
'jyr':['j', 'Fy', 'real'],
'jzr':['j', 'Fz', 'real'],
'jxi':['j', 'Fx', 'imag'],
'jyi':['j', 'Fy', 'imag'],
'jzi':['j', 'Fz', 'imag'],
'jxr':['j', 'x', 'real'],
'jyr':['j', 'y', 'real'],
'jzr':['j', 'z', 'real'],
'jxi':['j', 'x', 'imag'],
'jyi':['j', 'y', 'imag'],
'jzi':['j', 'z', 'imag'],
'hxr':['h', 'Ex', 'real'],
'hyr':['h', 'Ey', 'real'],
'hzr':['h', 'Ez', 'real'],
'hxi':['h', 'Ex', 'imag'],
'hyi':['h', 'Ey', 'imag'],
'hzi':['h', 'Ez', 'imag'],
'hxr':['h', 'x', 'real'],
'hyr':['h', 'y', 'real'],
'hzr':['h', 'z', 'real'],
'hxi':['h', 'x', 'imag'],
'hyi':['h', 'y', 'imag'],
'hzi':['h', 'z', 'imag'],
}
radius = None
@@ -50,10 +51,11 @@ class Rx(SimPEG.Survey.BaseRx):
"""Field Type projection (e.g. e b ...)"""
return self.knownRxTypes[self.rxType][0]
@property
def projGLoc(self):
"""Grid Location projection (e.g. Ex Fy ...)"""
return self.knownRxTypes[self.rxType][1]
# @property
# def projGLoc(self, u):
# """Grid Location projection (e.g. Ex Fy ...)"""
# return u._GLoc(self.rxType[0])
# return self.knownRxTypes[self.rxType][1]
@property
def projComp(self):
@@ -61,15 +63,46 @@ class Rx(SimPEG.Survey.BaseRx):
return self.knownRxTypes[self.rxType][2]
def projectFields(self, src, mesh, u):
P = self.getP(mesh)
u_part_complex = u[src, self.projField]
# get the real or imag component
real_or_imag = self.projComp
u_part = getattr(u_part_complex, real_or_imag)
projGLoc = u._GLoc(self.knownRxTypes[self.rxType][0])
if projGLoc == 'CC':
P = self.getP(mesh, projGLoc)
Z = 0.*P
if mesh.dim == 3:
if mesh._meshType == 'CYL' and mesh.isSymmetric and u_part.size > mesh.nC: # TODO: there must be a better way to do this!
if self.knownRxTypes[self.rxType][1] == 'x':
P = sp.hstack([P,Z])
elif self.knownRxTypes[self.rxType][1] == 'z':
P = sp.hstack([Z,P])
elif self.knownRxTypes[self.rxType][1] == 'y':
raise Exception('Symmetric CylMesh does not support y interpolation, as this variable does not exist.')
else:
if self.knownRxTypes[self.rxType][1] == 'x':
P = sp.hstack([P,Z,Z])
elif self.knownRxTypes[self.rxType][1] == 'y':
P = sp.hstack([Z,P,Z])
elif self.knownRxTypes[self.rxType][1] == 'z':
P = sp.hstack([Z,Z,P])
else:
projGLoc += self.knownRxTypes[self.rxType][1]
P = self.getP(mesh, projGLoc)
return P*u_part
def projectFieldsDeriv(self, src, mesh, u, v, adjoint=False):
P = self.getP(mesh)
projGLoc = u._GLoc(self.knownRxTypes[self.rxType][0])
if projGLoc != 'CC':
projGLoc += self.knownRxTypes[self.rxType][1]
P = self.getP(mesh, projGLoc)
if not adjoint:
Pv_complex = P * v
SimPEG/EM/Utils/testingUtils.py
+2 -2
View File
@@ -5,9 +5,9 @@ import sys
from scipy.constants import mu_0
def getFDEMProblem(fdemType, comp, SrcList, freq, verbose=False):
cs = 5.
cs = 10.
ncx, ncy, ncz = 6, 6, 6
npad = 3
npad = 5
hx = [(cs,npad,-1.3), (cs,ncx), (cs,npad,1.3)]
hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)]
hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)]
SimPEG/Mesh/TensorMesh.py
+10
View File
@@ -233,6 +233,10 @@ class BaseTensorMesh(BaseMesh):
'Fz' -> z-component of field defined on faces
'N' -> scalar field defined on nodes
'CC' -> scalar field defined on cell centers
# 'CCVx' -> x-component of a field defined on cell centers
# 'CCVy' -> y-component of a field defined on cell centers
# 'CCVz' -> z-component of a field defined on cell centers
"""
if self._meshType == 'CYL' and self.isSymmetric and locType in ['Ex','Ez','Fy']:
raise Exception('Symmetric CylMesh does not support %s interpolation, as this variable does not exist.' % locType)
@@ -256,6 +260,12 @@ class BaseTensorMesh(BaseMesh):
Q = sp.hstack(components)
elif locType in ['CC', 'N']:
Q = Utils.interpmat(loc, *self.getTensor(locType))
# elif locType in ['CCVx', 'CCVy', 'CCVz']:
# Q = Utils.interpmat(loc, 'CC')
# Zero = 0.*Q
else:
raise NotImplementedError('getInterpolationMat: locType=='+locType+' and mesh.dim=='+str(self.dim))
SimPEG/Survey.py
+5 -2
View File
@@ -35,7 +35,7 @@ class BaseRx(object):
"""Number of data in the receiver."""
return self.locs.shape[0]
def getP(self, mesh):
def getP(self, mesh, projGLoc=None):
"""
Returns the projection matrices as a
list for all components collected by
@@ -48,7 +48,10 @@ class BaseRx(object):
if mesh in self._Ps:
return self._Ps[mesh]
P = mesh.getInterpolationMat(self.locs, self.projGLoc)
if projGLoc is None:
projGLoc = self.projGLoc
P = mesh.getInterpolationMat(self.locs, projGLoc)
if self.storeProjections:
self._Ps[mesh] = P
return P
tests/em/fdem/forward/test_FDEM_forward.py
+61 -40
View File
@@ -7,26 +7,28 @@ from SimPEG.EM.Utils.testingUtils import getFDEMProblem
testEB = True
testHJ = True
testEJ = False
verbose = False
TOL = 1e-5
TOLEBHJ = 1e-5
TOLEJHB = 1e-1
FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order
CONDUCTIVITY = 1e1
MU = mu_0
freq = 1e-1
addrandoms = True
addrandoms = False
SrcList = ['RawVec', 'MagDipole_Bfield', 'MagDipole', 'CircularLoop']
def crossCheckTest(fdemType, comp):
def crossCheckTest(fdemType1, fdemType2, comp, TOL=TOLEBHJ):
l2norm = lambda r: np.sqrt(r.dot(r))
prb1 = getFDEMProblem(fdemType, comp, SrcList, freq, verbose)
prb1 = getFDEMProblem(fdemType1, comp, SrcList, freq, verbose)
mesh = prb1.mesh
print 'Cross Checking Forward: %s formulation - %s' % (fdemType, comp)
print 'Cross Checking Forward: %s formulation - %s' % (fdemType1, comp)
m = np.log(np.ones(mesh.nC)*CONDUCTIVITY)
mu = np.log(np.ones(mesh.nC)*MU)
@@ -42,16 +44,8 @@ def crossCheckTest(fdemType, comp):
if verbose:
print ' Problem 1 solved'
if fdemType == 'e':
prb2 = getFDEMProblem('b', comp, SrcList, freq, verbose)
elif fdemType == 'b':
prb2 = getFDEMProblem('e', comp, SrcList, freq, verbose)
elif fdemType == 'j':
prb2 = getFDEMProblem('h', comp, SrcList, freq, verbose)
elif fdemType == 'h':
prb2 = getFDEMProblem('j', comp, SrcList, freq, verbose)
else:
raise NotImplementedError()
prb2 = getFDEMProblem(fdemType2, comp, SrcList, freq, verbose)
# prb2.mu = mu
survey2 = prb2.survey
@@ -71,57 +65,84 @@ def crossCheckTest(fdemType, comp):
class FDEM_CrossCheck(unittest.TestCase):
if testEB:
def test_EB_CrossCheck_exr_Eform(self):
self.assertTrue(crossCheckTest('e', 'exr'))
self.assertTrue(crossCheckTest('e', 'b', 'exr'))
def test_EB_CrossCheck_eyr_Eform(self):
self.assertTrue(crossCheckTest('e', 'eyr'))
self.assertTrue(crossCheckTest('e', 'b', 'eyr'))
def test_EB_CrossCheck_ezr_Eform(self):
self.assertTrue(crossCheckTest('e', 'ezr'))
self.assertTrue(crossCheckTest('e', 'b', 'ezr'))
def test_EB_CrossCheck_exi_Eform(self):
self.assertTrue(crossCheckTest('e', 'exi'))
self.assertTrue(crossCheckTest('e', 'b', 'exi'))
def test_EB_CrossCheck_eyi_Eform(self):
self.assertTrue(crossCheckTest('e', 'eyi'))
self.assertTrue(crossCheckTest('e', 'b', 'eyi'))
def test_EB_CrossCheck_ezi_Eform(self):
self.assertTrue(crossCheckTest('e', 'ezi'))
self.assertTrue(crossCheckTest('e', 'b', 'ezi'))
def test_EB_CrossCheck_bxr_Eform(self):
self.assertTrue(crossCheckTest('e', 'bxr'))
self.assertTrue(crossCheckTest('e', 'b', 'bxr'))
def test_EB_CrossCheck_byr_Eform(self):
self.assertTrue(crossCheckTest('e', 'byr'))
self.assertTrue(crossCheckTest('e', 'b', 'byr'))
def test_EB_CrossCheck_bzr_Eform(self):
self.assertTrue(crossCheckTest('e', 'bzr'))
self.assertTrue(crossCheckTest('e', 'b', 'bzr'))
def test_EB_CrossCheck_bxi_Eform(self):
self.assertTrue(crossCheckTest('e', 'bxi'))
self.assertTrue(crossCheckTest('e', 'b', 'bxi'))
def test_EB_CrossCheck_byi_Eform(self):
self.assertTrue(crossCheckTest('e', 'byi'))
self.assertTrue(crossCheckTest('e', 'b', 'byi'))
def test_EB_CrossCheck_bzi_Eform(self):
self.assertTrue(crossCheckTest('e', 'bzi'))
self.assertTrue(crossCheckTest('e', 'b', 'bzi'))
if testHJ:
def test_HJ_CrossCheck_jxr_Jform(self):
self.assertTrue(crossCheckTest('j', 'jxr'))
self.assertTrue(crossCheckTest('j', 'h', 'jxr'))
def test_HJ_CrossCheck_jyr_Jform(self):
self.assertTrue(crossCheckTest('j', 'jyr'))
self.assertTrue(crossCheckTest('j', 'h', 'jyr'))
def test_HJ_CrossCheck_jzr_Jform(self):
self.assertTrue(crossCheckTest('j', 'jzr'))
self.assertTrue(crossCheckTest('j', 'h', 'jzr'))
def test_HJ_CrossCheck_jxi_Jform(self):
self.assertTrue(crossCheckTest('j', 'jxi'))
self.assertTrue(crossCheckTest('j', 'h', 'jxi'))
def test_HJ_CrossCheck_jyi_Jform(self):
self.assertTrue(crossCheckTest('j', 'jyi'))
self.assertTrue(crossCheckTest('j', 'h', 'jyi'))
def test_HJ_CrossCheck_jzi_Jform(self):
self.assertTrue(crossCheckTest('j', 'jzi'))
self.assertTrue(crossCheckTest('j', 'h', 'jzi'))
def test_HJ_CrossCheck_hxr_Jform(self):
self.assertTrue(crossCheckTest('j', 'hxr'))
self.assertTrue(crossCheckTest('j', 'h', 'hxr'))
def test_HJ_CrossCheck_hyr_Jform(self):
self.assertTrue(crossCheckTest('j', 'hyr'))
self.assertTrue(crossCheckTest('j', 'h', 'hyr'))
def test_HJ_CrossCheck_hzr_Jform(self):
self.assertTrue(crossCheckTest('j', 'hzr'))
self.assertTrue(crossCheckTest('j', 'h', 'hzr'))
def test_HJ_CrossCheck_hxi_Jform(self):
self.assertTrue(crossCheckTest('j', 'hxi'))
self.assertTrue(crossCheckTest('j', 'h', 'hxi'))
def test_HJ_CrossCheck_hyi_Jform(self):
self.assertTrue(crossCheckTest('j', 'hyi'))
self.assertTrue(crossCheckTest('j', 'h', 'hyi'))
def test_HJ_CrossCheck_hzi_Jform(self):
self.assertTrue(crossCheckTest('j', 'hzi'))
self.assertTrue(crossCheckTest('j', 'h', 'hzi'))
if testEJ:
# def test_EJ_CrossCheck_jxr_Jform(self):
# self.assertTrue(crossCheckTest('e', 'j', 'jxr'))
# def test_EJ_CrossCheck_jyr_Jform(self):
# self.assertTrue(crossCheckTest('e', 'j', 'jyr'))
# def test_EJ_CrossCheck_jzr_Jform(self):
# self.assertTrue(crossCheckTest('e', 'j', 'jzr'))
# def test_EJ_CrossCheck_jxi_Jform(self):
# self.assertTrue(crossCheckTest('e', 'j', 'jxi'))
# def test_EJ_CrossCheck_jyi_Jform(self):
# self.assertTrue(crossCheckTest('e', 'j', 'jyi'))
# def test_EJ_CrossCheck_jzi_Jform(self):
# self.assertTrue(crossCheckTest('e', 'j', 'jzi'))
def test_EJ_CrossCheck_jxr_Jform(self):
self.assertTrue(crossCheckTest('e', 'j', 'exr', TOL=TOLEJHB))
# def test_EJ_CrossCheck_jyr_Jform(self):
# self.assertTrue(crossCheckTest('e', 'j', 'eyr'))
# def test_EJ_CrossCheck_jzr_Jform(self):
# self.assertTrue(crossCheckTest('e', 'j', 'ezr'))
# def test_EJ_CrossCheck_jxi_Jform(self):
# self.assertTrue(crossCheckTest('e', 'j', 'exi'))
# def test_EJ_CrossCheck_jyi_Jform(self):
# self.assertTrue(crossCheckTest('e', 'j', 'eyi'))
# def test_EJ_CrossCheck_jzi_Jform(self):
# self.assertTrue(crossCheckTest('e', 'j', 'ezi'))
if __name__ == '__main__':
unittest.main()