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

Change Data to Survey

Browse Source
This commit is contained in:
rowanc1
2014-03-07 14:01:46 -08:00
parent e074cca07c
commit a14d771515
10 changed files with 60 additions and 60 deletions
Download Patch File Download Diff File Expand all files Collapse all files
simpegEM/FDEM/FDEM.py
+15 -15
View File
@@ -1,7 +1,7 @@
from SimPEG import Problem, Solver, Utils, np, sp
from scipy.constants import mu_0
from FieldsFDEM import FieldsFDEM
from DataFDEM import DataFDEM
from SurveyFDEM import SurveyFDEM
class ProblemFDEM_e(Problem.BaseProblem):
@@ -14,21 +14,21 @@ class ProblemFDEM_e(Problem.BaseProblem):
\dcurl E + i \omega B = 0 \\\\
\dcurl^\\top \MfMui B - \MeSig E = \Me \j_s
"""
def __init__(self, mesh, model, **kwargs):
Problem.BaseProblem.__init__(self, mesh, model, **kwargs)
def __init__(self, model, **kwargs):
Problem.BaseProblem.__init__(self, model, **kwargs)
solType = 'b'
storeTheseFields = 'e'
dataPair = DataFDEM
surveyPair = SurveyFDEM
solveOpts = {'factorize':False, 'backend':'scipy'}
j_s = None
def getFieldsObject(self):
return FieldsFDEM(self.mesh, self.data.nTx,
self.data.nFreq, store=self.storeTheseFields)
return FieldsFDEM(self.mesh, self.survey.nTx,
self.survey.nFreq, store=self.storeTheseFields)
####################################################
@@ -67,11 +67,11 @@ class ProblemFDEM_e(Problem.BaseProblem):
:rtype: scipy.sparse.csr_matrix
:return: A
"""
omega = self.data.omega[freqInd]
omega = self.survey.omega[freqInd]
return self.mesh.edgeCurl.T*self.MfMui*self.mesh.edgeCurl + 1j*omega*self.MeSigma
def getRHS(self, freqInd):
omega = self.data.omega[freqInd]
omega = self.survey.omega[freqInd]
#TODO: this needs to also depend on your transmitter!
return -1j*omega*self.Me*self.j_s
@@ -83,13 +83,13 @@ class ProblemFDEM_e(Problem.BaseProblem):
F = self.getFieldsObject()
for freqInd in range(self.data.nFreq):
for freqInd in range(self.survey.nFreq):
A = self.getA(freqInd)
b = self.getRHS(freqInd)
e = Solver(A, options=self.solveOpts).solve(b)
F.set_e(e, freqInd)
omega = self.data.omega[freqInd]
omega = self.survey.omega[freqInd]
#TODO: check if mass matrices needed:
b = -1./(1j*omega)*self.mesh.edgeCurl*e
F.set_b(b, freqInd)
@@ -103,13 +103,13 @@ class ProblemFDEM_e(Problem.BaseProblem):
if u is None:
u = self.fields(m)
Jvs = range(self.data.nFreq)
P = self.data.projectFieldsDeriv(u)
Jvs = range(self.survey.nFreq)
P = self.survey.projectFieldsDeriv(u)
for i, freqInd in enumerate(range(self.data.nFreq)):
for i, freqInd in enumerate(range(self.survey.nFreq)):
e = u.get_e(freqInd)
omega = self.data.omega[freqInd]
# for txInd in self.data.nTx
omega = self.survey.omega[freqInd]
# for txInd in self.survey.nTx
dMe_dsig = self.mesh.getEdgeInnerProductDeriv(m, v=e)
dsig_dm = self.model.transformDeriv(m)
b = 1j*omega * ( dMe_dsig * ( dsig_dm * v ) )
simpegEM/FDEM/DataFDEM.py → simpegEM/FDEM/SurveyFDEM.py
+5 -5
View File
@@ -1,10 +1,9 @@
from SimPEG import Utils, np, sp
from SimPEG.Data import BaseData
from SimPEG import Survey, Utils, np, sp
from FieldsFDEM import FieldsFDEM
class DataFDEM(BaseData):
class SurveyFDEM(Survey.BaseSurvey):
"""
docstring for DataFDEM
docstring for SurveyFDEM
"""
txLoc = None #: txLoc
@@ -25,7 +24,7 @@ class DataFDEM(BaseData):
return self.freq.size
def __init__(self, **kwargs):
BaseData.__init__(self, **kwargs)
Survey.BaseSurvey.__init__(self, **kwargs)
Utils.setKwargs(self, **kwargs)
@property
@@ -50,3 +49,4 @@ class DataFDEM(BaseData):
# Interpolation Matrices
####################################################
simpegEM/FDEM/__init__.py
+1 -1
View File
@@ -1,3 +1,3 @@
from FieldsFDEM import FieldsFDEM
from DataFDEM import DataFDEM
from SurveyFDEM import *
from FDEM import ProblemFDEM_e
simpegEM/TDEM/BaseTDEM.py
+10 -10
View File
@@ -13,21 +13,21 @@ class MixinInitialFieldCalc(object):
storeTheseFields = 'b'
def getInitialFields(self):
if self.data.txType == 'VMD_MVP':
if self.survey.txType == 'VMD_MVP':
# Vertical magnetic dipole, magnetic vector potential
F = self._getInitialFields_VMD_MVP()
else:
exStr = 'Invalid txType: ' + str(self.data.txType)
exStr = 'Invalid txType: ' + str(self.survey.txType)
raise Exception(exStr)
return F
def _getInitialFields_VMD_MVP(self):
if self.mesh._meshType is 'CYL1D':
MVP = Sources.MagneticDipoleVectorPotential(np.r_[0,0,self.data.txLoc], np.c_[np.zeros(self.mesh.nN), self.mesh.gridN], 'x')
MVP = Sources.MagneticDipoleVectorPotential(np.r_[0,0,self.survey.txLoc], np.c_[np.zeros(self.mesh.nN), self.mesh.gridN], 'x')
elif self.mesh._meshType is 'TENSOR':
MVPx = Sources.MagneticDipoleVectorPotential(self.data.txLoc, self.mesh.gridEx, 'x')
MVPy = Sources.MagneticDipoleVectorPotential(self.data.txLoc, self.mesh.gridEy, 'y')
MVPz = Sources.MagneticDipoleVectorPotential(self.data.txLoc, self.mesh.gridEz, 'z')
MVPx = Sources.MagneticDipoleVectorPotential(self.survey.txLoc, self.mesh.gridEx, 'x')
MVPy = Sources.MagneticDipoleVectorPotential(self.survey.txLoc, self.mesh.gridEy, 'y')
MVPz = Sources.MagneticDipoleVectorPotential(self.survey.txLoc, self.mesh.gridEz, 'z')
MVP = np.concatenate((MVPx, MVPy, MVPz))
# Initialize field object
@@ -86,8 +86,8 @@ class MixinTimeStuff(object):
class ProblemBaseTDEM(MixinTimeStuff, MixinInitialFieldCalc, BaseProblem):
"""docstring for ProblemTDEM1D"""
def __init__(self, mesh, model, **kwargs):
BaseProblem.__init__(self, mesh, model, **kwargs)
def __init__(self, model, **kwargs):
BaseProblem.__init__(self, model, **kwargs)
####################################################
@@ -144,7 +144,7 @@ class ProblemBaseTDEM(MixinTimeStuff, MixinInitialFieldCalc, BaseProblem):
def forward(self, m, RHS, CalcFields, F=None):
if F is None:
F = FieldsTDEM(self.mesh, self.data.nTx, self.nTimes, store=self.storeTheseFields)
F = FieldsTDEM(self.mesh, self.survey.nTx, self.nTimes, store=self.storeTheseFields)
dtFact = None
for tInd, t in enumerate(self.times):
@@ -165,7 +165,7 @@ class ProblemBaseTDEM(MixinTimeStuff, MixinInitialFieldCalc, BaseProblem):
def adjoint(self, m, RHS, CalcFields, F=None):
if F is None:
F = FieldsTDEM(self.mesh, self.data.nTx, self.nTimes, store=self.storeTheseFields)
F = FieldsTDEM(self.mesh, self.survey.nTx, self.nTimes, store=self.storeTheseFields)
dtFact = None
for tInd, t in reversed(list(enumerate(self.times))):
simpegEM/TDEM/DataTDEM.py → simpegEM/TDEM/SurveyTDEM.py
+4 -4
View File
@@ -1,10 +1,10 @@
from SimPEG import Utils, np
from SimPEG.Data import BaseData
from SimPEG.Survey import BaseSurvey
from FieldsTDEM import FieldsTDEM
class DataTDEM1D(BaseData):
class SurveyTDEM1D(BaseSurvey):
"""
docstring for DataTDEM1D
docstring for SurveyTDEM1D
"""
txLoc = None #: txLoc
@@ -20,7 +20,7 @@ class DataTDEM1D(BaseData):
return self.timeCh.size
def __init__(self, **kwargs):
BaseData.__init__(self, **kwargs)
BaseSurvey.__init__(self, **kwargs)
Utils.setKwargs(self, **kwargs)
def projectFields(self, u):
simpegEM/TDEM/TDEM_b.py
+7 -7
View File
@@ -2,7 +2,7 @@ from BaseTDEM import ProblemBaseTDEM
from FieldsTDEM import FieldsTDEM
from SimPEG.Utils import mkvc
import numpy as np
from DataTDEM import DataTDEM1D
from SurveyTDEM import SurveyTDEM1D
class ProblemTDEM_b(ProblemBaseTDEM):
"""
@@ -16,12 +16,12 @@ class ProblemTDEM_b(ProblemBaseTDEM):
with \\\(\\b\\\) defined on cell faces and \\\(\e\\\) defined on edges.
"""
def __init__(self, mesh, model, **kwargs):
ProblemBaseTDEM.__init__(self, mesh, model, **kwargs)
def __init__(self, model, **kwargs):
ProblemBaseTDEM.__init__(self, model, **kwargs)
solType = 'b'
dataPair = DataTDEM1D
surveyPair = SurveyTDEM1D
####################################################
# Internal Methods
@@ -51,12 +51,12 @@ class ProblemTDEM_b(ProblemBaseTDEM):
u = self.fields(m)
p = self.Gvec(m, v, u)
y = self.solveAh(m, p)
return self.data.dpred(m, u=y)
return self.survey.dpred(m, u=y)
def Jtvec(self, m, v, u=None):
if u is None:
u = self.fields(m)
p = self.data.projectFieldsAdjoint(v)
p = self.survey.projectFieldsAdjoint(v)
y = self.solveAht(m, p)
w = self.Gtvec(m, y, u)
return w
@@ -88,7 +88,7 @@ class ProblemTDEM_b(ProblemBaseTDEM):
def Gtvec(self, m, v, u=None):
if u is None:
u = self.fields(m)
tmp = np.zeros((self.mesh.nE,self.data.nTx))
tmp = np.zeros((self.mesh.nE,self.survey.nTx))
for i in range(self.nTimes):
tmp += v.get_e(i)*u.get_e(i)
p = -mkvc(self.model.transformDeriv(m).T*self.mesh.getEdgeMassDeriv().T*tmp)
simpegEM/TDEM/__init__.py
+2 -2
View File
@@ -1,4 +1,4 @@
from BaseTDEM import ProblemBaseTDEM
from DataTDEM import DataTDEM1D
from SurveyTDEM import SurveyTDEM1D
from FieldsTDEM import FieldsTDEM
from TDEM_b import ProblemTDEM_b
from TDEM_b import ProblemTDEM_b
simpegEM/Tests/test_FDEM.py
+2 -2
View File
@@ -27,9 +27,9 @@ class TDEM_bDerivTests(unittest.TestCase):
'rxType':'bz',
'freq': np.logspace(0,3,4),
}
dat = EM.FDEM.DataFDEM(**opts)
dat = EM.FDEM.SurveyFDEM(**opts)
prb = EM.FDEM.ProblemFDEM_e(mesh, model)
prb = EM.FDEM.ProblemFDEM_e(model)
prb.pair(dat)
sigma = np.log(np.ones(mesh.nC)*1e-3)
simpegEM/Tests/test_TDEM_b_DerivAdjoint.py
+12 -12
View File
@@ -13,10 +13,10 @@ class TDEM_bDerivTests(unittest.TestCase):
hx = Utils.meshTensors(((0,cs), (ncx,cs), (npad,cs)))
hy = Utils.meshTensors(((npad,cs), (ncy,cs), (npad,cs)))
mesh = Mesh.Cyl1DMesh([hx,hy], -hy.sum()/2)
active = mesh.vectorCCz<0.
model = Model.ActiveModel(mesh, active, -8, nC=mesh.nCz)
model = Model.ComboModel(mesh,
model = Model.ActiveModel(mesh, active, -8, nC=mesh.nCz)
model = Model.ComboModel(mesh,
[Model.LogModel, Model.Vertical1DModel, model])
@@ -26,9 +26,9 @@ class TDEM_bDerivTests(unittest.TestCase):
'rxType':'bz',
'timeCh':np.logspace(-4,-2,20),
}
self.dat = EM.TDEM.DataTDEM1D(**opts)
self.dat = EM.TDEM.SurveyTDEM1D(**opts)
self.prb = EM.TDEM.ProblemTDEM_b(mesh, model)
self.prb = EM.TDEM.ProblemTDEM_b(model)
self.prb.setTimes([1e-5, 5e-5, 2.5e-4], [10, 10, 10])
self.sigma = np.ones(mesh.nCz)*1e-8
@@ -45,21 +45,21 @@ class TDEM_bDerivTests(unittest.TestCase):
prb = self.prb
sigma = self.sigma
u = prb.fields(sigma)
Ahu = prb.AhVec(sigma, u)
V1 = Ahu.get_b(0)
V2 = 1/prb.getDt(0)*prb.MfMui*u.get_b(-1)
self.assertTrue(np.linalg.norm(V1-V2)/np.linalg.norm(V2) < 1.e-6)
V1 = Ahu.get_e(0)
self.assertTrue(np.linalg.norm(V1) < 1.e-6)
for i in range(1,u.nTimes):
dt = prb.getDt(i)
V1 = Ahu.get_b(i)
V2 = 1/dt*prb.MfMui*u.get_b(i-1)
self.assertTrue(np.linalg.norm(V1)/np.linalg.norm(V2) < 1.e-6)
@@ -175,7 +175,7 @@ class TDEM_bDerivTests(unittest.TestCase):
d_sig = 10*np.random.rand(prb.model.nP)
derChk = lambda m: [prb.data.dpred(m), lambda mx: -prb.Jvec(sigma, mx)]
derChk = lambda m: [prb.survey.dpred(m), lambda mx: -prb.Jvec(sigma, mx)]
print '\n'
print 'test_Deriv_J'
passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=6, eps=1e-20)
@@ -286,7 +286,7 @@ class TDEM_bDerivTests(unittest.TestCase):
V1 = d.dot(prb.Jvec(sigma, m))
V2 = m.dot(prb.Jtvec(sigma, d))
self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6)
if __name__ == '__main__':
simpegEM/Tests/test_forward_EMproblem.py
+2 -2
View File
@@ -27,9 +27,9 @@ class TDEM_bTests(unittest.TestCase):
'rxType':'bz',
'timeCh':np.logspace(-4,-2,20),
}
self.dat = EM.TDEM.DataTDEM1D(**opts)
self.dat = EM.TDEM.SurveyTDEM1D(**opts)
self.prb = EM.TDEM.ProblemTDEM_b(mesh, model)
self.prb = EM.TDEM.ProblemTDEM_b(model)
self.prb.setTimes([1e-5, 5e-5, 2.5e-4], [150, 150, 150])
self.sigma = np.ones(mesh.nCz)*1e-8